ANÁLISE NÃO LINEAR GEOMÉTRICA DE ESTRUTURAS

UTILIZANDO UMA FORMULAÇÃO CO-ROTACIONAL

RENATO CÉSAR GAVAZZA MENIN

EXAME DE QUALIFICAÇÃO DE DOUTORADO

DEPARTAMENTO DE ENGENHARIA CIVIL E AMBIENTAL

FACULDADE DE TECNOLOGIA

UNIVERSIDADE DE BRASÍLIA

UNIVERSIDADE DE BRASÍLIA FACULDADE DE TECNOLOGIA

DEPARTAMENTO DE ENGENHARIA CIVIL E AMBIENTAL PROGRAMA EM ESTRUTURAS E CONSTRUÇÃO CIVIL

ANÁLISE NÃO LINEAR GEOMÉTRICA DE ESTRUTURAS UTILIZANDO UMA FORMULAÇÃO CO-ROTACIONAL

RENATO CÉSAR GAVAZZA MENIN ORIENTADOR: WILLIAM TAYLOR MATIAS SILVA

EXAME DE QUALIFICAÇÃO DE DOUTORADO

BRASÍLIA / DF: MARÇO DE 2004.

ii

ÍNDICE Capítulo Página 1. INTRODUÇÃO 1

1.1 MOTIVAÇÃO 1 1.2 FORMULAÇÃO CO-ROTACIONAL 1 1.3 ALGORITMOS DE SOLUÇÃO NÃO LINEAR 3 1.4 OBJETIVOS 3 1.5 DESCRIÇÃO DO TRABALHO 4

2. FORMULAÇÃO CO-ROTACIONAL / TRELIÇA 2D E 3D 5 2.1 INTRODUÇÃO 5 2.2 DESCRIÇÃO CINEMÁTICA 6 2.2.1 Sistema de coordenadas 7 2.3 DESLOCAMENTOS DEFORMACIONAIS 8 2.3.1 Movimento deformacional em função dos deslocamentos nodais 9 2.4 ENERGIA DE DEFORMAÇÃO 13 2.4.1 Vetor de forças internas 14 2.4.2 Matriz de rigidez tangente 15 3. FORMULAÇÃO CO-ROTACIONAL / PÓRTICO PLANO 17 3.1 MODELOS MATEMÁTICOS DE ELEMENTOS DE VIGA 17 3.2 DESCRIÇÃO CINEMÁTICA 18 3.2.1 Sistema de coordenadas 19 3.3 DESLOCAMENTOS DEFORMACIONAIS 19 3.3.1 Movimento deformacional em função dos deslocamentos nodais 20 3.3.2 Derivadas parciais dos deslocamentos deformacionais 23 3.4 ESFORÇOS RESULTANTES 24 3.5 ENERGIA DE DEFORMAÇÃO DA VIGA 24 3.6 VETOR DE FORÇA INTERNAS 25 3.7 MATRIZ DE RIGIDEZ TANGENTE 26 3.8 ÂNGULOS DE ROTAÇÃO NA FORMULAÇÃO CO-ROTACIONAL 27 4. MATRIZ DE ROTAÇÃO NO ESPAÇO 31 4.1 ROTAÇÕES FINITAS NO ESPAÇO 31 4.2 MATRIZ DE ROTAÇÃO PARA PEQUENAS ROTAÇÕES 35 5. FORMULAÇÃO CO-ROTACIONAL / PÓRTICO ESPACIAL 37 5.1 DESCRIÇÃO CINEMÁTICA 37 5.2 EQUAÇÕES DE EQUILÍBRIO 40 5.3 MUDANÇA DA VARIÁVEL ITERATIVA DE ROTAÇÃO 42 5.4 OPERADOR DE PROJEÇÃO 46 5.5 VETOR DE FORÇAS INTERNAS E MATRIZ DE RIGIDEZ TANGENTE 51 5.6 ESFORÇOS RESULTANTES 55

iii

6. FORMULAÇÃO CO-ROTACIONAL / ELEMENTO DE CASCA TRIANGULAR 56 6.1 DESCRIÇÃO CINEMÁTICA 56 6.1.1 Vetor de deslocamentos deformacionais 59 6.2 VARIAÇÃO DOS DESLOCAMENTOS DEFORMACIONAIS 61 6.2.1 Variação da matriz de transformação Tn 61 6.2.2 Variação da matriz de rotação de corpo rígido R0n 61 6.2.3 Variação de ud com relação a v 62 6.2.4 Variação de ud co-rotacional com relação a v 63 6.2.5 Variação de θd co-rotacional com relação a wd co-rotacional 64 6.2.6 Variação de θd co-rotacional com relação a v 65 6.2.7 Variação de vd co-rotacional com relação a v 65 6.3 VARIAÇÕES DE ALTA ORDEM 66 6.4 VETOR DE FORÇAS INTERNAS 67 6.5 MATRIZ DE RIGIDEZ TANGENTE CONSISTENTE 69 6.5.1 Matriz de rigidez tangente material KM 69 6.5.2 Matriz de rigidez tangente geométrica KGR 69 6.5.3 Matriz de rigidez tangente geométrica KGM 70 6.5.4 Matriz de rigidez tangente geométrica KGP 71 6.5.5 Matriz de rigidez tangente 71 7. ELEMENTO DE CASCA TRIANGULAR 72 7.1 INTRODUÇÃO 72 7.2 MATRIZ DE RIGIDEZ BÁSICA DO ELEMENTO ANDES 72 7.3 MATRIZ DE RIGIDEZ DE ALTA ORDEM DO ELEMENTO ANDES 73 7.4 DEFINIÇÕES GEOMÉTRICAS PARA O ELEMENTO TRIANGULAR 75 7.5 ELEMENTO TRIANGULAR DE MEMBRANA 77 7.5.1 Matriz de rigidez básica 77 7.5.2 Matriz de rigidez de alta ordem 78 7.5.3 Matriz de rigidez 80 7.6 ELEMENTO TRIANGULAR À FLEXÃO 81 7.6.1 Matriz de rigidez básica 81 7.6.2 Matriz de rigidez de alta ordem 82 7.7 EXTENSÕES NÃO-LINEARES PARA O ELEMENTO TRIANGULAR 84 8. EXEMPLOS NUMÉRICOS 85 8.1 TRELIÇAS PLANAS E ESPACIAIS 85 8.1.1 Estrutura articulada 2D não simétrica 85 8.1.2 Estrutura articulada 3D em forma de torre 87 8.1.3 Estrutura articulada 3D em forma de cúpula / Papadrakakis [1987] 88 8.1.4 Estrutura articulada 3D em forma de cúpula / Choong & Hangai [1993] 89 8.2 PÓRTICOS PLANOS E ESPACIAIS 90 8.2.1 Viga em balanço com momento fletor aplicado na extremidade livre 91 8.2.2 Viga em balanço com carregamento transversal na extremidade livre 92 8.2.3 Viga em balanço com carregamento axial aplicado na extremidade livre 93 8.2.4 Pórtico de Lee 95 8.2.5 Arcos circulares de grande altura 96 8.2.6 Viga engastada espacial com curvatura de 45 graus no plano 98 8.3 ANÁLISE DE RESULTADOS 99 REFERÊNCIAS BIBLIOGRAFICAS 100

iv

LISTA DE FIGURAS Figura Página 2.1 – Elemento finito de barra articulado nas configurações inicial e atual 6 2.2 – Posição de uma partícula PR na configuração co-rotacionada CR 8 2.3 – Movimento do elemento de barra no plano 10 3.1 – Elemento de viga C1 nas configurações inicial e atual 18 3.2 – Posição de uma partícula PR na configuração co-rotacionada CR 20 3.3 – Deslocamentos globais e deslocamentos globais rotacionados 21 3.4 – Deslocamentos deformacionais no sistema local 22 3.5 – Esforços resultantes e convenções de sinais positivos 24 3.6 – Casos especiais 1, 2 e 3 de rotações deformacionais 29 3.7 – Casos especiais 4, 5, 6 e 7 de rotações deformacionais 30 4.1 – Movimento de corpo rígido no espaço 31 4.2 – Caráter não vetorial das rotações no espaço 33 5.1 – Elemento de pórtico no espaço 37 5.2 – Translações e rotações de um elemento genérico de pórtico 39 5.3 – Esforços resultantes e convenção de sinais positivos 55 6.1 – Configuração inicial, co-rotacionada e deformada do elemento de casca 56 6.2 – Translações nodais no elemento de casca triangular 57 6.3 – Rotações nodais no elemento de casca triangular 59 7.1 – Dimensões geométricas e vetores unitários no elemento triangular 75 7.2 – Sistema de coordenadas triangulares 76 8.1 – Estrutura articulada 2D não simétrica 85 8.2 – Trajetória não linear de equilíbrio para a estrutura articulada 2D 86 8.3 – CST e função sign para a estrutura articulada 2D 86 8.4 – Estrutura articulada 3D em forma de torre / Vista em planta e corte 87 8.5 – Trajetória de equilíbrio e curva CST x Carregamento 87 8.6 – Estrutura em forma de cúpula / Papadrakakis [1987] 88 8.7 – Trajetória de equilíbrio e curva CST x Deslocamento 88 8.8 – Estrutura articulada 3D em forma de cúpula / Choong & Hangai [1993] 89 8.9 – Trajetória não linear de equilíbrio para a estrutura 89 8.10 – CST e função sign para cúpula / Choong & Hangai [1993] 90 8.11 – Viga em balanço com momento aplicado na extremidade livre 91 8.12 – Configurações deformadas da estrutura 91 8.13 – Trajetórias de equilíbrio da extremidade livre da viga 92 8.14 – Viga engastada e livre, carregamento e propriedades geométricas e mecânicas 92 8.15 – Configurações deformadas da estrutura 92 8.16 – Trajetórias de equilíbrio da extremidade livre 93 8.17 – Esforços internos: força axial, cortante e momento fletor 93

v

8.18 – Viga em balanço submetida a carregamento axial na extremidade livre 94 8.19 – Configuração deformada da estrutura 94 8.20 – Trajetórias de equilíbrio da extremidade livre 94 8.21 – Pórtico de Lee 95 8.22 – (a) Trajetória de equilíbrio e (b) Configurações deformadas 95 8.23 – CST, função sign e número de pivots negativos 96 8.24 – (a) Arco biengastado e (b) Arco rotulado-engastado 96 8.25 – (a) Trajetórias de equilíbrio e (b) Curva CST x carregamento 97 8.26 – Configurações deformadas (a) Arco biengastado e (b) Arco rotulado-engastado 98 8.27 – Configuração inicial, deformada e propriedades mecânicas e geométricas 98 8.28 – Trajetória de equilíbrio (a) direção Y e (b) direções X e Z 99

LISTA DE TABELAS Tabela Página 8.1 – Programas computacionais 85 8.2 – Carregamentos em Newton correspondentes aos pontos limites 97 8.3 – Deslocamentos da extremidade da viga 99

1

CAPÍTULO 1 INTRODUÇÃO

1.1 MOTIVAÇÃO Hoje em dia se presencia o aumento da utilização de estruturas cada vez mais esbeltas em várias áreas da engenharia, tais como: edificações, pontes, cascos de navios, fuselagens de aviões, cúpulas de coberturas e estruturas off-shore. Porém, devido a esta esbeltes, que é possível graças à utilização de materiais com alta resistência e baixo peso próprio, estas estruturas podem estar sujeitas a fenômenos de instabilidade de equilíbrio, que podem ocorrer localmente ou de maneira global. Portanto, é necessário que o engenheiro tenha ferramentas que sejam capazes de realizar uma análise qualitativa e quantitativa do comportamento destas estruturas, tanto na fase pré-crítica, na qual estes fenômenos ainda não ocorreram, quanto na fase posterior a perda de estabilidade de equilíbrio, denominada fase pós-crítica. O fato de uma sistema estrutural apresentar instabilidade de equilíbrio não implica necessariamente que o mesmo tenha perdido a sua capacidade portante. A perda ou não desta capacidade portante está intimamente relacionada com a natureza da instabilidade de equilíbrio que possa ocorrer no sistema. Desta maneira, se torna necessário conhecer a natureza deste fenômeno para melhor avaliar o desempenho da capacidade portante de uma estrutura na fase pós-crítica. Nos estudos relacionados a estes tipos de fenômenos, observa-se que em um grande número de casos a estrutura ou componente estrutural se comporta elasticamente mesmo na fase pós-crítica de modo que ocorrem apenas não-linearidades geométricas (grandes deslocamentos e pequenas deformações) e com isso possibilitando que se adote como hipótese simplificada que as deformações sejam pequenas ou mesmo infinitesimais, usualmente dentro do regime elástico. Esta simplificação resulta em uma série de benefícios na construção dos modelos de elementos finitos para a análise de instabilidade, de modo a permitir o uso de modelos lineares para obter a resposta deformacional do sistema ao passo que as grandes translações e rotações de corpo rígido, que caracterizam a não linearidade geométrica nas fases pré e pós-crítica possam ser tratadas separadamente.

A formulação co-rotacional para a análise não linear geométrica de estruturas é portanto baseada na separação explícita dos movimentos de corpo rígido (translações e rotações) e dos movimentos deformacionais. Esta separação segrega a não-linearidade aos movimentos de corpo rígido, de modo a permitir a reutilização de modelos lineares de elementos finitos já existentes, respeitando-se certas limitações de modelagem (pequenas deformações), conforme comentado anteriormente. No presente trabalho, a formulação co-rotacional foi utilizada como descrição cinemática para avaliar o comportamento não linear de diversos tipos de estruturas planas e espaciais discretizadas por elementos finitos de treliças, vigas ou cascas planas triangulares. Portanto, este trabalho representa a fusão de dois grandes tópicos independentes: a formulação co-rotacional e a implementação computacional de algoritmos de continuação da resposta estrutural para solução de sistemas não lineares. Um pequeno relato histórico destes tópicos será apresentado nas próximas seções. 1.2 FORMULAÇÃO CO-ROTACIONAL

O conceito da formulação co-rotacional foi introduzido inicialmente por Wempner [1969] que desenvolveu uma formulação para o estudo de cascas sujeitas a pequenas deformações e grandes deslocamentos e por Belytschko & Hsieh [1973] que estudaram elementos de viga sujeitos a grandes rotações e propuseram um método baseado em um sistema de coordenadas curvilinear denominado “convected coordinates”.

2

Posteriormente, Fraeijs de Verbeke [1976] desenvolveu uma formulação co-rotacional para a análise dinâmica de estruturas, porém utilizando um único sistema de eixos co-rotacionais para a estrutura como um todo, estando mais voltada para uma solução analítica do problema do que para uma formulação de elementos finitos, sendo denominada “shadow element”. A determinação deste sistema de eixos único para a estrutura como um todo criava uma série de dificuldades de modo que o conceito da configuração fantasma ou shadow element, foi levado para o nível do elemento por vários pesquisadores, dentre os quais se destacam Bergan & Horrigmoe [1976] e Bergan & Nygard [1989]. Nos trabalhos de Bergan e Nygard, o conceito da configuração fantasma transformou-se em uma ferramenta de visualização muito útil que facilitou o entendimento da formulação co-rotacional. Este conceito foi usado pelos autores para eliminar os movimentos de corpo rígido de cada um dos elementos e obter apenas o movimento deformacional, a partir do qual pode ser computado o vetor de forças internas do elemento. Entretanto, as derivadas do vetor de forças internas não foram usadas diretamente na formação da matriz de rigidez tangente, fato que conduziu a uma perda de consistência.

Outra importante contribuição é atribuída à Rankin & Brogan [1986], que introduziram a chamada formulação EICR (Element Independent Corotational Formulation), que foi posteriormente refinada por Rankin & Nour-Omid [1988] e por Nour-Omid & Rankin [1991]. Esta formulação não utiliza explicitamente o conceito do “shadow element”, mas o caminho para a obtenção dos deslocamentos deformacionais, que se baseia no uso de operadores de projeção, também chamados projetores, é bastante similar ao processo utilizado por Bergan & Nygard [1989]. Estes autores usaram a formulação co-rotacional diretamente para formar a matriz de rigidez tangente, proporcionando uma matriz de rigidez consistente.

Entretanto, a formulação de Nour-Omid e Rankin [1991] ainda apresentava restrições no número de graus de liberdade que poderiam participar na rotação do sistema de coordenadas do elemento e ao mesmo tempo manter a consistência da matriz de rigidez tangente. Para resolver este problema, Haugen [1994] desenvolveu um trabalho aplicado para o estudo de cascas planas discretizadas por elementos triangulares e quadrangulares que continham o grau de liberdade de rotação torsional (drilling), combinando as principais características das duas formulações anteriores (shadow element e EICR), ou seja, mesclando a natureza invariável da formulação de Bergan e o equilíbrio e a consistência da formulação de Rankin.

Outras contribuições importantes são atribuídas a Cole [1990] e Crisfield [1990] que apresentaram formulações consistentes para a análise não linear geométrica de pórticos planos e espacias; Crisfield & Shi [1994] que apresentaram um procedimento para a análise dinâmica não-linear de treliças planas segundo a formulação co-rotacional; Pacoste & Eriksson [1996] que estudaram problemas de instabilidade para elementos de viga no plano e no espaço, comparando as descrições lagrangiana total e co-rotacional e posteriormente Pacoste [1998] que fez estudos de instabilidade de cascas utilizando elementos finitos planos e triangulares de casca contendo três nós e seis graus de liberdade por nó, seguindo basicamente a formulação descrita por Nour-Omid & Rankin [1991] através da utilização de projetores, porém implementando uma parametrização das rotações finitas no espaço que leva a uma mudança de variáveis adicional, de modo que as variáveis relacionadas à rotações no espaço se tornem aditivas e com isso tornando desnecessário eventuais procedimentos de atualização; Cardona [1989] que utilizou o conceito da formulação co-rotacional para o estudo de mecanismos; Menin & Taylor [2003a] que estudaram o comportamento pós-crítico de sistemas de barras articuladas no plano e no espaço utilizando distintas medidas de deformações; Menin & Taylor [2003b] que estudaram problemas de instabilidade de pórticos planos discretizados com elementos de viga de Euler-Bernoulli; e finalmente Menin & Taylor [2004] e Monteiro [2004] que estudaram problemas de não-linearidade geométrica de pórticos espaciais baseando-se no conceito de operadores de projeção.

3

1.3 ALGORITMOS DE SOLUÇÃO NÃO LINEAR Antes da metade dos anos 70, problemas estruturais envolvendo não-linearidade

geométrica eram geralmente tratados com métodos puramente incrementais sob controle de carga. Estes métodos têm a grande desvantagem de poder desviar a solução da trajetória de equilíbrio. O erro associado neste caso é dependente do passo de carga e freqüentemente é acumulativo durante a análise, tanto que requer um passo de carga muito pequeno para uma análise mais precisa. Somando-se a este fato, a possibilidade de resolução do sistema ir além de um ponto crítico utilizando controle de carga era muita pequena ou mesmo impossível. Esta dificuldade motivou o desenvolvimento de métodos incrementais-iterativos, nos quais os incrementos foram seguidos pelas iterações de correção do equilíbrio. Estas correções trazem a solução de volta para a trajetória de equilíbrio e o algoritmo é menos dependente do tamanho do passo de carga utilizado.

Em uma análise envolvendo não-linearidade geométrica, as estruturas usualmente alcançam um nível de carga máximo, no qual se tornam incapazes de resistir a mais incrementos de carga, até que uma mudança significativa na geometria ocorra. Estes níveis de carga são chamados de pontos críticos (limites ou de bifurcação) sendo caracterizados por uma matriz de rigidez tangente singular. Um método baseado em controle de carga pode ser capaz de detectar um ponto limite mas em geral não é capaz de ir além deste ponto. A necessidade de atravessar um ponto limite é motivada pelo fato de que em muitos casos a estrutura possui capacidade resistente de carga que pode ainda ser aproveitada. Existem na literatura diversas estratégias baseadas no controle de carga-deslocamento, nos quais tanto a carga quanto o deslocamento podem variar simultaneamente durante os passos incrementais, permitindo que os algoritmos sejam capazes de atravessar um ponto limite e obter a continuação da resposta, dentre os quais se destacam o controle de deslocamento e os métodos de comprimento de arco. Segundo Haugen [1994] nenhum algoritmo é aplicável a todos os tipos de problemas, porém os algoritmos do tipo comprimento de arco são geralmente considerados os mais versáteis em termos de alcance de problemas que eles podem resolver. No presente trabalho foi utilizado o método de comprimento de arco cilíndrico conforme descrito por Crisfield [1991]. 1.4 OBJETIVOS O presente trabalho tem por objetivo mostrar a evolução dos conhecimentos adquiridos pelo autor em relação à formulação co-rotacional para a análise de problemas envolvendo não linearidade geométrica de estruturas planas e espaciais discretizadas com elementos finitos de treliças, vigas ou cascas, durante o período que vai de Março/2003 à Fevereiro/2004.

Terminada esta primeira etapa do programa de doutorado realizada no Departamento de Engenharia Estrutural da Universidade de Brasília (UnB/Brasil), o próximo objetivo desta pesquisa será o estudo do comportamento pós-crítico de estruturas discretizadas com elementos finitos de cascas planas triangulares utilizando a descrição cinemática co-rotacional descrita por Haugen [1994] para uma análise não linear estática e estender estes conceitos para efetuar o desenvolvimento, aplicação e avaliação de elementos finitos de cascas de baixa ordem e alta performance para obter a resposta não linear deste tipo de estruturas em uma análise dinâmica, durante o período que vai de Março/2004 à Fevereiro/2005. Esta segunda etapa a ser realizada após a defesa do exame de qualificação de doutorado deverá ser completada no Departamento de Engenharia Aeroespacial da Universidade do Colorado (CAS Center for Aerospace Structures) em Boulder/USA, sob a orientação do Prof. Carlos Felippa, utilizando-se a chamada bolsa sanduíche (PDEE), durante um período de 12 meses.

4

Em relação à instituição de destino, ou seja, o CAS (Center for Aerospace Structures), este é pertencente ao Departamento de Engenharia Aeroespacial da Universidade do Colorado em Boulder nos Estados Unidos. O centro foi criado em 1986 para conduzir pesquisas sobre simulações computacionais para o projeto, análise e controle de estruturas aeroespaciais. Atualmente o CAS concentra os estudos em uma série de assuntos dentre os quais pode-se destacar: mecânica estrutural analítica e computacional, projeto de estruturas aeroespaciais, implementação de métodos de simulação em programação paralela, problemas de interação fluido-estrutura, aeroelasticidade, métodos variacionais aplicados a mecânica estrutural, termomecânica, dinâmica dos fluidos computacional, modelos de turbulência, métodos de decomposição de domínio e subestruturação, monitoramento e detecção de danos estruturais e dinâmica de estruturas aeroespaciais. Dentre os laboratórios que o CAS possui, destacam-se o Laboratório de Dinâmica Estrutural e Controle, o Laboratório de Computação de Alta Performance e o Laboratório de Simulação de Vôo. 1.5 DESCRIÇÃO DO TRABALHO O presente trabalho é composto de oito capítulos, procurando-se mostrar a evolução do autor em relação à formulação co-rotacional, em ordem crescente de complexidade. O capítulo 2 é destinado a descrição cinemática de treliças planas e espaciais utilizando-se diferentes medidas de deformações, com base nos trabalhos apresentados por Felippa [2001] e Menin & Taylor [2003a].

No capítulo 3 é desenvolvida a formulação co-rotacional para elementos de pórticos planos discretizados utilizando-se o modelo matemático de viga de Euler-Bernoulli (C1), sem o acoplamento dos efeitos do esforço axial e de flexão, conforme os trabalhos de Felippa [2001], Cortivo & Taylor [2002], e Menin & Taylor [2003b]. O capítulo 4 é destinado a obtenção da matriz de rotação que descreve o movimento de corpo rígido no espaço, cujo comportamento esta fundamentado no teorema de Euler, bem como na obtenção do pseudo-vetor de rotação no espaço pelo algoritmo atribuído à Spurrier, descrito em Monteiro [2004], sendo ambos utilizados nos capítulos seguintes para o estudo de grandes rotações envolvendo elementos de pórtico espacial e cascas. Uma vez conhecida a matriz de rotação no espaço de acordo com o capítulo anterior, no capítulo 5 é feita a descrição cinemática da formulação co-rotacional para o caso de elementos de pórtico espaciais baseando-se em especial nos trabalhos apresentados por Nour-Omid & Rankin [1991] e Monteiro [2004]. No capítulo 6 é apresentada a descrição cinemática co-rotacional para um elemento de casca plana e triangular, segundo a formulação CSSE (Consistent Symmetrizable Self-Equilibrated), descrita por Haugen [1994]. O capítulo 7 é destinado a descrição do elemento finito de casca linear que fornece a matriz de rigidez linear Ke e o vetor de forças internas fe, utilizados na formulação co-rotacional de cascas apresentada no capítulo anterior. O elemento finito linear de casca triangular utilizado no presente trabalho é baseado na formulação ANDES (Assumed Natural Deviatoric Strains), conforme apresentado nos trabalhos de Haugen [1994] e Felippa e Militello [1992], valendo ressaltar que o elemento possui o grau de liberdade de rotação torsional (drilling) como parte integrante do elemento de membrana. Finalmente, no capítulo 8 é feita uma pequena descrição dos programas computacionais desenvolvidos para realizar as análises e são apresentados uma série de exemplos numéricos para comprovar a eficiência da formulação co-rotacional para a análise não linear geométrica de estruturas.

5

CAPÍTULO 2 FORMULAÇÃO CO-ROTACIONAL / TRELIÇA 2D E 3D

2.1 INTRODUÇÃO

Atualmente, são geralmente usados três tipos de descrições cinemáticas na análise não linear geométrica através do método dos elementos finitos, que podem ser distinguidas entre si basicamente pela escolha da configuração de referência. A primeira delas é a descrição lagrangiana total, na qual a configuração de referência é raramente ou nunca mudada e em geral é igual à configuração inicial ao longo de toda a análise, sendo as tensões e deformações medidas com relação à esta configuração. A segunda descrição é a lagrangiana atualizada, para a qual a última configuração em equilíbrio, uma vez atingida passa a ser a próxima configuração de referência para os passos subsequentes, sendo as tensões e deformações redefinidas assim que a configuração de referência é atualizada. Já na chamada descrição co-rotacional, a configuração de referência é dividida em duas partes, sendo as tensões e deformações medidas à partir de uma configuração co-rotacionada, ao passo que a configuração inicial é mantida como configuração de referência para medir os deslocamentos de corpo rígido.

Segundo Felippa [2001], a descrição lagrangiana total permanece sendo a formulação mais utilizada hoje em dia, ao passo que o interesse pela descrição lagrangiana atualizada está diminuindo bastante e sendo gradualmente substituída pela descrição co-rotacional. No presente capítulo, será apresentada a descrição cinemática referente à formulação co-rotacional, e mais especificamente para o caso de elementos de barra bi-articulados (treliças) no plano e no espaço de acordo com Felippa [2001] e Menin & Taylor [2003a].

O principal conceito desta formulação é a divisão ou decomposição da configuração de referência em duas parcelas:

1. A configuração inicial (C0) que é mantida fixa ao longo de toda a análise, servindo como uma configuração de referência fixa. Usualmente se adota um sistema de coordenadas global para toda a estrutura.

2. A configuração co-rotacionada (CR) que varia de elemento para elemento. Para cada elemento a configuração CR pode ser obtida através do deslocamento de corpo rígido em relação à configuração C0. O sistema de coordenadas se move conjuntamente com o elemento, sendo a deformação do elemento medida em relação ao sistema de coordenadas local da configuração CR.

No presente trabalho, são utilizadas quatro medidas distintas de deformações, sendo duas delas descritas em coordenadas materiais (deformação de engenharia e Green-Lagrange) e duas descritas em coordenadas espaciais (deformação de Biot e Almansi). A descrição cinemática do elemento é feita em relação à configuração inicial ou indeformada, supondo uma relação linear entre o par conjugado de tensão e deformação do tipo σ = Eε, sendo E o módulo de elasticidade do material e adotando-se o mesmo valor de E para as distintas medidas de deformações. Será demonstrado através dos exemplos numéricos que no caso de deformações infinitesimais, as configurações inicial e atual se confundem e portanto se obtém unicidade na resposta de uma estrutura independentemente da configuração em que se escolhe o modelo constitutivo e do tipo de deformação que se utilize, porém no regime de deformações finitas esta hipótese implica na definição de materiais diferentes e em conseqüência não se obtém unicidade na resposta. Para obter unicidade na resposta em regime de deformações finitas é necessário recorrer a transformações tensoriais que fazem o mapeamento do tensor constitutivo entre as configurações inicial ou indeformada e a atual, tema que não será tratado neste trabalho.

6

2.2 DESCRIÇÃO CINEMÁTICA Considerando-se um elemento finito de barra articulado que se move no espaço de acordo com a Fig. (2.1), por simplicidade será admitido como hipótese que os eixos locais (x0

e,y0e,z0

e) do elemento na configuração inicial C0 estão alinhados com os sistemas de coordenadas globais material e espacial, designados por (X,Y,Z) e (x,y,z) respectivamente. Somando-se a esta primeira hipótese, também é admitido que a origem do sistema de eixos locais em C0 está situado na metade do comprimento inicial do elemento, designado por L0. O elemento de barra se move da configuração inicial C0 até a configuração atual C, cujos eixos locais podem ser definidos como (xe,ye,ze). A chamada configuração co-rotacionada, designada por CR se move conjuntamente com o elemento até a configuração C, posicionando-se simetricamente com respeito a configuração atual. Pode ser observado também pela Fig. (2.1) que os eixos locais co-rotacionados (xR

e,yRe,zR

e) coincidem com os eixos locais (xe,ye,ze) em C.

Figura 2.1 – Elemento finito de barra articulado nas configurações inicial e atual

Tomando-se uma partícula P0 de coordenadas (X,Y) em C0, que se move ao ponto PR de coordenadas (xR,yR) em CR, e em seguida se move para o ponto P de coordenadas (x,y) em C, então, o deslocamento total u da partícula, em coordenadas globais pode ser descrito por: u = x – X (2.1) Este deslocamento pode então ser decomposto em uma parte deformacional uD e outra que corresponde ao deslocamento de corpo rígido uR de modo que: u = uR + uD = (xR – X) + (x – xR) (2.2) Na formulação co-rotacional, as equações do movimento deformacional são escritas em função das coordenadas locais (xe,ye,ze) em C, conforme a seguinte equação: uDe = Q.uD (2.3) sendo Q uma matriz de rotação 2x2 (treliça plana) ou 3x3 (treliça espacial) para transformar do sistema global (X,Y,Z) ao sistema local (xe,ye,ze).

ψ (+)

Re yy ,

Re xx ,

Re xxxX ,,, 0

Re zzzZ ,,, 0

Re yyyY ,,, 0

0C0O X

0P

uRuR

e zz ,C

RC RO

RP DuP

RP Du P

RxRu

x

0O X0P

u0u

e

e

e

7

Os deslocamentos deformacionais uDe são utilizados para obter o vetor de forças internas e a matriz de rigidez tangente, conforme será comentado posteriormente. 2.2.1 Sistema de Coordenadas

Os sistemas de coordenadas local (xe,ye,ze) na configuração atual C e global (x,y,z) se relacionam através da seguinte equação: xe = Q.(x – u0) (2.4) Esta relação expressa na equação acima pode ser visualizada através da Fig. (2.1), sendo u0 o vetor que representa o deslocamento do ponto O0 em C0 ao ponto O em C. A matriz de rotação Q que aparece nas Eq. (2.3) e (2.4) pode ser definida, segundo Gere & Weaver [1981], no caso de treliça plana como:

Q =

− xy

yx

CCCC

(2.5)

Já no caso de treliça espacial, a matriz de rotação Q pode ser expressa, também segundo Gere & Weaver [1981] por:

Q =

++−

+−+

+−

2222

22

22

22

0zx

x

zx

z

zx

zyzx

zx

yx

zyx

CC

C

CC

CCC

CCCC

CC

CC

CCC

(2.6)

Deve-se ressaltar que a matriz Q apresentada na Eq. (2.6) é valida para todas as posições do elemento de barra no espaço, exceto quando o elemento está alinhado com o eixo Y, para o qual Cx = 0 e Cz = 0, e a matriz de rotação Q resulta em:

Q =

−

1000000

y

y

CC

(2.7)

sendo (Cx, Cy) no caso de treliça plana e (Cx, Cy, Cz) no caso de treliça espacial, os cosenos diretores do elemento de barra na configuração atual C (direção do eixo local xe), em relação ao sistema global de coordenadas, conforme será comentado posteriormente. Uma vez que a matriz Q é uma matriz ortogonal, ou seja QTQ = QQT = I, sendo I a matriz identidade, então a inversa da Eq. (2.4) pode ser escrita como: x = QTxe + u0 (2.8)

8

2.3 DESLOCAMENTOS DEFORMACIONAIS Nos passos seguintes, será apresentada a obtenção dos deslocamentos deformacionais em coordenadas locais, definido anteriormente na Eq. (2.3). No caso de treliças planas, as coordenadas das partículas PR em CR e P em C são definidas pelas seguintes equações:

xR =

+

−=

0

0.vu

YX

CCCC

yx

xy

yx

R

R = QTX + u0 (2.9)

x =

++

=

vYuX

yx

= X + u = IX + (2.10)

A interpretação geométrica da Eq. (2.9) para o caso de uma treliça plana pode ser vista na Fig. (2.2). Nesta figura os cosenos diretores (Cx, Cy) do elemento de barra na configuração atual C (direção do eixo local xe) são calculados em função do ângulo ψ entre os eixos locais x0

e e xe no sentido anti-horário, sendo designados respectivamente por (cosψ, senψ).

Figura 2.2 – Posição de uma partícula PR na configuração co-rotacionada CR Lembrando que uD = x – xR, conforme definido na Eq. (2.2) e substituindo-se os valores de x e xR definidos respectivamente pelas Eq.(2.9) e (2.10), obtém-se:

uD =

−−

+

−−

−=

−−

=

0

0.1

1vvuu

YX

CCCC

yyxx

vu

xy

yx

R

R

D

D = (I – QT)X + u – u0 (2.11)

Por último, pode-se obter o deslocamento deformacional em relação às coordenadas locais através da transformação de coordenadas dada pela Eq. (2.3): uDe = Q.uD = (Q – I)X + Q(u – u0) (2.12)

Re yyyY ,,, 0

ψ (+)

0u0P

ψsen.Xψcos.Y

ψsen.Y

ψcos.X

Re xx ,

Re yy ,

RP

ex0

Ry0v

RORC

0C

0O XY

Rx

e

e

9

No caso de treliças espaciais, as coordenadas das partículas PR em CR e P em C são definidas de forma análoga ao caso bidimensional pelas seguintes equações:

xR =

+

=

0

0

0

.wvu

ZYX

zyx

T

R

R

R

Q = QTX + u0 (2.13)

x =

+++

=

wZvYuX

zyx

= X + u = IX + u (2.14)

Lembrando que uD = x – xR, conforme definido na Eq. (2.2) e substituindo-se os valores

de x e xR definidos respectivamente pelas Eq.(2.13) e (2.14), obtém-se:

uD = =

−−−

=

R

R

R

D

D

D

zzyyxx

wvu

(I – QT).

−−−

+

0

0

0

wwvvuu

ZYX

= (I – QT)X + u – u0 (2.15)

sendo a matriz de rotação QT que aparece nas Eq. (2.13) e (2.15), definida pela transposta da matriz de rotação expressa na equação (2.6) ou (2.7), conforme a posição da barra em relação ao eixo global Y. Vale lembrar que o deslocamento deformacional em relação às coordenadas locais no caso tridimensional também pode ser expresso pela Eq. (2.12) definida para o caso bidimensional, bastando usar os vetores e matrizes correspondentes no espaço. 2.3.1 Movimento deformacional em função dos deslocamentos nodais Neste item, os deslocamentos deformacionais definidos anteriormente para um ponto genérico são particularizados para os nós das extremidades dos elementos de barra. No caso específico de treliças planas, as coordenadas nodais do elemento em C0 com relação aos eixos locais são X2 = – X1 = ½ L0 e Y2 = Y1 = 0, sendo L0 o comprimento do elemento nesta configuração. Os deslocamentos dos nós de extremidade podem ser então definidos por:

u =

( )( )( )( )

( )( )( )( )

−−

=

=

=

0,210,210,210,21

,,,,

0

0

0

0

22

22

11

11

2

2

1

1

LvLuLvLu

YXvYXuYXvYXu

vuvu

2

1

uu

(2.16)

De maneira similar, o movimento deformacional em função dos deslocamentos nodais pode ser expresso da seguinte forma:

uDe =

( )( )( )( )

( )( )( )( )

−−

=

=

=

0,210,210,210,21

,,,,

0

0

0

0

22

22

11

11

2

2

1

1

LvLuLvLu

YXvYXuYXvYXu

vuvu

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

2

1

uu (2.17)

10

Uma vez definida a Eq. (2.17), a Eq.(2.12) pode então ser re-escrita em função dos deslocamentos nodais:

uDe =

−−

−

+

−−−−

−

−=

y

x

y

x

xy

yx

xy

yx

eD

eD

eD

eD

CC

CC

L

vvuuvvuu

CCCC

CCCC

vuvu

1

1

21

0000

0000

0

02

02

01

01

2

2

1

1

(2.18)

Uma vez que o campo de deslocamentos do elemento é linear em X e em Y, o elemento permanece reto na configuração atual C e portanto se pode escrever que u0 = ½ .(u1 + u2), v0 = ½ .(v1 + v2) (2.19) O próximo passo é definir os valores dos cosenos diretores (Cx, Cy) em função dos deslocamentos nodais, e em seguida achar o comprimento do elemento (L) na configuração atual, conforme apresentado na Fig. (2.3). Vale ressaltar que nesta dedução, a configuração inicial já não se encontra mais alinhada com os eixos globais.

Figura 2.3 – Movimento do elemento de barra no plano

Inicialmente é feito a rotação dos deslocamentos nodais em relação ao sistema de eixos locais na configuração inicial (x0

e,y0e,z0

e), definido em função do ângulo α:

u21rot =

−

=

−−

=

21

21

021021

021021

12

12

21

21

vu

LXLYLYLX

vvuu

vu

rotrot

rotrot

rot

rot

(2.20)

( )( )

−=−=

−==−==

1221

1221

012021

012021

sencos

vvvuuu

LYYLYLXXLX

αα

(2.21)

L

X

Y

rotv21

rotu2

roto uL 21+

rotv2

exey

ψ (+)

rotu1

eoy e

ox

oL

rotv1α

11

Uma vez conhecidos os deslocamentos nodais rotacionados, segundo a Eq. (2.20), é possível definir as demais variáveis cinemáticas envolvidas na formulação co-rotacional em função das relações geométricas apresentadas na Fig. (2.3):

Cx = cosψ = LuL rot

210 +;

Lv

Crot

y21sen == ψ (2.22)

( ) ( )2

21

2

210rotrot vuLL ++= (2.23)

Nas próximas seções, serão obtidos o vetor de forças internas e a matriz de rigidez tangente dos elementos através das derivadas primeira e segunda do funcional da energia de deformação respectivamente. Assim, será calculada a derivada primeira de L em relação aos deslocamentos nodais u, de modo que:

ψcos12

=∂∂

−=∂∂

uL

uL , ψsen

12

=∂∂

−=∂∂

vL

vL (2.24)

cuja forma vetorial se expressa como

−−

=∂∂

ψψψψ

sencossencos

uL (2.25)

Utilizando as Eq. (2.22), (2.23) e (2.24) e a relação sen2ψ + cos2ψ = 1, pode-se obter a derivada segunda de L em relação aos deslocamentos nodais u:

−−−−

−−−−

=∂∂

ψψψψψψψψψψψψ

ψψψψψψψψψψψψ

22

22

22

22

2

2

coscossencoscossencossensencossensen

coscossencoscossencossensencossensen

1L

Lu

(2.26)

O procedimento utilizado para o caso bidimensional pode também ser estendido de forma análoga para o caso de uma treliça espacial, na qual as coordenadas nodais do elemento em C0 com relação aos eixos locais são X2 = – X1 = ½ L0, Y2 = Y1 = 0 e Z1 = Z2 = 0 sendo L0 o comprimento do elemento nesta configuração. Os deslocamentos dos nós de extremidade podem ser então definidos por:

u =

( )( )( )( )( )( )

( )( )( )( )( )( )

−−−

=

=

=

0,0,210,0,210,0,210,0,210,0,210,0,21

,,,,,,,,,,,,

0

0

0

0

0

0

222

222

222

111

111

111

2

2

2

1

1

1

2

1

LwLvLuLwLvLu

ZYXwZYXvZYXuZYXwZYXvZYXu

wvuwvu

uu

(2.27)

12

De maneira similar, o movimento deformacional em função dos deslocamentos nodais pode ser expresso da seguinte forma para o caso tridimensional:

uDe =

( )( )( )( )( )( )

( )( )( )( )( )( )

−−−

=

=

=

0,0,210,0,210,0,210,0,210,0,210,0,21

,,,,,,,,,,,,

0

0

0

0

0

0

222

222

222

111

111

111

2

2

2

1

1

1

2

1

LwLvLuLwLvLu

ZYXwZYXvZYXuZYXwZYXvZYXu

wvuwvu

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

uu (2.28)

Uma vez definida a Eq. (2.28), a Eq.(2.12) pode então ser re-escrita para o caso

tridimensional em função dos deslocamentos nodais:

uDe =

−

−+

−−−−−−

=

2

2

2

1

1

1

02

02

02

01

01

01

2

2

2

1

1

1

ZYXZYX

wwvvuuwwvvuu

wvuwvu

eD

eD

eD

eD

eD

eD

IQ00IQ

Q00Q

3x33x3

3x33x3

3x33x3

3x33x3 (2.29)

sendo Q3x3 a matriz de rotação dada pela Eq.(2.6) ou (2.7) conforme a posição do elemento de barra, 03x3 uma matriz quadrada nula de ordem 3 e I a matriz identidade de ordem 3. Uma vez que o campo de deslocamentos do elemento é linear em X, Y e Z o elemento permanece reto na configuração atual C e portanto se pode escrever que u0 = ½ .(u1 + u2), v0 = ½ .(v1 + v2), w0 = ½ .(w1 + w2) (2.30) O próximo passo é definir os valores dos cosenos diretores (Cx, Cy, Cz) em função dos deslocamentos nodais, e em seguida achar o comprimento do elemento (L) na configuração atual, lembrando que nesta dedução, a configuração inicial já não se encontra mais alinhada com os eixos globais. Assim como no caso bidimensional, inicialmente é feita a rotação dos deslocamentos nodais em relação ao sistema de eixos locais definidos na configuração inicial:

u21rot =

=

−−−

=

21

21

21

12

12

12

21

21

21

wvu

wwvvuu

wvu

rotrot

rotrot

rotrot

rot

rot

rot

Q (2.31)

sendo Q a matriz de rotação definida pelas equações (2.6) ou (2.7), empregando-se os seguintes cosenos diretores:

( )( )( )

−==−==−==

012021

012021

012021

LZZLZCLYYLYC

LXXLXC

rotz

roty

rotx

(2.32)

13

Uma vez conhecidos os deslocamentos nodais rotacionados, segundo a Eq. (2.31), é possível definir as demais variáveis cinemáticas envolvidas na formulação co-rotacional de elementos de treliças espaciais em função de relações geométricas análogas às apresentadas para o caso bidimensional na Fig. (2.3):

Cx = LuL rot

210 +;

Lv

Crot

y21= ;

Lw

Crot

z21= (2.33)

( ) ( ) ( )2

21

2

21

2

210rotrotrot wvuLL +++= (2.34)

Nas próximas seções, serão obtidos o vetor de forças internas e a matriz de rigidez

tangente dos elementos através das derivadas primeira e segunda do funcional da energia de deformação respectivamente. Assim, será calculada a derivada primeira de L em relação aos deslocamentos nodais u, utilizando-se as Eq. (2.33) e (2.34) de modo que:

xCuL

uL

=∂∂

−=∂∂

12

, yCvL

vL

=∂∂

−=∂∂

12

, zCwL

wL

=∂∂

−=∂∂

12

(2.35)

cuja forma vetorial se expressa como

−−−

=∂∂

z

y

x

z

y

x

CCCCCC

Lu

(2.36)

Utilizando as Eq. (2.33), (2.34) e (2.35) e a relação Cx

2+ Cy2+ Cz

2 = 1, pode-se obter a derivada segunda de L em relação aos deslocamentos nodais u:

( ) ( )( ) ( )

( ) ( )( ) ( )

( ) ( )( ) ( )

+−−+−

−+−+−

−−++−

+−+−−

+−−+−

+−−−+

=∂∂

2222

2222

2222

2222

2222

2222

2

2 1

yxzyzxyxzyzx

zyzxyxzyzxyx

zxyxzyzxyxzy

yxzyzxyxzyzx

zyzxyxzyzxyx

zxyxzyzxyxzy

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCC

LLu

(2.37)

2.4 ENERGIA DE DEFORMAÇÃO Admitindo uma relação linear entre o par conjugado de tensão e deformação, pode-se descrever a seguinte relação entre tensões e deformações nas configurações inicial e atual: σX = E.εX em C0 (2.38) σx = E.εx em C (2.39)

14

sendo E o módulo de elasticidade longitudinal do material. Considerando-se esta hipótese e definindo a área da seção transversal dos elementos em C0 e C como sendo respectivamente A0 e A, a energia de deformação de um elemento de treliça nas configurações inicial e atual podem ser definidas da seguinte forma:

U0 = ∫0

0

20.

21 L

eX dXEA ε em C0 (2.40)

U = ∫L

ex dxEA

0

2.21 ε em C (2.41)

No presente trabalho, foram utilizadas quatro medidas distintas de deformação, sendo duas delas descritas em coordenadas materiais (deformação de engenharia e Green-Lagrange) e duas em coordenadas espaciais (Biot e Almansi), sendo expressas respectivamente pelas seguintes equações:

( )

( )

−=−

=

−=−

=

−=−

=

−=−

=

−

−

22

20

2

10

22

0

20

20

0

121

2

1

121

2

1

λε

λε

λε

λε

LLL

LLL

LLL

LLL

almansix

biotx

greenX

engX

(2.42)

2.4.1 Vetor de forças internas

O vetor de forças internas descrito em coordenadas materiais pode ser obtido através da primeira derivada do funcional da energia de deformação em relação aos deslocamentos nodais, de modo que, derivando a Eq. (2.40) se obtém a seguinte expressão:

∫ ∫ ∂∂

=∂

∂=

∂∂

=0 0

0 000

0 221ˆ

L LeX

XeX

X dXEAdXEAU

uuuf e ε

εε

ε (2.43)

Já o vetor de forças internas descrito em coordenadas espaciais pode ser obtido de forma análoga ao anterior, derivando a Eq. (2.41) em relação aos deslocamentos nodais

∫ ∫ ∂∂

=∂∂

=∂∂

=L L

exx

exx dxEAdxEAU

0 0

221

uuuf e ε

εε

ε (2.44)

O esforço axial nos elementos de treliça nas configurações inicial e atual podem ser definidos respectivamente da seguinte forma: N0 = EA0εX, N = EAεx (2.45)

15

Levando-se em conta a Eq. (2.25) para o caso bidimensional ou Eq. (2.36) para o caso tridimensional, além das Eq. (2.42) e (2.45) e efetuando a integração das Eq. (2.43) e (2.44), se chega às expressões do vetor de forças internas em coordenadas materiais e espaciais dados respectivamente por:

u

f e

∂∂

=LN 00

ˆ β , u

f e

∂∂

=LNβ ,

⇒=⇒=

−⇒=⇒=

−

−

AlmansiBiotLagrangeGreen

engenhariadef

20

10

0

0 .1

λβλβ

λββ

(2.46)

sendo os coeficientes β0 e β resultantes da integração das Eq. (2.43) e (2.44), em função das diferentes medidas de deformação. É importante ressaltar que no item 2 foi suposto que o sistema de eixos locais do elemento na configuração inicial C0 esta alinhado com os sistemas de eixos globais material e espacial. Em uma formulação mais geral, se supõe que existe uma certa inclinação entre estes sistemas de eixos e portanto se deve escrever o vetor de forças internas em coordenadas materiais e espaciais, em relação ao sistema de eixos globais através da seguinte relação: eTe

g fRf ˆˆ = , eTeg fRf = (2.47)

= T

TT

Q00Q

R (2.48)

onde RT é a matriz de rotação que transforma do sistema de coordenadas local para o sistema de coordenadas global, 0 é uma matriz nula e QT, a transposta da matriz Q, expressa na Eq. (2.5) para o caso bidimensional ou Eq. (2.6) e (2.7) para o caso tridimensional, com os cosenos diretores definidos nas Eq. (2.21) ou (2.32), conforme se trate de treliça plana ou de treliça espacial respectivamente. 2.4.2 Matriz de rigidez tangente A matriz de rigidez tangente em relação às coordenadas materiais pode ser obtida pelo cálculo da segunda derivada do funcional da energia de deformação definida em C0 com relação aos deslocamentos nodais;

eL

XX

TXX dXEA

U∫

∂∂

+

∂∂

∂∂

=∂

∂=

0

02

2

020

2ˆ

uuuuK e ε

εεε (2.49)

Levando-se em conta as Eq. (2.25) e (2.26) para o caso bidimensional ou Eq. (2.36) e (2.37) para o caso tridimensional, além das Eq. (2.42) e (2.45) e efetuando a integração da Eq. (2.49), se chega à expressão da matriz de rigidez tangente em coordenadas materiais:

2

2

0000

0ˆuuu

K e

∂∂

+

∂∂

∂∂

=LNLL

LEA T

βα ( )

−⇒−=⇒=

.2/131

20

0

LagrGreenengenharia

λαα

(2.50)

16

Procedendo de forma análoga, a matriz de rigidez tangente em coordenadas espaciais pode ser obtida pela segunda derivada do funcional da energia de deformação definida em C:

eL

xx

Txx dxEAU

∫

∂∂

+

∂∂

∂∂

=∂∂

=0

2

2

2

2

uuuuK e ε

εεε

(2.51)

Levando-se em conta as Eq. (2.25) e (2.26) para o caso bidimensional ou Eq. (2.36) e (2.37) para o caso tridimensional, além das Eq. (2.42) e (2.45) e efetuando a integração da Eq. (2.51), se chega à expressão da matriz de rigidez tangente em coordenadas espaciais:

2

2

uuuK e

∂∂

+

∂∂

∂∂

=LNLL

LEA T

βα ( )

⇒−=⇒−=

−−

−−

AlmansiBiot

2/352324

12

λλαλλα

(2.52)

Conforme comentado anteriormente para o vetor de forças internas, é necessário transformar do sistema de coordenadas local para o sistema de coordenadas global através da seguinte operação: RKRK eTe

gˆˆ = , RKRK eTe

g = (2.53) sendo RT a matriz de rotação definida na Eq. (2.48).

17

CAPÍTULO 3 FORMULAÇÃO CO-ROTACIONAL / PÓRTICO PLANO

3.1 MODELOS MATEMÁTICOS DE ELEMENTOS DE VIGA Na mecânica estrutural são geralmente utilizados dois modelos distintos para discretizar os elementos de viga que compõem as estruturas. O primeiro deles é referente à teoria clássica de viga, também chamado modelo de Euler-Bernoulli (C1), no qual as forças cortantes transversais são obtidas a partir do equilíbrio do elemento, porém o seu efeito no cálculo das deformações é desprezado. A hipótese fundamental deste modelo é que as seções transversais permanecem planas e normais ao eixo longitudinal deformado, com a rotação ocorrendo em torno do eixo neutro que passa através do centróide da seção transversal. O segundo é o modelo de Timoshenko (C0) que corrige a teoria clássica de viga com o efeito das deformações cisalhantes de primeira ordem. Nesta teoria, as seções transversais permanecem planas e sofrem rotações em torno do eixo neutro da mesma forma que ocorre no modelo de Euler-Bernoulli, porém elas não permanecem normais ao eixo longitudinal deformado, sendo esta diferença angular em relação ao eixo normal produzida pelo esforço cortante transversal que é assumido constante ao longo da seção transversal. De acordo com Faria [1998] e Felippa [2001], os dois modelos matemáticos descritos anteriormente estão baseados na hipótese de pequenas deformações e comportamento elástico e isotrópico do material. Além disso, ambos não consideram mudanças nas dimensões da seção transversal à medida que a viga se deforma. Estes modelos podem considerar a não linearidade geométrica devido ao efeito de grandes deslocamentos e rotações desde que as outras hipóteses comentadas anteriormente sejam obedecidas.

No presente capítulo, será apresentada a descrição cinemática referente à formulação co-rotacional, e mais especificamente para o caso de pórticos planos discretizados através de elementos de viga utilizando-se o modelo matemático de Euler-Bernoulli, sem o acoplamento dos efeitos do esforço axial e da flexão, conforme os trabalhos apresentados por Felippa [2001], Cortivo & Taylor [2002] e Menin & Taylor [2003b], procurando-se seguir o mesmo formato do capítulo anterior, no qual a formulação co-rotacional foi desenvolvida para discretizar elementos de treliça, de modo a facilitar a visualização das semelhanças e diferenças existentes entre o elemento de viga e o elemento de barra articulado.

Na descrição cinemática co-rotacional desenvolvida para o elemento de viga de Euler-Bernoulli, a configuração de referência é dividida em duas partes, sendo as tensões e deformações medidas à partir de uma configuração co-rotacionada, ao passo que a configuração inicial é mantida como configuração de referência para medir os deslocamentos de corpo rígido. Desta forma, o deslocamento total de cada nó da viga analisada pode ser decomposto em duas partes, sendo a primeira referente ao deslocamento de corpo rígido, no qual ocorrem apenas translações e rotações; e a segunda que esta associada aos deslocamentos deformacionais, consistindo de deformações na direção axial do elemento (alongamentos ou encurtamentos) e na deformação angular do elemento em torno do eixo z. Vale lembrar que são considerados rotações e translações finitas para cada nó do elemento de viga em análise, porém suas deformações são admitidas infinitesimais. A energia de deformação total é obtida pelo somatório das energias de deformação axial e de flexão isoladamente, sem consideração do acoplamento entre ambas, conforme comentado anteriormente.

No final deste capítulo, são descritas algumas informações importantes para o cálculo das rotações de corpo rígido e deformacional, necessárias na implementação computacional do elemento de viga de Euler-Bernoulli segundo a formulação co-rotacional.

18

3.2 DESCRIÇÃO CINEMÁTICA Considerando-se um elemento finito de viga prismático e reto que se move no plano de acordo com a Fig. (3.1), por simplicidade será admitido como hipótese inicial que os eixos locais (x0

e,y0e) do elemento na sua configuração inicial C0 estão alinhados com os sistemas de

coordenadas globais material e espacial, designados por (X,Y) e (x,y) respectivamente, tendo a sua origem O0 localizada na metade do comprimento inicial do elemento, definido por L0. O elemento de viga se move da configuração inicial C0 até a configuração atual C de modo que o eixo longitudinal passe através dos nós de extremidade na configuração atual, definindo assim o eixo local xe. A origem do sistema de eixos locais (xe,ye) é então definida na metade da distância entre os nós de extremidade, formando um ângulo ψ com o eixo local x0

e. A chamada configuração co-rotacionada CR se move conjuntamente com o elemento até a configuração C, posicionando-se simetricamente com respeito à configuração atual, de modo que os eixos locais co-rotacionados (xR

e,yRe) coincidam com os eixos locais (xe,ye) em C.

1

ψ (+)

xe,xRe

X,x,x0e,xR

Y,y,y0e,yR

ye,yRe

XO0

C0

CR

PR

P0

uD

uRu1

2

2

O0

P0

u

uDPRP

u0

OR

uR

x

xR

X

C

P

Figura 3.1 – Elemento finito de viga C1 nas configurações inicial e atual

Tomando-se uma partícula P0 de coordenadas (X,Y) em C0, que se move ao ponto PR de

coordenadas (xR,yR) em CR, e em seguida se move ao ponto P de coordenadas (x,y) em C, então, o deslocamento total u da partícula , em coordenadas globais pode ser descrito por:

u = x – X =

−−

θYyXx

(3.1)

sendo θ definido como sendo a rotação total, obtida pelo somatório da rotação de corpo rígido ψ com a rotação deformacional θ , ou seja: θ = ψ + θ (3.2) Este deslocamento u pode então ser decomposto em uma parte deformacional uD e outra que corresponde ao deslocamento de corpo rígido uR de modo que: u = uR + uD = (xR – X) + (x – xR) (3.3)

19

Na formulação co-rotacional, as equações do movimento deformacional são escritas em função das coordenadas locais (xe,ye) em C, conforme a seguinte equação: uD

e = Q.uD (3.4) sendo Q uma matriz de rotação 3x3 utilizada para transformar do sistema global (X,Y) ao sistema local (xe,ye). Os deslocamentos deformacionais uD

e são utilizados para obter o vetor de forças internas e a matriz de rigidez tangente, conforme será comentado posteriormente. 3.2.1 Sistema de coordenadas Os sistemas de coordenadas local (xe,ye) na configuração atual C e global (x,y) se relacionam através da seguinte equação: xe = Q.(x – u0) (3.5) Esta relação expressa na equação acima pode ser visualizada através da Fig. (3.1), sendo u0 o vetor que representa o deslocamento do ponto O0 em C0 ao ponto O em C. A matriz de rotação Q que aparece nas equações (3.4) e (3.5) pode ser definida, segundo Gere & Weaver [1981], no caso de pórticos planos como:

−=

10000

xy

yx

CCCC

Q (3.6)

sendo (Cx,Cy) os cosenos diretores do elemento de viga na configuração atual C (direção do eixo local xe), em relação ao sistema global de coordenadas, conforme será comentado posteriormente. Uma vez que a matriz Q é uma matriz ortogonal, ou seja QTQ = QQT = I, sendo I a matriz identidade de ordem 3, então a inversa da Eq. (3.5) pode ser escrita como: x = QTxe + u0 (3.7) 3.3 DESLOCAMENTOS DEFORMACIONAIS Nos passos seguintes, será apresentada a obtenção dos deslocamentos deformacionais em coordenadas locais, definido anteriormente na Eq. (3.4) por uD

e. No caso de um elemento de viga no plano, as coordenadas das partículas PR em CR e P em C são definidas pelas equações:

xR =

+=

+

−=+

ψθψθ Po

R

R

Po

xy

yx

yx

vu

YX

CCCC

0

0

10000

0T uXQ (3.8)

x = X + u = IX + u =

+=

+

θθθθ PoPo

yx

vu

YX

100010001

(3.9)

20

É importante enfatizar que as Eq. (3.8) e (3.9) definidas anteriormente, foram deduzidas para um ponto genérico qualquer P0 e que para um ponto situado sobre o eixo local x0

e na configuração inicial C0, conforme será o caso dos nós de extremidade do elemento de viga, o termo θPo será nulo, ou seja, θPo = 0. A interpretação geométrica da Eq. (3.8) pode ser vista na Fig. (3.2). Nesta figura os cosenos diretores (Cx,Cy) do elemento de viga na configuração atual (direção do eixo local xe) são calculados em função do ângulo ψ entre os eixos locais x0

e e xe no sentido anti-horário, sendo designados respectivamente por (cosψ,senψ).

X.cosψ

Y.senψ

Y.cosψX.senψ

ye,yRe

xe,xRe

PR 2ψ(+)

C0

CR

O0

OR

P0

1

1

2θPo

u0

xR

v0

yR

x0e

Figura 3.2 – Posição de uma partícula PR na configuração co-rotacionada CR

Lembrando que uD = x – xR, conforme definido na Eq. (3.3) e substituindo-se os valores de x e xR definidos respectivamente pelas Eq. (3.8) e (3.9), obtém-se:

uD = x – xR = (I – QT)X + u – u0 =

−−−

ψθR

R

yyxx

(3.10)

Finalmente, pode-se obter o deslocamento deformacional em relação às coordenadas locais através da transformação de coordenadas dada pela Eq. (3.4): uD

e = Q.uD = (Q – I)X + Q(u – u0) (3.11) 3.3.1 Movimento deformacional em função dos deslocamentos nodais Neste item, os deslocamentos deformacionais definidos anteriormente para um ponto genérico são particularizados para os nós das extremidades dos elementos definidos na Fig. (3.2) com os índices 1 e 2. Para o caso específico de elementos de vigas prismáticos e retos no plano, as coordenadas nodais do elemento em C0 com relação ao sistema de eixos locais (X,Y) são X2 = – X1 = ½ L0 e Y2 = Y1 = 0, sendo L0 o comprimento de elemento nesta configuração. Os deslocamentos totais u e deformacionais uD

e dos nós de extremidade podem ser então definidos por:

21

u =

( )( )( )( )( )( )

−−−

=

=

0,2/10,2/10,2/10,2/10,2/10,2/1

0

0

0

0

0

0

2

2

2

1

1

1

LLvLuLLvLu

vu

vu

θ

θ

θ

θ

2

1

uu

, uDe =

( )( )( )( )( )( )

−−−

=

=

0,2/10,2/10,2/10,2/10,2/10,2/1

0

0

0

0

0

0

2

2

2

1

1

1

LLvLuLLvLu

vu

vu

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

eD

θ

θ

θ

θ

2

1

uu (3.12)

Uma vez conhecida a Eq. (3.12), a Eq. (3.11) definida para um ponto genérico pode então ser re-escrita em função dos deslocamentos nodais:

uDe =

−−−−−−

−

−

+

−−

−

=

ψθ

ψθ

θ

θ

ψψ

ψψ

ψψ

ψψ

ψ

ψ

ψ

ψ

2

02

02

1

01

01

0

2

2

2

1

1

1

1000000000000000010000000000

0

10

1

21

vvuu

vvuu

CSSC

CSSC

SC

SC

L

vu

vu

eD

eD

eD

eD

(3.13)

sendo Cψ = cosψ e Sψ = senψ os cosenos diretores definidos em função da rotação de corpo rígido ψ. Já os deslocamentos translacionais do ponto O0 em C0 ao ponto O em C podem ser expressos pelas seguintes equações: u0 = ½ .(u1 + u2), v0 = ½ .(v1 + v2) (3.14) O próximo passo é definir os valores dos cosenos diretores (Cψ, Sψ) em função dos deslocamentos nodais, e em seguida achar o comprimento do elemento (L) na configuração atual, conforme apresentado na Fig. (3.3). Vale ressaltar que nesta dedução, com o intuito de se representar um caso mais geral, a configuração inicial já não se encontra mais alinhada com os eixos globais.

L0

L

ψθ1

θ2

u1 rot

u2 rot

v1 rot

v2 rot

L0 + u21 rot

v21 rot

x0 ey0

e

xe

α

y e

L0

L

α

ψθ1

θ2

v1

v2

u1

u2

XY

Figura 3.3 – Deslocamentos globais e deslocamentos globais rotacionados

22

Na Fig. (3.3) é mostrado o movimento do elemento de viga no plano, sendo que do lado esquerdo são indicados os deslocamentos globais, ou seja em relação ao sistema de eixos globais (X,Y) e do lado direito os deslocamentos globais rotacionados. Portanto, o primeiro passo é fazer a rotação dos deslocamentos nodais globais em relação ao sistema de eixos locais na configuração inicial (x0

e,y0e), definido em função do ângulo α, lembrando que não é

necessário transformar os ângulos de rotação:

u21rot =

−

=

−−

=

21

21

12

12

21

21

vu

CSSC

vvuu

vu

rotrot

rotrot

rot

rot

αα

αα (3.15)

( )( )

−=−=

−===−===

1221

1221

012021

012021

//sen//cos

vvvuuu

LYYLYSLXXLXC

αα

α

α

(3.16)

Uma vez conhecidos os deslocamentos nodais rotacionados, segundo a Eq. (3.15), é possível definir as demais variáveis cinemáticas envolvidas na formulação co-rotacional em função das relações geométricas apresentadas na Fig. (3.3):

Cψ = cos ψ = LuL rot

210 + (3.17)

Sψ = sen ψ = L

v rot21 (3.18)

L = ( ) ( )2

21

2

210rotrot vuL ++ (3.19)

Definidas as principais variáveis cinemáticas envolvidas na formulação co-rotacional, serão apresentados na Fig. (3.4), à seguir, os deslocamentos deformacionais em relação ao sistema de eixos locais na configuração atual C, de modo a possibilitar uma melhor visualização e entendimento dos mesmos. É importante salientar que são válidas as seguintes relações: u2D

e = – u1De = d/2, v2D

e = v1De = 0, ψθθ −= 11 e ψθθ −= 22 .

xe,xRe

ye,yRe

+ d/2

- d/2

21

1θ

2θ

Figura 3.4 – Deslocamentos deformacionais no sistema local

23

3.3.2 Derivadas parciais dos deslocamentos deformacionais

Nas próximas seções, serão obtidos o vetor de forças internas e a matriz de rigidez tangente dos elementos através das derivadas primeira e segunda do funcional da energia de deformação. À partir das Eq. (3.17), (3.18) e (3.19) pode-se definir as seguintes expressões:

ψCuL

uL

=∂∂

−=∂∂

12

; ψSvL

vL

=∂∂

−=∂∂

12

; 012

=∂∂

=∂∂

θθLL

L

SuC

uC 2

12

ψψψ =∂

∂−=

∂

∂;

LCS

vC

vC ψψψψ −=

∂

∂−=

∂

∂

12

; 012

=∂

∂=

∂

∂

θθψψ CC

(3.20)

LCS

uS

uS ψψψψ −=

∂

∂−=

∂

∂

12

; L

CvS

vS 2

12

ψψψ =∂

∂−=

∂

∂; 0

12

=∂

∂=

∂

∂

θθψψ SS

L

Suu

ψψψ−=

∂∂

−=∂∂

12

; L

Cvv

ψψψ=

∂∂

−=∂∂

12

; 012

=∂∂

=∂∂

θψ

θψ

Derivando-se cada um dos elementos do vetor uD

e de deslocamentos deformacionais, definido na Eq.(3.13) em relação aos deslocamentos globais u, obtém-se a seguinte expressão:

∂∂∂∂∂∂

−−

−−−−

−−

=

∂∂∂∂

∂∂

2

2

2

1

1

1

2

2

2

1

1

1

1//0//00000002/2/02/2/0//1//00000002/2/02/2/

θ

θ

θ

θ

ψψψψ

ψψψψ

ψψψψ

ψψψψ

vu

vu

LCLSLCLS

SCSCLCLSLCLS

SCSC

vu

vu

eD

eD

eD

eD

(3.21)

Lembrando que u2D

e = – u1De = ½ .d, a Eq. (3.21) pode ser re-escrita da seguinte forma:

∂∂∂∂∂∂

−−−−

−−=

∂∂∂

2

2

2

1

1

1

2

1

1//0//0//1//00

θ

θ

θθ

ψψψψ

ψψψψ

ψψψψ

vu

vu

LCLSLCLSLCLSLCLS

SCSCd (3.22)

Levando em conta as Eq. (3.20) e (3.22) definidas acima, as segundas derivadas dos

deslocamentos deformacionais podem ser definidas por:

−−−−

−−−−

=∂∂

00000000000000000000

1

22

22

22

22

2

2

ψψψψψψ

ψψψψψψ

ψψψψψψ

ψψψψψψ

CCSCCSCSSCSS

CCSCCSCSSCSS

Ld

u (3.23)

24

−−−−−−

−−−−−−

=∂∂

=∂∂

0000000202020200000002020202

1

2222

2222

2222

2222

222

2

21

2

ψψψψψψψψ

ψψψψψψψψ

ψψψψψψψψ

ψψψψψψψψ

θθ

CSSCCSCSSCCSCSCS

CSCSCSSCCSCSSCCS

Luu (3.24)

3.4 ESFORÇOS RESULTANTES Os esforços resultantes em cada um dos elementos de viga na configuração atual C são N, V, M1 e M2 , sendo N o esforço normal, V o esforço cortante e M1 e M2 os momentos fletores nas extremidades inicial e final. Os esforços N e V são constantes ao longo de todo o elemento ao passo que o momento fletor M = M(xe) varia linearmente ao longo do elemento por se tratar de um modelo hermitiano. Estes esforços, com as respectivas convenções de sinais são mostrados na Fig. (3.5) apresentada à seguir, sendo obtidos a partir das deformações, de acordo com as seguintes equações apresentadas por Harrison [1973]:

N = EA0ε = dL

EA

0

0 ; V = ( )210

21 6 θθ +=+

LLEI

LMM

(3.25)

M1 = ( )210

22 θθ +LEI ; M2 = ( )21

0

22 θθ +LEI

sendo E o módulo de elasticidade do material; A0 a área da seção transversal; I o momento de inércia da seção; L e L0 os comprimentos dos elementos nas configurações inicial e atual respectivamente; e ε a deformação de engenharia, ou seja, ε = (L – L0) / L0 = d / L0.

M1

M2

N

N

VV

C

Figura 3.5 – Esforços resultantes e convenções de sinais positivos

3.5 ENERGIA DE DEFORMAÇÃO DA VIGA A energia de deformação da viga, considerando-se apenas deformações infinitesimais pode ser expressa pela seguinte equação: U = U A + U F (3.26) sendo U A e U F as energias de deformação axial e de flexão respectivamente. Desta forma, adotam-se as seguintes expressões:

25

U A = 2

0

0200 2

121 d

LEA

LEA =ε (3.27)

U F = ( )2221

21

0

2 θθθθ ++LEI (3.28)

3.6 VETOR DE FORÇAS INTERNAS O vetor de forças internas f é obtido pela derivada primeira do funcional da energia de deformação U em relação aos deslocamentos globais u:

f = u∂

∂U (3.29)

Substituindo-se os valores U A e U F definidos nas Eq.(3.27) e (3.28) na Eq.(3.29) e usando-se as derivadas primeiras de d, 1θ e 2θ , chega-se às seguintes expressões:

f A = uuu ∂

∂=

∂∂

=

∂∂ dNdd

LEA

dL

EA

0

02

0

0

21 = N [ ]TSCSC 00 ψψψψ −− (3.30)

f F = ( )

∂

∂+

∂∂

+∂∂

+∂∂

=

++

∂∂

uuuuu2

22

1211

10

2221

21

0

2222 θθ

θθθ

θθθθθθθ

LEI

LEI

(3.31) f F = [ ]TMVCVSMVCVS 21 ψψψψ −− Uma vez definidos os vetores f A e f F, através das Eq.(3.30) e (3.31), pode-se obter o vetor de forças internas somando-se as contribuições devidas aos esforços axial e de flexão: f = f A + f F (3.32)

É importante ressaltar que no item 3.2 foi suposto que o sistema de eixos locais do elemento na configuração inicial C0 esta alinhado com os sistemas de eixos globais material e espacial. Em uma formulação mais geral, se supõe que existe uma certa inclinação entre estes sistemas de eixos definida em função do ângulo α, conforme apresentado na Fig. (3.3) e portanto se deve re-escrever o vetor de forças internas, em relação ao sistema de eixos globais através da seguinte relação: fRf T=g (3.33)

= T

TT

Q00Q

R (3.34)

onde RT é a matriz de rotação que transforma do sistema de coordenadas local para o sistema de coordenadas global, 0 é uma matriz nula de ordem 3 e QT, a transposta da matriz Q, expressa na Eq. (3.6), com os cosenos diretores definidos nas Eq. (3.16).

26

3.7 MATRIZ DE RIGIDEZ TANGENTE A matriz de rigidez tangente é obtida pela segunda derivada do funcional da energia de deformação U em relação ao vetor de deslocamentos globais u, ou de forma similar através da primeira derivada do vetor de forças internas f, podendo ser decomposta em duas partes: a matriz de rigidez material (KM) e a matriz de rigidez geométrica (KG).

K = GM KKuf

u+=

∂∂

=∂∂

2

2U (3.35)

Substituindo-se os valores de f A e f F apresentados nas Eq. (3.30) e (3.31) na Eq. (3.35) e usando-se as derivadas primeiras e segundas de d, 1θ e 2θ definidas nas Eq. (3.22), (3.23) e (3.24) chega-se às seguintes expressões:

∂∂

+

∂∂

∂∂

=∂

∂2

2

0

0

uuuuf A dddd

LEA T

= KM axial + KG

axial (3.36)

( ) ( )

∂

∂++

∂∂

+∂∂

=∂

∂uuuu

f F2

211

210

222 θθθ

θθθ

LEI = KM

flexão + KG flexão (3.37)

flexãoaxial

MMM KKK += ; flexãoG

axialGG KKK += (3.38)

−−−−

−−−−

=

00000000000000000000

22

22

22

22

0

0

ψψψψψψ

ψψψψψψ

ψψψψψψ

ψψψψψψ

SCSSCSCSCCSC

SCSSCSCSCCSC

LEAaxial

MK (3.39)

−−

−−−−

−−

−−

−−

−−−−

=

233

133

366366

366366

133

233

366366

366366

2

2

2

22

2

2

22

2

22

2

2

2

22

2

2

22

2

22

2

0

LC

LS

LC

LS

LC

LC

LCS

LC

LC

LCS

LS

LCS

LS

LS

LCS

LS

LC

LS

LC

LS

LC

LC

LCS

LC

LC

LCS

LS

LCS

LS

LS

LCS

LS

LEIflexão

ψψψψ

ψψψψψψψψ

ψψψψψψψψ

ψψψψ

ψψψψψψψψ

ψψψψψψψψ

MK (3.40)

27

−−−−

−−−−

=

00000000000000000000

22

22

22

22

ψψψψψψ

ψψψψψψ

ψψψψψψ

ψψψψψψ

CCSCCSCSSCSS

CCSCCSCSSCSS

LNaxial

GK (3.41)

−−−−−−

−−−−−−

=

0000000202020200000002020202

2222

2222

2222

2222

yy

yy

yy

yy

flexãoG

CSSCCSCSSCCSCSCS

CSCSCSSCCSCSSCCS

LV

ψψψψψψ

ψψψψψψ

ψψψψψψ

ψψψψψψ

K (3.42)

Conforme comentado anteriormente para o vetor de forças internas, é necessário transformar do sistema de coordenadas local para o sistema de coordenadas global através da operação apresentada à seguir, sendo RT a matriz de rotação definida na Eq. (3.34): KRRK T

g = (3.43) 3.8 ÂNGULOS DE ROTAÇÃO NA FORMULAÇÃO CO-ROTACIONAL Neste item, serão descritas algumas informações importantes para o cálculo das rotações de corpo rígido e deformacional, necessárias na implementação computacional do elemento de viga de Euler-Bernoulli segundo a formulação co-rotacional. Vale lembrar que as rotações totais θ podem ser obtidas diretamente a partir do vetor de deslocamentos globais u. Inicialmente é feita a determinação da rotação de corpo rígido (ψ), em função dos cosenos diretores Cψ e Sψ, descritos nas equações (3.17) e (3.18), definidos pela posição da viga na configuração deformada. A escolha de qual coseno diretor usar na determinação da rotação de corpo rígido é definida em função do quadrante no qual a viga se encontra na sua configuração deformada, lembrando que na maioria dos programas computacionais, dentre os quais está o Fortran, no qual foi feita a implementação computacional do elemento de viga, os valores do arco-seno (sen-1) são fornecidos no primeiro e quarto quadrantes, de modo que correspondem a ângulos entre –π/2 e π/2 ao passo que os valores de arco-coseno (cos-1) são fornecidos no primeiro e segundo quadrantes, compreendendo valores entre 0 e π. Em função destes fatores e seguindo a recomendação de Crisfield [1991], as rotações de corpo rígido foram calculadas conforme as equações (3.44) descritas a seguir: Cψ>0 e Sψ>0 ⇒ (quad = 1) e ψ = sen-1(Sψ) : 0 ≤ ψ ≤ π/2 Cψ<0 e Sψ>0 ⇒ (quad = 2) e ψ = cos-1(Cψ) : π/2 ≤ ψ ≤ π (3.44) Cψ<0 e Sψ<0 ⇒ (quad = 3) e ψ = -cos-1(Cψ) : -π ≤ ψ ≤ -π/2 Cψ>0 e Sψ<0 ⇒ (quad = 4) e ψ = sen-1(Sψ) : - π/2 ≤ ψ ≤ 0

28

Uma vez conhecida a rotação total (θ ) e a rotação de corpo rígido (ψ), o próximo passo é a determinação da rotação deformacional (θ ). É importante ressaltar que no presente trabalho, foram admitidas rotações totais de até 720 graus tanto no sentido anti-horário (rotações positivas), quanto no sentido horário (rotações negativas) e que as rotações deformacionais são pequenas e em geral não ultrapassam 10 graus, de modo que os valores das rotações de corpo rígido são bem próximos dos respectivos valores das rotações totais. Em função da maneira pela qual são calculadas as rotações de corpo rígido pelo Fortran através da equação (3.44) e devido ao fato das rotações poderem ser tanto positivas quanto negativas, são necessárias algumas transformações de ângulos de modo que a rotação deformacional seja calculada de forma correta. No caso de rotações totais de até 2π, ou seja 360 graus, são apresentadas nas equações (3.45) e (3.46) as fórmulas para o cálculo das rotações deformacionais em função do quadrante em que se encontra a rotação de corpo rígido, para o caso de rotações positivas ou negativas, respectivamente: (0 ≤ θ ≤ 2π) e (quad = 1) ⇒ ψθθ −= (0 ≤ θ ≤ 2π) e (quad = 2) ⇒ ψθθ −= (3.45) (0 ≤ θ ≤ 2π) e (quad = 3) ⇒ ( )πψθθ 2+−= (0 ≤ θ ≤ 2π) e (quad = 4) ⇒ ( )πψθθ 2+−= (-2π ≤ θ ≤ 0) e (quad = 1) ⇒ ( ) ψπθθ −+= 2 (-2π ≤ θ ≤ 0) e (quad = 2) ⇒ ( ) ψπθθ −+= 2 (3.46) (-2π ≤ θ ≤ 0) e (quad = 3) ⇒ ψθθ −= (-2π ≤ θ ≤ 0) e (quad = 4) ⇒ ψθθ −= Já para o caso de rotações totais entre 2π e 4π, ou seja entre 360 e 720 graus, são apresentadas nas equações (3.47) e (3.48) as fórmulas para o cálculo das rotações deformacionais em função do quadrante em que se encontra a rotação de corpo rígido, para o caso de rotações positivas ou negativas, respectivamente: (2π ≤ θ ≤ 4π) e (quad = 1) ⇒ ( )πψθθ 2+−= (2π ≤ θ ≤ 4π) e (quad = 2) ⇒ ( )πψθθ 2+−= (3.47) (2π ≤ θ ≤ 4π) e (quad = 3) ⇒ ( )πψθθ 4+−= (2π ≤ θ ≤ 4π) e (quad = 4) ⇒ ( )πψθθ 4+−= (-4π ≤ θ ≤ -2π) e (quad = 1) ⇒ ( ) ψπθθ −+= 4 (-4π ≤ θ ≤ -2π) e (quad = 2) ⇒ ( ) ψπθθ −+= 4 (3.48) (-4π ≤ θ ≤ -2π) e (quad = 3) ⇒ ( ) ψπθθ −+= 2 (-4π ≤ θ ≤ -2π) e (quad = 4) ⇒ ( ) ψπθθ −+= 2 Na determinação da rotação deformacional admitindo-se rotações totais de até 720 graus, conforme comentado acima, ainda existem alguns casos especiais de rotações que devem ser considerados na implementação computacional. Estes casos especiais envolvem basicamente a ocorrência simultânea de ângulos referentes a rotações totais e de corpo rígido no primeiro e quarto quadrantes respectivamente ou vice-versa.

29

Existem na realidade sete casos especiais que serão tratados a seguir. Os três primeiros casos, que podem ser visualizados através da Fig. (3.6), estão associados a ocorrência de ângulos de rotação de corpo rígido ou totais iguais a zero, sendo as rotações deformacionais calculadas segundo as equações (3.49) a (3.51) apresentada a seguir: Caso 1: (ψ = 0.0) ⇒ θθ = (3.49) Caso 2: (θ = 0.0) ⇒ ψθ −= (3.50) Caso 3: (ψ = 0.0 e θ = 0.0) ⇒ 0.0=θ (3.51)

negθ

negθθ

θ

0.0=θ

Caso 1 Caso 2

Caso 3

Figura 3.6 – Casos especiais 1, 2 e 3 de rotações deformacionais

Na Fig.(3.6), os ângulos referentes a rotações totais são indicados na cor verde, os referentes a rotações de corpo rígido em azul e as rotações deformacionais em vermelho, podendo estas últimas serem positivas (sentido anti-horário) ou negativos (sentido horário). Os demais 4 casos especiais envolvem a ocorrência simultânea de ângulos associados a rotações totais e de corpo rígido no primeiro e quarto quadrantes, para os quais as transformações propostas nas equações (3.45) a (3.48) conduzem a erros. Para resolver estes problemas, deve-se utilizar as seguintes equações na implementação computacional: Caso 4: (0 ≤ θ ≤ π/2) ⇒ ψθθ −= (3.52) (2π ≤ θ ≤ 2.5π) ⇒ ( )πψθθ 2+−= Caso 5: (1.5π ≤ θ ≤ 2π) ⇒ ( )πψθθ 2+−= (3.53) (3.5π ≤ θ ≤ 4π) ⇒ ( )πψθθ 4+−= Caso 6: (-2π ≤ θ ≤ -1.5π) ⇒ ( ) ψπθθ −+= 2 (3.54) (-4π ≤ θ ≤ -3.5π) ⇒ ( ) ψπθθ −+= 4 Caso 7: (-π/2 ≤ θ ≤0) ⇒ ψθθ −= (3.55) (-2.5π ≤ θ ≤-2π) ⇒ ( ) ψπθθ −+= 2

30

Os casos especiais 4 e 5 são referentes a rotações totais positivas ao passo que os casos especiais 6 e 7 estão associados a rotações totais negativas. A primeira linha de cada um destes casos se refere a rotações totais em valor absoluto menores que 2π e a segunda linha de cada caso se refere a rotações totais em valor absoluto entre 2π e 4π. Estes casos especiais podem ser visualizados através da Fig. (3.7), utilizando-se a mesma convenção de cores definida para a Fig. (3.6) e lembrando que as diferenças angulares entre as rotações totais e de corpo rígido apresentadas em ambas as figuras, encontram-se exageradas apenas para facilitar a visualização.

Caso 4

Caso 5

Caso 6

Caso 7

θ θ

negθ negθ

)(+θ

)(+θ

)(−θ

)(−θ

Figura 3.7 – Casos especiais 4, 5, 6 e 7 de rotações deformacionais

31

CAPÍTULO 4 MATRIZ DE ROTAÇÃO NO ESPAÇO

4.1 ROTAÇÕES FINITAS NO ESPAÇO As rotações no plano (2D), por exemplo no plano xy, podem ser definidas por um único escalar, ou seja, um ângulo de rotação θ em torno do eixo z, sendo válida a propriedade comutativa: θ1 + θ2 = θ2 + θ1, uma vez que os ângulos de rotação são apenas escalares. Por outro lado, o estudo das rotações no espaço é um pouco mais complicado, sendo o assunto fundamentado no teorema de Euler, apresentado por Felippa [2001]: O movimento genérico de um corpo rígido fixo por um ponto pode ser descrito por uma única rotação do corpo em torno de um eixo passando por tal ponto. Consequentemente, para se descrever uma rotação no espaço são necessários não só o ângulo mas também a direção ou eixo de rotação. Estes são normalmente os mesmos atributos que caracterizam um vetor, entretanto as rotações finitas no espaço não obedecem as leis do cálculo vetorial, em especial não sendo mais válida a propriedade comutativa, uma vez que a alteração da ordem de aplicação de duas rotações sucessivas não resulta na mesma resposta a menos que o eixo de rotação seja mantido fixo. A seguir será apresentada a obtenção da matriz de rotação no espaço conforme descrito por Monteiro [2004]. O movimento de dois pontos genéricos P e Q de um corpo rígido é mostrado na Fig. (4.1), sendo o movimento completo separado em duas etapas. Inicialmente, todos os pontos transladam da mesma quantidade dtrans, de modo que P e Q se deslocam para P’ e Q’ respectivamente. Em seguida, o corpo é submetido a uma rotação representado pelo ângulo θ em torno do eixo direcionado segundo a direção unitária e que passa pelo ponto P’, fazendo com que o ponto Q’ mude para Q’’.

e

drot

dtrans

rnro

P

P’

Q

Q’

Q’’

O

Q’

Q’’

O

a b

θ

R

x

y

z

Figura 4.1 – Movimento de corpo rígido no espaço

Adotando-se a notação r0 = P’Q’, rn = P’Q’’ e drot = Q’Q’’, conforme apresentados na Fig. (4.1) acima, a matriz de rotação Rθ e o vetor drot podem ser definidos respectivamente através das seguintes equações: rn = Rθ . r0 (4.1) drot = rn – r0 (4.2)

32

Por se tratar de um corpo rígido, a trajetória descrita pelo ponto Q entre as posições Q’ e Q’’, define um arco de circunferência de raio R que pertence ao plano perpendicular ao eixo de rotação cujo centro está localizado no ponto O. Através da Fig. (4.1), pode-se definir que: rn = r0 + drot = r0 + a + b (4.3) sendo os vetores a e b descritos pelas equações a seguir, onde α é o ângulo definido entre o eixo de rotação (e) e o vetor r0:

αsen0

0

0

0

rre

arere

aa×

=××

= (4.4)

a

aebaeaebb ×

=××

= (4.5)

Observando-se a Fig. (4.1), pode-se definir que: R=αsen0r θsenR=a ( )θcos1−= Rb (4.6) Substituindo-se os valores definidos pela Eq. (4.6) nas equações (4.4) e (4.5) os vetores a e b podem ser re-escritos da seguinte forma: a = senθ .e × r0 b = (1-cosθ).e × e × r0 (4.7) Inserindo os vetores a e b, definidos acima, na Eq. (4.3) e utilizando-se a seguinte propriedade associada ao produto vetorial entre dois vetores m e n definida em Argyris [1982] e em Felippa [2001]:

m × n = Sm . n Sm =

−−

−

00

0

12

13

23

mmmm

mm (4.8)

sendo m1, m2 e m3 os componentes do vetor m, obtém-se a seguinte expressão: rn ( ) 000 reerer ××−+×+= .cos1.sen θθ (4.9) rn ( ) 0ee0e0 rSSrSr ...cos1..sen θθ −++= (4.10) rn ( )( )2.cos1.sen ee SSI θθ −++= .r0 (4.11) que ao ser comparada com a Eq. (4.1), fornece a seguinte equação para definir a matriz de rotação no espaço, conhecida na literatura como fórmula de Rodrigues: ( ) 2.cos1.sen eeθ SSIR θθ −++= (4.12)

33

Conforme comentado anteriormente, ao contrário do caso plano, as rotações no espaço não podem ser tratadas vetorialmente. De modo a mostrar o caráter não vetorial das rotações, será aplicado uma seqüência de três rotações sucessivas de 90o ao corpo rígido representado pelo retângulo A0 contido inicialmente no plano xy da Fig. (4.2). Na figura do lado esquerdo é aplicado inicialmente uma rotação em torno do eixo x, levando o corpo da posição A0 para a posição A1. Em seguida, aplica-se uma rotação em torno do eixo y, levando o corpo da posição A1 para A2 e por último, aplica-se uma rotação em torno do eixo z que leva o corpo da posição A2 para a posição final A3. Já na figura apresentada do lado direito, são também aplicadas três rotações sucessivas de 90o ao corpo rígido representado pelo retângulo A0 contido inicialmente no plano xy, porém as rotações foram aplicadas em ordem inversa, ou seja, inicialmente aplica-se uma rotação em torno do eixo z, em seguida em torno do eixo y e por último uma rotação em torno do eixo x. Pode-se ver claramente através da Fig. (4.2) que a alteração da ordem na qual foram aplicadas as rotações levou a resultados completamente diferentes, elucidando o caráter não vetorial das rotações no espaço. No entanto, vale lembrar que quando efetuadas em torno de um único eixo, as rotações podem ser tratadas vetorialmente como ocorre nos casos bidimensionais.

A0 A0

A1 A1

A2

A2

A3

A3

x x

y y

z z

Figura 4.2 – Caráter não vetorial das rotações no espaço

Vale ressaltar que por definição, o vetor r0 não sofre alterações no seu comprimento com uma rotação de corpo rígido e portanto pode-se dizer que: r0

T.r0 = rnT.rn

r0

T.r0 =(Rθ.r0)T(Rθ.r0) r0

T.r0 =r0T.Rθ

T. Rθ.r0 r0

T.(I – RθT. Rθ ). r0 = 0 (4.13)

cuja validade para um vetor r0 qualquer implica em:

RθT. Rθ = I (4.14)

34

Consequentemente, pode-se afirmar que a matriz de rotação é ortogonal e portanto, segundo Cole [1990], Rθ em um ponto pode ser caracterizada por apenas três componentes (θ1, θ2 e θ3) em vez de nove. Esta representação da matriz de rotação na forma de um vetor θ, geralmente conhecido na literatura por pseudovetor de rotação, conforme apresentada abaixo:

θ e.

3

2

1

θθθθ

=

= 2

32

22

1 θθθθ ++= (4.15)

permite re-escrever a Eq. (4.12) em função das componentes de rotação:

( )( )

22

2

2/2/sen.

21sen

θθθ SSIRθ

θθ

θ++= (4.16)

O problema inverso, ou seja, a obtenção das componentes de rotação θ a partir de Rθ, simbolizado pela seguinte operação: θ = rot (Rθ), pode ser obtido facilmente em função da parte anti-simétrica de Rθ expressa por:

e

Tθθ SRR .sen

2θ=

− (4.17)

que pode ser re-escrita expandindo-se as matrizes como

−−

−=

−−−−−−

00

0sen

00

0

21

12

13

23

23321331

32231221

31132112

eeee

ee

RRRRRRRRRRRR

θ (4.18)

sendo Rij os componentes de Rθ. Em função da equação acima pode-se dizer que:

re =

−−−

=

=

1221

3113

2332

3

2

1

21sen.sen

RRRRRR

eee

θθ (4.19)

portanto, lembrando que θ = θ.e, sendo e um vetor unitário e θ ≥ 0, tem-se que:

0θ = se 0=r (4.20)

( )rr.rθ 1sen −= se 0≠r

Observa-se que os valores de θ estão restritos ao intervalo (0,π/2). Entretanto, pode-se

facilmente ampliar o intervalo de amplitudes para (0,π) empregando-se um algoritmo muito conhecido atribuído a Spurrier, descrito por vários autores tais como: Monteiro [2004], Crisfield [1997], Cole [1990], sendo apresentado a seguir:

35

tr(Rθ) = R11 + R22 + R33 m = max[tr(Rθ), R11, R22, R33] n = )(21 θRtrm −+ se m = tr(Rθ)

q = 2n q =

n21 [R32 – R23 R13 – R31 R21 – R12]T

se m = R11

q = ( )233221 RRn

− q = n2

1 [n2 R21 + R12 R31 + R13]T

(4.21) se m = R22

q = ( )311321 RRn

− q = n2

1 [R12 + R21 n2 R32 + R23]T

se m = R33

q = ( )122121 RRn

− q = n2

1 [R13 + R31 R23 + R32 n2]T

O algoritmo descrito acima obtém, a partir da matriz de rotação, o escalar q e o vetor q empregados no cálculo do vetor w, definido por Crisfield [1997] como

qqew 2

2tan.2 =

=θ (4.22)

que define o pseudovetor de rotação através das seguintes expressões: 0θ = se 0=w (4.23a)

ww.

wθ

= −

2tan.2 1 se 0≠w (4.23b)

Vale ressaltar que, as quantidades q e q são as principais componentes da álgebra dos quatérnios, comumente utilizada no estudo de rotações finitas no espaço por um grande número de autores com o intuito de reduzir o esforço computacional. 4.2 MATRIZ DE ROTAÇÃO PARA PEQUENAS ROTAÇÕES Admitindo-se que o corpo rígido apresentado na Fig. (4.1) gire de um pequeno ângulo ∆θ, tal que (∆θ < 5o), em torno do eixo de rotação, tem-se que o comprimento do arco descrito pela trajetória de Q entre as posições Q’ e Q’’, pode ser definido de maneira aproximada pela corda que une os pontos Q’ e Q’’, representada por drot: R.∆θ ≅ rotd (4.24) ao passo que o deslocamento drot pode ser considerado perpendicular ao plano formado pelo eixo de rotação e pelo vetor r0:

36

0red

rere

dd rot

0

0rotrot ×=

××

=R

(4.25)

substituindo-se a equação (4.24) em (4.25), obtém-se: 00rot r∆θred ×=×∆= .θ (4.26)

A equação acima pode então ser substituída na Eq. (4.3) para se obter a matriz de rotação no caso de pequenas rotações:

rn = r0 + ∆θ × r0

rn = r0 + S∆θ . r0

rn = ( I + S∆θ ).r0 = R∆θ . r0 ⇒ R∆θ = ( I + S∆θ ) (4.27)

que é a linearização da fórmula de Rodrigues. Em decorrência da linearização, pequenas rotações podem ser tratadas vetorialmente. Se o corpo sofre uma rotação ∆θ1, seguida de uma rotação ∆θ2, tem-se: rn = R∆θ2.R∆θ1 .r0 = R∆θ12 .r0 (4.28) sendo a matriz de rotação resultante R∆θ12 = ( I + S∆θ2 ). ( I + S∆θ1 ) = I + S∆θ1 + S∆θ2 + S∆θ2 .S∆θ1 ≅ I + S∆θ1 + S∆θ2 (4.29) onde a parcela S∆θ2.S∆θ1 pode ser desprezada por se tratar de um infinitésimo de ordem superior. Nessas condições, R∆θ12 = R∆θ21 e portanto, ∆θ1 + ∆θ2 = ∆θ2 + ∆θ1.

37

CAPÍTULO 5 FORMULAÇÃO CO-ROTACIONAL / PÓRTICO ESPACIAL

5.1 DESCRIÇÃO CINEMÁTICA Uma vez conhecida a matriz de rotação que descreve o movimento de um corpo rígido no espaço conforme definido anteriormente, no presente capítulo será descrita a formulação co-rotacional referente a pórticos espaciais baseando-se em especial nos trabalhos apresentados por Monteiro [2004] e Nour-Omid & Rankin [1991]. O movimento de um elemento de pórtico espacial desde a sua configuração inicial C0 até a configuração atual C é mostrado na Fig. (5.1), sendo utilizados três sistemas distintos de eixos cartesianos ortogonais da mesma forma que a empregada por Rankin & Brogan [1986]:

• Sistema global utilizado para definir a conectividade entre os elementos da estrutura. • Sistema local E que sofre translações e giros ao acompanhar o elemento, sendo em

geral conhecido na literatura como sistema co-rotacional. • Sistemas nodais A e B que estão atrelados respectivamente aos nós inicial (1) e final

(2) de cada elemento.

C0

x

y

z

a1

a2 a3

b1

b2b3

A

A0

BB0

x0y0

z0e10

e20 e30

u1

u2

Rθ1

Rθ2

1

1

2 2

Figura 5.1 – Elemento de pórtico no espaço

O sistema local E, definido na configuração inicial C0 como E0, é dado por: E0 = [ e10 e20 e30 ] = Rθ0 (5.1)

e10 = 21

21

XX e30 =

veve

10

10

××

e20 = e30 × e10 (5.2)

sendo X21 = X2 – X1, tal que o vetor Xi representa a posição inicial do nó i na configuração C0 e v é um vetor contido no plano local x0y0 do elemento. No presente trabalho, a origem do sistema local é o nó 1 e o posicionamento de seus eixos locais coincide com as direções principais de inércia do elemento.

38

A orientação inicial dos sistemas nodais A e B na configuração C0, denominadas A0 e B0 respectivamente, é escolhida igual a dos eixos locais:

A0 = B0 = E0 (5.3)

O movimento do elemento de pórtico entre as configurações inicial e atual pode ser subdividido em duas etapas. Observando-se a Fig. (5.1) novamente, inicialmente os nós 1 e 2 sofrem translações representadas por u1 e u2, resultando em deformações axiais do elemento e em seguida, os sistemas nodais A e B sofrem rotações de θ1 e θ2 com relação a referência inicial E0, gerando curvaturas no elemento. A posição dos nós em C pode então ser dada por: xi = Xi + ui (5.4) ao passo que a orientação dos eixos nodais pode ser obtida por: A = [ a1 a2 a3 ] = Rθ1 . A0 B = [ b1 b2 b3 ] = Rθ2 . B0 (5.5) O sistema co-rotacional E na configuração atual C, por sua vez pode ser definido da seguinte forma, segundo Rankin & Nour-Omid [1988]: E = [ e1 e2 e3 ] = RθE (5.6)

e1 = 21

21

xx e3 =

21

21

aeae

×× e2 = e3 × e1 (5.7)

observando-se que o eixo e1 aponta no sentido da corda que une os nós 1 e 2.

Uma vez definido o sistema de eixos co-rotacionais, estabelece-se a mudança de coordenadas entre o sistema global xyz e o local xeyeze através das seguintes operações descritas por Monteiro [2004], sendo x e X quantidades tensoriais de primeira e segunda ordem, respectivamente:

xe = ET. x Xe = ET.X.E (5.8)

O próximo passo é a determinação dos deslocamentos e rotações deformacionais medidos no sistema local, também denominados co-rotacionais, obtidos a partir dos deslocamentos generalizados globais do elemento. Inicialmente, será calculada a posição do i-ésimo nó do elemento genérico apresentado na Fig. (5.2) em sua configuração atual C, designada pelo vetor xi

e. xi

e = ET. ( xi – x1) = ET. xi1 (5.9) O deslocamento co-rotacional ui

c pode então ser obtido, escrevendo-se localmente a relação definida na Eq. (5.4) e substituindo-se o valor de xi

e definido na equação acima: ui

c = xie – Xi

e ui

c = ET. xi1 – Xie

ui

c = ET. (Xi1 + ui1) – Xie (5.10)

39

sendo a posição inicial local do nó i em C0, definida através da seguinte equação: Xi

e = E0T. (Xi – X1) = E0

T. Xi1 (5.11)

C0

C

Ri

Rdi

E0

E

ui

u1

Xi X1

xi x11

1

i

i

Figura 5.2 – Translações e rotações de um elemento genérico de pórtico

Lembrando que o deslocamento corrente de rotação θi compõe-se de uma parcela rígida θr seguida de uma parcela deformacional θd,i e que a rotação de corpo rígido do elemento é por definição a rotação entre os eixos E0 e E, então: Rθi = Rθd,i . Rθr (5.12) Rθr = E . E0

T (5.13) substituindo-se a Eq.(5.13) em (5.12), obtém-se: Rθd,i = Rθi . Rθr

T = Rθi .( E . E0T)T

Rθd,i = Rθi . E0 . E T (5.14) A rotação co-rotacional θi

c pode então ser obtida efetuando-se a mudança de coordenadas proposta na Eq. (5.8): Rθi

c = ET. Rθi . E0 θic = rot(Rθi

c) (5.15) A mudança de variáveis, inerente à formulação co-rotacional fica portanto definida pelas equações (5.10) e (5.15). Os deslocamentos co-rotacionais nodais d1

c e d2c do elemento de

pórtico podem então ser definidos através das seguintes expressões:

d1c =

c1

c1

θu

=

000

c1u θ1

c = rot(ET.A) (5.16)

40

d2c =

c2

c2

θu

−=

00

0

2

LLcu θ2

c = rot(ET.B) (5.17)

sendo os comprimentos inicial e final do elemento definidos por: 21X=0L 21x=L (5.18) 5.2 EQUAÇÕES DE EQUILÍBRIO Uma vez atingida uma configuração de equilíbrio, o trabalho virtual externo (δWe), realizado pelas forças nodais generalizadas externas Fext sobre os deslocamentos virtuais nodais generalizados de toda a estrutura δD, será igual ao trabalho virtual interno (δWi), proveniente do somatório do trabalho virtual realizado pela distribuição de tensão que atua em cada um dos elementos, de modo que: δWe + δWi = 0 (5.19) δWe = δDT. Fext (5.20) δWi = ( )∑ ∫−

Vo

dV0.σε Tδ (5.21)

Em virtude do princípio da invariância da energia interna para o caso de movimentos de corpo rígido, o trabalho realizado pelas forças internas durante deslocamentos virtuais rígidos será nulo e consequentemente, este trabalho é proveniente somente de δdd que representa a parcela deformacional de δd. Além disso, por se tratar de uma quantidade escalar, o trabalho virtual interno proveniente dos deslocamentos deformacionais medidos no sistema global e local devem ser iguais e portanto: δWi (δd) = δWi (δdd) = δWi (δdc) (5.22) sendo δdc os deslocamentos virtuais co-rotacionais generalizados do elemento e portanto representa a parcela de deslocamentos deformacionais δdd escrito no sistema local. Vale ressaltar que é muito mais vantajoso definir as tensões e deformações conjugadas no sistema local, uma vez que devido ao fato do campo de deslocamentos co-rotacionais ser exclusivamente deformacional, pode-se linearizar a relação deslocamento-deformação desde que localmente as deformações sejam pequenas. Esta hipótese torna-se mais realista a medida que a malha é refinada.

Supondo uma relação deformação-deslocamento local linear, é possível escrever a seguinte expressão para o campo de deformações ε em função dos deslocamentos co-rotacionais dc do elemento: ε = Bd

c. dc ε = ε(xe,ye,ze) (5.23) lembrando que neste caso Bd

c independe de dc. Substituindo-se:

41

δε = ccd

cc dBd

dε δδ .. =

∂∂ (5.24)

na Eq. (5.21) obtém-se a seguinte expressão para o trabalho virtual interno:

δWi = ( ) ( ) ( ) ccc

dc fdσBd .... 0

TTT

Vo

dV ∑∑ ∫ −=− δδ (5.25)

sendo fc as forças co-rotacionais generalizadas do elemento que realizam trabalho com os deslocamentos virtuais δdc, tal que: ( )∫=

Vo

TdV0..σBf c

dc σ = σ(xe,ye,ze) (5.26)

Admitindo-se as seguintes expressões que relacionam os deslocamentos virtuais (δdc), (δde) e (δd), definidas por Monteiro [2004] e Nour-Omid & Rankin [1991]: δdc = P. δde δde = GT. δd (5.27) então:

δdc = P. GT. δd (5.28) e substituindo-se a equação acima na Eq. (5.25), obtém-se:

δWi = FDfPGd T .... TcT δδ −=− ∑ (5.29) sendo F o vetor de forças nodais generalizadas de toda a estrutura obtido a partir da contribuição f de cada elemento: F = ∑ ∑ ∑== ffGfPG ecT ... (5.30) É importante ressaltar que a transformação P que aparece na Eq. (5.27) atua como um filtro para os deslocamentos deformacionais, uma vez que ao ser aplicada sobre o vetor de deslocamentos virtuais nodais δde extrai-lhe a sua parcela deformacional δdc. Portanto, se P elimina a parcela de deslocamento rígido de δde, então ao ser aplicada sobre a sua parcela deformacional, resulta nela mesma: δdc = P. δdc (5.31)

δdc = P. P. δde = P2. δde (5.32)

e consequentemente:

P = P2 (5.33)

caracterizando-o como um operador de projeção segundo Rankin & Nour-Omid [1988] e Rao & Mitra [1971].

42

Substituindo-se as equações (5.29) e (5.20) em (5.19) e considerando a arbitrariedade de δD, resulta no seguinte sistema não linear, cuja solução pode ser obtida pelo processo incremental-iterativo de Newton-Raphson: Ψ = Fext – F = 0 (5.34) Vale enfatizar o fato de que no método de Newton-Raphson, os deslocamentos nodais generalizados são atualizados da seguinte forma: Dk+1 = Dk + δDk (5.35) sendo invariavelmente inconsistente para os deslocamentos de rotação, representando apenas uma estimativa arbitrária para a posição de equilíbrio e portanto desprovida de sentido físico. Uma nova estimativa para o deslocamento translacional ui é obtida a partir de seu valor corrente adicionando-lhe o deslocamento iterativo δui: ui

k+1 = ui k + δui

k (5.36a) ao passo que a atualização para a rotação nodal, não podendo ser feita da mesma maneira, é obtida atualizando-se a matriz de rotação associada. Dada a estimativa inconsistente natural do método: θi

k+1 = θi k + δθi

k (5.36b) a matriz de rotação associada ao deslocamento θi

é atualizada consistentemente na forma: Rθi

k+1 = Rθi k + δRθi

k + ... (5.37) sendo a variação δRθ calculada diferenciando a Eq. (4.12) em relação a θ e considerando a identidade definida na Eq. (5.38) que pode ser facilmente obtida após manipulação algébrica: Se. Rθ = Rθ. Se = cosθ. Se + senθ. Se

2 (5.38) δRθ = Se. Rθ. δθ + senθ. δSe + (1– cosθ).(δSe. Se + Se. δSe) (5.39) 5.3 MUDANÇA DA VARIÁVEL ITERATIVA DE ROTAÇÃO Nesta seção será efetuada uma mudança na variável iterativa de rotação de acordo com as formulações co-rotacionais apresentadas por Monteiro [2004] e Nour-Omid & Rankin [1991], segundo a qual é importante diferenciar o deslocamento virtual de rotação δθi da rotação do sistema nodal δ

iθ , uma vez que esta última não está relacionada com a variação do campo de deslocamento de rotação do elemento e sim com o giro instantâneo da tríade nodal. Portanto, será empregada a rotação do eixo nodal δ iθ como variável iterativa de rotação no lugar do deslocamento virtual δθi. Com esta mudança, a variável de rotação passa a ser θ , estando associada ao mesmo estado de rotação de θ: θRR =θ (5.40)

43

cuja matriz de rotação pode ser atualizada na forma:

ki

ki

ki θθδθ RRR .1 =+ (5.41)

porém, é importante salientar que δ ≠θ δθ, conforme será mostrado posteriormente. O sistema não linear definido na Eq. (5.34) pode então ser re-escrito através da seguinte equação, onde a barra sobre as quantidades indica dependência da variável θ : 0FFψ ext =−= (5.42) Conforme comentado por Monteiro [2004], é importante mencionar que no sistema local a matriz de rigidez constitutiva ck é não simétrica, uma vez que a energia potencial do elemento é definida em função dos deslocamentos co-rotacionais dc, ou seja, os graus de liberdade do elemento na formulação elástico-linear, e não em termos dos deslocamentos

cd que dependem das variáveis associadas a rotação da tríade nodal. Considerando a Eq. (5.40), substituindo-se a matriz de rotação linearizada apresentada na Eq. (4.27) em (5.41) e lembrando que δ ≠θ δθ, então: ( ) k

ik

ik

ik

ik

ik

ik

ik

i θθδθθθδθθδθ RSRRSIRRR ...1 +=+==+ (5.43) que ao ser comparada com a Eq. (5.37) fornece: T

θθθδ δ RRS .= (5.44) A expressão definida na equação acima estabelece a relação existente entre a rotação δ θ dos eixos cartesianos associados a matriz Rθ com a variação consistente δRθ no sistema global de coordenadas. O próximo passo é a determinação da relação entre as rotações δθ e δ θ obtida aplicando-se as equações (5.39) e (4.12) em (5.44), que após algebrismos que foram apresentados por Monteiro [2004] resulta na seguinte expressão:

( )θθθθ23 ...cos1sen...sen SSSSSSθθS δδθ

θδθ

θδθ

θθθθθδ −

−++

−= T (5.45)

a partir da qual, pode-se obter: δ θ = Λθ

-1. δθ (5.46) sendo:

Λθ -1 = TθθSI

θθ

3θ2

sencos1senθ

θθθ

θθ

θ −+

−+=

∂∂ (5.47)

A relação inversa é expressa por: δθ = Λθ . δ θ (5.48) onde:

44

2θ .

21

θθ SSIθθΛ ξ+−=

∂∂

= (5.49)

( )θθ

θθθξsen2

cos1sen22

+−= (5.50)

Uma vez feita a atualização dos deslocamentos nodais ui através da Eq. (5.36) e da matriz de rotação através da Eq. (5.41), pode-se obter os deslocamentos co-rotacionais generalizados do elemento utilizando-se as expressões (5.10) e (5.15). Conhecido o vetor de deslocamentos co-rotacionais dc e a matriz de rigidez constitutiva clássica do elemento kc, o vetor de foças co-rotacionais f c pode ser calculado na forma:

f c = kc. dc f c =

c

c

2

1

ff

= ci

cic

i mnf (5.51)

Porém, como foi feita uma mudança da variável de rotação, torna-se necessário determinar as quantidades cf e ck do elemento, sendo as forças cf obtidas diretamente a partir de f c aplicando-se a regra da cadeia sobre a energia potencial do elemento φ :

[ ] cTcT

c

T

c

cT

cc fH

ddd

df .. =

∂∂

∂∂

=

∂∂

=φφ (5.52)

sendo Hc uma matriz quadrada 12x12 definida por:

=

∂∂

= cc

cc

c

cc

2221

1211

HHHH

ddH (5.53)

e Hij

c sub-matrizes quadradas 6x6 que, considerando a Eq. (5.49), podem ser expressas por:

=

∂

∂=ci

ijc

i

ci

ij

ijc

ijθ

δδ

δ

Λ00I

θθ

0

0IH ..

. (5.54)

onde δij é o delta de Kronecker e 0 são matrizes quadradas 3x3 nulas. Já a matriz de rigidez

ck do elemento pode ser obtida da seguinte forma:

ck = ( )[ ]cTccc

c

fHdd

f .∂

∂=

∂∂ (5.55)

= ( ) ( )[ ] cccTccc

cTc

21.. kkfHdd

fH +=∂

∂+

∂∂ (5.56)

A parcela simétrica c

1k é expressa por:

45

( ) ( ) ccTcc

c

c

cTcc HkH

dd

dfHk ...1 =

∂∂

∂∂

= (5.57)

enquanto a parcela não simétrica c

2k pode ser expressa pela seguinte equação:

( )( )

∂∂

∂∂

= c

ccTc

cc

ddfH

dk ..2 = ( )

( )c

cTc

cTc

c HfHfH

d.

.

.

222

111

∂∂ =

( )( )

( )( )

∂∂

∂∂

∂∂

∂∂

cTc

c

c

cTc

c

c

c

c

2

2

1

1

.

.

2

1

mΛd

dn

mΛd

dn

θ

θ (5.58)

lembrando que o vetor f c entra na derivada da equação acima como uma constante, tem-se:

[ ] 12321

xc

c

c

c

0dn

dn

=

∂∂

=

∂∂

(5.59)

ao passo que as demais sub-matrizes podem ser representadas após algebrismos simples por:

( )( ) [ ]00Ω0mΛd

ccTc c 11.

1=

∂∂

θ ( )( )cT

cc

c 11

1 .1

mΛθ

Ωθ∂

∂= (5.60)

( )( ) [ ]ccTc c 22.

2Ω000mΛ

d=

∂∂

θ ( )( )cT

cc

c 22

2 .2

mΛθ

Ωθ∂

∂= (5.61)

sendo 0 matrizes quadradas 3x3 nulas, ao passo que as derivadas que aparecem nas equações (5.60) e (5.61) podem ser expressas pelas seguinte equação, conforme demostrado em detalhes por Monteiro [2004] e Nour-Omid & Rankin [1991]:

( )( ) ( ) TTTTm

T θmSθmmθImθSmΛθ

.....2...21 2

θθ µξ +−++−=∂∂ (5.62)

( ) ( ))2/sen(4

2/sen8sen.12

2

θθθθθθ

θξ

θµ −+

==dd (5.63)

Finalmente, substituindo-se as equações (5.59), (5.60) e (5.61) na Eq. (5.58), pode-se obter:

ccc HΩk .2 =

12122

1

xc

cc

=

Ω0

Ω0

Ω (5.64)

46

cuja influência segundo Nour-Omid & Rankin [1991] é desprezível quando as rotações co-rotacionais são moderadas ou quando a malha é refinada, para as quais a matriz ck do elemento é simétrica para efeito prático. 5.4 OPERADOR DE PROJEÇÃO Nesta seção será apresentada a determinação do operador de projeção P do elemento, definido em função da mudança da variável iterativa de rotação proposta na seção anterior. Partindo da Eq. (5.27) apresentada anteriormente, tem-se:

12122221

1211

xe

c

=

∂∂

=PPPP

ddP

66x

ej

ci

ej

ci

ej

ci

ej

ci

ej

ci

ij

∂

∂

∂

∂

∂

∂

∂

∂

=

∂

∂=

θθ

uθ

θu

uu

dd

P (5.65)

Inicialmente será feita a determinação da variação do deslocamento co-rotacional ui

c, definido na Eq. (5.10): 11 .. i

Ti

Tci uExEu δδδ += (5.66)

Aplicando a propriedade (5.44) à tríade local do elemento e efetuando-se a mudança do sistema de coordenadas apresentado em (5.8), obtém-se: 11 ... i

Ti

Tci E

uExSEu δδ θδ +−= 11 .... i

Ti

TTE

uExEESE δθδ +−= ( ) ( )ee

iee

ieE

11. uuxxS δδθδ

++−−= (5.67)

lembrando que o nó inicial 1 foi definido como a origem do sistema local de coordenadas, tem-se que a sua posição e seus deslocamentos locais são nulos e portanto: e

ie

Ex

ei

ei

ci e

ie

EuθSuxSu δδδδ

θδ+=+−= .. (5.68)

Uma vez conhecida a variação do deslocamento co-rotacional segundo a equação acima, o próximo passo é a determinação da variação da matriz de rotação associada ao deslocamento co-rotacional c

iθ , definida pela Eq. (5.15), de modo que: ( ) 0.. ERERER i

Ti

Tic θθθ

δδδ += (5.69)

Procedendo de forma análoga ao caso do deslocamento co-rotacional e aplicando-se a propriedade (5.44) à tríade local do elemento e efetuando-se a mudança do sistema de coordenadas apresentado em (5.8), obtém-se:

47

( ) 0..... ERSERSER iiT

iT

i Ec θθδθθδθδ +−=

( ) 0....... EREESEEESE i

Ti

TTTE θθδθδ +−=

( ) 0...... EREESEESE i

Ti

TTE θθδθδ +−= ( ) ce

Ee ii θθδθδ

RS .−

= (5.70) isolando-se o termo entre parênteses na equação acima: ( ) T

iii cceE

e θθθδθδδ RRS .=

− (5.71)

que resulta, considerando a Eq. (5.44), na seguinte expressão: e

Ee

ic

i θθθ δδδ −= (5.72) Em função das equações (5.68) e (5.72) as derivadas dos deslocamentos co-rotacionais

ciu e c

iθ com relação aos deslocamentos eiu e e

iθ podem ser definidas por:

IuθS

uu

uθS

uu

... ijej

eE

xej

ei

ej

eE

xej

ci

ei

ei

δ+

∂∂

=

∂

∂+

∂∂

=

∂

∂

∂∂

=

∂

∂+

∂∂

=

∂

∂e

j

eE

xej

ei

ej

eE

xej

ci

ei

ei θ

θSθu

θθS

θu

..

(5.73)

∂∂

−=

∂∂

−

∂

∂=

∂

∂e

j

eE

ej

eE

ej

ei

ej

ci

uθ

uθ

uθ

uθ

∂∂

−=

∂∂

−

∂

∂=

∂

∂e

j

eE

ijej

eE

ej

ei

ej

ci

θθI

θθ

θθ

θθ

.δ

Substituindo-se as derivadas definidas na Eq. (5.73) em (5.65), pode-se obter a seguinte

expressão para a sub-matriz ijP :

66

..

x

ej

eE

ej

eE

ej

eE

xej

eE

x

ijij

ei

ei

∂∂

−

∂∂

−

∂∂

∂∂

+

=

θθ

uθ

θθS

uθS

II

P δ (5.74)

ou de forma compacta T

jiijij ΓΨIP .−= δ (5.75)

48

sendo:

[ ]Tx

xi e

i

ei IS

IS

=

−=Ψ

T

ej

eE

T

ej

eE

ej

eE

j

∂∂

=

∂∂

∂∂

=dθ

θθ

uθΓ (5.76)

Consequentemente, o operador de projeção do elemento pode ser expresso por:

12122212

2111

1212 ....

xTT

TT

x

−

=

ΓΨΓΨΓΨΓΨ

II

P (5.77)

ou de forma compacta =P I – Ψ Γ T (5.78) sendo

[ ]Txx

xee

iISIS

ΨΨ

Ψ2

3122

1 =

=

T

LL

−

=

100010001

0000

000

100010001

000000000

(5.79)

T

e

eE

T

e

eE

e

eE

x

∂∂

=

∂∂

∂∂

=

=

dθ

dθ

dθ

ΓΓ

Γ213122

1 (5.80)

O próximo passo será a obtenção da matriz Γ definida acima que depende da geometria atual e da orientação dos eixos locais do elemento. Inicialmente, a propriedade (5.44) é aplicada sobre a tríade local do elemento E e em seguida é efetuada a mudança de coordenadas apresentada na Eq. (5.8), tal que: T

EEES .δθδ =

EEEEESE ..... TTT

Eδθδ =

EES δ

θδ.T

eE= (5.81)

que ao ser expandida fornece no sistema local de coordenadas:

( )( )( )

[ ][ ][ ]

−−

=

−

−

=e

e

e

eTe

eTe

eTe

eE

1

1

3

12

13

32

.010.100.010

...

eee

eeeeee

θδδδ

δδδ

δ (5.82)

49

À seguir serão então determinadas as variações δe1e e δe3

e que aparecem na equação anterior. A primeira variação pode ser obtida a partir das equações (5.7) e (5.18):

2121 xxe .1121 L

LLδδδ −= = ( ) 2111 ueeI δ..1 T

L− (5.83)

que localmente pode se escrito por:

( ) eeTeee

LL 2121111 uueeIe δδδ .100010000

1..1

=−= (5.84)

A variação δe3

e pode também ser obtida a partir da Eq. (5.7), no sistema local de coordenadas, da seguinte forma:

( ) ( ) eeeee

eeee

e31

1

1

1

3 eaeae

aeae

e ..1.12

2

2

2

××

−××

= δδδ (5.85)

Definindo o vetor a2

e como:

=

3

2

1

qqq

e2a (5.86)

e lembrando que, conforme definido na Eq. (5.7), os vetores e3 e a2 são ortogonais então q3 = 0 (5.87)

=

×

=×

2

2

1

00

0001

qqq

ee21 ae (5.88)

Logo,

( ) eeeeee

qq

q 3113 eaeaee ..1

2

222

2

δδδδ −×+×= (5.89)

A partir das equações (5.3) e (5.5) tem-se:

=

010

..2 01 ERa θ (5.90)

cuja variação, considerando-se a Eq. (5.44), pode ser expressa por:

50

2101012 aSERSERa .010

...010

.. 1 θδθθδθδδ =

=

= (5.91)

que escrita localmente: eee

ee 1a22 θSaSa21δδ

θδ.. −== (5.92)

Definindo

=

3

2

1

1

ααα

eθ (5.93)

e calculando e

2aδ a partir da Eq. (5.92) pode obter: 312 .δαδ qq = (5.94) Substituindo-se as equações (5.84), (5.91) e (5.94) em (5.89) e definindo η = q1 /q2, tem-se a seguinte expressão para a variação e

3eδ :

( ) eeeeee

qq

q 321123 eaeeae2

31

2

.1 δαδδδ −×+×−= (5.95)

eee

L 321 eθ

u..

0001000

.00

00100

31 δαηδη

ηδ

ηη −

−+

−

−= (5.96)

Uma vez conhecidas as variações de e

1eδ e e3eδ definidas anteriormente pelas equações

(5.84) e (5.96) respectivamente, pode-se então calcular as componentes de eEθδ em (5.82):

ee

ee

ee

ee

eE L 1

3

2

1

00000001

010100

00θuθ 21 δ

ηδ

η

θδθδθδ

δ

−+

−−

=

= (5.97)

Em função das equações (5.80) e (5.97), finalmente obtém-se a matriz Γ na forma:

TT

e

eE

LL

L

−−

−

−=

∂∂

=000000000

010100

00

0000000

010100

001

ηηη

dθ

Γ (5.98)

51

que pode ser decomposta na seguinte maneira: ΞΥΥΓ 21 .+= (5.99) sendo

T

=

000000000

000000000

000000001

000000000

1Υ (5.100)

T

−

−

−=

000000000

010100

000

000000010

010100000

2Υ (5.101)

TL

L

=

1000100

1ηη

Ξ (5.102)

5.5 VETOR DE FORÇAS INTERNAS E MATRIZ DE RIGIDEZ TANGENTE O vetor de forças internas e a matriz de rigidez tangente da estrutura após a mudança da variável iterativa de rotação podem ser expressos pelas seguintes equações: ∑∑∑ === ffGfPGF ecT ... (5.103)

∑∑ =

∂∂

=

∂∂

= tt kdf

DFK (5.104)

sendo o vetor de forças internas do elemento no sistema global e a matriz G definidos por:

cT fPGf ..= ( )E

EE

EE

G 4

1212

diag

x

=

= (5.105)

ao passo que a matriz de rigidez tangente pode ser obtida pela variação das forças internas f do elemento:

cT fPGf δδ ..1 = ( ) 321.. ffffPGf δδδδδ ++== cT cT fPGf ..2 δδ = (5.106)

cT fPGf ..3 δδ =

52

Inicialmente, será calculada a variação das forças co-rotacionais cfδ , que em função das equações (5.27) e (5.55), pode ser expressa por:

dGPkddd

dd

dfd

dff δδδδ ...... Tc

e

e

c

c

cc

c

cc =

∂∂

∂∂

∂∂

=

∂∂

= (5.107)

Portanto

( ) dkdGPkPGf δδδ ...... 11 t

TcT == (5.108) sendo

Tett GkG

df

k .. 11

1 =

∂∂

= PkPk ..1cTe

t = (5.109)

O próximo passo é a determinação da variação do operador de projeção Pδ que pode ser obtido calculando-se a variação do operador P definido na Eq. (5.78) e lembrando que as componentes 1Υ e 2Υ definidas na equação (5.99) são constantes, tal que: TT

2.... ΥΞΨΓΨP T δδδ −−= (5.110) Observando-se que as matrizes Ψ e Γ são bi-ortonormais, ou seja, Γ T.Ψ = I , e considerando a Eq. (5.99), a condição de bi-ortonormalidade pode ser escrita na forma: IΞΥΨΥΨ =+ ... 21

TT (5.111) Diferenciando-se a equação acima e sabendo que ( 2.ΥΨT ) admite inversa: ( ) ΓΨΥΨΞ ... -1

2TT δδ −= (5.112)

Assim, a variação do operador de projeção pode ser escrito como

( )( ) TTTTT

21

2 ....... ΥΥΨΨΓΨΓΨP −+−= δδδ

( )( ) TTTTT2

1

2 ....... ΥΥΨΨΓΨΓΨ−

+−= δδ 1...... −+−= ΨΨΓΨΓΨ δδ TT (5.113) que considerando-se a condição de bi-ortonormalidade e a Eq. (5.78) ( )TTTT ΓΓΨΨΓΨΓΨP ....... 1 −−+−= δδδ ( ) TT ΓΨΓΨI .... δ−−= TΓΨP ...δ−= (5.114)

53

Portanto ( ) cTT fΓΨPGf ....2 δδ −= cTT fPΨΓG ..... δ−= eT fΨΓG .... δ−= (5.115) e levando-se em conta a derivada das equações (5.79) e (5.10) ( )ee

ee 2nn2 ...21

xSxSΓGf 1 δδδ +=

( )cc

ee 2nn...

21uSuSΓG 1 δδ += (5.116)

podendo ser compactado na seguinte forma

cT dFΓGf δδ .~..2 −=

312

2

1

~

x

n

n

e

e

=

0S

0S

F (5.117)

Portanto, cT dFΓGf δδ .~..2 −=

ddd

ddFΓG δ...~..

∂∂

∂∂

−=e

e

cT

dkdGPFΓG δδ ....~.. 2t

TT =−= (5.118) sendo

Tett GkG

dfk .. 2

22 =

∂∂

= PFΓk .~.2Te

t = (5.119)

Por último, será calculada a variação δG no sistema local de coordenadas. Considerando a Eq. (5.81) e a matriz G definida na Eq. (5.105):

( )eE

diag

x

θδ

δδ

δδ

δ SG

EE

EE

G 4

1212

.=

= (5.120)

54

que fornece para 3fδ ecT fGfPGf 3 ... δδδ == ( ) ( ) e

Ee

eeE

diagdiag θSGfSGf

δθδ

.... 44 −== (5.121)

podendo ser re-escrita da seguinte forma:

eEθHGf δδ ..3 −=

3122

2

1

1

xm

n

m

n

e

e

e

e

=

SSSS

H (5.122)

e considerando a identidade

dGΓddd

dθ

θ δδδ ... TTe

e

eEe

E =

∂∂

∂∂

= (5.123)

escreve-se na forma dkdGΓHGf δδδ ..... 33 t

TT =−= (5.124) sendo:

Tett GkG

df

k .. 33

3 −=

∂∂

= Tet ΓHk ..3 −= (5.125)

Finalmente, agrupando-se as equações (5.109), (5.119) e (5.125) obtém-se a matriz de rigidez tangente do elemento de pórtico espacial no sistema global de eixos Te

ttttt GkGkkkk ..321 =++= (5.126) sendo TTcTe

te

te

te

t ΓHPFΓPkPkkkk ...~...321 −−=++= (5.127) A parcela 1tk é conhecida por rigidez material ou constitutiva e as parcelas 2tk e 3tk juntas formam a chamada rigidez geométrica do elemento. Vale ressaltar que a matriz de rigidez tangente do elemento é não simétrica, entretanto testes numéricos apresentados por Cole [1990] e Crisfield [1990] evidenciam que no equilíbrio de sistemas conservativos a matriz de rigidez tangente da estrutura é simétrica e consequentemente, na implementação computacional foi feita uma simetrização da matriz de rigidez tangente, calculando-a como a média entre a parte simétrica e a parte anti-simétrica segundo a equação abaixo:

( )Tttt kkk += .

21 (5.128)

55

5.6 ESFORÇOS RESULTANTES Os elementos de pórtico espacial estão submetidos a seis esforços resultantes em cada um

dos nós, representados na configuração atual por N, Vy, Vz, T, My e Mz, totalizando doze esforços por elemento. O esforço normal N, os cortantes Vy e Vz e o momento torçor T são constantes ao longo de todo o elemento, ao passo que os momentos fletores My = My (xe) e Mz = Mz (xe) variam linearmente ao longo do elemento, admitindo-se tratar de um modelo hermitiano. Estes esforços, com as respectivas convenções de sinais são mostrados a seguir na Fig. (5.3), sendo obtidos a partir das deformações de acordo com as seguintes equações apresentadas por Harrison [1973]:

dLAE

AEN0

00

... == ε c

xx

LIG

T θ.

=

( )cy

cy

yy L

IEM 21

01 2

..2θθ += ( )c

yc

yy

y LIE

M 210

2 2..2

θθ +=

( )cz

cz

zz L

IEM 21

01 2..2

θθ += ( )cz

cz

zz L

IEM 21

02 2..2

θθ += (5.129)

LMM

V zzy

211

+=

LMM

V zzy

212

+−=

LMM

V yyz

211

+−=

LMM

V yyz

212

+=

sendo E e G os módulos de elasticidade longitudinal e transversal; A0 a área da seção transversal; Ix o momento de inércia associado à torção de Saint Venant e Iy e Iz os momentos de inércia associados à flexão nas direções y e z; ε a deformação de engenharia e os ângulos ( )c

zc

yc

x 111 ,, θθθ e ( )cz

cy

cx 222 ,, θθθ definidos respectivamente pelas equações (5.16) e (5.17).

N

N

V1y

V2y

V1z

V2z

M1z

M2z

M1y

M2y

T

T

x0

y0

z0

Figura 5.3 – Esforços resultantes e convenção de sinais positivos

56

CAPÍTULO 6 FORMULAÇÃO CO-ROTACIONAL / ELEMENTO DE CASCA TRIANGULAR

6.1 DESCRIÇÃO CINEMÁTICA Neste capítulo será apresentada a descrição do movimento referente a formulação co-rotacional para um elemento de casca triangular, submetido à deslocamentos e rotações finitas, porém com deformações infinitesimais, seguindo-se em especial os trabalhos de Haugen [1994] e Felippa & Militello [1992]. Vale ressaltar que no presente capítulo serão feitas algumas mudanças de notação em relação aos capítulos anteriores de modo a facilitar o entendimento da formulação do elemento de casca, procurando-se manter ao máximo a notação empregada nos dois trabalhos comentados anteriormente. Considerando-se um elemento de casca triangular, conforme apresentado na Fig.(6.1), a sua configuração inicial é representada por C0, estando a mesma associada a um sistema de eixos locais definido em função da tríade de vetores ortogonais (i1

0, i20, i3

0), em relação a um sistema de eixos globais definido por (I1, I2, I3). Na configuração inicial, a transformação de um vetor referido em relação ao sistema global, definido por x, para o sistema local, definido por x~ , assim como a operação inversa, podem ser expressas pelas seguintes equações:

xTx .~0= xTx ~.0

T= ( )( )( )

=T

T

T

03

02

01

0

iii

T (6.1)

sendo T0 a matriz de transformação que representa o sistema local de eixos em C0.

I1I2

I3

i1 o

i2 o

i3 o

i1 n

i2 n

i3 n

C0

C0nCn

Figura 6.1 – Configurações inicial, co-rotacionada e deformada do elemento de casca

O movimento do elemento é caracterizado pelo deslocamento da configuração inicial C0

para a configuração atual ou deformada Cn, sendo a configuração co-rotacionada C0n uma configuração intermediária entre C0 e Cn, obtida por um movimento de corpo rígido a partir da configuração inicial. As configurações deformada e co-rotacional estão associadas a um mesmo sistema local de eixos definido em função da tríade de vetores ortogonais (i1

n, i2n, i3

n).

57

A origem do sistema local de eixos é obtida pelo deslocamento de corpo rígido do centróide do elemento, estando o mesmo associado a média aritmética das coordenadas nodais do elemento e em geral não coincide com o centro de massa do mesmo. Os sistemas de coordenadas nas configurações inicial a atual podem ser obtidos pelas rotações de corpo rígido R0n, a partir da configuração inicial:

0

0 . inn

i iRi = ⇔ Tn

Tn 00 .TRT = (6.2)

Uma vez que T0, Tn e R0n são matrizes ortogonais, a rotação de corpo rígido R0n pode ser

então obtida multiplicando-se a equação anterior por T0:

00 .TTR Tnn = ⇔ n

TTn TTR .00 = (6.3)

Considerando-se um nó qualquer do elemento de casca que se move da posição inicial r0 até a posição deformada rn, segundo a Fig.(6.2) apresentada a seguir, o deslocamento nodal pode ser expresso por: urr += 0n ⇒ 0rru −= n (6.4) De acordo com a formulação co-rotacional, o deslocamento pode então ser decomposto em duas parcelas, sendo a primeira parcela referente ao movimento de corpo rígido e a outra referente ao movimento deformacional. dr uuu += (6.5)

I1 I2

I3

i1 0

i2 0

i3 0

i1n

i2n

i3n

r 0

rc 0

r n

r 0n

ur

uud

uc

x0n

x0

C0

C0n

Cn

Figura 6.2 – Translações nodais no elemento de casca triangular

O deslocamento de corpo rígido ur que aparece na equação anterior pode ser obtido pela diferença na posição nodal entre as configurações inicial C0 e co-rotacional C0n. Analogamente, o deslocamento deformacional ud pode ser obtido pela diferença na posição nodal entre as configurações co-rotacional C0n e deformada Cn e portanto:

58

00 rru −= nr ⇔ 000 xrr += c (6.6)

nn

d0rru −= ⇔ 0

00000 .xRurxurr ncc

ncc

n ++=++= (6.7) sendo x0, o vetor que contém as coordenadas nodais na configuração inicial em relação ao centróide do elemento, conforme apresentado na Fig.(6.2). Combinando estas duas equações e após algebrismos simples, obtém-se o deslocamento de corpo rígido: ( ) 0

0 xIRuu −+= ncr (6.8)

Os deslocamentos deformacionais podem então ser obtidos substituindo a equação acima na Eq. (6.5), resultando em:

( )( )0

0 xIRuuuuu −+−=−= ncrd ( ) 0

0 xIRuu −−−= nc (6.9) que no sistema local pode ser expresso por: ( ) 0

0~.~~.~ xTIRTuuuTu T

nnncdnd −−−== ( ) 0

0~.~~~ xIRuu −−−= nc (6.10)

onde as quantidades medidas no sistema local de coordenadas (co-rotacional) são indicadas com um sinal de tio (tilde) sobreposto sobre o vetor ou matriz. Na equação acima, o vetor 0~x representa as coordenadas nodais da configuração co-rotacionada com relação ao sistema local, podendo ser expresso por: n

n0000 ..~ xTxTx == (6.11)

Já a rotação de um nó qualquer, associada ao movimento entre a configuração inicial C0 e a configuração deformada Cn pode ser descrita pelo tensor ou matriz de rotação R, conforme mostrado na Fig.(6.3). Seguindo-se os trabalhos de Rankin & Brogan [1986] e Nour-Omid & Rankin [1991] e procedendo de forma análoga ao estudo das translações, assume-se que o tensor R pode ser decomposto em uma rotação de corpo rígido R0n e uma rotação deformacional Rd, tal que: nd 0.RRR = (6.12) Substituindo-se a Eq.(6.3) na equação acima, a rotação deformacional no sistema global e no sistema local (co-rotacional) de coordenadas podem ser expressas respectivamente por: n

TTnd TTRRRR ... 00 == (6.13)

T

nT

nnT

nT

ndnd 00 ........~ TRTTTTRTTRTR === (6.14)

59

Conforme comentado anteriormente no capítulo 4, o tensor de rotação R pode ser expresso pela conhecida fórmula de Rodrigues, mostrada novamente na Eq.(6.15) para uma rápida visualização, a partir da qual pode-se obter o pseudo-vetor de rotação θ , utilizando-se um algoritmo atribuído à Spurrier, também comentado nos capítulos anteriores, sendo esta operação apresentada a seguir na Eq.(6.16):

( )( )

22

2

.2/

2/sen21sen

θθ θθ

θθ SSIR ++= , sendo ( )θS Spin=θ (6.15)

( )

==

3

2

1

θθθ

Rθ rot (6.16)

lembrando que a rotação deformacional expressa pelas equações (6.13) ou (6.14) é assumida pequena, porém finita.

I1 I2

I3

i1 0

i2 0

i3 0

i1n

i2n

i3n

R0n

R0n = TnTT0

RRd = RR0n

T

C0

C0n

Cn

Figura 6.3 – Rotações nodais no elemento de casca triangular

6.1.1 Vetor de deslocamentos deformacionais Um elemento de casca com N nós tem o seu vetor posição nas configurações inicial C0 e deformada Cn indicados respectivamente por xa

0 e xan, sendo “a” usado como índice que

representa o número do nó em questão. O campo de deslocamentos do elemento pode ser descrito em função dos deslocamentos nodais de translação ua e da matriz de rotação nodal Ra. Portanto, o conjunto (ua, Ra) para a = 1,...,N define o vetor de deslocamentos nodais v em relação ao sistema global de coordenadas. Na formulação co-rotacional, porém, é necessário estabelecer o vetor de deslocamentos nodais deformacionais vd com base no campo de deslocamentos nodais expresso por v . Este vetor vd contém as translações e rotações deformacionais ordenados da seguinte forma:

TdNdNddd θuθuv ......11= (6.17)

60

Na determinação do vetor de deslocamentos deformacionais, inicialmente é definido o sistema local de coordenadas T0 com base no vetor posição xa

0 na configuração inicial C0 e em seguida determina-se o vetor de coordenadas nodais no sistema local em C0:

( )000

0~caa xxTx −= ⇔ ∑

=

=N

aac N 1

00 1 xx (6.18)

a seguir, define-se o sistema local de coordenadas Tn com base no vetor posição xa

n = xa0 + ua

na configuração deformada C e de forma análoga determina-se o vetor de coordenadas nodais no sistema local em C:

( )nc

nan

na xxTx −=~ ⇔ ∑

=

=N

a

na

nc N 1

1 xx (6.19)

Uma vez conhecidos os vetores posição em C e C0 no sistema local, conforme as

equações (6.18) e (6.19), calcula-se para cada um dos nós a = 1,...,N as translações e rotações deformacionais de acordo com as seguintes equações: 0~~~

an

ada xxu −= (6.20) ( ) ( )T

andda rotrot 0..~~ TRTRθ == (6.21) e finalmente armazenam-se as translações e rotações deformacionais descritas pelas equações (6.20) e (6.21) no vetor de deslocamentos deformacionais no sistema local de eixos:

TdNdNddd θuθuv ~~......~~~

11= (6.22) Vale a pena ressaltar que, de forma similar ao caso de pórticos espaciais apresentado no capítulo anterior, durante o processo incremental-iterativo de resolução do sistema, efetua-se uma mudança na variável iterativa de rotação. Neste processo, os graus de liberdade são definidos pelo vetor δ v, no qual as translações são vetores de deslocamentos incrementais δ ua ao passo que as rotações são representadas por giros instantâneos da tríade nodal δ wa. O vetor de deslocamentos incrementais pode então ser representado por: T

NN wuwuv δδδδδ ......11= (6.23) Portanto, durante o processo de resolução do sistema e atualização dos deslocamentos, os deslocamentos translacionais são atualizados de forma consistente da forma tradicional ao passo que para os graus de liberdade rotacionais é feita uma atualização consistente da matriz de rotação associada ao giro instantâneo da tríade nodal, conforme as equações abaixo: aaa uuu δ+= (6.24) ( ) aaa RwRR .δ= (6.25)

61

6.2 VARIAÇÃO DOS DESLOCAMENTOS DEFORMACIONAIS O objetivo da presente seção é definir a variação do vetor de deslocamentos deformacionais do elemento no sistema co-rotacional vd

R, com relação aos graus de liberdade do elemento v, podendo esta variação ser expressa através da seguinte relação:

vv

vv δδ .

∂∂

=R

dRd (6.26)

sendo o índice R, utilizado para indicar que uma variável é expressa no sistema co-rotacional. Serão apresentados a seguir uma série de variações necessárias para a determinação da relação definida acima na Eq.(6.26). 6.2.1 Variação da matriz de transformação Tn A variação da matriz de transformação Tn pode ser definida em função da variação com relação a rotação instantânea da tríade nodal w~ no sistema local. Esta relação pode ser expressa conforme a seguinte equação apresentada por Nour-Omid & Rankin [1991]:

( ) nii

nn Spinw

wTw

TT .~~.~ δδδ −=

∂∂

= (6.27)

lembrando que Spin( w~ ) = Tn.Spin(w).Tn

T e que Tn é uma matriz ortogonal, então: ( ) ( )( ) ( )wTTTwTTwT δδδδ SpinSpinSpin nn

Tnnnn .....~ −=−=−= (6.28)

e uma vez que Spin(δ w) é uma matriz anti-simétrica, pode-se obter: ( )( ) ( )( ) ( ) T

nT

nT

nT

n SpinSpinSpin TwTwwTT ... δδδδ =−−=−= (6.29) 6.2.2 Variação da matriz de rotação de corpo rígido R0n A Eq.(6.29) apresentada acima pode ser utilizada na determinação da variação da matriz de rotação de corpo rígido R0n, definida anteriormente na Eq. (6.3), tal que: 00000 .).(... TTwTTTTTTR T

nT

nT

nT

nn Spin δδδδδ ==+= nSpin 0).( Rwδ= (6.30) que ao ser referenciada em relação ao sistema local fornece a seguinte equação: ( ) T

nnnT

nnnn Spin TRwTTRTR .).(...~000 δδδ ==

nSpin 0

~).~( Rwδ= (6.31)

62

6.2.3 Variação de ud com relação a v Para um elemento contendo N nós, o deslocamento do centróide uc pode ser definido como a média dos deslocamentos nodais:

∑=

=N

aac N 1

.1 uu (6.32)

Aplicando-se a Eq.(6.9) que contém os deslocamentos de translação deformacionais para o caso de um nó específico “a”, obtém-se a seguinte equação:

( ) 00

1..1

an

N

bbada N

xIRuuu −−

−= ∑

=

( ) 00

1.. an

N

bbab xIRuP −−

= ∑

=

, sendo IP .1

−=

Nabab δ (6.33)

cuja variação pode ser definida por

00

1.. an

N

bbabda xRuPu δδδ −

= ∑

=

(6.34)

lembrando que δPab = 0, δI = 0, δxa

0 = 0. Substituindo-se a Eq.(6.30) na segunda parcela do lado direito da Eq.(6.34), obtém-se: ( ) ( ) ( ) wxxwxRwxR δδδδ ..... 000

00

0n

an

aanan SpinSpinSpin −=== ( ) vGx δ..0n

aSpin−= (6.35) onde a matriz G relaciona a variação instantânea da tríade nodal δw com a variação dos deslocamentos nodais δv, através da seguinte expressão:

∑=

==N

aa

1.. avGvGw δδδ (6.36)

Esta matriz G pode ser decomposta em submatrizes nodais Ga e será comentada com

maiores detalhes para o caso de um elemento de casca triangular posteriormente. Substituindo-se a Eq.(6.35) em (6.34) obtém-se a seguinte variação:

[ ] ( )( )∑=

+=n

bbb

naabda Spin

1

0 .. vGx0Pu δδ (6.37)

O próximo passo será transformar a equação acima do sistema global para o sistema local

de coordenadas. Com este propósito, inicialmente aplica-se a Eq.(6.11) para o caso de um nó específico “a”, podendo obter-se a seguinte equação:

63

nanaa

000

0 ..~ xTxTx == (6.38) que simplesmente estabelece que as coordenadas da configuração co-rotacional no sistema local são constantes. Consequentemente, a transformação da Eq.(6.37) para o sistema local de coordenadas fornece a seguinte equação:

[ ] ( )( )∑=

+=n

bbb

nanabda Spin

1

0 ~.~..~~ vGxT0Pu δδ

[ ] ( )( )∑=

+=n

bbbaab Spin

1

0 ~.~.~~ vGx0P δ (6.39)

na qual o vetor bv~δ contém os graus de liberdade de translação e rotação de acordo com: [ ]T

bbb wuv ~~~ δδδ = (6.40) 6.2.4 Variação de ud co-rotacional com relação a v Para o caso de um nó específico “a” do elemento, a relação entre as coordenadas nodais na configuração co-rotacional C0n e na configuração deformada Cn pode ser expressa por: da

na

na uxx += 0 (6.41)

cuja variação, substituindo-se a Eq.(6.34), pode ser definida através da seguinte equação:

00

1

00 . an

N

bab

nada

na

na xRuPxuxx b δδδδδδ −

−+=+= ∑

=

(6.42)

Lembrando a Eq.(6.35), tem-se: 0

00 . an

na xRx = (6.43)

cuja variação pode ser expressa por 0

00

00

00 ... ananan

na xRxRxRx δδδδ =+= (6.44)

que ao ser substituída na Eq.(6.42) fornece a seguinte equação:

−=−

−+= ∑∑

==

N

bab

na

N

bab

na

na

1

0

1

0bb uPxuPxx δδδδδ (6.45)

Calculando-se a variação da Eq.(6.41) no sistema co-rotacional: R

daR

danR

aRn

a uuxx δδδδ =+= 0 (6.46)

64

uma vez que as coordenadas nodais na configuração co-rotacional xa0n são constantes. A

variação de um vetor arbitrário pode ser definida segundo Haugen [1994] por: xwxx ×+= δδδ R ⇒ xwxx ×−= δδδ R (6.47) sendo δxR a variação do vetor no sistema co-rotacional e δw a rotação deste sistema em relação ao sistema global. Esta expressão pode então ser utilizada para calcular a variação do vetor de deslocamentos deformacionais ud no sistema co-rotacional, ou seja, δud

R. Portanto, aplicando-se a Eq.(6.47) para um nó específico “a” do elemento na configuração deformada Cn e substituindo-se as Eq.(6.45) e (6.46), pode-se obter: n

aan

aRn

a xwxx ×−= δδδ

( ) naa

N

bbab

naa

N

bbab

Rda Spin xwuPxwuPu ...

11δδδδδ −

=×−

= ∑∑

==

( ) [ ] ( )( ) bvGx0PwxuP δδδ ∑∑==

+=+

=

N

bb

naaba

na

N

bbab SpinSpin

11.... (6.48)

que transformando-se para o sistema local de coordenadas fornece:

[ ] ( )( ) bvGx0Pu ~.~.~~~1

δδ ∑=

+=N

bb

naab

Rda Spin (6.49)

6.2.5 Variação de θd co-rotacional com relação a wd co-rotacional A variação do pseudovetor de rotação deformacional θda para o caso específico de um nó “a” do elemento requer a sua variação com relação ao giro instantâneo da tríade nodal δwa, uma vez que foi efetuada uma mudança na variável iterativa de rotação, conforme comentado anteriormente para o caso de pórticos espaciais, tal que:

dadada

Rda wΛw

wθθ δδδ .. =

∂∂

= (6.50)

sendo a matriz Λ calculada anteriormente nas Eq.(5.49) e (5.50), sendo apresentadas novamente nas Eq.(6.51a) e (6.51b), para facilitar a sua visualização, podendo-se observar uma pequena mudança na notação em relação à formulação apresentada para o caso de elementos de pórticos espaciais.

( ) ( )θθIΛ 2..21 SpinSpin ξ+−= (6.51a)

( )θθ

θθθξsen..2

cos1.sen.22

+−= (6.51b)

65

6.2.6 Variação de θd co-rotacional com relação a v A parcela deformacional do giro instantâneo da tríade nodal wda pode ser obtido pela diferença entre o giro instantâneo total e a parcela referente ao giro de corpo rígido: rada www −= (6.52) cuja variação fornece a seguinte expressão:

vGwvvwwwww δδδδδδδ .. −=

∂∂

−=−= aii

rarada (6.53)

sendo

vGvvww δδδ .. =

∂∂

= ii

rr (6.54)

que fornece a seguinte equação para a variação da parcela deformacional do giro instantâneo da tríade nodal:

[ ]( )∑=

−=N

bbbabda

1.. vGI0w δδδ (6.55)

Finalmente, combinando-se as Eq.(6.50) e (6.55), pode-se então obter a variação do pseudovetor de rotação θda em relação aos graus de liberdade v no sistema co-rotacional, resultando na seguinte expressão:

[ ]( )∑=

−==N

bbbabda

Rda

1.... vGI0ΛwΛθ δδδδ (6.56)

6.2.7 Variação de vd co-rotacional com relação a v Uma vez conhecidas as variações das translações e rotações deformacionais em relação aos graus de liberdade v no sistema co-rotacional, segundo as Eq.(6.48) e (6.56), para o caso de um nó específico “a” do elemento e agrupando-as em um único vetor que contém translações e rotações, obtém-se:

( ) ( )

b

N

bb

na

Rda

RdaR

daSpinN

vGIx

000I

I00I

Λ00I

θuv δ

δδδ ..

./1.

1∑

=

−+

−+

=

=

( ) ( )

b

N

bb

naSpinN

vGI

x000I

I00I

Λ00I

δ.../1

.1

∑=

−−

−

=

( ) b

N

bRabTabaa vPPIH δ..

1∑

=

−−= (6.57)

66

Uma vez definida a Eq.(6.57) para um nó específico “a” do elemento, a Eq.(6.26) do início da seção, pode então ser finalmente expressa através da seguinte igualdade:

vPHvv

vv δδδ ... =

∂∂

=R

dRd (6.58)

onde

=

NNH0

0HH

K

MOM

L11

, sendo

=

Λ00I

H aa (6.59)

sendo Λ definido anteriormente na Eq.(6.51), ao passo que P é um operador de projeção não-linear expresso pela seguinte expressão: ( )RT PPIP −−= (6.60) na qual a parcela PT está associada às translações do centróide do elemento e a parcela PR está associada à rotação de corpo rígido em torno do centróide do mesmo, sendo expressos respectivamente, por:

=

TNNTN

NTT

T

PP

PPP

L

MOM

L

1

111

, sendo ( )

=

000I

P./1 N

Tab (6.61)

( )

( )G

Ix

Ix

GSP ..

1

−

−

==n

N

n

R

Spin

Spin

M (6.62)

6.3 VARIAÇÕES DE ALTA ORDEM Para o desenvolvimento do vetor de forças internas e da matriz de rigidez tangente nas próximas seções, é necessário conhecer a variação do operador de projeção não-linear PT pós multiplicado por um vetor de forças f, onde f é considerado constante em relação a variações. De forma similar, é também necessário o conhecimento da variação da matriz ΛT pós multiplicada por um vetor de momentos m, sendo m constante em relação a variações. Estas operações são em geral conhecidas na álgebra tensorial como contrações. Inicialmente, será determinada a variação do vetor ΛT definido na Eq.(6.51) após contração com um vetor de momentos m, tal que:

wMmΛmΛ δδδ .... =∂∂

= ii

TT w

w (6.63)

67

Esta relação foi apresentada de forma detalhada por Nour-Omid & Rankin [1991] e Monteiro [2004], podendo ser expressa por:

( ) ( ) ( ) ΛθmθθmmθImθmM ......2.....21 2

+−++−= TTTT SpinSpin µξ (6.64)

( ) ( )( )2/sen4

2/sen.8sen.2

2

θθθθθθµ −+

= (6.65)

O próximo passo é calcular a variação do operador de projeção PT contraído com o vetor de forças internas. Inicialmente, calcula-se a variação de P: GSGSPPPIP δδδδδδδ .. −−=−=−−= RRT (6.66) O vetor de forças internas f pode ser decomposto em um vetor de forças internas auto-equilibrado fb e um vetor de forças internas não-equilibrado fu: ub fff += ⇔ fPf .T

b = e ( ) fPIf .Tu −= (6.67)

de modo que a variação de PT contraída com o vetor de forças internas pode ser expresso por: ( )( )ub

TTTT ffSGSGfP T ++−= .... δδδ ( ) u

TTTTb

TTb

TT fSGSGfSGfSG ....... δδδδ +−−−= ( ) ub

TTu

TTTTb

TT fPfSGfSGSGfSG T ........ δδδδδ −−=+−−= (6.68) onde ST.fb = 0, uma vez que S simplesmente esta associada à rotações de corpo rígido e não gera nenhum trabalho ao ser multiplicada por um vetor de forças auto-equilibradas. Segundo Haugen [1994] o termo δPT.fu, em geral tende a zero quando C0n e Cn estão próximas uma vez que fu tende a zero. 6.4 VETOR DE FORÇAS INTERNAS Nesta seção será apresentada a determinação do vetor de forças internas segundo a formulação CSSE (Consistent Symmetrizable Self-Equilibrated) descrita por Haugen [1994], na qual é obtida um vetor de forças internas auto-equilibradas, que por sua vez gera uma matriz de rigidez simetrizável através da variação consistente do vetor de forças internas. O termo simetrizável significa que embora a matriz de rigidez tangente não seja simétrica, a razão de convergência nas iterações próximas do equilíbrio não é prejudicada quando é feita uma simetrização da matriz de rigidez tangente conforme comentado no capítulo 5. Na determinação do vetor de forças internas será admitido uma relação deformação-deslocamento linear, em relação ao sistema local de coordenadas das configurações C0n e Cn:

dvBε ~.~~ = ⇒ ddd

vBvvεε ~.~~.~~~ δδδ =

∂∂

= (6.69)

68

na qual dv~ é o vetor de deslocamentos deformacionais escrito em relação ao sistema local. Já as relações constitutivas no sistema local de coordenadas por sua vez, serão expressas por: εCσ ~.~~ = (6.70) sendoC~ uma matriz constante para o caso de uma material linear elástico. Aplicando-se a matriz de transformação T sobre a Eq.(6.58), pode-se obter esta relação, no sistema local de coordenadas, tal que: vTPHvPHv δδδ ..~.~~.~.~~ ==d (6.71) A variação do trabalho interno (δWi) proveniente do somatório do trabalho virtual realizado pela distribuição de tensões que atuam em cada um dos elementos pode ser obtida através da seguinte expressão: ( )∑∫−=

V

Ti dVW σε ~.~δδ (6.72)

Finalmente, substituindo-se as Eq.(6.69), (6.70) e (6.71) em (6.72), obtém-se: ( )( )∑∫−=

V

TTdi dVW εCBv ~.~.~.~δδ

( )∑∫−=

Vd

TTTTT dVvBCBHPTv ~.~.~.~.~.~..δ

∑−= fv .Tδ (6.73) sendo o vetor de forças internas no sistema global de coordenadas expresso por: ( ) ( )∫=

Vd

TTTT dV vBCBHPTvf ~.~.~.~.~.~.ˆ (6.74)

( ) fTfHPTvKHPTvf ~.~.~.~.~.~.~.~.ˆ T

eTTT

deTTT === (6.75)

onde eK~ é a matriz de rigidez linear do elemento; ef~ é o vetor de forças internas associado à

matriz de rigidez linear e ao vetor de deslocamentos deformacionais no sistema local; e f~ o vetor de forças internas no sistema local de coordenadas. Para o caso de um elemento com N nós e seis graus de liberdade por nó, sendo três graus de liberdade de translação e três graus de liberdade de rotação, a matriz de transformação T que aparece nas equações acima consiste de matrizes Tn repetidas 2N vezes ao longo da diagonal, segundo a expressão abaixo:

=

n

n

T0

0TT O (6.76)

69

6.5 MATRIZ DE RIGIDEZ TANGENTE CONSISTENTE A matriz de rigidez tangente K é dita consistente em relação ao vetor de forças internas para o caso específico desta matriz ser obtida pela variação de f com relação à variações de v:

vKvvfvff δδδδ ... =

∂∂

=∂∂

= iiv

(6.77)

Substituindo-se a Eq.(6.75) em (6.77), a variação de f com relação a v pode ser expressa por: e

TTe

TTe

TTe

TT fHPTfHPTfHPTfHPTf TTTT ~.~.~.~.~.~.~.~.~.~.~.~. δδδδδ +++= ( ) vKvKKKK δδ .. =+++= MGMGPGR (6.78) sendo KM a matriz de rigidez tangente material da estrutura ao passo que a matriz de rigidez tangente geométrica é obtida pela contribuição de três parcelas: KGR, KGP e KGM. A seguir, serão obtidas cada uma destas matrizes que compõem a matriz de rigidez tangente. 6.5.1 Matriz de rigidez tangente material KM

O último termo do lado direito da Eq.(6.78) gera a chamada matriz de rigidez tangente material KM:

deTT

deTT

eTT vKHPTvKHPTfHPT TTT ~.~.~.~.~.~.~.~.~.~.~. δδδ +=

vPHKHPT T ~.~.~.~.~.~. δe

TT= vTPHKHPT T δ..~.~.~.~.~. e

TT= vK δ.M= (6.79) sendo TKTTPHKHPTK TT .~..~.~.~.~.~. Me

TTM == (6.80)

PHKHPK ~.~.~.~.~~

eTT

M = (6.81) 6.5.2 Matriz de rigidez tangente geométrica KGR A matriz de rigidez tangente KGR está associada à variação da matriz de transformação T:

( )( )

( )( )

=

==

NrT

n

NrT

n

rT

n

rT

n

NT

n

NT

n

Tn

Tn

Te

TTT

SpinSpin

SpinSpin

mwTnwT

mwTnwT

mTnT

mTnT

fTfHPT

~.~.

~.~.

~.~.

~.~.

~.

~.

~.

~.

~.~.~.~.1

1

1

1

δδ

δδ

δδ

δδ

δδ MM (6.82)

70

eTTT fHPT ~.~.~.δ

( )( )

( )( )

( )( )

( )( )

r

N

N

T

rN

rN

r

r

T

SpinSpin

SpinSpin

SpinSpin

SpinSpin

w

mn

mn

T

wmwn

wmwn

T ~.

.~.~

.~.~

.

~.~~.~

~.~~.~

.1

1

1

1

δ

δδ

δδ

−=

−= MM

vTGFTwFT δδ ....~.. nm

Trnm

T −=−= vK δ.GR= (6.83) sendo TKTTGFTK .~..... GR

Tnm

TGR =−= ⇒ GFK .~

nmGR −= (6.84)

( )( )

( )( )

=

N

N

nm

SpinSpin

SpinSpin

mn

mn

F

~~

~~

~ 1

1

M eTT

N

N

fHP

mn

mn

f ~.~.~

~~

~~

~ 1

1

=

= M (6.85)

6.5.3 Matriz de rigidez tangente geométrica KGM A matriz de rigidez tangente geométrica KGM corresponde ao terceiro termo do lado direito da Eq.(6.78) e está associada a variação da matriz TH~δ :

=

=

dNN

dTT

NT

N

T

TTe

TTT

wM0

wM0

PT

mΛ0

mΛ0

PTfHPT

~.~

~.~

.~.

~.~

~.~

.~.~.~.~.1111

δ

δ

δ

δδ MM

d

N

TT v

M0

M0

PT ~.

~

~

.~.1

δ

= O

vTPMPT δ..~.~.~. TT= vK δ.GM= (6.86) com aM~ definido conforme a Eq.(6.63), tal que: TKTTPMPTK .~..~.~.~. GM

TTTGM == ⇒ PMPK ~.~.~~ T

GM = (6.87)

71

6.5.4 Matriz de rigidez tangente geométrica KGP A matriz de rigidez tangente geométrica KGP está associada à variação do operador de projeção não linear. Efetuando-se a decomposição do vetor de forças contraído e

T fH ~.~ em

uma parcela de forças auto-equilibradas bf~ e outra não equilibrada uf~ , conforme comentado nas seções anteriores, o segundo termo do lado direito da Eq.(6.78) pode ser transformado em: ( ) e

TTTTe

TTTTTe

TTT fHPIPTfHPSGTfHPT ~.~..~.~.~.~.~.~.~.~.~. −+−= δδδ u

TTb

TTT fPTfSGT ~.~.~.~.~. δδ +−= (6.88) De acordo com Haugen [1994], o termo u

TT fPT ~.~.δ , associado às forças não-equilibradas pode ser desprezado uma vez que produz uma contribuição muito pequena em relação à contribuição das forças auto-equilibradas, especialmente se as configurações co-rotacional C0n e deformada Cn estiverem próximas. Portanto desprezando a segunda parcela da Eq. (6.88):

( )[ ]

−= ∑

= i

iN

i

ni

TTe

TTT Spinmn

0xGTfHPT ~~

.~.~.~.~.~.1

δδ

( )[ ]

−−= ∑

=n

i

ni

N

ii

TT Spinθx0nGT ~~

.~.~.1 δ

δ

( )[ ]

−−= ∑

= di

diN

ii

TT Spinθu

0nGT ~~

.~.~.1 δ

δ

vKvTPFGT δδ ...~.~.~. GP

Tn

TT =−= (6.89) sendo TKTTPFGTK .~..~.~.~. GP

TTn

TTGP =−= ⇒ PFGK ~.~.~~ T

nT

GP −= (6.90) ( ) ( )[ ]T

Nn SpinSpin 0n0nF ~~~1 L= (6.91)

6.5.5 Matriz de rigidez tangente Finalmente, substituindo-se as Eq.(6.80), (6.84), (6.87) e (6.90) em (6.78), pode-se obter a matriz de rigidez tangente, conforme a seguinte equação:

( )TKKKKTK .~~~~. GPGRGMMT +++=

( )TPFGGFPMPPHKHPT T .~.~.~~.~~.~.~~.~.~.~.~ T

nT

nmT

eTT −−+= (6.92)

72

CAPÍTULO 7 ELEMENTO DE CASCA TRIANGULAR

7.1 INTRODUÇÃO Uma das grandes vantagens da formulação co-rotacional é a possibilidade de re-utilização de elementos finitos lineares para a análise de problemas envolvendo grandes deslocamentos e rotações, desde que seja admitido um regime de pequenas deformações. Neste capítulo será apresentado um elemento finito de casca linear que fornece a matriz de rigidez linear Ke e o vetor de forças internas fe, utilizados na análise não linear geométrica segundo a formulação co-rotacional de cascas apresentada no capítulo anterior. O elemento finito linear de casca triangular utilizado no presente capítulo é baseado na formulação ANDES (Assumed Natural Deviatoric Strains), conforme apresentado nos trabalhos de Haugen [1994], Alvin et al. [1992] e Felippa & Militello [1992]. Vale ressaltar que este elemento possui o grau de liberdade de rotação torsional (drilling) como parte dos componentes do elemento de membrana. Uma das características marcantes da formulação ANDES é a decomposição da matriz de rigidez linear Ke em duas parcelas, de acordo com a seguinte equação: ( ) ehbe fvKKvK =+= .. (7.1) sendo Kb a matriz de rigidez básica, construída de modo a garantir a convergência e Kh a matriz de rigidez de alta ordem, construída com o intuito de garantir a estabilidade e exatidão dos resultados. A matriz Kb independe da formulação utilizada uma vez que a mesma é inteiramente definida em função de uma hipótese de tensão constante e campo de deslocamentos no contorno pré-definido, sendo desenvolvida pela primeira vez, de uma forma mais elaborada, segundo a formulação FF (Free Formulation) por Bergan & Nygard [1989]. Já a matriz Kh, por outro lado, pode ser obtida utilizando-se diferentes formulações, tais como FF (Free Formulation), EFF (Extended Free Formulation) e ANDES. 7.2 MATRIZ DE RIGIDEZ BÁSICA DO ELEMENTO ANDES De acordo com os trabalhos apresentados por Alvin el at. [1992] e Haugen [1994], o procedimento para a construção da matriz de rigidez básica Kb do elemento pode ser descrito em quatro etapas: Etapa 1. Deve ser assumido um estado de tensão σ constante no interior do elemento, ao qual estão associadas as tensões no contorno σn: σTnσσ .. normaln == (7.2) sendo n o vetor unitário normal ao contorno do elemento e Tnormal uma matriz de transformação que simplesmente substitui o tensor produto σni = σij.nj por uma multiplicação matricial equivalente. Etapa 2. Conecta-se o campo de deslocamentos d do contorno aos graus de liberdade v: vNd .d= (7.3)

73

onde a matriz Nd contém as funções dos deslocamentos do contorno que devem satisfazer a continuidade entre os elementos e incluir os movimentos associados à deslocamentos de corpo rígido e à deformações constantes. Vale ressaltar que o campo de deslocamentos d do elemento não é construído de forma explícita na formulação ANDES. Etapa 3. Constrói-se a matriz L (force-lumping matrix), que de forma consistente relaciona as tensões do contorno σn às forças nodais do elemento, que por sua vez estão conjugadas aos graus de liberdade v, mediante a aplicação do princípio dos trabalhos virtuais: ∫ ∫ ∫ ====

S Se

TT

S

Tdn

Tnormal

Td

Tn

T dSdSdS fvσLvσNvσTNvσd .......... δδδδδ (7.4)

que fornece a seguinte equação para a matriz L: ∫∫ ==

S

Tdn

Snormal

Td dSdS NTNL . (7.5)

Etapa 4. Finalmente, a matriz de rigidez básica do elemento pode ser construída:

Tb V

LCLK ...1= (7.6)

sendo C a matriz tensão-deformação que representa as relações constitutivas do material e V o volume de um elemento tridimensional, o qual é substituído pela área ou pelo comprimento, caso o elemento seja bidimensional ou unidimensional, respectivamente. 7.3 MATRIZ DE RIGIDEZ DE ALTA ORDEM DO ELEMENTO ANDES Segundo os trabalhos apresentados por Felippa & Militello [1992] e Haugen [1994], o comportamento de alta ordem do elemento está fundamentado no uso da distribuição de tensões assumidas, ao invés de utilizar os modos de deslocamento. Existem basicamente, cinco etapas que devem ser seguidas para a construção da matriz de rigidez de alta ordem do elemento ANDES: Etapa 1. Selecionam-se locais no elemento onde as medidas de deformações em coordenadas naturais devem ser determinadas (natural straingage locations). Para a maioria dos elementos ANDES, estes locais para medição de deformações são escolhidos sobre as chamadas linhas de referência, mas isto não é uma regra geral. No caso do presente trabalho os locais para determinação das deformações serão os nós de canto do elemento triangular, ao passo que as linhas de referência serão os lados do elemento. Através de interpolações apropriadas, expressam-se as deformações naturais do elemento ε em função das leituras de medidas de deformação nestes locais: gAε .ε= (7.7) onde ε é um campo de deformações em coordenadas naturais, que deve incluir todos os estados de deformações constantes.

74

Etapa2. Fazer a transformação das deformações naturais ε em deformações cartesianas e: gAgATεTe .... === εdefdef (7.8) para cada ponto do elemento. Caso e ≡ ε ou caso seja possível trabalhar em coordenadas naturais esta etapa pode ser desprezada. Este é geralmente o caso quando Tdef é constante ao longo de todo o elemento, conforme ocorre para o elemento de casca triangular apresentado no presente trabalho. Etapa 3. Relacionam-se as leituras de medidas de deformações naturais g com os graus de liberdade v através da seguinte operação: vQg .= (7.9) onde Q é uma matriz de transformação que relaciona os deslocamentos com as deformações nodais. As técnicas para se fazer isto variam de elemento para elemento, de modo que não se tem uma regra padrão para o procedimento. Freqüentemente, por simplicidade, a análise é fracionada em subproblemas: ..... 2211 ++= vQvQg (7.10) onde v1, v2, ... são subconjuntos de v convenientemente escolhidos. Etapa 4. Separa-se o campo de deformações cartesianas em deformações principais e deformações desviadoras: ( )gAAeee .dd +=+= (7.11) sendo:

∫=V

def dVV εATA ..1 ∫ ∫ ==

V Vdd dVdV 0.gAe (7.12)

lembrando que para o caso da matriz Tdef ser constante, esta etapa pode ser feita empregando-se deformações naturais. Etapa 5. Finalmente, a matriz de rigidez de alta ordem pode ser expressa por: QKQK ... d

Th β= ⇔ ∫=

Vd

Tdd dVACAK .. (7.13)

onde β > 0 é um fator de escala utilizado apenas para melhorar a performance no caso de malhas pouco refinadas. Normalmente, combina-se o produto AQ em uma matriz de deslocamento-deformação única B, a qual é dividida em duas parcelas de forma análoga à divisão feita anteriormente com a matriz A: ( ) ( ) vBvBBvQAAvQAe ...... =+=+== dd (7.14)

75

de modo que ∫=

Vd

Tdh dVBCBK ...β (7.15)

O próximo passo será aplicar estas regras para a construção de um elemento finito de casca triangular contendo três nós. Uma vez que o elemento está contido em um plano local xy, conforme comentado a seguir, os esforços de membrana e de flexão podem ser desenvolvidos separadamente. 7.4 DEFINIÇÕES GEOMÉTRICAS PARA O ELEMENTO TRIANGULAR A geometria do elemento finito triangular utilizado no presente trabalho é apresentada na Fig.(7.1). O elemento apresenta lados retos e a sua geometria é completamente definida pela posição dos seus três nós de canto definidos por 1, 2 e 3, numerados no sentido anti-horário. O elemento é referenciado em relação à um sistema local de eixos cartesianos (x,y). Diferenças nodais são abreviadas escrevendo-se xij = xi – xj ou yij = yi – yj . Definindo-se li como sendo o comprimento do lado oposto ao nó de canto i e hi como sendo a altura do nó i até o lado i, conforme apresentados na Fig.(7.1), pode-se definir as seguintes relações:

22jkjki yxl += e

ii l

Ah .2= (7.16)

onde A é a área do elemento, podendo ser calculada por:

( ) ( ) ( )122131132332

321

321

111det.

21 yxyxyxyxyxyx

yyyxxxA −+−+−=

= (7.17)

valendo lembrar que área do elemento medida segundo a equação acima é sempre positiva, desde que os nós sejam numerados no sentido anti-horário.

(x1,y1) (x2,y2)

(x3,y3)

l1

l2

l3

1 2

3

h1

h2

h3 n1

n2

n3

s1

s2

s3

x

y

1 2

3

Figura 7.1 – Dimensões geométricas e vetores unitários no elemento triangular

76

Os vetores unitários ni e si, que definem respectivamente as direções normais e paralelas à direção do lado li podem ser expressas pelas seguintes equações:

=

=

0.1

0kj

kj

iiy

ix

i yx

lss

s e

−=

=

00ix

iy

iy

ix

i ss

nn

n (7.18)

Os pontos de um elemento triangular podem ser também definidos em relação à um sistema de coordenadas paramétrico (ζ1, ζ2, ζ3), sendo este sistema usualmente conhecido por sistema de coordenadas triangulares. Este sistema de coordenadas pode ser representado por um conjunto de linhas retas paralelas ao lado oposto ao nó de canto i em questão, conforme mostrado na Fig.(7.2). Neste sistema, as equações dos lados 23, 31 e 12 são representadas respectivamente por: ζ1 = 0, ζ2 = 0 e ζ3 =0. Já as coordenadas dos três nós de canto 1, 2 e 3 podem ser representadas por (1,0,0), (0,1,0) e (0,0,1). Os três pontos médios dos lados 12, 23 e 31 tem coordenadas (1/2, 1/2,0), (0, 1/2, 1/2) e (1/2,0, 1/2) e o centróide (1/3,1/3,1/3).

ζ1= 1

ζ1= 0

3/4 1/2

1/4 ζ2= 1

ζ3= 1

ζ2= 0

ζ3= 03/4

1/2

1/2

1/4

1/4

1

2 2

2

3 3 3 3/4

1 1

Figura 7.2 – Sistema de coordenadas triangulares

Vale a pena ressaltar que o sistema de coordenadas triangulares está associado à seguinte propriedade ou restrição: 1321 =++ ςςς (7.19) sendo a relação entre o sistema de coordenadas triangulares e cartesiano expresso por:

( ) ( )[ ]kjjkjkkii xyyyxxyxyxA 00.

21

−+−+−=ς (7.20)

onde x0 e y0 representam as coordenadas do centróide em relação ao sistema de coordenadas cartesiano, tal que:

( )3210 .31 xxxx ++= e ( )3210 .

31 yyyy ++= (7.21)

77

7.5 ELEMENTO TRIANGULAR DE MEMBRANA O elemento triangular de membrana utilizado no presente trabalho é baseado no elemento do tipo ANDES desenvolvido por Felippa & Militello [1992] e utilizado por Haugen [1994]. Os graus de liberdade vi do elemento de membrana consistem de translações planas u e v, respectivamente nas direções locais x e y, e na rotação torsional (drilling) θz:

=

zi

i

i

i vu

θv sendo

∂∂

−∂∂

=yu

xv

z .21θ (7.22)

7.5.1 Matriz de rigidez básica Ordenando-se os graus de liberdade em cada um dos nós como translações nas direções x e y e rotações torsionais em torno do eixo z, conforme apresentado na Eq.(7.22), a matriz L (force-lumping matrix) que associa as tensões de superfície com as forças nodais, definida na Eq.(7.4) pode ser expressa por:

σLf .=e sendo

=

xy

yy

xx

τσσ

σ (7.23)

=

3

2

1

LLL

L

=

3

2

1

e

e

e

e

fff

f

=

zi

yi

xi

ei

mff

f (7.24)

Segundo Felippa & Militello [1992], as tensões do contorno no elemento de membrana em um nó j serão somente função dos lados adjacentes ij e jk de modo que a matrix L total de cada um dos elementos pode ser subdividida em função da contribuição de cada um dos nós separadamente, conforme apresentando na Eq.(7.24), tal que:

( ) ( ) ( )

−−−−

−=

3/...6/.6/.0

0.

21

2222ijijkjkjkjijkjij

kiki

kiki

j

yxyxxxyyyxxy

αααL (7.25)

onde os índices nodais (i,j,k) estão submetidos à permutações cíclicas de (1,2,3) e o termo α é apenas um fator de peso utilizado para melhorar a performance do elemento. A matriz de rigidez básica pode então ser computada por:

Tb A

t LCLK ...= (7.26)

sendo t a espessura do elemento de casca e A a área do elemento.

78

7.5.2 Matriz de rigidez de alta ordem Felippa & Militello [1992] obtiveram o comportamento de alta ordem do elemento triangular do tipo ANDES em função dos graus de liberdade de alta ordem θhi, definidos como sendo os graus de liberdade de rotação torsional (drilling) θi menos as rotações associadas aos deslocamentos de corpo rígido e às deformações constantes θ0 do elemento CST (Constant Strain Triangle): 0θθθ −= ihi (7.27) sendo

[ ]v.000.41

2121131332320 yxyxyxA

−−−−−−=θ (7.28)

de modo que substituindo-se a Eq.(7.28) em (7.27), obtém-se:

v.400040004

..41

212113133232

212113133232

212113133232

3

2

1

=

AyxyxyxyxAyxyxyxyxAyx

Ah

h

h

θθθ

(7.29)

O próximo passo será a subdivisão das rotações de alta ordem em rotações principais ( ) 3/321 hhh θθθθ ++= e rotação desviadoras θθθ −= hii

' :

v.

31

4431

4431

44

3200

3100

3100

3100

3200

3100

3100

3100

3200

212113133232

'3

'2

'1

−−

−−

−−

=

Ay

Ax

Ay

Ax

Ay

Axθ

θθθ

(7.30)

sendo [ ]T

zzz vuvuvu 333222111 θθθ=v (7.31) que em forma matricial pode ser expresso pela seguinte equação: vHθ .vh θ= (7.32) Vale a pena enfatizar que a decomposição das rotações de alta ordem está associada com um decomposição similar das deformações de alta ordem que em coordenadas naturais podem ser expressas por: tb εεε += (7.33)

79

Campo de Deformações Flexional:

O vetor εb que aparece na equação anterior representa o campo de deformações flexionais de alta ordem em coordenadas naturais. Este campo de deformações é representado por deformações medidas em relação à direções paralelas aos lados do triângulo, no sentido anti-horário, sendo armazenadas para cada um dos nós de canto da seguinte forma: T

bbbb 133221 εεε=ε (7.34) estando estas deformações relacionadas às rotações de alta ordem desviadoras '

iθ através das seguintes operações:

'1

'3

'2

'1

1132

1134

1131

1323

1323

1325

1214

1212

1211

113

132

121

1 .. θQε b

b

b

b

b =

−−

−=

=θθθ

χρχρχρχρχρχρ

χρχρχρ

εεε

'2

'3

'2

'1

2133

2135

2133

2322

2321

2324

2214

2211

2212

313

232

221

2 .. θQε b

b

b

b

b =

−−

−=

=θθθ

χρχρχρχρχρχρ

χρχρχρ

εεε

(7.35)

'3

'3

'2

'1

3131

3134

3132

3321

3322

3324

3215

3213

3213

313

332

321

3 .. θQε b

b

b

b

b =

−−

−=

=θθθ

χρχρχρχρχρχρ

χρχρχρ

εεε

onde

2.34

ij

kij l

A=χ e 2.3

2

ij

jij

iij l

A−== χχ (7.36)

de modo que as deformações naturais ε bjk no nó de canto i foram definidas como ε bjk

i e o vetor ε b no nó de canto i definido por ε bi. As constantes ρi são apenas coeficientes numéricos a serem escolhidos de modo a melhorar a performance do elemento e segundo recomendação de Haugen [1994] devem ser utilizados os seguintes valores: ρ2 = 1, ρ3 = 1/2, ρ1 = ρ4 = ρ5 = 0 (7.37)

Uma vez conhecidas as matrizes Qbi, segundo a Eq.(7.35), as deformações de flexão de alta ordem que estão associadas às rotações de alta ordem desviadoras podem ser então interpoladas linearmente entre os nós do elemento triangular conforme a seguinte operação: ( ) ''

332211 ..... θBθQQQε bbbbb =++= ςςς (7.38)

80

Campo de Deformações Torsional: O campo de deformações torsional, definido em coordenadas naturais pelo vetor εt na Eq.(7.33), por sua vez está associado as rotações de alta ordem principais θ . Este campo pode ser obtido aplicando-se a cada um dos cantos do elemento triangular as mesmas rotações e ao mesmo tempo mantendo-se os deslocamentos dos cantos impedidos. As deformações torsionais estão relacionadas às rotações desviadoras θ através da seguinte relação:

θ.θ.3

133

13

322

32

211

21

13

32

21

t

t

t

t

t Bε =

=

=

ςχςχςχ

εεε

(7.39)

sendo ζij = ζi - ζj. De modo a facilitar a combinação com o campo de deformações flexional, é conveniente definir o campo de deformações torsional expandido, no qual cada uma das colunas recebem 1/3 da deformação:

=

=

θθθ

.

213

13213

13213

13

212

32212

32212

32

211

21211

21211

21

13

32

21

ςχςχςχςχςχςχςχςχςχ

εεε

t

t

t

tε (7.40)

Campo de Deformações Total: O campo de deformações total em coordenadas naturais pode então ser obtido pela combinação dos campos de deformação flexional e torsional, sendo expresso em relação aos graus de liberdade v através da seguinte equação: tb εεε +=

θ.. 'tb BθB +=

[ ] htb θBB .= [ ] vHBB .. vθtb= vB.= (7.41) 7.5.3 Matriz de Rigidez: Uma vez conhecido o campo de deformações total, a matriz de rigidez de alta ordem do elemento de membrana pode então ser calculada da seguinte forma: ∫=

A

Th dABCBK .. ε sendo defdef TCTC ..1−=ε (7.42)

onde a matriz Tdef transforma as deformações do sistema de coordenadas naturais para o sistema de coordenadas cartesiano e a matriz Tdef

-1 realiza a operação inversa.

81

Segundo Felippa & Militello [1992], as deformações naturais podem ser transformadas para deformações no sistema cartesiano através da seguinte operação:

eTε ..2

. 1

13132

132

13

32322

322

32

21212

212

21

13

32

21−=

=

= def

xy

yy

xx

eee

cssccssccssc

εεε

(7.43)

sendo c21 = x21 / l21, s21 = y21 / l21, c32 = x32 / l32, s32 = y32 / l32, c13 = x13 / l13, s13 = y13 / l13. Já a relação inversa pode ser definida por:

( ) ( ) ( )

+++=

31

23

12

21332212312

23221131231

22113323123

2133212

2322131

2211323

2133212

2322131

2211323

2 ..4

1

εεε

lyxxylyxxylyxxylxxlxxlxxlyylyylyy

Ae (7.44)

ou de forma compacta εTe .def= (7.45) 7.6 ELEMENTO TRIANGULAR À FLEXÃO O componente à flexão do elemento de casca empregado no presente trabalho é baseado no elemento finito linear de placa triangular AQR desenvolvido por Militello [1991] e utilizado por Haugen [1994]. Os graus de liberdade nodais vi associados à flexão do elemento são representados por duas rotações θx e θy no plano que contém o elemento e uma translação para fora do plano w:

=

i

yi

xi

i

wθθ

v sendo yw

x ∂∂

=θ e xw

y ∂∂

−=θ (7.46)

7.6.1 Matriz de rigidez básica

Procedendo de forma análoga ao caso do elemento de membrana apresentado na seção anterior e ordenando os graus de liberdade em cada um dos nós conforme a equação acima, a matriz L (force-lumping matrix) que associa as tensões de contorno com as forças nodais, definida na Eq.(7.4) pode ser expressa por:

σLf .=e sendo

=

xy

yy

xx

τσσ

σ (7.47)

=

3

2

1

LLL

L

=

3

2

1

e

e

e

e

fff

f

=

zi

yi

xi

ei

fmm

f (7.48)

82

Conforme comentado anteriormente para o elemento de membrana, as tensões do contorno do elemento em um nó j serão somente função dos lados adjacentes ij e jk de modo que a matrix L total de cada um dos elementos pode ser subdividida em função da contribuição de cada um dos nós separadamente, conforme apresentando na Eq.(7.48), podendo ser expressa segundo Haugen [1994] por:

−−−=

kiki

kikij

xyyx

00

000.

21L (7.49)

onde os índices nodais (i,j,k) estão submetidos à permutações cíclicas de (1,2,3). A matriz de rigidez básica pode então ser computada por:

Tb A

t LCLK ...= (7.50)

sendo t e A a espessura e a área do elemento de casca, respectivamente. 7.6.2 Matriz de rigidez de alta ordem Na determinação da matriz de rigidez de alta ordem para o elemento ANDES triangular desenvolvido por Militello [1991] empregam interpolações diretas das curvaturas em coordenadas naturais. Como linhas de referência, Militello [1991] escolheu os três lados do elemento, que funcionam como vigas hermitianas. As leituras de deformações nodais expressas em função dos graus de liberdade v, podem ser definidas por: vFQvQg .. ∗== F (7.51) sendo 3

313

232

232

121

121

31 kkkkkkT =g (7.52) 333222111 www yxyxyx

T θθθθθθ=v (7.53)

−−−−−−

−−−−−−

−−−−−−

=

446000226446226000

226446000000446226

000226446226000446

FQ (7.54)

83

=

313131

232323

232323

121212

121212

313131

FFFFFF

FFFFFF

FFFFFF

F sendo

=

=

=

31

31

31

312

3131

23

23

23

232

2323

12

12

12

122

1212

1

1

1

ln

ln

l

ln

ln

l

ln

ln

l

yx

yx

yx

F

F

F

(7.55)

Pode-se observar através do vetor g, apresentado na Eq.(7.52), que para cada um dos nós de canto do elemento triangular, existem duas curvaturas em coordenadas naturais, associadas às direções dos lados adjacentes ao nó em questão, totalizando seis curvaturas naturais por elemento. Porém, são necessárias três curvaturas em coordenadas naturais para que seja feita a transformação para o sistema de coordenadas cartesiano. A terceira medida é então determinada utilizando-se a seguinte regra de projeção: assume-se que a curvatura natural kij varia linearmente ao longo do lado ij e seja constante ao longo de direções perpendiculares ao lado ij. Desta forma é atribuída ao nó k uma terceira curvatura kij de acordo com a projeção do nó k no lado ij. Esta hipótese pode ser expressa em forma matricial por: gAk .k= (7.56)

++++

++=

23132131

13221232

32113121

000000000000

ςλςςλςςλςςλς

ςλςςλς

kA (7.57)

ij

kiijT

kiij l

l..ss−=λ (7.58)

A parte desviadora da deformação pode então ser obtida subtraindo-se a parcela referente à deformação principal:

dAA

kkkd ∫−= AAA (7.59)

++++

++=

23132131

13221232

32113121

000000000000

dddd

dddd

dddd

ςλςςλςςλςςλς

ςλςςλς

sendo ζdi = ζi – 1/3 . As curvaturas desviadoras no sistema cartesiano podem então ser expressas por: vBvQATgATkTk ....... dkddefkddefdefd ==== (7.60)

84

Finalmente, a matriz de rigidez de alta ordem do elemento triangular submetido à esforços de flexão pode então ser calculada a partir das deformações desviadoras por:

∫=

Ad

Tdh dABCBK .. (7.61)

7.7 EXTENSÕES NÃO-LINEARES PARA O ELEMENTO TRIANGULAR A matriz de rigidez linear Ke do elemento de casca triangular pode ser obtida somando-se os efeitos dos esforços de membrana e de flexão definidos nas duas seções anteriores, de modo que o elemento finito linear de casca triangular está praticamente pronto para ser incorporado na formulação co-rotacional comentada no capítulo 6. A extensão para a análise não linear geométrica de um elemento consiste basicamente em definir um procedimento que alinhe a configuração co-rotacional C0n o mais próximo possível da configuração deformada Cn, lembrando que a diferença entre estas duas configurações define o vetor de deslocamentos deformacionais vd. Em função do procedimento escolhido para realizar esta tarefa, serão produzidas distintas sub-matrizes Gi, definidas anteriormente na Eq.(6.36), que no sistema de coordenadas locais podem ser expressas por: vGw ~.~~ δδ =r ⇔ [ ]321

~~~~ GGGG = (7.62) Existem várias técnicas na literatura que podem ser utilizadas para fazer este alinhamento entre as configurações co-rotacional e deformada. No presente trabalho foi escolhida a metodologia que alinha um dos lados do elemento triangular, conforme apresentado por Haugen [1994]. Este procedimento é semelhante ao utilizado por Nour-Omid & Rankin [1991], diferindo apenas pelo fato de que enquanto Nour-Omid & Rankin [1991] escolheram o lado 13 para definir o vetor unitário e2 e o nó 1 como origem do sistema de coordenadas, no presente trabalho alinha-se a direção do lado 12 com o eixo unitário e1 e emprega-se a média das coordenadas nodais (centróide) como a origem do sistema de coordenadas. Calculando-se a variação consistente do sistema de eixos locais, as submatrizes Gi no sistema local podem ser definidas como:

−=

0000/200000000000

.21~

12

32

32

1

lAyx

AG

=

0000/200000000000

.21~

12

13

13

2

lAyx

AG (7.63)

=

0000000000000000

.21~

32

32

3 yx

AG

85

CAPÍTULO 8 EXEMPLOS NUMÉRICOS

Neste capítulo será feita a análise não linear geométrica de algumas estruturas planas e espaciais, constituídas por elementos finitos de treliça ou de viga, segundo as formulações cinemáticas co-rotacionais (CR) apresentadas nos capítulos anteriores. Para fazer a análise foram desenvolvidos quatro programas computacionais em linguagem Fortran, conforme apresentados na Tabela (8.1), observando-se que para o caso específico dos programas envolvendo elementos de treliça, foi também implementado nos mesmos programas a formulação lagrangiana total (LT), conforme apresentado por Taylor [2001]. É importante ressaltar que não foi feita a implementação computacional referente à formulação co-rotacional para elementos de casca triangular, uma vez que pretende-se utilizar os programas já implementados na Universidade do Colorado / EUA, durante o doutorado sandwich. As trajetórias de equilíbrio foram obtidas utilizando o método de comprimento de arco cilíndrico combinado com o método de Newton-Raphson descritos por Crisfield [1991].

Tabela 8.1 – Programas computacionais Programa Elemento Dimensão Cinemática Capítulo

Truss_AL2D Treliça 2D CR e LT 2Truss_AL3D Treliça 3D CR e LT 2Beam_AL2D Viga 2D CR 3Beam_AL3D Viga 3D CR 4 e 5

8.1 TRELIÇAS PLANAS E ESPACIAIS Nesta seção serão apresentados quatro exemplos de estruturas constituídas por elementos finitos de treliças planas ou espaciais, sendo as mesmas analisadas segundo as formulações co-rotacional e lagrangiana total. 8.1.1 Estrutura articulada 2D não simétrica A primeira estrutura analisada corresponde a uma cobertura articulada plana, abatida e não simétrica, cujas características geométricas são apresentadas na Fig.(8.1). Esta estrutura foi estudada por Powel & Simons [1981] e apresenta 18 nós e 33 elementos de barra, com rigidez axial EA = 9.0 x 106 e esta submetida à 3 cargas nodais P de igual magnitude.

PP

P

35.0 45.0

9.445Nó 9

Figura 8.1 – Estrutura articulada 2D não simétrica

86

Na Fig.(8.2), são apresentadas as trajetórias de equilíbrio não linear para a direção vertical do nó central da estrutura, representado na Fig.(8.1) com o número 9, utilizando-se as deformações de engenharia e Biot. Pode-se constatar que quando se adota a mesma medida de deformação para as duas formulações (co-rotacional e lagrangiana total), as trajetórias de equilíbrio coincidem, porém no caso de deformações medianas as respostas são discrepantes quando se adotam distintas medidas de deformações em uma mesma formulação.

0500010000150002000025000300003500040000

-20-15-10-50Deslocamento Vertical - Nó 9

Car

ga

engenharia biot

Figura 8.2 – Trajetória não linear de equilíbrio para estrutura articulada 2D Conforme se pode ver na Fig.(8.2), a trajetória de equilíbrio apresenta pontos limites (tangente horizontal) e turning points (tangente vertical). Para detectar os pontos limites e turning points foi calculado o parâmetro de rigidez CST- current stiffness parameter, segundo Crisfield [1991], que se anula quando um destes tipos de pontos é alcançado. Nas Fig. (8.3a) e (8.3c), o CST é mostrado em função do deslocamento vertical do nó central e do número de passos de carga, respectivamente. Através destas figuras constata-se que o CST se anula 6 vezes, sendo 4 correspondentes a pontos limites e 2 referentes a turning points.

Figura 8.3 – CST e função sign para a estrutura articulada 2D

-0.25

0

0.25

0.5

0.75

1

1.25

-20-15-10-50

Deslocamento Vertical - Nó 9

CST

-1.5

-1

-0.5

0

0.5

1

1.5

-20-15-10-50

Deslocamento Vertical - Nó 9

Funç

ão S

ign

-0.4-0.2

00.20.40.60.8

11.2

0 10 20 30 40 50

Passos de Carga

CST

-1.5

-1

-0.5

0

0.5

1

1.5

0 10 20 30 40 50

Passos de Carga

Funç

ão S

ign

a) b)

c) d)

87

Uma outra forma de detectar pontos críticos (pontos limites e de bifurcação) é através da chamada função sinal (sign) que toma o valor +1 se o número de pivots negativos da matriz de rigidez tangente triangularizada é igual entre os passos n–1 e n e toma o valor –1 quando o número de pivots é diferente entre os passos n–1 e n. Nas Fig. (8.3b) e (8.3d), a função sinal é mostrada em função do deslocamento vertical do nó central e do número de passos de carga, respectivamente. Através destas figuras constata-se que a função sign troca de sinal 4 vezes, correspondendo aos 4 pontos limites. Vale ressaltar que esta função não detecta turning points uma vez que estes tipos de pontos não são considerados pontos críticos. 8.1.2 Estrutura articulada 3D em forma de torre Neste exemplo é estudada uma estrutura articulada espacial em forma de torre com imperfeição geométrica, conforme mostrado na Fig.(8.4). Esta estrutura foi também analisada por Onate [1986] e possui 13 nós e 36 elementos, estando submetida à um carregamento vertical no seu topo (nó 1). A estrutura possui 27 graus de liberdade livres e 12 graus de liberdade restringidos em 4 apoios de segundo gênero, representados por círculos hachuriados na vista em planta. Todos os elementos de barra possuem a mesma seção, com área A = 1.00 e módulo de elasticidade E = 0.1 x 105.

11

12

111119 19 19

15

15

15

250

250

250

P

x

y

x

z

1

Figura 8.4 – Estrutura articulada 3D em forma de torre / Vista em planta e corte

Na Fig.(8.5a) é mostrada a trajetória de equilíbrio não linear para o deslocamento horizontal do nó 1 na direção x para as formulações co-rotacional e lagrangiana total utilizando-se a deformação de engenharia, podendo-se observar uma excelente concordância nos resultados entre as duas formulações estudadas.

0

20

40

60

80

100

0 50 100 150 200 250

Deslocamento

Car

rega

men

to

____ Co-Rotacional ♦ Lagrangiana

-0.4-0.2

00.20.40.60.8

11.2

0 10 20 30 40 50 60 70 80 90 100

Carregamento

CST Pcrit = 90.75

(b)(a)

Figura 8.5 – Trajetória de Equilíbrio e curva CST x Carregamento

88

Já na Fig.(8.5b) é apresentada a variação do parâmetro CST (Current Stiffness Parameter) em relação ao carregamento da estrutura, podendo-se constatar, conforme o esperado, que o CST anula-se apenas uma vez, correspondendo a um ponto limite que pode ser observado em função da tangente horizontal na trajetória de equilíbrio apresentada na Fig.(8.5a). Através deste parâmetro foi possível detectar a carga crítica da estrutura: Pcrit = 90.75. 8.1.3 Estrutura articulada 3D em forma de cúpula / Papadrakakis [1987] Este exemplo é referente à uma estrutura articulada espacial em forma de cúpula analisada por Papadrakakis [1987], cujas características geométricas são apresentadas na Fig.(8.6), sendo r1 = 25.00 e r2 = 50.00. A estrutura possui 13 nós e 24 elementos de barra bi-articulados de igual seção transversal com rigidez axial EA = 1.0 x 105, estando submetida a um carregamento vertical aplicado sobre o nó central no topo da estrutura. A cúpula possui 21 graus de liberdade livres e 18 graus de liberdade restringidos em 6 apoios de segundo gênero, representados na vista em planta da Fig.(8.6) com círculos hachuriados.

xy

z

x

6.2162.00

Pr1

r2

Figura 8.6 – Estrutura em forma de cúpula / Papadrakakis [1987]

A trajetória de equilíbrio primaria do nó central na direção vertical (eixo z) é apresentada na Fig.(8.7a) para as formulações co-rotacional e lagrangiana total, utilizando-se a deformação de engenharia, podendo-se comprovar a excelente concordância entre os resultados. Observando-se a trajetória de equilíbrio, pode-se constatar que existem 8 pontos limites (tangente horizontal) e 2 turning points (tangente vertical), sendo estes pontos detectados pelas 10 anulações do parâmetro CST que é apresentado na Fig.(8.7b) em função dos deslocamentos do nó central da estrutura na direção vertical.

-1000

-500

0

500

1000

-20-15-10-50

Deslocamentos

Car

rega

men

to

____ Co-Rotacional ♦ Lagrangiana -2

-1

0

1

2

3

4

-20 -15 -10 -5 0

Deslocamentos

CST

(a) (b)

Figura 8.7 – Trajetória de equilíbrio e Curva CST x Deslocamento

89

8.1.4 Estrutura articulada 3D em forma de Cúpula / Choong & Hangai [1993]

A quarta estrutura estudada nesta seção corresponde a uma cúpula articulada espacial analisada por Choong & Hangai [1993], cujas características geométricas são apresentadas na Fig.(8.8). Esta estrutura apresenta 25 nós e 60 elementos de barra bi-articulados, com rigidez axial EA = 1.0x104 estando submetida à 7 carregamentos verticais P de igual magnitude aplicados respectivamente dos nós 1 ao 7. Esta estrutura apresenta 6 nós restringidos nas direções x, y e z, representados na vista em planta por círculos hachuriados.

P

x

y

x

z

1

2

3

4

5

6

7

2550

50

50

Vista em Planta Vista em Corte

PP

2.06.2212.3314.16

Figura 8.8 – Estrutura articulada 3D em forma de cúpula / Choong & Hangai [1993]

Na Fig.(8.9), são apresentadas as trajetórias de equilíbrio não linear para a direção vertical do nó central da estrutura, representado na Fig.(8.8) com o número 1, utilizando-se as deformações de Green-Lagrange e Almansi. Pode-se constatar que quando se adota a mesma medida de deformação para as duas formulações (co-rotacional e lagrangiana total), as trajetórias de equilíbrio coincidem, porém no caso de deformações medianas as respostas são discrepantes quando se adotam distintas medidas de deformações. Conforme se pode ver na Fig.(8.9), a trajetória de equilíbrio apresenta pontos limites (tangente horizontal) e turning points (tangente vertical). Nas Fig.(8.10a) e (8.10c), o parâmetro CST é mostrado em função do deslocamento vertical do nó central e do número de passos de carga, respectivamente. Através destas figuras constata-se que o CST se anula 10 vezes, sendo 8 correspondentes a pontos limites e 2 referentes a turning points.

-80-60-40-20020406080

-60-40-200

Deslocamento Vertical - Nó 1

Car

ga

green almansi

Figura 8.9 – Trajetória não linear de equilíbrio para estrutura

90

Conforme comentado anteriormente, uma outra forma de detectar pontos críticos (pontos limites e pontos de bifurcação) é através da chamada função sinal (sign) que toma o valor +1 se o número de pivots negativos da matriz de rigidez tangente triangularizada é igual entre os passos n–1 e n e toma o valor –1 quando o número de pivots é diferente entre os passos n–1 e n. Nas Fig.(8.10b) e (8.10d), a função sinal é mostrada em função do deslocamento vertical do nó central e do número de passos de carga, respectivamente. Através destas figuras constata-se que a função sign troca de sinal 28 vezes, correspondendo aos 8 pontos limites e as demais 20 referentes a pontos de bifurcação. Novamente, vale ressaltar que esta função não detecta turning points uma vez que estes tipos de pontos não são considerados pontos críticos e que no presente trabalho não serão obtidas as trajetórias secundárias, restringindo-nos somente às trajetórias primárias.

-1

-0.5

0

0.5

1

1.5

-60-50-40-30-20-100

Deslocamento Vertical - Nó 1

CST

-1.5

-1

-0.5

0

0.5

1

1.5

-60-50-40-30-20-100

Deslocamento Vertical - Nó 1

Funç

ão S

ign

-1

-0.5

0

0.5

1

1.5

0 40 80 120 160

Passos de Carga

CST

-1.5

-1

-0.5

0

0.5

1

1.5

0 40 80 120 160

Passos de Carga

Funç

ão S

ign

a) b)

c) d)

Figura 8.10 – CST e função sign para Cúpula / Choong & Hangai [1993]

8.2 PÓRTICOS PLANOS E ESPACIAIS Nesta seção serão apresentados seis exemplos clássicos de estruturas discretizadas por elementos finitos de vigas planas ou espaciais, sendo as mesmas analisadas segundo as descrições cinemáticas co-rotacionais apresentadas nos capítulos 3, 4 e 5. Vale a pena enfatizar que os problemas envolvendo estruturas no plano podem ser resolvidos utilizando-se o programa computacional Beam_AL2D (Pórticos Planos) ou o programa Beam_AL3D (Pórticos Espaciais), conforme comentado na Tabela (8.1). Para o caso das estruturas planas estudadas no presente trabalho, todas apresentaram os mesmos resultados quando analisadas com os programas plano e espacial, conforme esperado. Ao contrário dos programas com elementos de treliça, para o caso dos elementos de viga, não foi feita a implementação computacional da formulação lagrangiana total.

91

8.2.1 Viga em balanço com momento fletor aplicado na extremidade livre Este exemplo corresponde à uma viga em balanço sujeita à flexão pura através da aplicação de um momento fletor positivo (sentido anti-horário) na sua extremidade livre, conforme mostrado na Fig.(8.11). Esta estrutura, com comprimento total L = 1000.0 mm, foi estudada anteriormente por Faria [1998] segundo a formulação lagrangiana total, sendo discretizada com 10 elementos de viga de mesmo comprimento e seção transversal quadrada com as seguintes propriedades mecânicas: E =30.0x106 N/mm2, A =1.0mm2 e Iz =8.333x102.

M (+)

X

Y

L

Figura 8.11 – Viga em balanço com momento aplicada na extremidade livre

No trabalho de Faria [1998] a viga foi submetida ao momento fletor M=15708.λ, sendo λ

um fator de carga que variou de 0 até 1. Quando o momento fletor M atinge o seu valor máximo de 15708.0 (λ=1), a viga se torna um círculo de diâmetro d = L/2π =318.31mm, conforme pode ser visto na configuração deformada da estrutura na Fig. (8.12a). No caso do presente trabalho, a viga foi submetida a um momento fletor máximo de M=31416.0 que corresponde a um fator de carga λ=2, para o qual a viga completa uma segunda volta em torno do eixo z. Para facilitar a visualização, é mostrada na Fig.(8.12a) a configuração deformada da estrutura para a primeira volta completa em torno do eixo z (λ=0.00,...,1.00) ao passo que na Fig.(8.12b) é mostrada a deformada para a segunda volta completa em torno do eixo z (λ=1.00,...,2.00). Este exemplo constitui um teste severo para verificar a eficiência da formulação cinemática e da implementação computacional para o caso da análise não linear de estruturas envolvendo grandes rotações.

0

100

200

300

400

500

600

700

-400 -200 0 200 400 600 800 1000

Coordenadas X (mm)

Coo

rden

adas

Y (m

m)

λ=0.00

λ=0.25λ=0.50

λ=0.75

λ=1.00

-50

0

50

100

150

200

250

300

350

-200 -150 -100 -50 0 50 100 150 200

Coordenadas X (mm)

Coo

rden

adas

Y (m

m)

λ=1.00

λ=1.47

λ=2.0

(b)(a)

Figura 8.12 – Configurações deformadas da estrutura

As trajetórias de equilíbrio para a extremidade da viga nas direções vertical e horizontal são mostradas nas Fig.(8.13a) e (8.13b) respectivamente. Nestas trajetórias, os deslocamentos vertical (V) e horizontal (U) são parametrizados em função do comprimento total L da viga. Vale lembrar que estas trajetórias de equilíbrio são referentes ao carregamento total λ=2.

92

05000

100001500020000250003000035000

0 0.2 0.4 0.6 0.8

[ V / L ]

Mom

ento

(N.m

m)

← 1 Volta

← 2 Voltas

♦

♦

05000

100001500020000250003000035000

0 0.2 0.4 0.6 0.8 1 1.2 1.4

[ U / L ]

Mom

ento

(N.m

m)

1 Volta → ♦

♦2 Voltas →

(a) (b)

Figura 8.13 – Trajetórias de Equilíbrio da extremidade livre da viga 8.2.2 Viga em balanço com carregamento transversal aplicado na extremidade livre Neste exemplo, estudado por Cole [1990], uma viga horizontal engastada e livre de seção transversal quadrada foi submetida a um carregamento vertical P na sua extremidade livre, conforme mostrado na Fig.(8.14), juntamente com as suas propriedades geométricas e mecânicas. A medida que a carga P aumenta, a extremidade livre da viga sofre deslocamentos para baixo e em direção ao suporte, sendo algumas configurações deformadas apresentadas na Fig.(8.15), para um carregamento máximo de PL2 / EI = 10. A resposta da estrutura exibe um comportamento não-linear porque o momento da viga não aumenta na mesma proporção que o carregamento P, uma vez que ocorre uma diminuição progressiva do braço de alavanca com as deformações da viga.

P

L

L = 3.2 mÁrea = 0.01 m2

Ix = 16.667 x 106

Iy = Iz = 8.333 x 106

E = 210 x 109

Poisson (ν) = 0.0

Figura 8.14 – Viga engastada e livre, carregamento e propriedades geométricas e mecânicas

P

P

P

P

P

Figura 8.15 – Configurações deformadas da estrutura

93

As trajetórias de equilíbrio para os deslocamentos horizontais (U) e verticais (V) da extremidade livre são apresentadas nas Fig.(8.16a) e (8.16b) respectivamente, juntamente com os valores encontrados por Cole [1990], podendo-se observar uma excelente concordância entre os resultados.

0123456789

10

0 0.1 0.2 0.3 0.4 0.5 0.6

[ U / L ]

[ PL2 /

EI ]

----- Menin / Taylor ♦ Cole

0123456789

10

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

[ V / L ]

[ PL2 /

EI ]

----- Menin / Taylor ♦ Cole

(a) (b)

Figura 8.16 – Trajetórias de equilíbrio da extremidade livre (a) deslocamento horizontal e (b) deslocamento vertical

Esta estrutura foi discretizada utilizando-se 8 elementos, sendo os esforços resultantes

ao longo da viga apresentados na Fig.(8.17) para um nível de carga de PL2 / EI = 7. Através desta figura, pode-se observar um aumento da rigidez da estrutura com proporções crescentes da carga sendo resistidas pelas forças axiais, acompanhada por uma diminuição da resistência a flexão.

0.001.002.003.004.005.006.007.008.00

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

[X / L]

[PL2 /

EI]

FL2 / EI

VL2 / EI

ML / 2EI

Figura 8.17 – Esforços internos: força axial, cortante e momento fletor

8.2.3 Viga em balanço com carregamento axial aplicado na extremidade livre Uma viga em balanço é submetida a um carregamento axial P em sua extremidade livre conforme indicado na Fig.(8.18), juntamente com as suas propriedades mecânicas e geométricas. Quando P é maior que a carga crítica elástica, a viga sofre flambagem e as deformações aumentam consideravelmente, conforme apresentado na Fig.(8.19). Para iniciar a flambagem, aplica-se uma pequena carga de perturbação transversal Q = P/1000 na direção vertical sobre a extremidade livre. Este exemplo, também estudado por Cole [1990] tem por objetivo mostrar a performance da formulação não-linear na modelagem do comportamento pós-flambagem de estruturas.

94

P

Q = P/1000

L

L = 3.2 mÁrea = 0.01 m2

Ix = 16.667 x 106

Iy = Iz = 8.333 x 106

E = 210 x 109

Poisson (ν) = 0.0

Figura 8.18 – Viga em balanço submetida a carregamento axial na extremidade livre

P = 0.0

P = 4.35 x 105

P = 4.69 x 105

P = 7.73 x 105

P = 3.86 x 106

Figura 8.19 – Configuração deformada da estrutura Na Fig.(8.19) são apresentadas as configurações deformadas da estrutura para os passos de carga 3, 5, 9 e 13, aos quais estão associados carregamentos axiais (P) de 4.35x105, 4.69x105, 7.73x105 e 3.86x106, respectivamente.

A estrutura foi modelada utilizando-se 8 elementos de mesmo comprimento, sendo o carregamento aplicado em 13 passos de carga utilizando-se o método de comprimento de arco cilíndrico. As trajetórias de equilíbrio para o deslocamento vertical (V) e horizontal (U) da extremidade livre da viga, parametrizados em função do comprimento da viga (L) são apresentadas nas Fig.(8.20a) e (8.20b), respectivamente. Nestas mesmas figuras são também indicados os valores encontrados por Cole [1990], podendo-se observar uma excelente concordância nos resultados. Em função da trajetória de equilíbrio pode-se determinar a carga crítica de flambagem da estrutura que parametrizada em função das propriedades mecânicas e geométricas pode ser expressa por Pcrit.L2 / E.I = 2.5.

0.02.55.07.5

10.012.515.017.520.022.5

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

[ V / L ]

[ PL2 /

EI ]

0.02.55.07.5

10.012.515.017.520.022.5

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

[ U / L]

[ PL2 /

EI ]

____ Menin / Taylor♦ Cole

____ Menin / Taylor♦ Cole

(a) (b)

Figura 8.20 – Trajetórias de equilíbrio da extremidade livre

95

8.2.4 Pórtico de Lee Esta estrutura bastante conhecida na literatura, corresponde a um pórtico plano formado pela união entre uma viga e uma coluna, cujas características geométricas são apresentadas na Fig.(8.21). Este exemplo foi estudado por Schweizerhof & Wriggers [1986] sendo cada uma das duas barras discretizada através de 10 elementos finitos de viga de mesmo comprimento. As propriedades mecânicas são: A0 = 6 cm2, E = 720 KN/cm2 e I = 2.0 cm4.

24 96 cm

P

120 cm

Nó 13

Figura 8.21 – Pórtico de Lee

Na Fig. (8.22a) são apresentadas as trajetórias de equilíbrio não linear para as direções

vertical e horizontal do nó 13, sobre o qual é aplicado o carregamento; e de modo a facilitar o entendimento das trajetórias de equilíbrio e obter uma melhor visualização do comportamento pós-crítico da estrutura são apresentadas na Fig. (8.22b), as configurações deformadas do sistema para os carregamentos referentes aos pontos limites PL1 e PL2, conforme comentado à seguir.

-1.5-1

-0.50

0.51

1.52

2.5

0 10 20 30 40 50 60 70 80 90 100

Deslocamento (cm)

Car

rega

men

to (K

N)

0

40

80

120

-40 0 40 80 120 160

Coordenadas X (cm)

Coo

rden

adas

Y (c

m) PL1

PL2

Inicial

----- vertical----- horizontal

(a) (b)

Figura 8.22 – (a) Trajetórias de equilíbrio e (b) Configurações deformadas

Conforme se pode ver na Fig.(8.22a), a trajetória de equilíbrio do deslocamento vertical apresenta pontos limites (tangente horizontal) e turning points (tangente vertical). Para detectar os pontos limites e turning points foi calculado o parâmetro de rigidez CST – current stiffness parameter, segundo Crisfield [1991], que se anula quando um destes tipos de pontos é alcançado. Nas Fig. (8.23a) e (8.23b), o CST é mostrado em função do deslocamento vertical do nó 13 e do carregamento atuante. Através destas figuras constata-se que o CST se anula 4 vezes, sendo 2 correspondentes a pontos limites e 2 referentes a turning points. Na Fig.(8.23b), verifica-se que o primeiro ponto limite (PL1) ocorre para um carregamento de 1.857 KN e o segundo (PL2) para -0.950 KN. O carregamento referente à PL1 é muito próximo ao valor apresentado por Schweizerhof & Wriggers [1986] que é de 1.851 KN, que entretanto não forneceu o valor referente à PL2.

96

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

0 10 20 30 40 50 60 70 80 90 100

Deslocamento Vertical (cm)

CST

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

-1 -0.5 0 0.5 1 1.5 2

Carregamento (KN)

CST

(a) (b)

-1.5

-1

-0.5

0

0.5

1

1.5

-90-80-70-60-50-40-30-20-100

Deslocamento Vertical (cm)

-1.5

-1

-0.5

0

0.5

1

1.5

0 10 20 30 40 50 60 70

Passos de Carga

----- Função sign----- Neg. Pivots

----- Função sign----- Neg. Pivots

(c) (d)

Figura 8.23 – CST, função sign e número de pivos negativos

Uma outra forma de detectar pontos críticos (pontos limites e pontos de bifurcação) é através da chamada função sinal (sign) que toma o valor +1 se o número de pivots negativos da matriz de rigidez triangularizada é igual entre os passos n–1 e n e toma o valor –1 quando o numero de pivots é diferente. Nas Fig. (8.23c) e (8.23d), a função sinal e o número de pivots negativos são mostrados em função do deslocamento vertical e do número de passos de carga. Através destas figuras constata-se que a função sign troca de sinal 2 vezes, correspondendo aos 2 pontos limites. Vale ressaltar que esta função não detecta turning points. 8.2.5 Arcos circulares de grande altura Esta estrutura analisada por Faria [1998] corresponde a um arco circular de grande altura submetido a ação de uma força vertical no ápice conforme apresentado na Fig.(8.24) com R = 100 cm, assumindo-se duas condições de contorno distintas (arco biengastado ou arco rotulado-engastado), devido as quais o arco apresenta diferentes comportamentos pós-críticos. A estrutura foi discretizada utilizando-se 20 elementos de mesmo comprimento, cujas propriedades mecânicas são expressas por: EA0 = 106 N e EI= 106 N.cm2.

P P

R R145o

R R145o

(a) (b)

Figura 8.24 – (a) Arco biengastado e (b) Arco rotulado-engastado

97

Na Fig. (8.25a), são apresentadas as trajetórias de equilíbrio não linear para a direção vertical do nó central da estrutura, para os arcos biengastado e rotulado-engastado, podendo se verificar que para o caso do arco biengastado a trajetória de equilíbrio apresenta apenas pontos limites (tangente horizontal) ao passo que para o arco rotulado-engastado ocorrem pontos limites e turning points (tangente vertical), entre os quais ocorre uma acentuada queda da rigidez da estrutura na fase pós-crítica.

-200

0

200

400

600

800

1000

-250-200-150-100-500

Deslococamento Vertical (cm)C

arre

gam

ento

(N)

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

-200 0 200 400 600 800 1000

Carregamento (N)

CST

BiengastadoRot - engastado

BiengastadoRot - engastado

(a) (b)

Figura 8.25 – (a) Trajetórias de equilíbrio e (b) Curva CST x carregamento

Conforme comentado nos exemplos anteriores, para detectar os pontos limites e turning points foi calculado o parâmetro de rigidez CST – current stiffness parameter, que se anula quando um destes tipos de pontos é alcançado. Na Fig.(8.25b), o CST é mostrado em função do carregamento no topo da estrutura, podendo-se constatar que o CST se anula 2 vezes para o arco biengastado, correspondendo aos dois pontos limites PL1 e PL2 ao passo que para o arco rotulado-engastado, o CST se anula 4 vezes, correspondendo aos 2 pontos limites e 2 turning points. Com o auxílio da Fig.(8.25b) é possível calcular os valores dos carregamentos nos quais ocorrem os pontos limites PL1 e PL2 para os dois arcos sendo estes valores apresentados na Tabela (8.2):

Tabela 8.2 – Carregamentos em Newton correspondentes aos pontos limites

Tipo de Arco PL 1 PL 2Biengastado 941.65 432.65

Rotulado-engastado 910.22 -81.60

Para facilitar o entendimento das trajetórias de equilíbrio e obter uma melhor visualização das diferenças no comportamento pós-crítico da estrutura ao se variar as condições de contorno, são apresentadas nas Fig.(8.26a) e (8.26b) as configurações deformadas do sistema para os carregamentos referentes aos pontos limites PL1 e PL2, conforme indicados anteriormente na Tabela (8.2). Ao se observar as deformadas das estruturas, verifica-se que no caso do arco biengastado a estrutura sofre deformações simétricas ao passo que o arco rotulado-engastado, sofre uma grande diminuição da rigidez acompanhada de deformações não simétricas. Analisando-se as trajetórias de equilíbrio na Fig.(8.25a) e os carregamentos correspondentes aos pontos limites na Tabela (8.2), pode-se ver claramente que o fato de substituir um engaste por uma rótula no apoio esquerdo leva a uma grande diminuição da rigidez na fase pré-crítica, ou seja, anterior ao ponto limite PL1, conduzindo à uma diminuição da capacidade resistente da estrutura em cerca de 31.43 N e um grande aumento nos deslocamentos, acompanhada por uma forte não linearidade na fase pós-crítica entre os 2 turning points, além de deformações não simétricas desde o início do carregamento.

98

-80

-40

0

40

80

120

-120 -60 0 60 120

Coordenadas X (cm)

Coo

rden

adas

Y (c

m)

-80

-40

0

40

80

120

-180 -120 -60 0 60 120

Coordenadas X (cm)

Coo

rden

adas

Y (c

m)

inicial

PL1

PL2

inicial

PL1

PL2

(a) (b)

Figura 8.26 – Configurações deformadas (a) Arco biengastado e (b) Arco rotulado-engastado

8.2.6 Viga engastada espacial com curvatura de 45 graus no plano Neste exemplo, uma viga engastada e livre inicialmente curva com seção transversal quadrada é submetida a um carregamento vertical P = 600 aplicado na sua extremidade livre. A estrutura foi discretizada utilizando-se 8 elementos de igual comprimento, sendo a configuração inicial e as propriedades mecânicas e geométricas mostradas na Fig.(8.27).

X

Y

Z

45o

R = 100

P = 0

P = 300

P = 600

Área = 1.0 x 1.0Ix = 1.6667 x 10-1

Iy = Iz = 8.3333 x 10-2

E = 1.0 x 107

Poisson (ν) = 0.0

Figura 8.27 – Configuração inicial, deformada e propriedades mecânicas e geométricas Neste exemplo, a estrutura foi submetida a um carregamento máximo de P = 600, conforme comentado anteriormente, utilizando-se 21 passos de carga no método de comprimento de arco cilíndrico definido em Crisfield [1991]. Na tabela (8.3) são apresentados os deslocamentos da extremidade livre da viga para carregamentos P de 300, 450 e 600, nas direções X, Y e Z. Nesta tabela são também apresentados valores encontrados por outros autores, tais como: Simo & Vu-Quoc [1986] e Cole [1990], podendo-se observar uma grande concordância nos resultados entre todos os autores.

99

Tabela 8.3 – Deslocamentos da extremidade livre da viga

Autor X Y Z X Y Z X Y ZSimo / Vu-Quoc -11.87 40.08 -6.96 -18.39 48.39 -10.67 -23.48 53.37 -13.50

Cole -11.95 40.25 -7.01 -18.50 48.58 -10.73 -23.61 53.58 -13.55Menin / Taylor -11.90 40.19 -7.02 -18.43 48.53 -10.74 -23.52 53.50 -13.56

P = 300 P = 450 P = 600

A seguir, serão mostradas as trajetórias de equilíbrio primárias para a extremidade livre da viga nas direções X, Y e Z. Na Fig.(8.28a) é apresentada a trajetória de equilíbrio para o deslocamento vertical, ou seja, na direção Y ao passo que na Fig.(8.28b) são mostradas as trajetórias de equilíbrio nas direções horizontais X e Z.

0100200300400500600700

0 10 20 30 40 50 60

Deslocamento Vertical - Extremidade Livre

Car

rega

men

to -

P

0100200300400500600700

-25 -20 -15 -10 -5 0

Deslocamentos Horizontais - Extremidade Livre

Car

rega

men

to -

P

Direção X

Direção ZDireção Y

(a) (b)

Figura 8.28 – Trajetória de equilíbrio (a) direção Y e (b) direções X e Z

8.3 ANÁLISE DE RESULTADOS Em relação aos exemplos analisados envolvendo estruturas constituídas por elementos de treliças, pode-se constatar que no caso de deformações infinitesimais, as configurações inicial e deformada se confundem e portanto, se obtém unicidade na resposta estrutural, independente do tipo de deformação utilizado, porém para o regime de deformações finitas, ao se adotar uma relação linear entre o par conjugado de tensão e deformação, não se obtém unicidade na resposta, implicando na utilização de materiais diferentes. Para se obter unicidade na resposta neste caso, seria necessário fazer uso de transformações tensoriais que fazem o mapeamento do tensor constitutivo entre as configurações inicial e atual. Pode também ser observado que para o caso de se adotar o mesmo tipo de deformação para as formulações co-rotacional e lagrangiana total, as trajetórias de equilíbrio são coincidentes. O parâmetro de rigidez CST e a função sinal (sign) permitiram detectar e classificar com grande precisão as singularidades encontradas na matriz de rigidez triangularizada bem como possibilitaram a obtenção das cargas críticas com boa aproximação em relação aos valores encontrados na literatura técnica especializada. Com a crescente utilização de estruturas cada vez mais esbeltas, aumenta a possibilidade de ocorrência de fenômenos de instabilidade de equilíbrio tanto na fase pré-crítica quanto na fase posterior a perda de estabilidade (fase pós-crítica), estando a perda ou não da capacidade portante intimamente ligada com a natureza da instabilidade de equilíbrio, de modo que no presente trabalho procurou-se obter uma melhor avaliação do desempenho da capacidade portante das estruturas nas fases pré e pós-crítica. O método de comprimento de arco cilíndrico utilizado para obter as trajetórias de equilíbrio possibilitou a simulação numérica com boa aproximação e de maneira geral todos os resultados encontrados apresentaram grande concordância em relação aos encontrados por outros autores.

100

REFERÊNCIAS BIBLIOGRÁFICAS

• Alvin, K., de la Fuente, H.M., Haugen, B., Felippa, C.A., 2002, “Membrane triangles with corner drilling freedoms – I. The EFF element”, Finite Elements in Analysis and Design, Vol.12, pp163-187.

• Argyris, J.H.,1982, “An excursion into large rotations”, Computer Meth. Appl. Mech.

Engrg., Vol. 32, pp. 85-155.

• Bathe, K.J. & Wilson, E., 1976, “Numerical Methods in Finite Element Analysis”, Prentice Hall Inc., New Jersey.

• Belytschko, T. & Hsieh, B.J., 1973, “Non-linear transient finite element analysis with

convected co-coordinates”, Int. J. Numer. Methods in Engineering, Vol.7, pp 255-271.

• Bergan, P.G. & Horrigmoe, G., 1976, “Incremental variational principles and finite element models for nonlinear problems”, Computer Methods in Applied Mech. Engineering, Vol. 7, pp 201-217.

• Bergan P.G. & Nygard, M.K., 1989, “Nonlinear shell analysis using Free Formulation

finite elements”, Finite Element Methods for Nonlinear Problems, Springer Verlag, Berlin, pp 317-338.

• Brebbia, C.A. & Ferrante, J., 1978, “Computational Methods for the Solution of

Engineering Problems”, 3rd Edition, Pentech Press, London.

• Cardona, A., 1989, “An integrated approach to mechanism analysis”, Ph.D thesis, University of Liege, Belgium.

• Choong, K.K. & Hangai, Y., 1993, “Review on methods of bifurcation analysis for

geometrically nonlinear structures”, Bulletin of IASS, Vol. 34, No 112, pp. 13-149.

• Clough, R.W. & Penzien, J., 1993, “Dynamics of Structures”, 2nd Edition, McGraw – Hill Incorporation.

• Cole, G., 1990, “Consistent co-rotational formulation for geometrically nonlinear

beam elements with special reference to large rotations”, Ph.D thesis – Kingston Polytechnic, UK.

• Cortivo, N. & Taylor, W.M.S., 2002, “Análise não-linear de pórticos planos,

utilizando uma formulação co-rotacional”, XXX Jornadas Sul-Americanas de Engenharia Estrutural, Brasília/DF, Brasil.

• Crisfield, M. A., 1990, “A consistent co-rotational formulation for non-linear three-

dimensional beam elements”, Computer Methods Appl. Mech. Engineering, Vol. 81, pp131-150.

101

• Crisfield, M. A., 1991, “Non-linear finite element analysis of solids and structures”, Volume 1: Essential, J. Wiley, Chichester.

• Crisfield, M. A., 1997, “Non-linear finite element analysis of solids and structures”,

Volume 2: Advanced Topics, J. Wiley, Chichester.

• Crisfield, M. A. & Shi, J., 1994, “A co-rotational element/time-integration strategy for non-linear dynamics”, Int. Journal for Numerical Methods in Engineering, Vol. 37, pp 1897-1913.

• Crisfield, M. A. & Moita, G. F., 1996, “A unified co-rotational framework for solids,

shells and beams”, International Journal of Solids and Structures, Vol. 33, No 20-22, pp. 2969-2992.

• Faria, H.P., 1998, “Análise não-linear de instabilidade elástica de pórticos planos”,

Dissertação de Mestrado, UnB / DF, Brasil.

• Felippa, C.A. & Militello, C., 1992, “Membrane triangles with corner drilling freedoms – II. The ANDES element”, Finite Elements in Analysis and Design, Vol.12, pp 189-201.

• Felippa, C.A. , 2000, “A systematic approach to the element independent co-rotational

dynamics of finite elements”, Report CU-CAS-00-03, Center for Aerospace Structures University of Colorado, Boulder/USA.

• Felippa, C.A., 2001, “Non-linear finite element methods”, Lecture notes for the course

non-linear finite element methods, Center for Aerospace Structures, University of Colorado, Boulder/USA.

• Fraeijs de Verbeke, B.M., 1976, “The dynamics of flexible bodies”, Int. J. Engineering

Science, Pergamon Press, 895-913.

• Gere, J. & Weaver, W., 1981, “Análise de Estruturas Reticuladas”, Editora Guanabara, Rio de Janeiro.

• Harrison, H.B., 1973, “Computer methods in structural analysis”, Prentice Hall Inc.,

New Jersey.

• Haugen, B., 1994, “Buckling and Stability Problems for Thin Shell Structures Using High Performance Finite Elements”, Ph.D Thesis, University of Colorado, USA.

• Horrigmoe, G. & Bergan, P.G., 1978, “Instability analysis of free-form shells by flat

finite elements”, Comput. Methods Applied Mech. Engrg., Vol. 16, pp 11-35.

• Luenberger, D.G. , 1984, “Linear and Nonlinear Programming”, 2nd Edition, Addison Wesley Publishing Company, Massachusetts.

102

• Menin, R.C.G. & Taylor, W.M.S, 2003a, “Resposta pós-crítica de sistemas articulados com diferentes deformações utilizando uma formulação co-rotacional”, XXIV Cilamce, Ouro Preto / MG, Brasil.

• Menin, R.C.G. & Taylor, W.M.S, 2003b, “Resposta pós-crítica de pórticos planos

discretizados com elementos de viga de Euler-Bernoulli utilizando uma formulação co-rotacional”, XXIV Cilamce, Ouro Preto / MG, Brasil.

• Menin, R.C.G. & Taylor, W.M.S, 2004, “Análise não-linear geométrica de pórticos

espaciais utilizando uma formulação co-rotacional”, XXXI Jornadas Sul Americanas de Engenharia Estrutural, Mendoza, Argentina.

• Militello, C. & Felippa, C.A., 1990, “Variational formulation of high performance

finite elements: Parametrized variational principles”, Computers & Structures, Vol. 36, No 1, pp 1-11.

• Militello, C., 1991, “Application of parametrized variational principles to the finite

element method”, Ph.D thesis, Department of Aerospace Engineering Science, University of Colorado, Boulder/USA.

• Monteiro, F.A.C., 2004, “Estudo de uma Formulação Co-Rotacional Geral: Aplicação

a Pórticos Espaciais”, Dissertação de Mestrado, ITA / SP.

• Neto, E.L. & Yshii, Y., 2002, “Formulação co-rotacional para pórticos planos”, XXX Jornadas Sul-Americanas de Engenharia Estrutural, Brasília/DF, Brasil.

• Newmark, N.M., 1959, “A Method of Computation for Structural Dynamics”, Journal

of the Engineering Division – ASCE Vol. 85 EM3 (July): 67-94.

• Nour-Omid, B. & Rankin, C.C., 1991, “Finite rotation analysis and consistent linearization using projectors”, Comp. Meth. in Applied Mechanics and Engineering, Vol. 93, pp. 353-384.

• Nygard, M.K., 1986, “The free formulation for non-linear finite elements with

application to shells”, Dr. Ing. Thesis, Div. of Structural Mechanics, Norwegian Institute of Technology, Trondhein, Norway.

• Onate, E., 1986, “Una formulación incremental para problemas de no linealidad

geométrica por el metodo de los elementos finitos”, UPC / Barcelona, Espanha.

• Pacoste, C. & Eriksson, A., 1996, “Beam element in instability problems”, Comp. Meth. in Applied Mechanics and Engineering, Vol. 144, pp. 163-197.

• Pacoste, C., 1998, “Co-rotational flat facet triangular elements for shell instability

analyses”, Comp. Meth. in Applied Mechanics and Engineering, Vol 156, pp. 75-110.

• Papadrakakis, M., 1987, “Analysis methods for spacial structures”, in G.de Roeck, A.Quiroga, M. van Laethen and E. Back, eds., “Shell and Spacial Structures”, Computational Aspects, Springer, Berlin, pp. 121-148.

103

• Park, K.C. & Stanley, G.M., 1986, “A curved C0 shell element based on assumed natural-coordinate strains”, J. Appl. Mech. Engrg., Vol. 53, pp 278-290.

• Paz, M., 1997, “Structural Dynamics – Theory and Computation”, 4th Edition. Kluwer

Academic Publishers. Boston.

• Powel, G. & Simons, J., 1981, “Improved iteration strategy for nonlinear structures”, International Journal for Numerical Methods in Engineering, Vol. 17, pp. 1455-1467.

• Rankin, C.C. & Brogan, F.A., 1986, “An element independent corotational procedure

for the treatment of large rotations”, ASME J. Pressure Vessel Technology, Vol. 108, pp 165-174.

• Rankin, C.C. & Nour-Omid, B., 1988, “The use of projectors to improve finite

element performance”, Computers & Structures, Vol. 30, pp 257-267.

• Rao, C.R. & Mitra, S.K., 1971, “Generalized inverse of matrices and its applications”, John Wiley and Sons Inc., New York.

• Schweizerhof, K. & Wriggers, P., 1986, “Consistent linearization for path following

methods in nonlinear F.E. Analysis”, Computer Meth. Appl. Mech. Engrg., Vol. 59, pp 261-279, North-Holland.

• Simo, J.C. & Vu-Quoc, L, 1986, “A three-dimensional finite strain rod model. Part 2:

Computational Aspects”, Computer Meth. Appl. Mech. Engrg., Vol. 58, pp 79-116.

• Taylor, W.M.S, 2001, “Análisis no lineal elástico de estructuras de barras articuladas com diferentes medidas de deformacion”, XXII Cilamce, Campinas / SP, Brasil.

• Timoshenko, S. & Woinowsky-Krieger,S., 1959, “Theory of plates and shells”,

Second Edition, McGraw-Hill Book Company, New York.

• Wempner, G.A., 1969, “Finite elements, finite rotations and small strains of flexible shells”, Int. J. Solids and Structures, Vol. 5, pp 117-153.