Método dos Elementos Finitos

MÉTODO DOS ELEMENTOS FINITOS Álvaro F M Azevedo httpwwwfeupptalvaro Portugal Faculdade de Engenharia da Universidade do Porto 1ª Edição Abril 2003 Método dos Elementos Finitos Álvaro F M Azevedo httpwwwfeupptalvaro Faculdade de Engenharia da Universidade do Porto Portugal 1ª Edição Abril 2003 iii PREFÁCIO O Método dos Elementos Finitos MEF apresenta actualmente um nível de desenvolvimento que permite a sua utilização pela generalidade dos projectistas de estruturas Enquanto que no passado muitos dos utilizadores do MEF estavam também envolvidos na respectiva programação em computador verificase hoje em dia que a quase totalidade dos projectistas de estruturas apenas se preocupa com a utilização do correspondente software e com a interpretação dos resultados obtidos Devido à grande complexidade associada ao desenvolvimento de modernos programas de computador dispondo de uma interface gráfica intuitiva o desenvolvimento de software tem sido cada vez mais restringido às empresas especializadas Por este motivo o utilizador programador quase desapareceu dando lugar ao mero utilizador Perante um problema de análise de estruturas e dispondo de um software intuitivo é perfeitamente acessível a um projectista a obtenção de resultados credíveis mesmo quando não tem acesso à fonte do código computacional ou quando desconhece as características do modelo que está a utilizar Será então necessário exigir que um estudante de Engenharia atribua parte do seu tempo à aprendizagem de formulações e metodologias que na vida profissional vai certamente ignorar Antecedendo a resposta a esta questão apresentamse algumas considerações Para que possa dar resposta em tempo útil à necessidade de justificação da segurança de uma estrutura um projectista que não conheça as técnicas correspondentes à formulação do MEF será tentado pela simples utilização de um qualquer software de cálculo Uma vez que não tem acesso aos modelos que estão programados nem tem bases para a sua compreensão procederá à utilização do software de acordo com o treino que recebeu ou com base em sucessivas improvisações A tentação para aceitar os resultados provenientes do programa é grande quaisquer que sejam esses resultados uma vez que considera que o software escolhido tem elevada qualidade Os potenciais perigos de uma utilização nestas condições são a não percepção de eventuais erros na introdução dos dados a ausência de correspondência entre o modelo seleccionado e a estrutura que está a ser analisada o facto de serem desprezadas importantes condicionantes etc Na ausência de uma comparação dos resultados provenientes do MEF com os oriundos de outros modelos existe o sério risco de a segurança de uma estrutura ser justificada com Método dos Elementos Finitos Prefácio iv base em cálculos completamente inadequados Este facto tem sido confirmado pelo elevado número de acidentes em estruturas acabadas de construir bem como pela grande quantidade de reparações que tem sido necessário efectuar em construções recentes A transmissão aos alunos dos fundamentos do MEF e também de uma introdução à correspondente programação em computador constituem certamente factores que conduzirão os futuros projectistas a uma utilização mais segura dos softwares de análise de estruturas Existe uma outra motivação para continuar a ser necessário ensinar as bases teóricas do MEF que consiste no facto de ser fundamental preparar hoje os inovadores de amanhã Uma vez que as ferramentas relativas à aplicação do MEF se encontram intimamente ligadas ao mundo da informática e uma vez que este apresenta uma constante e rápida evolução é garantido que dentro de alguns anos será necessário adaptar as técnicas de análise de estruturas às plataformas de computação que nessa altura existirem Se a actual base de conhecimentos ficar limitada a um reduzido número de pessoas certamente que será difícil encontrar no futuro investigadores que garantam o progresso da ciência Por todos estes motivos se conclui ser fundamental prosseguir com o ensino das técnicas em que se baseia a generalidade dos programas de elementos finitos A principal motivação para a escrita desta publicação foi a de organizar de um modo coerente algumas das formulações em que se baseou o desenvolvimento do programa FEMIX 40 Apesar de já existirem versões anteriores a actual versão do programa foi totalmente rescrita de modo a ser possível explorar uma muito mais versátil estruturação do código computacional Esperase com este empreendimento produzir um software em que seja simples desenvolver e testar novas formulações Por último desejo agradecer às pessoas que se têm empenhado no desenvolvimento do projecto FEMIX e que muito contribuíram para que todos os conceitos aqui expostos apresentem uma indispensável clareza e coerência Em particular um agradecimento àquele que esteve presente desde o início Joaquim Barros bem como aos entusiastas mais recentes José Sena Cruz e António Ventura Gouveia Agradeço também ao Luís Brás o trabalho que teve na preparação do modelo da ponte que figura na capa Álvaro F M Azevedo Abril 2003 v ÍNDICE 1 INTRODUÇÃO1 11 Tipo de análise2 12 Fundamentos do MEF4 13 Breve história do MEF5 14 Exemplo de aplicação do MEF6 2 TRANSFORMAÇÃO LINEAR DE COORDENADAS 13 21 Simbologia13 22 Caso geral14 23 Caso particular com S e S coincidentes18 24 Matriz de transformação de uma barra rectilínea no espaço19 25 Considerações finais 27 3 MÉTODO DOS DESLOCAMENTOS EM TRELIÇAS E PÓRTICOS29 31 Simbologia29 32 Referenciais31 33 Graus de liberdade 32 34 Matriz de transformação 34 35 Matriz de rigidez e vector solicitação 35 36 Assemblagem da matriz de rigidez global e do vector solicitação37 37 Introdução das condições de apoio 41 38 Faseamento da análise de um pórtico 3D 44 39 Matriz de rigidez de uma barra de treliça 3D no referencial local45 310 Matriz de rigidez de uma barra de pórtico 3D no referencial local 46 311 Considerações finais 47 4 ELEMENTOS FINITOS UNIDIMENSIONAIS 49 41 Simbologia49 42 Funções interpoladoras ou funções de forma50 43 Campo de deformações54 Método dos Elementos Finitos Índice vi 44 Princípio dos trabalhos virtuais56 45 Matriz de rigidez e vector solicitação 57 46 Elemento finito unidimensional com três nós60 47 Elemento finito unidimensional com substituição de variável 64 48 Considerações finais 70 5 QUADRATURA DE GAUSS73 51 Simbologia73 52 Integração de uma função polinomial73 53 Integrais múltiplos 79 54 Considerações finais 81 6 ESTADO PLANO DE TENSÃO 83 61 Simbologia83 62 Funções interpoladoras ou funções de forma85 63 Campo de deformações90 64 Princípio dos trabalhos virtuais92 65 Matriz de rigidez e vector solicitação 92 651 Cálculo de um elemento da matriz de rigidez95 652 Cálculo do vector solicitação correspondente a uma carga distribuída 97 66 Caso geral com substituição de variáveis 99 67 Algoritmo de cálculo da matriz de rigidez de um elemento isoparamétrico 108 68 Cálculo das tensões e deformações finais112 69 Considerações finais 113 7 FUNÇÕES INTERPOLADORAS 115 71 Simbologia115 72 Caso unidimensional116 73 Caso bidimensional118 74 Procedimento genérico para determinar as funções de forma 121 75 Elementos bidimensionais famílias Lagrangeana e serendipity 126 76 Propriedades das funções interpoladoras130 77 Interpolação Hermitiana132 Método dos Elementos Finitos Índice vii 78 Considerações finais 142 8 ASSEMBLAGEM DE ELEMENTOS FINITOS145 81 Simbologia145 82 Assemblagem da matriz de rigidez global e do vector solicitação146 83 Considerações finais 152 9 FORÇAS NODAIS EQUIVALENTES153 91 Simbologia153 92 Expressões genéricas das forças nodais equivalentes155 93 Força concentrada num ponto interior160 94 Carga distribuída por unidade de comprimento163 95 Carga distribuída por unidade de superfície 170 96 Carga distribuída por unidade de volume170 97 Considerações finais 172 10 SÓLIDOS ESTADO PLANO DE DEFORMAÇÃO E AXISSIMETRIA 175 101 Simbologia175 102 Elementos sólidos tridimensionais bricks 176 103 Estado plano de deformação 184 104 Estado axissimétrico 187 105 Considerações finais 192 11 FLEXÃO DE VIGAS 193 111 Simbologia193 112 Flexão composta plana194 113 Considerações finais 200 12 VIGA DE EULERBERNOULLI 203 121 Simbologia203 122 Viga de dois nós sem substituição de variável204 123 Viga de três nós sem substituição de variável 212 124 Viga de dois nós com substituição de variável212 Método dos Elementos Finitos Índice viii 125 Considerações finais 220 13 VIGA DE TIMOSHENKO223 131 Simbologia223 132 Viga de dois nós com substituição de variável224 133 Considerações finais 237 ANEXO A UTILIZAÇÃO DO PROGRAMA FEMIX 31 239 A1 Instalação 239 A2 Preparação dos dados240 A3 Execução do programa245 A4 Visualização gráfica246 A5 Considerações finais 248 1 CAPÍTULO 1 INTRODUÇÃO No âmbito da Engenharia de Estruturas o Método dos Elementos Finitos MEF tem como objectivo a determinação do estado de tensão e de deformação de um sólido de geometria arbitrária sujeito a acções exteriores Este tipo de cálculo tem a designação genérica de análise de estruturas e surge por exemplo no estudo de edifícios pontes barragens etc Quando existe a necessidade de projectar uma estrutura é habitual procederse a uma sucessão de análises e modificações das suas características com o objectivo de se alcançar uma solução satisfatória quer em termos económicos quer na verificação dos prérequisitos funcionais e regulamentares As técnicas descritas nesta publicação apenas correspondem à fase de análise do comportamento de uma estrutura cuja geometria materiais e acções são a priori conhecidos Nos cursos de Engenharia Civil e de Engenharia Mecânica é tradicional começarse por ensinar a análise de estruturas limitada às vigas pórticos treliças e grelhas As estruturas deste tipo recebem a designação de reticuladas por serem constituídas por barras prismáticas cuja secção transversal apresenta dimensões muito inferiores ao comprimento do seu eixo As estruturas não reticuladas são em geral estudadas como meios contínuos eg paredes lajes cascas sólidos Nas estruturas reticuladas surgem já muitos conceitos que são comuns à generalidade das estruturas tais como o de equilíbrio compatibilidade tensão deformação relação entre tensão e deformação etc No âmbito das estruturas reticuladas tornase particularmente simples explicar o método das forças e o método dos deslocamentos bem como outras técnicas que em geral são difíceis de estender aos meios contínuos Antes do aparecimento do MEF a análise dos meios contínuos era efectuada por resolução directa dos sistemas de equações de derivadas parciais que regem o fenómeno tendo em consideração as necessárias condições fronteira Para facilitar a aplicação desta técnica a problemas não elementares era comum recorrer a séries de Fourier 11 Devido à sua complexidade estes procedimentos só eram aplicáveis a meios contínuos homogéneos e de geometria simples Para tentar ultrapassar algumas destas limitações era frequente a substituição de derivadas exactas por derivadas Introdução Álvaro F M Azevedo 2 aproximadas calculadas com base em grelhas de pontos Da aplicação desta técnica resulta o método das diferenças finitas que antes do aparecimento dos computadores apresentava o inconveniente de requerer a resolução de grandes sistemas de equações lineares Para evitar este inconveniente foram propostos diversos métodos de relaxação baseados na sucessiva diminuição de um conjunto de resíduos 11 Devido à morosidade associada à aplicação de qualquer um destes métodos tornavase muito atractiva a substituição do problema real por outro semelhante de modo a se poder recorrer a resultados publicados em tabelas ou ábacos Com o grande desenvolvimento que o MEF teve na década de 60 12 e com a banalização do recurso ao computador passou a ser prática corrente a análise de estruturas de geometria arbitrária constituídas por múltiplos materiais e sujeitas a qualquer tipo de carregamento Este avanço é tão significativo que os outros métodos atrás referidos deixaram praticamente de ser utilizados Actualmente o seu interesse restringese ao de fornecer soluções teóricas de problemas simples para validar métodos aproximados A formulação do MEF pode ser baseada no método dos deslocamentos em modelos de equilíbrio ou em métodos híbridos e mistos 13 De todos estes métodos aquele que apresenta uma maior simplicidade e consequentemente uma maior versatilidade é o método dos deslocamentos sendo este o único que é abordado nesta publicação Associado ao método dos deslocamentos surgem muitos conceitos que se supõe que o leitor já domina no âmbito das estruturas reticuladas como por exemplo as noções de grau de liberdade deslocamento generalizado força generalizada equilíbrio matriz de rigidez vector solicitação assemblagem introdução de condições de apoio etc Nesta publicação alguns destes conceitos são de novo abordados sendo dada particular ênfase à sua generalização aos meios contínuos bidimensionais e tridimensionais 11 Tipo de análise Quando surge a necessidade de resolver um problema de análise de uma estrutura a primeira questão que se coloca é a sua classificação quanto à geometria modelo do material constituinte e acções aplicadas O modo como o MEF é formulado e aplicado depende em parte das simplificações inerentes a cada tipo de problema Referemse em seguida alguns aspectos que é necessário ter em consideração na fase que antecede a análise de uma estrutura Introdução Álvaro F M Azevedo 3 Análise dinâmica ou estática As acções sobre as estruturas são em geral dinâmicas devendo ser consideradas as forças de inércia associadas às acelerações a que cada um dos seus componentes fica sujeito Por este motivo seria de esperar que a análise de uma estrutura teria obrigatoriamente de ter em consideração os efeitos dinâmicos Contudo em muitas situações é razoável considerar que as acções são aplicadas de um modo suficientemente lento tornando desprezáveis as forças de inércia Nestes casos a análise designase estática Nesta publicação apenas são considerados problemas em que se supõem válidas as simplificações inerentes a uma análise estática Análise não linear ou linear Na análise de uma estrutura sólida é habitual considerar que os deslocamentos provocados pelas acções exteriores são muito pequenos quando comparados com as dimensões dos componentes da estrutura Nestas circunstâncias admitese que não existe influência da modificação da geometria da estrutura na distribuição dos esforços e das tensões ie todo o estudo é feito com base na geometria inicial indeformada Se esta hipótese não for considerada a análise é designada não linear geométrica É também frequente considerar que ao nível do material que constitui a estrutura a relação entre tensões e deformações é linear Nos casos em que esta simplificação não é considerada é necessário recorrer a algoritmos específicos de análise não linear material Nesta publicação apenas se aborda o caso da análise linear quer geométrica quer material Tipo de estrutura As estruturas podem ser classificadas quanto à sua geometria como reticuladas laminares ou sólidas Estas últimas são as mais genéricas sendo classificadas como sólidas as que não apresentarem características que as permitam enquadrar no grupo das laminares ou das reticuladas As estruturas laminares são as que se desenvolvem para ambos os lados de uma superfície média mantendose na sua vizinhança É o caso de uma lâmina cuja Introdução Álvaro F M Azevedo 4 espessura é muito inferior às restantes dimensões Quando a superfície média é plana a estrutura laminar pode ser classificada como parede laje ou casca plana Uma parede apenas se encontra sujeita a acções paralelas ao seu plano médio Uma laje pode ter aplicadas forças perpendiculares ao plano médio e momentos cujo vector está contido no plano médio Uma estrutura laminar plana sujeita a outros tipos de acções é designada casca plana Quando a superfície média não é plana temse uma casca tridimensional As estruturas reticuladas são as constituídas por barras prismáticas cujas dimensões transversais são muito menores do que o comprimento do respectivo eixo Neste tipo de estruturas é habitual distinguir os pórticos das treliças conforme é ou não considerada a compatibilidade de rotações nas extremidades de barras adjacentes É possível tratar com grande eficiência uma classe de problemas de análise de estruturas designados axissimétricos Estes ocorrem quando a estrutura é um sólido de revolução e as acções são todas axissimétricas em relação ao mesmo eixo Neste tipo de problemas é ainda possível distinguir o caso do sólido de revolução do caso da lâmina de revolução Será também tratado como um caso particular a análise de uma estrutura que consiste num sólido cuja geometria a acções se repetem indefinidamente ao longo de um eixo rectilíneo Tratase do estado plano de deformação que pode ser estudado com base numa geometria bidimensional 12 Fundamentos do MEF A formulação do MEF requer a existência de uma equação integral de modo que seja possível substituir o integral sobre um domínio complexo de volume V por um somatório de integrais estendidos a sub domínios de geometria simples de volume Vi Esta técnica é ilustrada com o seguinte exemplo que corresponde ao integral de volume de uma função f n i V V i f dV dV f 1 1 Introdução Álvaro F M Azevedo 5 Em 1 pressupõese que n i iV V 1 2 Se for possível calcular todos os integrais estendidos aos sub domínios Vi basta efectuar o somatório correspondente ao segundo membro de 1 para se obter o integral estendido a todo o domínio Cada sub domínio Vi corresponde a um elemento finito de geometria simples eg segmento de recta triângulo quadrilátero tetraedro paralelepípedo O somatório indicado em 1 vai dar origem à operação designada assemblagem que apresenta muitas semelhanças com a que é efectuada nas estruturas reticuladas A equação integral referida no início desta secção é proveniente da aplicação do método dos resíduos pesados ou de um princípio variacional 13 No caso da aplicação do MEF à análise de estruturas a formulação mais intuitiva é a que se baseia no Princípio dos Trabalhos Virtuais PTV sendo a única que é abordada nesta publicação 13 Breve história do MEF Em 12 encontrase uma descrição detalhada da evolução do método dos elementos finitos ao longo do tempo Em 13 é efectuado o seu enquadramento com outros métodos da mesma família Apresentase em seguida apenas uma breve referência às principais fases do desenvolvimento do MEF Na generalidade dos casos é muito difícil definir a data em que determinado avanço do conhecimento foi efectuado No caso particular do MEF é referido por vários autores que a publicação mais antiga em que é utilizada a designação elemento finito é o artigo 14 que data de 1960 e tem como autor Ray Clough Anteriormente eram já conhecidas algumas técnicas que vieram a ser incorporadas no MEF sem este aparecer ainda com as principais características que hoje em dia possui Os grandes passos do desenvolvimento do MEF que o conduziram ao formato que actualmente apresenta maior aceitação foram dados na década de 60 e início da de 70 Inicialmente os elementos finitos mais comuns eram os triangulares e os tetraédricos passandose mais tarde a dar preferência aos quadriláteros e aos hexaedros Introdução Álvaro F M Azevedo 6 Ao contrário de outros métodos que eram utilizados no passado o MEF só tem utilidade prática se se dispuser de um computador digital Este requisito é devido à grande quantidade de cálculos que é necessário realizar nomeadamente na resolução de grandes sistemas de equações lineares Assim se compreende que o rápido desenvolvimento do MEF tenha praticamente coincidido com a generalização da utilização de computadores nos centros de investigação Com a proliferação de microcomputadores ocorrida no final da década de 80 e na década de 90 o MEF chega finalmente às mãos da generalidade dos projectistas de estruturas 14 Exemplo de aplicação do MEF Apresentase em seguida um exemplo de aplicação do MEF que consiste na análise de uma estrutura do tipo consola curta de pequena espessura sujeita às acções indicadas na Figura 11 Nestas condições podese admitir que se trata de um meio contínuo sujeito a um estado plano de tensão 15 Na Figura 11 está representada a malha utilizada que é constituída por 92 elementos finitos quadriláteros sendo cada um destes elementos definido por 8 nós Encontramse também assinalados os 10 nós que estão ligados ao meio exterior Depois de completada a análise da estrutura pelo MEF ficase a conhecer os valores aproximados dos deslocamentos e das tensões instaladas Na Figura 12 está representada a malha deformada pela acção das forças aplicadas à estrutura Para permitir uma melhor visualização dos deslocamentos estes são multiplicados por um factor de ampliação Como referência é também representada a malha original indeformada Com o tipo de visualização utilizado na Figura 13 é possível ter uma percepção imediata dos locais em que as tensões principais apresentam maiores valores bem como da trajectória das tensões dentro da estrutura Neste tipo de representação cada segmento de recta está orientado segundo uma direcção principal de tensão e a sua grandeza é proporcional ao valor da correspondente tensão normal A cor verde indica que se trata de uma tracção e à cor vermelha está associada uma compressão Na Figura 14 o valor da componente vertical do vector deslocamento é representado em cada ponto por intermédio de uma codificação por cores Consultando a escala Introdução Álvaro F M Azevedo 7 lateral ficase a conhecer a ordem de grandeza do deslocamento vertical em qualquer ponto da estrutura Na Figura 15 o tipo de visualização gráfica coincide com o da Figura 14 tratandose também da representação de um campo escalar por intermédio de uma codificação por cores O campo representado na Figura 15 é o das tensões normais σy sendo y o eixo vertical Esta componente do tensor das tensões é sempre perpendicular a facetas horizontais Fig 11 Consola curta malha de elementos finitos e acção exterior Introdução Álvaro F M Azevedo 8 Fig 12 Consola curta malha deformada representada sobre a estrutura indeformada Fig 13 Consola curta tensões principais e respectivas direcções Introdução Álvaro F M Azevedo 9 Fig 14 Consola curta campo de deslocamentos verticais Introdução Álvaro F M Azevedo 10 Fig 15 Consola curta campo de tensões normais segundo um eixo vertical Introdução Álvaro F M Azevedo 11 BIBLIOGRAFIA 11 Timoshenko S P Goodier J N Theory of Elasticity Third Edition McGrawHill 1988 12 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 13 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 14 Clough R W The Finite Element in Plane Stress Analysis Proc 2nd ASCE Conf on Electronic Computation Pittsburgh Pa September 1960 15 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 Introdução Álvaro F M Azevedo 12 13 CAPÍTULO 2 TRANSFORMAÇÃO LINEAR DE COORDENADAS Neste capítulo é apresentada a dedução da expressão que permite transformar as coordenadas de um ponto no espaço de um referencial S para outro S Quer os eixos de S quer os de S são definidos por versores cujas componentes se encontram no referencial geral S Estes três referenciais apresentam origem comum ponto O Sendo P um ponto genérico no espaço a transformação das componentes do vector OP coincide com a transformação das coordenadas do ponto P 21 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada neste capítulo Tabela 21 Simbologia relativa à transformação linear de coordenadas S Sistema de coordenadas referencial O Origem do sistema de coordenadas P Ponto genérico p Vector posição do ponto P x Eixo do sistema de coordenadas e Versor de um eixo do sistema de coordenadas A Matriz de transformação de S em S B Matriz de transformação de S em S g Referencial geral a Referencial auxiliar l Referencial local Transformação Linear de Coordenadas Álvaro F M Azevedo 14 α Ângulo entre eixos dos referenciais auxiliar e local Τ Matriz de transformação i Primeiro nó de uma barra j Segundo nó de uma barra L Comprimento de uma barra 22 Caso geral Na Figura 21 encontramse representados os três referenciais S S e S um ponto genérico P e o vector p OP 1x 2x 1x 2x 3x 3x 1x 2x 3x 1ˆe 2ˆe 3ˆe 1ˆe 2ˆe 3ˆe 1ˆe 2ˆe 3ˆe p O P Fig 21 Referenciais e ponto genérico P Os três referenciais que se supõem directos e ortonormados são definidos do seguinte modo Transformação Linear de Coordenadas Álvaro F M Azevedo 15 3 2 1 3 2 1 3 2 1 O x x x S O x x x S O x x x S 1 Versores de cada referencial 3 2 1 3 2 1 3 2 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e e e Versores de S e e e de S Versores e e e de S Versores 2 Ponto genérico x x x S P 3 2 1 3 Vector posição do ponto P 3 2 1 x x x OP p 4 Nota todos os versores e vectores apresentam as suas componentes no referencial S Versores do referencial S 100 ˆ 010 ˆ 001 ˆ 3 2 1 e e e 5 Vector p 3 2 1 p x x x 6 3 3 2 2 1 1 3 3 2 2 1 1 3 3 2 2 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ x e x e x e p x e x e x e p x e x e x e p 7 As coordenadas do ponto P no referencial 3 2 1 x x x S obtêmse projectando o vector p sobre os versores do referencial S Transformação Linear de Coordenadas Álvaro F M Azevedo 16 3 3 3 2 2 1 1 3 3 2 3 3 2 2 1 1 2 2 1 3 3 2 2 1 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e x e x e x e p e x e x e x e x e p e x e x e x e x e p e x 8 3 3 3 3 2 2 3 1 1 3 2 3 3 2 2 2 2 1 1 2 1 3 3 1 2 2 1 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ x e e x e e x e e x x e e x e e x e e x x e e x e e x e e x 9 Matricialmente temse 3 2 1 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 3 2 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ x x x e e e e e e e e e e e e e e e e e e x x x 10 A x x 11 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e e e e e e e e e e e e e e e e e e A 12 Nesta expressão x são as coordenadas de P no referencial S x são as coordenadas de P no referencial S e A é a matriz de transformação de S em S De um modo semelhante temse 3 3 3 2 2 1 1 3 3 2 3 3 2 2 1 1 2 2 1 3 3 2 2 1 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e x e x e x e p e x e x e x e x e p e x e x e x e x e p e x 13 3 3 3 3 2 2 3 1 1 3 2 3 3 2 2 2 2 1 1 2 1 3 3 1 2 2 1 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ x e e x e e x e e x x e e x e e x e e x x e e x e e x e e x 14 3 2 1 3 3 3 2 3 1 2 3 2 2 2 1 1 3 1 2 1 1 3 2 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ x x x e e e e e e e e e e e e e e e e e e x x x 15 Transformação Linear de Coordenadas Álvaro F M Azevedo 17 B x x 16 3 3 3 2 3 1 2 3 2 2 2 1 1 3 1 2 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e e e e e e e e e e e e e e e e e e B 17 Comparando 12 com 17 verificase que B AT 18 A expressão 16 pode escreverse da seguinte forma A x x T 19 Substituindo 11 em 19 temse A A x x T 20 Concluindose que I A A T 21 sendo I a matriz identidade Multiplicando ambos os membros de 21 por A1 à direita obtémse A1 A T 22 Quando a inversa de uma matriz coincide com a sua transposta dizse que a matriz é ortogonal Assim se conclui que a matriz de transformação A é uma matriz ortogonal Vaise agora proceder à análise do significado de cada um dos elementos de A A expressão 11 pode escreverse do seguinte modo 3 j 1 j ij i a x x 23 Transformação Linear de Coordenadas Álvaro F M Azevedo 18 sendo aij o elemento genérico da matriz A Em 12 verificase que j i ij e e a ˆ ˆ 24 Recorrendo à definição de produto escalar temse j i j i ij e e e e a cos ˆ ˆ ˆ ˆ 25 Uma vez que os versores dos referenciais possuem norma unitária j i ij e e a cos ˆ ˆ 26 e a matriz de transformação A pode ser obtida a partir dos cosenos dos ângulos entre versores dos referenciais S e S 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 ˆ cos ˆ ˆ cos ˆ ˆ ˆ cos ˆ cos ˆ ˆ cos ˆ ˆ ˆ cos ˆ cos ˆ cos ˆ ˆ ˆ ˆ cos e e e e e e e e e e e e e e e e e e A 27 23 Caso particular com S e S coincidentes Reproduzemse em seguida as expressões 5 11 e 12 100 ˆ 010 ˆ 001 ˆ 3 2 1 e e e 28 A x x 29 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e e e e e e e e e e e e e e e e e e A 30 No caso de os referenciais S e S serem coincidentes verificase que Transformação Linear de Coordenadas Álvaro F M Azevedo 19 i i e e ˆ ˆ 31 x A x 32 Substituindo 31 em 30 obtémse 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e e e e e e e e e e e e e e e e e e A 33 Atendendo a 28 verificase em 33 que a primeira linha da matriz A contém as componentes do versor 1ˆe no referencial S A segunda e terceira linhas contêm as componentes em S dos versores 2ˆe e 3ˆe em S Componentes de e em S s de e Componente em S s de e Componente A 3 2 1 3 3 ˆ ˆ ˆ 34 24 Matriz de transformação de uma barra rectilínea no espaço Nesta secção são utilizadas as expressões deduzidas nas secções anteriores com o objectivo de chegar à matriz de transformação de uma barra de treliça 3D e de pórtico 3D No âmbito da análise de estruturas pelo método dos deslocamentos admitemse as seguintes hipóteses é conhecida a geometria da estrutura que é constituída por barras prismáticas de eixo rectilíneo e de secção constante para cada barra são conhecidas as coordenadas dos dois nós extremos ficando assim definida a localização do seu eixo baricêntrico é conhecida a posição dos eixos principais centrais de inércia da secção transversal da barra 21 Transformação Linear de Coordenadas Álvaro F M Azevedo 20 Considerese um ângulo α que será definido adiante e que posiciona o referencial local principal central de inércia PCI em relação a um referencial auxiliar Assim vão ser considerados os seguintes referenciais local PCI l S auxiliar a S geral g S 0 α 35 O referencial geral g é aquele em relação ao qual todos os pontos e todos os vectores estão definidos sendo os seus versores definidos por 28 O referencial auxiliar a ao qual corresponde um ângulo α nulo tem o primeiro eixo coincidente com o eixo da barra e o segundo eixo perpendicular ao plano vertical que contem a barra O terceiro eixo é aquele que faz com que o referencial seja directo e ortonormado Este referencial será adiante definido com mais rigor O referencial local l tem como primeiro eixo o eixo da barra sendo os restantes eixos os eixos principais centrais de inércia da secção transversal da barra O ângulo α define a posição do referencial local l em relação ao referencial auxiliar a Vão ser em seguida definidas duas transformações transformação de g para a transformação de a para l A primeira transformação é realizada com a seguinte expressão que é semelhante a 32 g ag a x x T 36 sendo T ag a matriz que transforma as coordenadas de um ponto do referencial g para o referencial a A segunda transformação permite obter as coordenadas de um ponto no referencial l a partir das suas coordenadas no referencial a sendo semelhante à definida por 11 Transformação Linear de Coordenadas Álvaro F M Azevedo 21 a la l x x T 37 Substituindo 36 em 37 chegase a g ag la l x T x T 38 Uma vez que se pretende uma matriz de transformação de g para l g l x T x 39 comparando 38 com 39 concluise que T Tla T ag 40 Na Figura 22 é definida a posição do referencial auxiliar a em relação ao referencial geral g e à barra g1 g2 g 3 a1 a2 a 3 i j i j Fig 22 Posição do referencial a em relação ao referencial g Em relação à Figura 22 considerase ainda o seguinte o eixo g3 é vertical e orientado para cima o eixo baricêntrico da barra é definido pelos nós i e j Transformação Linear de Coordenadas Álvaro F M Azevedo 22 é em geral vantajoso considerar a convenção de ser sempre i j Assim o primeiro nó da barra é o nó i e o segundo é o nó j Esta convenção clarifica todo o processo de estudo da barra sem lhe introduzir qualquer limitação o eixo a1 coincide com o eixo baricêntrico da barra ie o eixo que é definido pelos centros de gravidade de todas as secções transversais da barra o eixo a1 encontrase orientado do nó i para o nó j o eixo a2 é perpendicular ao plano g3a1 e está orientado de acordo com o sentido do produto vectorial entre os versores de g3 e a1 o eixo a3 está contido no plano g3a1 e resulta do produto vectorial entre os versores de a1 e a2 desta forma o referencial a1a2a3 é sempre directo e ortonormado Para se calcular a matriz de transformação de g para a 36 vaise recorrer à expressão 34 Assim a primeira linha de T ag é constituída pelas componentes do versor a1 no referencial g e assim sucessivamente O cálculo das componentes do versor a1 é feito com base nas coordenadas dos nós i e j Coordenadas do nó i no referencial g i i i x x x 3 2 1 Coordenadas do nó j no referencial g j j j x x x 3 2 1 O comprimento da barra é calculado com a seguinte expressão 2 3 3 2 2 2 2 1 1 i j i j i j x x x x x x L 41 O vector 1a que em geral não tem norma unitária obtémse por subtracção das coordenadas dos nós i e j i j i j i j x x x x x x a 3 3 2 2 1 1 1 42 O versor 1ˆa obtémse dividindo o vector 1a pela respectiva norma Transformação Linear de Coordenadas Álvaro F M Azevedo 23 L a a 1ˆ 1 43 Para posterior referência designamse as componentes do versor 1ˆa por A1 A2 e A3 3 2 1 1 ˆ a A A A 44 Tal como foi atrás referido o eixo a2 é definido pelo produto vectorial entre os versores dos eixos g3 e a1 sendo 100 ˆ3 g 1 3 2 ˆ ˆ a g a 45 Uma vez que deste produto vectorial não resulta um versor é necessário dividir o vector 2 a pela respectiva norma 2 2 2ˆ a a a 46 Para posterior referência designamse as componentes do versor 2ˆa por B1 B2 e B3 3 2 1 2 ˆ a B B B 47 Para que o referencial a seja directo e ortonormado calculase o versor 3ˆa como sendo o resultado do produto vectorial entre 1ˆa e 2ˆa Do produto vectorial entre versores perpendiculares entre si resulta sempre um versor 2 1 3 ˆ ˆ ˆ a a a 48 Para posterior referência designamse as componentes do versor 3ˆa por C1 C2 e C3 3 2 1 3 ˆ a C C C 47 De acordo com o que foi deduzido os elementos da matriz de transformação do referencial g para o referencial a 36 são os seguintes 3 2 1 3 2 1 3 2 1 C C C B B B A A A T ag 48 Transformação Linear de Coordenadas Álvaro F M Azevedo 24 O resultado do produto vectorial expresso em 45 é um vector nulo sempre que o versor 1ˆa seja paralelo ao versor 3ˆg Supondo que o eixo 3ˆg é sempre vertical hipótese considerada atrás esta situação singular ocorre sempre que a barra é vertical Para estes casos é então necessário definir a matriz de transformação T ag com outro critério Na Figura 23 e na Figura 24 encontrase a posição do referencial a em relação ao referencial g para os casos da barra vertical orientada para cima e orientada para baixo g1 g2 g3 a1 a2 a3 i j i j 1 00 ˆ 0 1 0 ˆ 00 1 ˆ 3 2 1 a a a Fig 23 Posição do referencial a em relação ao referencial g para o caso da barra vertical orientada para cima g1 g2 g 3 a1 a2 a3 i j i j 1 00 ˆ 0 1 0 ˆ 00 1 ˆ 3 2 1 a a a Fig 24 Posição do referencial a em relação ao referencial g para o caso da barra vertical orientada para baixo Considerando as seguintes expressões para os versores do referencial a ficam cobertas as duas situações esquematizadas nas Figuras 23 e 24 3 2 1 3 3 1 00 ˆ A A A L x x a i j 49 Transformação Linear de Coordenadas Álvaro F M Azevedo 25 3 2 1 2 0 1 0 ˆ B B B a 50 3 2 1 3 3 3 00 ˆ C C C L x x a i j 51 Tal como em 48 a matriz de transformação T ag é constituída por 3 2 1 3 2 1 3 2 1 C C C B B B A A A T ag 52 Procedese em seguida à definição da matriz T la que foi referida em 37 Esta matriz de transformação relaciona as coordenadas de um ponto no referencial auxiliar a com as suas coordenadas no referencial local l As considerações que se seguem baseiamse na Figura 25 em que estão representados os referenciais a e l O referencial l é constituído pelo eixo da barra e pelos eixos principais centrais de inércia da secção transversal l1 l 2 l3 a1 a2 a3 i j i j α α Fig 25 Posição do referencial l em relação ao referencial a De acordo com a Figura 25 podese constatar o seguinte os eixos a1 e l1 coincidem os eixos l2 e l3 estão rodados de um ângulo α em relação aos eixos a2 e a3 Transformação Linear de Coordenadas Álvaro F M Azevedo 26 A transformação entre os referenciais a e l é um caso de transformação entre dois referenciais distintos do geral Nesta situação podese recorrer à matriz definida em 27 que corresponde a uma transformação entre os referenciais S e S Neste caso o referencial S é o referencial a e o referencial S é o referencial l A matriz de transformação é neste caso calculada com base nos cosenos dos ângulos formados pelos eixos dos dois referenciais 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 cos cos cos cos cos cos cos cos cos l a l a l a a l a l a l l a l a a l T la 53 De acordo com a Figura 25 temse α α α α cos cos 90 cos 90 cos 90 cos 90 cos cos 90 cos 90 0 cos T la 54 α α α α cos sin 0 sin cos 0 0 0 1 T la 55 As matrizes de transformação T ag e T la encontramse já definidas De acordo com 40 a matriz de transformação T do referencial geral para o local é definida do seguinte modo T Tla T ag 56 Tal como foi indicado em 39 a correspondente transformação é efectuada com a seguinte expressão g l x T x 57 As expressões aqui deduzidas e que permitem calcular a matriz T foram baseadas na informação de que é habitual dispor numa análise de um pórtico 3D pelo método dos deslocamentos ie das coordenadas dos nós e do ângulo α Transformação Linear de Coordenadas Álvaro F M Azevedo 27 Uma vez que a matriz T é ortogonal a transformação do referencial local para o geral é efectuada com a seguinte relação l T g x T x 58 25 Considerações finais As expressões da matriz de transformação deduzidas neste capítulo podem ser directamente utilizadas na formulação da matriz de rigidez de elementos de treliça ou de pórtico 3D bem como na formulação dos respectivos vectores de forças nodais equivalentes BIBLIOGRAFIA 21 Brazão Farinha J S Correia dos Reis A Tabelas Técnicas Edições Técnicas E T L 1998 Transformação Linear de Coordenadas Álvaro F M Azevedo 28 29 CAPÍTULO 3 MÉTODO DOS DESLOCAMENTOS EM TRELIÇAS E PÓRTICOS Com o objectivo de apresentar alguns conceitos como o de assemblagem e introdução de condições de apoio fazse aqui uma sucinta descrição do método dos deslocamentos aplicado à análise de treliças e pórticos tridimensionais 31 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação do método dos deslocamentos em treliças e pórticos Tabela 31 Simbologia relativa ao método dos deslocamentos em estruturas reticuladas g Referencial geral a Referencial auxiliar l Referencial local i Primeiro nó de uma barra j Segundo nó de uma barra α Ângulo entre eixos dos referenciais auxiliar e local xg Coordenadas de um ponto no referencial geral xl Coordenadas de um ponto no referencial local T Matriz de transformação a Deslocamento ou deslocamento generalizado θ Rotação F Força ou força generalizada M Momento Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 30 a Deslocamentos nodais nos graus de liberdade da estrutura no referencial geral ag Deslocamentos nodais nos graus de liberdade da barra no referencial geral al Deslocamentos nodais nos graus de liberdade da barra no referencial local K Matriz de rigidez da estrutura no referencial geral Kg Matriz de rigidez da barra no referencial geral Kl Matriz de rigidez da barra no referencial local F Forças nodais equivalentes à acção exterior nos graus de liberdade da estrutura no referencial geral Fg Forças nodais equivalentes à acção exterior nos graus de liberdade da barra no referencial geral Fl Forças nodais equivalentes à acção exterior nos graus de liberdade da barra no referencial local L Índice correspondente a um grau de liberdade não prescrito livre P Índice correspondente a um grau de liberdade prescrito R Reacção num apoio da estrutura n Número de graus de liberdade não prescritos livres p Número de graus de liberdade prescritos E Módulo de Young de um material A Área da secção transversal de uma barra L Comprimento de uma barra G Módulo de distorção de um material I Momento de inércia da secção transversal de uma barra It Momento de inércia de torção da secção transversal de uma barra Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 31 32 Referenciais De acordo com o que foi descrito no Capítulo 2 na formulação da matriz de rigidez de uma barra de eixo rectilíneo e de secção constante são considerados dois referenciais directos e ortonormados o geral g1g2g3 e o local l1l2l3 O referencial geral é aquele em que se encontram expressas as coordenadas de todos os nós que depois são utilizados para definir a posição das barras O referencial local é definido pelos seguintes eixos l1 é o eixo da barra e l2 e l3 são os eixos principais centrais de inércia da secção transversal da barra ver a Figura 31 g 1 g 2 g3 l 1 l2 l 3 i j i j Fig 31 Barra i j referencial geral g e referencial local l Considerase habitualmente sem perda de generalidade que a barra definida pelos nós i e j tem o nó i coincidente com a origem dos dois referenciais e o nó j sobre o semieixo positivo l1 É também habitual considerar que o número do nó i é inferior ao número do nó j i j Os eixos l2 e l3 podem ser trocados entre si tendo em atenção que o referencial local deve ser sempre directo A troca de l2 com l3 obriga a trocar entre si os valores dos momentos de inércia em relação a l2 e l3 Em qualquer dos casos é necessário definir criteriosamente o ângulo α ver o Capítulo 2 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 32 A transformação de coordenadas entre os referenciais g e l é efectuada com a seguinte expressão em que T é a matriz de transformação 3x3 definida também no Capítulo 2 g l x T x 1 Nesta expressão xg são as coordenadas de um ponto no referencial g e lx são as coordenadas desse mesmo ponto no referencial l A equação 1 também pode ser utilizada para transformar as componentes de um vector do referencial g para o referencial l 33 Graus de liberdade Num ponto do espaço pertencente a um corpo sujeito a deslocamentos e deformações podem ser considerados seis graus de liberdade três de deslocamento e três de rotação 6 5 4 3 2 1 3 2 1 3 2 1 a a a a a a a a a a θ θ θ 2 Designase por deslocamentos generalizados o agrupamento dos três deslocamentos e das três rotações num só vector com seis componentes ver a Figura 32 a 1 a 2 a3 a 4 θ 1 a 5 θ 2 a 6 θ 3 Fig 32 Deslocamentos generalizados Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 33 No estudo de um pórtico 3D são considerados os seis deslocamentos generalizados em cada ponto nodal da barra ou da estrutura O caso da treliça 3D em que apenas são considerados três deslocamentos em cada ponto nodal a1 a2 e a3 pode ser adaptado do pórtico 3D bastando eliminar tudo o que diz respeito a rotações e momentos Para se passar da treliça 3D para a treliça 2D basta suprimir tudo o que diz respeito a um dos três graus de liberdade Os pórticos 2D grelhas e vigas contínuas são também simplificações do caso do pórtico 3D Por ser o caso mais genérico de aqui em diante apenas se desenvolve a formulação da barra de pórtico 3D Em correspondência com os seis deslocamentos generalizados são consideradas seis forças generalizadas 3 forças e 3 momentos que se representam na Figura 33 F 1 F 2 F3 F4 Μ 1 F5 Μ2 F6 Μ3 Fig 33 Forças generalizadas Na Figura 34 encontrase representada uma barra de dois nós i e j Em cada nó são considerados seis graus de liberdade em correspondência com os seis deslocamentos generalizados 2 Assim o número de graus de liberdade da barra é doze Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 34 3 1 4 2 5 6 9 7 10 8 11 12 i j i j l 1 l 2 l 3 Fig 34 Graus de liberdade da barra i j no referencial local Em correspondência com os doze graus liberdade representados na Figura 34 têmse também as forças e os momentos que actuam nas extremidades da barra 34 Matriz de transformação A matriz de transformação T referida em 1 é uma matriz 3x3 cujos componentes são 33 32 31 23 22 21 13 12 11 T T T T T T T T T T 3 A transformação dos doze deslocamentos generalizados representados na Figura 34 pode ser efectuada com a seguinte relação desde que a matriz de transformação T passe a ser uma matriz 12x12 constituída pela repetição de 3 quatro vezes 12 1 12 12 12 1 g l a T a 4 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 35 g g g g g g g g g g g g l l l l l l l l l l l l a a a a a a a a a a a a T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T T a a a a a a a a a a a a 12 11 10 9 8 7 6 5 4 3 2 1 33 32 31 23 22 21 13 12 11 33 32 31 23 22 21 13 12 11 33 32 31 23 22 21 13 12 11 33 32 31 23 22 21 13 12 11 12 11 10 9 8 7 6 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 35 Matriz de rigidez e vector solicitação Supondo o caso de uma barra de eixo rectilíneo e secção constante a respectiva matriz de rigidez no referencial local Kl bem como o vector de forças nodais equivalentes a diversos tipos de acções Fl podem ser directamente obtidos com base num formulário de estruturas 31 ver também as Secções 39 e 310 Assim partese do princípio que se dispõe da matriz l K e do vector l F que se relacionam com a habitual equação 12 1 12 1 12 12 l l l F a K 6 sendo la o vector dos deslocamentos generalizados da barra no referencial local As equações 4 e 5 são válidas quer para os deslocamentos generalizados quer para as forças generalizadas tendose também 12 1 12 12 12 1 g l F T F 7 Uma vez que a matriz de transformação é ortogonal ie T T T 1 8 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 36 multiplicamse ambos os membros de 7 por T T e obtémse 12 1 12 12 12 1 l T g F T F 9 Substituindo em 9 a equação 6 12 1 12 12 12 1 l l l a K F 10 resulta 12 1 12 12 12 12 12 1 l l T g a K T F 11 Substituindo 4 em 11 chegase a 12 1 12 12 12 12 12 12 12 1 g l T g a T K T F 12 Uma vez que a relação de rigidez da barra no referencial geral é 12 1 12 1 12 12 g g g F a K 13 Da comparação de 12 com 13 concluise que a matriz de rigidez da barra de pórtico 3D no referencial geral é dada por 12 12 12 12 12 12 12 12 T K T K l T g 14 O vector solicitação F g pode ser calculado com a expressão 9 Depois de serem conhecidos os deslocamentos g a é possível calcular as acções nas extremidades das barras no referencial local recorrendo à seguinte expressão que resulta da substituição de 4 em 10 12 1 12 12 12 12 12 1 g l l a T K F 15 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 37 36 Assemblagem da matriz de rigidez global e do vector solicitação Depois de calculadas todas as matrizes de rigidez das barras no referencial geral com recurso à expressão 14 é necessário proceder ao cálculo da matriz de rigidez global da estrutura Uma operação semelhante tem de ser efectuada com os vectores solicitação das diversas barras A assemblagem na matriz de rigidez global das matrizes de rigidez das diversas barras é em seguida apresentada com base no exemplo da Figura 35 1 A a 1 2 3 4 B C D a 2 a 3 a 4 F 1 F 2 F3 F4 Fig 35 Assemblagem num exemplo unidimensional A estrutura representada na Figura 35 é unidimensional tem quatro nós 1 a 4 e quatro barras A a D Cada barra tem as suas características nomeadamente o módulo de Young E a área da secção transversal A e o comprimento L Em cada nó existe um único grau de liberdade Em correspondência com os quatro graus de liberdade existem quatro deslocamentos nodais a e quatro forças nodais equivalentes à acção exterior F Cada barra tem dois graus de liberdade um em cada extremidade Para cada barra é conhecida a matriz de rigidez 2x2 no referencial geral cuja designação se simplifica de acordo com 22 21 12 11 22 21 12 11 A A A A K K K K K A Barra A A A A A 16 22 21 12 11 22 21 12 11 B B B B K K K K K B Barra B B B B B 17 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 38 22 21 12 11 22 21 12 11 C C C C K K K K K Barra C C C C C C 18 22 21 12 11 22 21 12 11 D D D D K K K K K D Barra D D D D D 19 Atendendo à numeração global dos graus de liberdade 1 a 4 as matrizes de rigidez das barras passam a ser 0 0 0 0 0 0 0 0 0 0 0 0 2 1 22 21 12 11 A A A A K A Barra A 20 0 0 0 0 0 0 0 0 0 0 0 0 3 2 22 21 12 11 B B B B K B Barra B 21 22 21 12 11 0 0 0 0 0 0 0 0 0 0 0 0 4 3 C C C C K Barra C C 22 22 21 12 11 0 0 0 0 0 0 0 0 0 0 0 0 4 2 D D D D K D Barra D 23 O vector dos deslocamentos em todos os graus de liberdade da estrutura é 4 3 2 1 a a a a a 24 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 39 A a 1 B C D a 2 a 3 a 4 F 1 A F 2 A F1 B F2 B F1 C F2 C F1 D F2 D Fig 36 Vectores das forças nodais equivalentes a acções exteriores Atendendo à numeração global dos graus de liberdade os vectores das forças nodais equivalentes às acções nas diversas barras são ver a Figura 36 0 0 2 1 2 1 A A A F F F Barra A 25 0 0 3 2 2 1 B B B F F F Barra B 26 C C C F F F C Barra 2 1 0 0 4 3 27 D D D F F F D Barra 2 1 0 0 4 2 28 Os vectores e matrizes indicados em 2028 relacionamse entre si de acordo com as seguintes equações Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 40 A A F a K 29 B B F a K 30 C C F a K 31 D D F a K 32 A soma dos primeiros membros das equações 2932 é igual à soma dos seus segundos membros resultando D C B A D C B A F F F F a K a K a K a K 33 D C B A D C B A F F F F a K K K K 34 Uma vez que a relação de rigidez envolvendo todos os graus de liberdade da estrutura é F K a 35 concluise que D C B A K K K K K 36 e D C B A F F F F F 37 Adicionando as matrizes 2023 de acordo com 36 chegase a 22 22 21 21 12 11 22 21 12 12 11 11 22 21 12 11 0 0 0 0 D C C D C C B B D B D B A A A A K 38 Adicionando os vectores solicitação 2528 de acordo com 37 chegase a Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 41 D C C B D B A A F F F F F F F F F 2 2 1 2 1 1 2 1 39 O procedimento de assemblagem aqui exposto é generalizável ao caso em que existem seis graus de liberdade em cada nó Para esse fim é suficiente considerar que por exemplo C12 em vez de ser um escalar é uma matriz 6x6 contendo os elementos da matriz KC que relacionam os graus de liberdade do nó 1 com os graus de liberdade do nó 2 37 Introdução das condições de apoio O sistema de equações 35 ainda não pode ser resolvido porque falta entrar em linha de conta com as condições de apoio da estrutura Estas condições fronteira correspondem a apoios fixos ou assentamentos de apoio Os apoios fixos podem sempre ser tratados como assentamentos de apoio de valor nulo Por este motivo no desenvolvimento que se segue apenas são referidos os assentamentos de apoio O sistema de equações 35 relaciona forças e deslocamentos que se encontram no referencial geral englobando todos os graus de liberdade da estrutura Tendo em vista a consideração das condições de apoio os graus de liberdade da estrutura são divididos em dois grupos L graus de liberdade não prescritos livres P graus de liberdade prescritos Assim o sistema de equações 35 passa a ter a seguinte organização por blocos P P L P L PP PL LP LL R F F a a K K K K F K a 0 40 Em 40 aL é o vector que engloba os deslocamentos segundo os graus de liberdade não prescritos e aP engloba os prescritos O mesmo tipo de subdivisão é efectuado com o vector das forças nodais equivalentes à acção exterior F O vector adicional em que Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 42 figura RP contém as reacções de apoio que consistem nas forças ainda desconhecidas que fazem com que os deslocamentos em apoios assumam os valores prescritos Designando por n o número de graus de liberdade não prescritos e por p o número de graus de liberdade prescritos são especificadas na Tabela 32 as dimensões das submatrizes que figuram em 40 Tabela 32 Dimensões das submatrizes presentes em 40 K LL n x n K LP n x p K PL p x n K PP p x p L FL a n x 1 P P P R F a p x 1 Esta divisão em submatrizes obriga a fazer uma reorganização das linhas e das colunas da matriz K que figura em 35 bem como das componentes dos vectores a e F Na Tabela 33 é apresentado o significado dos elementos das quatro submatrizes de K indicadas em 40 Tabela 33 Significado dos elementos das submatrizes de K indicadas em 40 Deslocamento unitário imposto segundo um grau de liberdade Forças de fixação num grau de liberdade K LL Livre Livre K LP Livre Prescrito K PL Prescrito Livre K PP Prescrito Prescrito No novo sistema de equações indicado em 40 as incógnitas são aL e P R Os elementos de K P a FL e FP têm valores conhecidos Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 43 O sistema de equações 40 pode ser escrito do seguinte modo L P LP L LL F a K a K 41 P P P PP L PL R F a K a K 42 A equação 41 pode ser rescrita do seguinte modo P LP L L LL a K F a K 43 Em 43 K LL é uma matriz quadrada que em geral é não singular aL é o vector das incógnitas e os valores dos vectores e matrizes que estão no segundo membro são conhecidos Por este motivo 43 constitui um sistema de equações lineares que depois de resolvido fornece os valores dos deslocamentos L a A equação 42 pode ser rescrita do seguinte modo P P PP L PL P F a K a K R 44 Uma vez que os deslocamentos aL já são conhecidos esta expressão fornece os valores das reacções em graus de liberdade prescritos P R O modo de introdução das condições de apoio aqui descrito tem as seguintes vantagens na fase do processo que requer um maior volume de cálculos e uma grande quantidade de memória de armazenamento ie na fase de resolução do sistema de equações 43 o número de equações e incógnitas é n em vez de ser np em comparação com o método em que é adicionado à diagonal principal de K um número elevado o método aqui proposto apresenta menos problemas numéricos principalmente quando se utilizam métodos iterativos para resolver o sistema de equações A principal desvantagem do método aqui proposto é a necessidade de agrupar os elementos de K em diversas submatrizes Esta nova arrumação causa algumas dificuldades principalmente quando se utilizam técnicas de armazenamento esparso em banda ou em perfil Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 44 38 Faseamento da análise de um pórtico 3D Tendo em vista a análise de uma estrutura do tipo pórtico 3D pelo método dos deslocamentos sugerese o seguinte algoritmo Para cada barra Calcular a matriz de transformação T 3 e em seguida calcular 5 Calcular a matriz de rigidez da barra no referencial local l K Calcular a matriz de rigidez da barra no referencial geral g K com 14 Assemblar g K em K ver a Secção 36 Calcular o vector das forças nodais equivalentes à acção exterior na barra no referencial local l F Calcular g F com 9 Assemblar g F em F ver a Secção 36 Introduzir as condições de apoio ver a Secção 37 Resolver o sistema de equações lineares 43 determinando assim os deslocamentos Calcular as reacções nos apoios com 44 Para cada barra Passar os deslocamentos relativos à barra corrente do vector a para o vector g a Calcular l F com 15 Fim Embora seja possível utilizar o procedimento sugerido sem recursos informáticos é hoje em dia preferível implementálo por intermédio de um programa de computador Neste domínio surgem muitas alternativas tais como a selecção da linguagem de Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 45 programação o modo de criar os dados do problema o modo de armazenamento da informação as técnicas numéricas utilizadas o recurso ou não a bibliotecas de operações matriciais etc 39 Matriz de rigidez de uma barra de treliça 3D no referencial local Na Figura 37 encontrase representada uma barra de treliça espacial de eixo rectilíneo e secção constante A sua matriz de rigidez 45 expressa no referencial local l depende das seguintes grandezas E módulo de Young constante em todos os pontos da barra A área da secção transversal da barra considerada constante L comprimento da barra 3 1 4 2 5 6 i j i j l 1 l 2 l3 Fig 37 Treliça 3D graus de liberdade da barra i j no referencial local 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 EA L L EA EA L L EA K l 45 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 46 310 Matriz de rigidez de uma barra de pórtico 3D no referencial local Na Figura 38 encontrase representada uma barra de pórtico espacial de eixo rectilíneo e secção constante A sua matriz de rigidez 4650 expressa no referencial local l depende das seguintes grandezas E módulo de Young constante em todos os pontos da barra A área da secção transversal da barra considerada constante L comprimento da barra G módulo de distorção 32 I2 momento de inércia da secção transversal da barra em relação ao eixo l2 I3 momento de inércia da secção transversal da barra em relação ao eixo l3 It momento de inércia de torção da secção transversal da barra 33 34 Nota l2 e l3 são eixos principais centrais de inércia da secção transversal da barra 3 1 4 2 5 6 9 7 10 8 11 12 i j i j l1 l2 l 3 Fig 38 Pórtico 3D graus de liberdade da barra i j no referencial local j j l ij l i j l ii l l K K K K K 46 Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 47 L EI L EI L EI L EI L GI L EI L EI L EI L EI L EA K t ii l 3 2 3 2 2 2 2 2 3 2 2 3 3 3 4 0 0 0 6 0 0 4 0 6 0 0 0 0 0 0 0 0 6 0 12 0 0 6 0 0 0 12 0 0 0 0 0 0 47 L EI L EI L EI L EI L GI L EI L EI L EI L EI L EA K t i j l 3 2 3 2 2 2 2 2 3 2 2 3 3 3 2 0 0 0 6 0 0 2 0 6 0 0 0 0 0 0 0 0 6 0 12 0 0 6 0 0 0 12 0 0 0 0 0 0 48 T i j l l ij K K 49 L EI L EI L EI L EI L GI L EI L EI L EI L EI L EA K t j j l 3 2 3 2 2 2 2 2 3 2 2 3 3 3 4 0 0 0 6 0 0 4 0 6 0 0 0 0 0 0 0 0 6 0 12 0 0 6 0 0 0 12 0 0 0 0 0 0 50 311 Considerações finais Neste capítulo não foi considerada a possibilidade da a barra apresentar eixo não rectilíneo nem o facto de a secção transversal ser variável ao longo do eixo da barra Não foi também considerada a contribuição das tensões tangenciais para a deformação habitualmente designada deformação por esforço transverso A inclusão destas características faz com que a formulação apresentada neste capítulo perca a simplicidade atrás evidenciada Mais adiante serão apresentadas formulações da matriz de rigidez de uma barra recorrendo a técnicas específicas do Método dos Elementos Finitos em particular a formulação de viga de Timoshenko Com este tipo de elementos de barra é possível ter em consideração a deformação por esforço transverso o eixo curvilíneo e a secção variável Método dos Deslocamentos em Treliças e Pórticos Álvaro F M Azevedo 48 BIBLIOGRAFIA 31 Brazão Farinha J S Correia dos Reis A Tabelas Técnicas Edições Técnicas E T L 1998 32 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 33 Segadães Tavares A Análise Matricial de Estruturas Laboratório Nacional de Engenharia Civil Curso 129 Lisboa 1973 34 Massonnet C Résistance des Matériaux Dunod Paris 1968 49 CAPÍTULO 4 ELEMENTOS FINITOS UNIDIMENSIONAIS Antes de expor o método dos elementos finitos MEF de um modo aplicável a meios contínuos bidimensionais e tridimensionais apresentase com algum detalhe o caso unidimensional Quando apenas se considera uma dimensão o método resultante não tem grande interesse prático mas serve como introdução às técnicas que mais adiante serão expostas para os casos mais genéricos O método dos elementos finitos que adiante será exposto baseiase no método dos deslocamentos e na discretização de uma estrutura em subestruturas Cada uma dessas subestruturas designase por elemento finito e tem comportamento conhecido sendo o comportamento do todo considerado como a soma das partes Cada elemento finito tem n nós sendo apenas considerados explicitamente os deslocamentos generalizados nesses nós Os deslocamentos nos restantes pontos do elemento finito obtêmse por interpolação dos deslocamentos dos nós 41 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação do método dos elementos finitos Tabela 41 Simbologia relativa ao método dos elementos finitos n Número de nós do elemento finito L Comprimento da barra prismática x Coordenada cartesiana u Campo de deslocamentos a Deslocamento nodal N Função interpoladora ou função de forma Elementos Finitos Unidimensionais Álvaro F M Azevedo 50 ε Deformação B Matriz de deformação L 1 Operador diferencial L 1 d d x V Volume da barra prismática σ Tensão normal p Acção exterior distribuída por unidade de comprimento F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial local A Área da secção transversal da barra prismática E Módulo de elasticidade ou módulo de Young D Matriz de elasticidade ε σ D K Matriz de rigidez do elemento finito no referencial local c Coeficiente de um termo de um polinómio x Coordenada cartesiana de um nó de um elemento finito s Coordenada local E Módulo de elasticidade num nó do elemento finito A Área da secção transversal num nó do elemento finito J Jacobiano da transformação J d x d s 42 Funções interpoladoras ou funções de forma Na Figura 41 encontrase representado um elemento finito unidimensional com dois nós e com comprimento L 2 Elementos Finitos Unidimensionais Álvaro F M Azevedo 51 u x x 1a 2 a 1 2 L 2 x 1 x 1 Fig 41 Elemento finito unidimensional de dois nós O único eixo coordenado que é considerado é o eixo x ocorrendo todos os deslocamentos paralelamente a x A função u x corresponde ao campo de deslocamentos verificandose o seguinte 2 1 1 1 a u a u 1 sendo portanto a1 e a2 os deslocamentos dos nós Considerese agora como aproximação que a lei de variação do deslocamento entre os nós 1 e 2 é linear Nestas circunstâncias a seguinte função u x representa o campo de deslocamentos porque é linear em x e respeita 1 a x a a a u x 2 2 1 2 2 1 2 Os valores numéricos dos parâmetros a1 e a2 passarão a ser conhecidos depois de analisada a estrutura Colocando a1 e a2 em evidência em 2 chegase à seguinte expressão 2 1 2 1 2 1 2 1 2 1 a x a x u x 3 Em 3 temse uma soma de produtos de funções lineares de x pelos deslocamentos nodais a1 e a2 Elementos Finitos Unidimensionais Álvaro F M Azevedo 52 A equação 3 pode ser escrita em forma matricial 2 1 2 1 2 1 2 1 2 1 a a x x u x 4 ou 2 1 2 1 a a N x N x u x 5 sendo x x N x x N 2 1 2 1 2 1 2 1 2 1 6 e u N a 7 com 2 1 2 1 N N N x N x N 8 e 2 1 a a a 9 O gráfico das funções lineares N1 e N2 indicadas em 6 encontrase representado na Figura 42 Elementos Finitos Unidimensionais Álvaro F M Azevedo 53 N1x x 1 1 1 N2x 1 1 1 x Fig 42 Gráfico das funções N1x e N2x A principal característica dos gráficos das funções N1x e N2x é salientada na Tabela 42 e consiste no facto de a função N1x assumir o valor unitário no nó 1 e nulo nos restantes nós A função N2x assume o valor unitário no nó 2 e nulo nos restantes nós Esta característica será clarificada adiante quando se apresentarem exemplos de elementos finitos com mais do que dois nós Tabela 42 Características das funções N1x e N2x x 1 1 N1x 1 0 N2x 0 1 Apresentamse em seguida as funções de forma N1x e N2x para o caso da barra de dois nós de comprimento L ver a Figura 43 u x x 1a 2 a 1 2 L x L 2 x L 2 Fig 43 Elemento finito unidimensional de dois nós com comprimento L Elementos Finitos Unidimensionais Álvaro F M Azevedo 54 De um modo semelhante ao que foi descrito para o elemento de comprimento L 2 temse sucessivamente x L a a a a x u 1 2 2 1 2 10 2 1 1 2 1 1 2 1 a L x a L x u x 11 2 1 1 2 1 1 2 1 a a L x L x u x 12 L x x N L x x N 1 2 1 1 2 1 2 1 13 43 Campo de deformações O campo de deformações na barra é definido do seguinte modo d x ε d u 14 Atendendo a 5 temse 2 2 1 1 x a N d x N x a d ε 15 Uma vez que os deslocamentos nodais a1 e a2 não dependem de x da derivação resulta 2 2 1 1 d x a d N d x a d N ε 16 que em notação matricial fica 2 1 2 1 a a d x d N d x d N ε 17 Elementos Finitos Unidimensionais Álvaro F M Azevedo 55 Designando por B a matriz d x d N d x d N B 2 1 18 e atendendo a 9 temse ε B a 19 Designando por L 1 o seguinte operador diferencial d x d L 1 20 a equação 14 escrevese ε L1 u 21 Atendendo a 7 temse ε L1 N a 22 Comparando 22 com 19 concluise que L N B 1 23 De acordo com 18 e com 6 para o caso da barra de comprimento L 2 os elementos da matriz B são os seguintes 2 1 2 1 2 1 d x N d d x N d 24 2 1 2 1 B 25 No caso da barra de comprimento L de 18 e 13 chegase a Elementos Finitos Unidimensionais Álvaro F M Azevedo 56 L d x N d L d x N d 1 1 2 1 26 L L B 1 1 27 De 9 19 e 27 concluise que no caso da barra de comprimento L se tem L L L a a a a L L B a 1 2 2 1 1 1 ε 28 Neste exemplo simples a expressão do campo de deformações corresponde ao que se considera habitualmente para uma barra sujeita a um esforço axial Uma vez que ε não depende da coordenada x este elemento finito apresenta deformação constante 44 Princípio dos trabalhos virtuais Considerese um corpo sujeito a um conjunto de forças de volume e de superfície que lhe provocam uma deformação Com base no seu estado de equilíbrio estático a configuração do corpo é modificada por um conjunto de deslocamentos muito pequenos e compatíveis com as condições fronteira que se designam deslocamentos virtuais O princípio dos trabalhos virtuais ou princípio dos deslocamentos virtuais estabelece que o trabalho realizado pelas tensões internas na deformação virtual do corpo é igual ao trabalho realizado pelas forças exteriores nos deslocamentos virtuais dos seus pontos de aplicação 41 42 De um modo mais simplista é comum afirmar que o trabalho interno de deformação é igual ao trabalho externo das forças aplicadas Trabalho Interno Trabalho Externo 29 Apresentase em seguida uma versão simplificada do princípio dos trabalhos virtuais PTV adaptada ao caso das barras sujeitas a deslocamentos e forças apenas axiais Nas expressões que se seguem o prefixo δ indica que os deslocamentos ou deformações são virtuais Elementos Finitos Unidimensionais Álvaro F M Azevedo 57 L T V T u p d L dV δ δε σ 30 Nesta expressão o vector δε apenas tem a componente correspondente à extensão segundo o eixo da barra o vector σ apenas contem a tensão normal na secção transversal da barra o campo de deslocamentos δ u e a acção exterior distribuída p apenas referem a componente segundo o eixo da barra ver a Figura 44 u x x 1 2 L x L 2 x L 2 p F2 F1 Fig 44 Elemento finito unidimensional sujeito a uma acção axial uniformemente distribuída Neste caso a expressão do PTV 30 passa a ser a seguinte L T V T u p d L dV δ δ ε σ 31 45 Matriz de rigidez e vector solicitação Com base no princípio dos trabalhos virtuais apresentado na secção anterior vaise em seguida proceder à dedução das expressões da matriz de rigidez e do vector solicitação que são utilizados no método dos deslocamentos Designando por A a área da secção transversal da barra temse A d x dV 32 Uma vez que o eixo da barra coincide com o eixo x temse d x d L 33 Elementos Finitos Unidimensionais Álvaro F M Azevedo 58 A equação 19 referida à deformação virtual é a seguinte B δ a δ ε 34 que é equivalente a T T T B δ a δ ε 35 A relação constitutiva ou relação tensãodeformação é neste caso ε σ D 36 apresentando a matriz de elasticidade D apenas um elemento que consiste no módulo de Young E Substituindo 19 em 36 temse σ D B a 37 A equação 7 referida à deformação virtual é a seguinte a N u δ δ 38 que é equivalente a T T T N a u δ δ 39 Substituindo todas estas equações em 31 passa a terse o PTV expresso por 2 2 2 2 L L T T L L T T p d x a N a B D B a A d x δ δ 40 Uma vez que os deslocamentos nodais não dependem de x podem passar para fora do integral 2 2 2 2 L L T T L L T T p d x N a B D B A d x a a δ δ 41 Elementos Finitos Unidimensionais Álvaro F M Azevedo 59 De acordo com o PTV a equação 41 é verdadeira para qualquer conjunto de deslocamentos virtuais concluindose assim que 2 2 2 2 L L T L L T p d x N B D B A d x a 42 Comparando esta equação com a relação de rigidez que é utilizada no método dos deslocamentos F K a 43 temse no caso da barra unidimensional 2 2 L L T B D B A d x K 44 2 2 L L T p d x N F 45 As expressões 4245 são aplicáveis quando as seguintes grandezas são variáveis ao longo da barra módulo de Young E área da secção transversal A e carga distribuída p Apresentase em seguida o desenvolvimento das expressões 44 e 45 para o caso de E A e p serem constantes 2 2 L L T B B d x E A K 46 Atendendo a 27 2 2 1 1 1 1 L L L d x L L L E A K 47 E A L E A L E A L E A L K 48 Elementos Finitos Unidimensionais Álvaro F M Azevedo 60 Neste caso simples os elementos da matriz de rigidez coincidem com os que se obtêm directamente pelo método dos deslocamentos Partindo de 45 temse neste caso em que p é constante 2 2 L L T N d x p F 49 Atendendo a 8 e a 13 temse 2 2 1 2 1 1 2 1 L L x d x L L x p F 50 2 2 L p L p F 51 Esta expressão também coincide com a que se obtém por processos mais simples 46 Elemento finito unidimensional com três nós Considerese o elemento finito unidimensional com três nós representado na Figura 45 cujo comprimento é L 2 u x x 1a 1 2 L 2 x 1 x 1 3 a2 3a x 0 Fig 45 Elemento finito unidimensional de três nós Elementos Finitos Unidimensionais Álvaro F M Azevedo 61 De um modo semelhante ao que foi apresentado na Secção 42 considerase que a função u x é aproximada pelo seguinte polinómio de segundo grau 2 2 1 0 c x c x c u x 52 Pretendese que a função 52 respeite nos nós os valores dos respectivos deslocamentos sendo 3 2 1 1 0 1 a u a u a u 53 Atendendo a 52 temse 3 2 2 1 0 2 2 2 1 0 1 2 2 1 0 1 1 0 0 1 1 a c c c a c c c a c c c 54 que é equivalente a 3 2 1 2 1 0 1 1 1 0 0 1 1 1 1 a a a c c c 55 Explicitando c0 c1 e c2 temse 3 2 1 2 1 0 50 1 50 50 0 50 0 1 0 a a a c c c 56 Substituindo as expressões de c0 c1 e c2 em 52 chegase a 2 3 2 1 3 1 2 50 50 50 50 x a a a x a a a u x 57 que é equivalente a 3 2 2 2 1 2 50 50 1 50 50 x a x a x x a x u x 58 Elementos Finitos Unidimensionais Álvaro F M Azevedo 62 Em notação matricial temse 3 2 1 2 2 2 50 50 1 50 50 a a a x x x x x u x 59 Considerando 3 2 1 3 2 1 a a a N x x N N x u x 60 temse x x x N x x N x x x N 2 1 2 1 1 2 1 2 1 2 3 2 2 2 1 61 Neste caso 3 2 1 3 2 1 N N N N x x N N x N 62 u N a 63 3 2 1 a a a a 64 Na Figura 46 estão representados os gráficos das funções N1x N2x e N3x indicadas em 61 Elementos Finitos Unidimensionais Álvaro F M Azevedo 63 1 1 1 1 1 1 1 1 1 x x x N1x N2x N3x Fig 46 Gráfico das funções N1x N2x e N3x Na Tabela 43 encontramse algumas características das funções de forma representadas na Figura 46 comparar com a Tabela 42 Tabela 43 Características das funções N1x N2x e N3x x 1 0 1 N1x 1 0 0 N2x 0 1 0 N3x 0 0 1 Elementos Finitos Unidimensionais Álvaro F M Azevedo 64 Generalizando a expressão 18 para o caso do elemento de três nós resulta d x d N d x d N d x d N B 3 2 1 65 Atendendo a 61 os elementos da matriz B são neste caso os seguintes 2 1 2 2 1 x x x B 66 O cálculo da matriz de rigidez K e do vector solicitação F pode ser efectuado por um processo semelhante ao indicado na Secção 45 não sendo aqui desenvolvido 47 Elemento finito unidimensional com substituição de variável Na Figura 47 encontrase representado um elemento finito unidimensional com três nós e geometria qualquer u x x 1a 1 2 1x x 3 2a 3a x 2x 3x x Fig 47 Elemento finito unidimensional de três nós com geometria arbitrária As coordenadas dos nós são 1x 2x e 3x Tal como nos casos descritos anteriormente E representa o módulo de Young A é a área da secção transversal e p é a acção axial distribuída Todas estas grandezas podem eventualmente depender de x É possível calcular a matriz de rigidez K e o vector solicitação F com 44 e 45 utilizando como variável a coordenada x Contudo e tendo em vista a generalização deste estudo aos casos bidimensionais e tridimensionais vai ser efectuada uma substituição de variável do tipo Elementos Finitos Unidimensionais Álvaro F M Azevedo 65 x s x 67 A função x s neste caso seleccionada corresponde a uma interpolação coincidente com a que foi efectuada na Secção 46 para a função deslocamento u x em que foi utilizada a interpolação 60 conjuntamente com as funções de forma 61 3 2 1 3 2 1 x x x N s N s N s x s 68 3 3 2 2 1 1 N s x N s x N s x x s 69 s s s N s s N s s s N 2 1 2 1 1 2 1 2 1 2 3 2 2 2 1 70 De um modo semelhante ao que se verificou em 53 temse 3 2 1 1 0 1 x x x x x x 71 A substituição de variável 67 encontrase esquematizada na Figura 48 Elementos Finitos Unidimensionais Álvaro F M Azevedo 66 x 1 2 1x x 3 x 2x 3x x s 1 2 s 1 3 s 0 s 1 Fig 48 Substituição da variável x Após a substituição da variável x o integral 44 passa a ser 1 1 d s d s B D B A d x K T 72 com D B A e dxds dependentes da nova variável s Se não forem constantes D que coincide com E e A são interpolados com as mesmas funções de forma que foram utilizadas para interpolar as coordenadas dos nós ie a interpolação é efectuada tal como em 69 3 3 2 2 1 1 N s E N s E N s E E s 73 3 3 2 2 1 1 N s A N s A N s A A s 74 Nestas funções iE e iA são os valores no nó i do módulo de Young e da área da secção transversal A expressão de dxds que se passa a designar por J obtémse por derivação de 69 resultando 3 3 2 2 1 1 d s x d N x d s d N d s x d N d s d x J 75 Por derivação de 70 em ordem a s obtémse Elementos Finitos Unidimensionais Álvaro F M Azevedo 67 2 1 2 2 1 3 2 1 s d s N d s d s N d s d s N d 76 ficando 3 2 1 2 1 2 2 1 x s s x x s d s d x J 77 Para avaliar o integral 72 é ainda necessário definir a matriz B em função de s Atendendo à adaptação de 18 ao elemento de três nós que foi também utilizada em 65 existe a necessidade de calcular as derivadas das funções de forma em ordem a x mas expressas em função de s Com este objectivo e uma vez que as funções de forma Ni 61 dependem de x que por sua vez depende de s 69 temse recorrendo à regra da cadeia d s d x d x d N d s d N d s N x s d i i i 78 d x d N d s d x d s d N i i 79 Multiplicando ambos os membros de 79 pela inversa de dxds resulta d s d N d s d x d x d N i i 1 80 Uma vez que dxds é um escalar pode escreverse d s d x d s N d d x N d i i 81 Elementos Finitos Unidimensionais Álvaro F M Azevedo 68 sendo de acordo com 76 e 77 3 2 1 1 2 1 2 2 1 2 1 x s s x x s s d x N d 82 3 2 1 2 2 1 2 2 1 2 x s s x x s s d x N d 83 3 2 1 3 2 1 2 2 1 2 1 x s s x x s s d x N d 84 A matriz B apresenta os seguintes componentes 2 1 2 2 1 1 s s s J B 85 Depois de definidos todos os componentes da função integranda de 72 é possível efectuar as seguintes simplificações 1 1 E A J B B d s K T 86 sendo d x d N d x d N d x N d d x N d d x N d d x N d B B T 3 2 1 3 2 1 1 3 3 1 87 Elementos Finitos Unidimensionais Álvaro F M Azevedo 69 d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x d N d x N d B B T 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 3 3 88 Atendendo a 81 e ao facto de ser J dxds temse d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s d N d s N d J B B T 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 2 1 89 A expressão genérica do elemento Kij da matriz K é 1 1 d s d s d N d s d N J E A K j i ij 90 Como exemplo apresentase em seguida a expressão do elemento K13 da matriz de rigidez do elemento finito de acordo com 90 e 76 1 1 13 2 1 2 1 d s s s J E A K 91 Considerese agora um caso particular de uma barra de comprimento total L e nó 2 centrado ver a Figura 47 com 2 0 2 3 2 1 L x x L x 92 Neste caso particular a expressão de J calculada com 77 não depende de s sendo Elementos Finitos Unidimensionais Álvaro F M Azevedo 70 2 L d s d x J 93 Se além de J ser constante E e A também forem constantes é simples calcular o integral 91 resultando L E A K 3 13 1 94 Apresentase em seguida um exemplo numérico em que o nó 2 não se encontra centrado no elemento finito de três nós ver a Figura 47 05 03 02 3 2 1 x x x 95 Neste caso concreto a expressão de J calculada com 77 é 2 3 s d s d x J 96 Supondo E e A constantes temse de acordo com 91 1 1 2 13 2 3 4 1 d s s s E A K 97 Na prática é conveniente resolver os integrais 90 e 97 recorrendo a uma técnica de integração numérica que será descrita no Capítulo 5 48 Considerações finais A formulação pelo MEF aqui efectuada no âmbito de um problema muito simples serve como introdução às técnicas que se aplicam em meios contínuos com duas ou três dimensões de que são exemplo os estados planos de tensão as cascas e os sólidos Muitas das expressões matriciais que aqui foram apresentadas coincidem com as que Elementos Finitos Unidimensionais Álvaro F M Azevedo 71 surgem nos casos mais genéricos sendo apenas necessário redefinir as dimensões e os elementos dos vectores e das matrizes BIBLIOGRAFIA 41 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 42 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 Elementos Finitos Unidimensionais Álvaro F M Azevedo 72 73 CAPÍTULO 5 QUADRATURA DE GAUSS Muitos dos integrais que é necessário calcular no âmbito da aplicação do Método dos Elementos Finitos MEF não são triviais ie ou a primitiva da função integranda não existe explicitamente ou é demasiado complicada para viabilizar a sua utilização prática Por este motivo é essencial recorrer a técnicas de integração numérica que também recebem a designação de quadratura Neste capítulo é descrita e justificada a quadratura de Gauss por ser a mais utilizada no âmbito do MEF 51 51 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada no estudo da quadratura de Gauss Tabela 51 Simbologia relativa à quadratura de Gauss c Coeficiente de um termo de um polinómio I Valor exacto do integral J Valor do integral calculado de acordo com a quadratura de Gauss P Posição de um ponto de Gauss ou ponto de amostragem W Peso weight associado a um ponto de Gauss ou ponto de amostragem n Número de pontos de Gauss utilizados numa direcção p Grau de um polinómio 52 Integração de uma função polinomial Na Figura 51 encontrase representada uma função polinomial de grau 5 cuja expressão genérica é a seguinte Quadratura de Gauss Álvaro F M Azevedo 74 5 5 4 4 3 3 2 2 1 0 c x c x c x c x c x c f x 1 f x x 1 1 1P P2 3 P Fig 51 Função polinomial de grau 5 O integral exacto do polinómio 1 no intervalo 11 é 1 1 f x d x I 2 1 1 5 5 4 4 3 3 2 2 1 0 d x c x c x c x c x c x c I 3 4 2 0 5 2 3 2 2 c c c I 4 Para facilitar a sua comparação com uma expressão que vai ser em seguida apresentada o segundo membro de 4 é rescrito da seguinte forma 5 4 3 2 1 0 0 5 2 0 3 2 0 1 2 c c c c c c I 5 Suponhase agora que se pretende avaliar o integral de f x por intermédio do somatório de avaliações da função f x em determinados locais multiplicadas por adequados pesos No caso do polinómio de grau 5 indicado em 1 será adiante mostrado que para se obter um resultado exacto se deve avaliar a função f x em três pontos de amostragem Pi e multiplicar cada um desses valores por pesos Wi ver a Figura 51 O integral avaliado desta forma é designado por J sendo Quadratura de Gauss Álvaro F M Azevedo 75 3 3 2 2 1 1 W f P f P W W f P J 6 Mais adiante será deduzido o valor adequado para os seguintes parâmetros posição dos pontos de amostragem P1 P2 e P3 em que a função f x deve ser avaliada ver a Figura 51 valores dos pesos W1 W2 e W3 Uma vez que f x é um polinómio do tipo 1 a expressão 6 passa a ser 5 3 5 4 3 4 3 3 3 2 3 2 3 1 0 3 5 2 5 4 2 4 3 2 3 2 2 2 2 1 0 2 5 1 5 4 1 4 3 1 3 2 1 2 1 1 0 1 c P c P c P c P c P c W c P c P c P c P c P c W c P c P c P c P c P W c J 7 No segundo membro de 7 podemse colocar em evidência os coeficientes ci resultando 5 5 3 3 5 2 2 5 1 1 4 4 3 3 4 2 2 4 1 1 3 3 3 3 3 2 2 3 1 1 2 2 3 3 2 2 2 2 1 1 1 3 3 2 2 1 1 0 3 2 1 c W P W P P W c W P W P P W c W P W P P W c W P W P P W W P c W P P W c W W W J 8 Neste exemplo relativo ao polinómio de grau 5 indicado em 1 pretendese que a expressão de J 8 seja exactamente igual à de I 5 I J 9 Igualando os segundos membros de 5 e de 8 resulta Quadratura de Gauss Álvaro F M Azevedo 76 5 5 3 3 5 2 2 5 1 1 4 4 3 3 4 2 2 4 1 1 3 3 3 3 3 2 2 3 1 1 2 2 3 3 2 2 2 2 1 1 1 3 3 2 2 1 1 0 3 2 1 5 4 3 2 1 0 0 5 2 0 3 2 0 1 2 c W P W P P W c W P W P P W c W P W P P W c W P W P P W W P c W P P W c W W W c c c c c c 10 Uma vez que os coeficientes ci são arbitrários para que a igualdade 10 se verifique sempre é suficiente que 0 2 5 0 2 3 0 1 2 5 3 3 5 2 2 5 1 1 4 3 3 4 2 2 4 1 1 3 3 3 3 2 2 3 1 1 2 3 3 2 2 2 2 1 1 3 3 2 2 1 1 3 2 1 W P W P W P W P W P W P W P W P W P W P W P W P W P W P W P W W W 11 Para obter os valores de P1 P2 P3 W1 W2 e W3 resolvese o sistema de seis equações não lineares a seis incógnitas 11 A respectiva solução é 0 55555 55556 9 5 0 88888 88889 9 8 0 55555 55556 9 5 0 77459 66692 5 3 0 0 0 77459 66692 5 3 3 2 1 3 2 1 W W W P P P 12 O valor exacto do integral de um polinómio de grau 5 no intervalo 11 pode ser obtido com 5 3 9 5 0 9 8 5 3 9 5 f f f J I 13 No caso de a função f x ser genérica ie não polinomial ou polinomial de grau superior a 5 a expressão 13 fornece um valor aproximado do integral I 2 Quadratura de Gauss Álvaro F M Azevedo 77 5 3 9 5 0 9 8 5 3 9 5 1 1 f f f f x d x 14 O valor do integral calculado com o segundo membro de 14 é tanto mais correcto quanto mais a função f x se aproximar de um polinómio do tipo 1 Se se desejar um valor mais correcto para o integral existe a possibilidade de se utilizar mais pontos de amostragem Pi e correspondentes pesos Wi Os pontos de amostragem também são designados por pontos de Gauss O estudo que foi aqui realizado com um polinómio de grau 5 pode ser feito de um modo semelhante com polinómios de qualquer grau Na Tabela 52 apresentase os resultados que se obtêm quando se faz o estudo com polinómios de grau 1 grau 3 grau 5 e grau 7 Em 52 encontrase uma tabela que fornece os valores das posições dos pontos de amostragem e dos pesos para um número de pontos de Gauss no intervalo 110 Com base na Tabela 52 podemse extrair as seguintes conclusões com n pontos de Gauss obtémse o valor exacto do integral de um polinómio de grau p 2 n 1 ou inferior quando se pretende a solução exacta do integral de um polinómio de grau p o número de pontos de Gauss que se tem de utilizar é n p 1 2 ou superior Nota quando p é par devese substituir o seu valor pelo número ímpar imediatamente superior Nota o intervalo de integração de todos os integrais referidos no âmbito da quadratura de Gauss é o intervalo 11 Quadratura de Gauss Álvaro F M Azevedo 78 Tabela 52 Posições dos pontos de amostragem e respectivos pesos Número de pontos de Gauss Grau do polinómio que é possível integrar de um modo exacto Posições dos pontos de Gauss e respectivos pesos n p 2 n 1 Pi Wi 1 1 2 0 1 1 W P 2 3 1 1 3 1 3 1 2 1 2 1 W W P P 3 5 5 9 8 9 9 5 5 3 0 5 3 3 2 1 3 2 1 W W W P P P 4 7 034785 48451 065214 51549 065214 51549 034785 48451 086113 63116 0 33998 10436 0 33998 10436 63116 086113 4 3 2 1 4 3 2 1 W W W W P P P P Para justificar a expressão p 2 n 1 ver a Tabela 52 é suficiente considerar o seguinte sugerese que se acompanhem as seguintes considerações com o exemplo do polinómio de grau p 5 atrás descrito suponhase que se pretende integrar de um modo exacto um polinómio de grau p sendo p um número ímpar Quadratura de Gauss Álvaro F M Azevedo 79 o número de coeficientes ci no polinómio de grau p é igual a p 1 uma vez que existem p 1 coeficientes ci o sistema de equações não lineares 11 vai ter p 1 equações para que o sistema de equações 11 possa ser resolvido o número de incógnitas deve ser também p 1 uma vez que as incógnitas são as posições dos pontos de Gauss e respectivos pesos P1 P2 P3 W1 W2 W3 o número de pontos de Gauss n tem de ser metade do número de incógnitas p 1 ie n p 1 2 nesta expressão podese explicitar p resultando p 2 n 1 que é o resultado que se pretendia demonstrar Qualquer que seja o valor de n o valor de p que se obtém é sempre um número ímpar É por este motivo que conforme foi atrás referido se deve passar p para o valor ímpar imediatamente superior quando se utiliza a expressão n p 1 2 e o valor de p é par A expressão genérica da quadratura de Gauss com n pontos é n i i i W f P J 1 15 53 Integrais múltiplos Apresentase em seguida a adaptação da integração numérica descrita na secção anterior ao caso do integral duplo 1 1 1 1 f x y d x d y I 16 Considerando em primeiro lugar o integral em ordem a x temse de acordo com 15 1 1 1 d y W f P y J x n i i i 17 sendo nx o número de pontos de Gauss utilizados na direcção x Quadratura de Gauss Álvaro F M Azevedo 80 Considerando que a função integranda de 17 é uma função gy temse 1 1 g y d y J 18 com x n i i i W f P y y g 1 19 Substituindo agora o integral em ordem a y em 18 por um somatório do tipo 15 resulta y n j j j W g P J 1 20 sendo ny o número de pontos de Gauss utilizados na direcção y Atendendo a 19 a expressão 20 passa a ser y x n j n i j i i j W f P P W J 1 1 21 que é equivalente a x y n i n j j i j i f P P W W J 1 1 22 O número de pontos de Gauss associados à direcção x nx pode ser diferente do número de pontos de Gauss associados à direcção y ny A selecção destes números deve atender ao modo como a função f xy varia com x e com y Assim se na direcção x a função f xy se assemelhar a um polinómio de grau 5 e na direcção y a um de grau 7 deve ser nx 3 e ny 4 ver a Tabela 52 No caso do integral triplo podese generalizar 22 resultando x y z n i n j n k k j i k j i f P P P W W W x y z d x d y d z f 1 1 1 1 1 1 1 1 1 23 Quadratura de Gauss Álvaro F M Azevedo 81 No caso do integral do produto das funções f e g temse x y z n i n j n k k j i k j i k j i g P P P f P P P W W W x y z g x y z d x d y d z f 1 1 1 1 1 1 1 1 1 24 o que permite uma avaliação sequencial de f e g no ponto de Gauss Pi Pj Pk Esta consideração é extensiva a qualquer combinação de funções eg adição divisão etc Quando se tem por exemplo o integral de um produto de matrizes podese avaliar cada uma das matrizes em cada ponto de Gauss e só em seguida fazer o produto matricial Assim se evita ter de explicitar a função que resulta do produto matricial de diversas funções 54 Considerações finais O procedimento de integração numérica genericamente designado quadratura de Gauss tem como principal vantagem o facto de poder ser facilmente incluído num programa de computador destinado à análise de estruturas pelo MEF A principal dificuldade associada à sua utilização reside na necessidade de escolher um número de pontos de Gauss adequado à precisão pretendida BIBLIOGRAFIA 51 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 52 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 Quadratura de Gauss Álvaro F M Azevedo 82 83 CAPÍTULO 6 ESTADO PLANO DE TENSÃO Neste capítulo é descrita com pormenor a formulação de elementos finitos destinados à discretização de problemas de análise de estruturas que se enquadram no caso particular designado Estado Plano de Tensão 61 Apresentase em primeiro lugar o caso do elemento finito quadrado de dimensões fixas seguindose o elemento finito rectangular L1xL2 e por último o caso mais geral de geometria arbitrária A formulação aqui descrita baseiase no método dos deslocamentos e na discretização do domínio em elementos finitos de n nós apresentando algumas semelhanças com o que foi descrito no Capítulo 4 61 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação do método dos elementos finitos Tabela 61 Simbologia relativa ao método dos elementos finitos L Dimensão do elemento finito n Número de nós do elemento finito x Coordenada cartesiana u Campo de deslocamentos a Deslocamento nodal h Espessura do elemento finito laminar x Coordenada cartesiana de um nó de um elemento finito m Número de direcções consideradas no estado plano de tensão m 2 N Função interpoladora ou função de forma Estado Plano de Tensão Álvaro F M Azevedo 84 c Coeficiente de um termo de um polinómio p Número de graus de liberdade do elemento finito p n x m ε Extensão γ Distorção L Operador diferencial q Número de componentes do vector ε e do vector σ B Matriz de deformação V Volume do elemento finito laminar σ Tensão normal τ Tensão tangencial p Acção exterior distribuída por unidade de comprimento S Superfície do elemento finito laminar E Módulo de elasticidade ou módulo de Young ν Coeficiente de Poisson D Matriz de elasticidade ε σ D K Matriz de rigidez do elemento finito no referencial geral F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial geral s Coordenada local curvilínea s Coordenada local de um nó de um elemento finito NV Vector das funções interpoladoras ou funções de forma J Jacobiano da transformação E Módulo de elasticidade num nó do elemento finito ν Coeficiente de Poisson num nó do elemento finito h Espessura do elemento finito num nó Estado Plano de Tensão Álvaro F M Azevedo 85 P Posição de um ponto de Gauss ou ponto de amostragem W Peso weight associado a um ponto de Gauss ou ponto de amostragem nGPi Número de pontos de Gauss associado à direcção si J Valor do integral calculado de acordo com a quadratura de Gauss 62 Funções interpoladoras ou funções de forma Na Figura 61 encontrase representado um elemento finito quadrado com quatro nós e com dimensões L1xL2 2x2 u1 x1 x2 x1 1 2 L1 2 L2 2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 h x1 x2 Fig 61 Elemento finito quadrado de quatro nós As coordenadas dos nós são armazenadas na matriz x cujo elemento genérico ijx corresponde à coordenada cartesiana do nó i segundo a direcção xj 1 1 1 1 1 1 1 1 42 41 32 31 22 21 12 11 x x x x x x x x x 1 De acordo com a simbologia atrás apresentada a matriz x tem dimensões nxm A espessura do elemento finito laminar representado na Figura 61 é designada por h que pode também ser uma função de x1 e de x2 Estado Plano de Tensão Álvaro F M Azevedo 86 A função u x corresponde ao campo de deslocamentos verificandose o seguinte 2 1 2 2 1 1 x x u u x x u x 2 Cada uma das componentes de u x é interpolada separadamente com base em funções de forma Ni x1 x2 e nos deslocamentos dos nós ver a Figura 61 41 2 1 4 31 2 1 3 21 2 1 2 11 2 1 1 2 1 1 a x x N a x x N a x x N a x x N x u x 3 42 2 1 4 32 2 1 3 22 2 1 2 12 2 1 1 2 1 2 a x x N a x x N a x x N a x x N x x u 4 Em 3 e 4 bem como na Figura 61 aij corresponde ao deslocamento do nó i segundo a direcção xj Notese que o número de funções de forma Ni coincide com o número de nós do elemento finito n As considerações que se seguem serão apenas efectuadas com a componente u1 do campo de deslocamentos A sua extensão à componente u2 seria trivial A função u1x1 x2 deve assumir nos nós os valores nodais do campo de deslocamentos Atendendo às coordenadas dos nós indicadas em 1 pretendese que 41 1 31 1 21 1 11 1 1 1 1 1 1 1 1 1 a u a u a u a u 5 Para que as condições expressas em 5 sejam respeitadas as funções de forma a utilizar em 3 devem possuir as características indicadas na Tabela 62 Estado Plano de Tensão Álvaro F M Azevedo 87 Tabela 62 Características das funções Ni x1 x2 Nó 1 2 3 4 x1 x2 1 1 1 1 1 1 1 1 N1 x1 x2 1 0 0 0 N2 x1 x2 0 1 0 0 N3 x1 x2 0 0 1 0 N4 x1 x2 0 0 0 1 As seguintes funções polinomiais respeitam as condições indicadas na Tabela 62 que consistem no facto da a função Ni ter de assumir um valor unitário no nó i e um valor nulo nos restantes nós 4 1 1 4 1 1 4 1 1 4 1 1 2 1 2 1 4 2 1 2 1 3 2 1 2 1 2 2 1 2 1 1 x x x x N x x x x N x x x x N x x x x N 6 A função N1 x1 x2 pode tomar a seguinte forma 2 1 2 1 2 1 1 0 25 0 25 0 25 0 25 x x x x x x N 7 Um polinómio de segundo grau completo tem a seguinte expressão genérica 2 2 5 2 1 4 2 1 3 2 2 1 1 0 1 2 c x c x x c x c x c x c f x x 8 Comparando 7 com 8 verificase que a função de forma N1 x1 x2 é um polinómio de segundo grau incompleto porque lhe faltam os termos que em 8 se encontram sublinhados Considerações idênticas poderiam ser feitas em relação às restantes funções de forma Estado Plano de Tensão Álvaro F M Azevedo 88 Armazenando os deslocamentos nodais da seguinte forma 42 41 32 31 22 21 12 11 a a a a a a a a a 9 temse atendendo a 3 e a 4 42 41 32 31 22 21 12 11 4 3 2 1 4 3 2 1 2 1 0 0 0 0 0 0 0 0 a a a a a a a a N N N N N N N N u u 10 que em notação matricial se reduz a 1 1 p p m m a N u 11 sendo p n x m no caso da Figura 61 p 4 x 2 8 A matriz N é 4 3 2 1 4 3 2 1 0 0 0 0 0 0 0 0 N N N N N N N N N 12 Os gráficos das funções Ni x1 x2 definidas em 6 encontramse representados na Figura 62 ver também a Figura 61 No caso do elemento finito rectangular de dimensões L1 x L2 representado na Figura 63 as funções de forma seriam Estado Plano de Tensão Álvaro F M Azevedo 89 2 2 1 1 2 1 2 1 4 2 2 1 1 2 1 2 1 3 2 2 1 1 2 1 2 1 2 2 2 1 1 2 1 2 1 1 2 2 1 2 2 1 2 2 1 2 2 1 x L x L L L x x N x L x L L L x x N x L x L L L x x N x L x L L L x x N 13 x1 x2 N1 x1 x2 N4 x1 x2 N3 x1 x2 N2 4 1 2 3 Fig 62 Gráficos das funções Ni x1 x2 para um elemento de dimensões L1xL2 2x2 Estado Plano de Tensão Álvaro F M Azevedo 90 u1 x1 x2 x1 1 2 L1 L2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 h x1 x2 Fig 63 Elemento finito rectangular de quatro nós 63 Campo de deformações O campo de deformações num estado plano de tensão é definido do seguinte modo 61 2 1 1 2 2 1 12 2 1 0 0 u u x x x x γ ε ε 14 ou de um modo mais compacto 1 1 m q m q u L ε 15 Em 15 q é o número de componentes do vector ε que são neste caso três e L é o seguinte operador diferencial 1 2 2 1 0 0 x x x x L 16 Estado Plano de Tensão Álvaro F M Azevedo 91 Substituindo 11 em 15 temse 1 1 p p m q m q a N L ε 17 Designando por B o produto L N p m q m q p N L B 18 a expressão 17 passa a 1 1 p p q q a B ε 19 sendo de acordo com 18 16 e 12 4 3 2 1 4 3 2 1 1 2 2 1 0 0 0 0 0 0 0 0 0 0 N N N N N N N N x x x x B 20 1 4 2 4 1 3 2 3 1 2 2 2 1 1 2 1 2 4 2 3 2 2 2 1 1 4 1 3 1 2 1 1 0 0 0 0 0 0 0 0 x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N B 21 No caso do elemento com dimensões L1xL2 2x2 a matriz B é constituída pelas derivadas de 6 de acordo com 21 2 1 2 1 2 1 2 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 4 1 x x x x x x x x x x x x x x x x B 22 No caso do elemento de dimensões L1xL2 a matriz B é constituída pelas derivadas de 13 de acordo com 21 Estado Plano de Tensão Álvaro F M Azevedo 92 2 2 1 1 2 2 1 1 1 1 1 1 2 2 2 2 2 1 2 2 2 2 2 0 0 0 2 0 0 2 0 0 0 2 1 x L x L x L x L x L x L x L x L L L B L L L L L L L L 23 64 Princípio dos trabalhos virtuais Considerese um estado plano de tensão constituído por um elemento finito formulado de acordo com o que foi exposto nas secções anteriores Supondo que apenas existem acções distribuídas por unidade de comprimento na periferia do elemento finito do Princípio dos Trabalhos Virtuais PTV que foi exposto no Capítulo 4 resulta a seguinte equação L T V T u p d L dV δ δε σ 24 Nesta expressão o vector δε apresenta componentes em correspondência com o vector ε definido em 14 e 15 O vector σ é o seguinte 12 2 1 τ σ σ σ 25 65 Matriz de rigidez e vector solicitação Com base no princípio dos trabalhos virtuais referido na secção anterior vaise em seguida proceder à dedução das expressões da matriz de rigidez e do vector solicitação que são utilizados no método dos deslocamentos aplicado à análise de um estado plano de tensão Designando por h a espessura do elemento finito temse em 24 h d S dV 26 Estado Plano de Tensão Álvaro F M Azevedo 93 em que dS representa o elemento de superfície A equação 19 referida à deformação virtual é a seguinte B δ a δ ε 27 que é equivalente a T T T B δ a δ ε 28 A relação entre tensões e deformações é para um estado plano de tensão e no caso dos materiais isotrópicos 61 12 2 1 2 2 2 2 12 2 1 2 1 0 0 0 1 1 0 1 1 γ ε ε ν ν ν ν ν ν ν τ σ σ E E E E E 29 ou de um modo mais compacto ε σ D 30 sendo a matriz de elasticidade D a seguinte ν ν ν ν ν ν ν 2 1 0 0 0 1 1 0 1 1 2 2 2 2 E E E E E D 31 A matriz de elasticidade D depende do módulo de Young E e do coeficiente de Poisson ν Substituindo 19 em 30 temse σ D B a 32 A equação 11 referida à deformação virtual é a seguinte Estado Plano de Tensão Álvaro F M Azevedo 94 a N u δ δ 33 que é equivalente a T T T N a u δ δ 34 Substituindo todas estas equações em 24 passa a terse o PTV expresso por L T T S T T p d L a N a B D B a h d S δ δ 35 Uma vez que dS dx1 dx2 e os deslocamentos nodais não dependem das variáveis x1 e x2 os vectores δ aT e a podem passar para fora do integral L T T S T T p d L N a B D B h d S a a δ δ 36 De acordo com o PTV a equação 36 é verdadeira para qualquer conjunto de deslocamentos virtuais concluindose assim que L T S T p d L N B D B h d S a 37 Comparando esta equação com a relação de rigidez que é utilizada no método dos deslocamentos F K a 38 temse no caso do estado plano de tensão S T B D B h d S K 39 L T p d L N F 40 O vector a encontrase definido em 9 Nas expressões 3740 admitese que as seguintes grandezas podem não ser constantes no domínio de integração módulo de Young E coeficiente de Poisson ν espessura h e carga distribuída p Estado Plano de Tensão Álvaro F M Azevedo 95 No caso do elemento finito rectangular representado na Figura 63 a expressão da matriz de rigidez 39 passa a ser 2 2 2 2 2 1 2 2 1 1 L L L L T B D B h d x d x K 41 A matriz B corresponde à expressão 23 e no caso dos materiais isotrópicos a matriz D é dada por 31 Uma vez que h é um escalar as dimensões da matriz K coincidem com as do produto D B B T L L p q q q q p T p p B D B K 42 No caso do elemento finito de quatro nós temse L L 3 8 3 3 8 3 8 8 B D B K T 43 Atendendo a 40 as dimensões do vector solicitação F coincidem com as do produto p N T L L 1 1 m m p T p p N F 44 No caso do elemento finito de quatro nós temse L L 2 1 8 2 8 1 p N F T 45 651 Cálculo de um elemento da matriz de rigidez Apresentase em seguida o cálculo do elemento K58 da matriz de rigidez do elemento finito representado na Figura 61 com E 200 000 MPa ν 0 e h 03 m De acordo com 41 temse Estado Plano de Tensão Álvaro F M Azevedo 96 1 1 1 1 2 1 B D B h d x d x K T 46 O cálculo de D B B T pode ser efectuado com base nos somatórios correspondentes aos produtos matriciais q k q p pj kp ki ij T B D B D B B 1 1 47 sendo q 3 Para calcular K58 é suficiente desenvolver os somatórios para o caso i 5 j 8 3 1 3 1 8 5 58 k p p kp k T D B B B D B 48 3 1 38 3 5 28 2 5 18 1 5 k k k k k k k B D B B D B D B B 49 38 33 35 38 23 25 38 13 15 28 32 35 28 22 25 28 12 15 18 31 35 18 21 25 18 11 15 B D B B D B D B B B D B B D B D B B B D B B D B D B B 50 Consultando as matrizes B 22 e D 31 com ν 0 verificase facilmente que neste exemplo só o último monómio de 50 é não nulo Assim temse 1 4 2 3 38 33 35 58 2 x N E x N B D B D B B T 51 Atendendo a 22 e 31 e ao facto de ser E 200 000 chegase a 4 1 100 000 4 1 2 1 58 x x D B B T 52 2 1 58 1 6 250 1 x x D B B T 53 Atendendo a 46 e ao facto de ser h 03 m temse Estado Plano de Tensão Álvaro F M Azevedo 97 1 1 1 1 2 1 58 58 h d x d x B D B K T 54 1 1 1 1 2 1 2 1 58 1 1 875 1 d x d x x x K 55 MN m K 7 500000 000 58 56 652 Cálculo do vector solicitação correspondente a uma carga distribuída Na Figura 64 encontrase representado o elemento finito da Figura 61 sujeito a uma carga distribuída no bordo 23 x1 1 2 L1 2 L2 2 3 4 x2 2 MNm 3 MNm 9 MNm 4 MNm p1 p2 dL Fig 64 Elemento finito sujeito a uma carga distribuída As forças nodais equivalentes à acção distribuída no bordo calculamse com a expressão 40 que aqui se reproduz L T p d L N F 57 sendo N a matriz 12 e p o seguinte vector 2 1 p p p 58 Estado Plano de Tensão Álvaro F M Azevedo 98 Neste exemplo dL coincide com dx2 e todos os pontos do domínio de integração apresentam coordenada x1 1 Assim o integral de linha 57 passa a ser 1 1 p d 2x N F T 59 e nos elementos da matriz N que são as funções de forma 6 devese substituir x1 por 1 obtendose 0 1 2 1 1 2 1 1 0 1 2 4 2 4 2 1 4 2 2 3 2 3 2 1 3 2 2 2 2 2 2 1 2 2 1 2 1 2 1 1 x N x N x x N x x N x N x x N x x N x N x x N x N x N x x N 60 Atendendo a 12 e a 58 o produto p N T que figura em 59 é o seguinte 2 1 4 4 3 3 2 2 1 1 0 0 0 0 0 0 0 0 p p N N N N N N N N N T p 61 Para os valores das acções indicados na Figura 64 são as seguintes as expressões das funções p1 e p2 2 2 2 2 2 1 3 6 3 x x p x p x 62 Com base em 61 60 e 62 temse Estado Plano de Tensão Álvaro F M Azevedo 99 2 2 2 2 2 2 6 3 3 0 0 0 0 2 1 0 0 2 1 2 1 0 0 2 1 0 0 0 0 x x x x x x N T p 63 resultando de 59 0 0 7 000 000 3 333 333 5 000 000 2 666 667 0 0 42 41 32 31 22 21 12 11 F F F F F F F F F 64 Em 64 Fij representa a componente de F que está associada ao nó i e que actua na direcção xj Nos nós 1 e 4 são nulas as componentes da força nodal equivalente à carga distribuída no bordo 23 Neste exemplo simples os valores indicados em 64 coincidem com as reacções que se obteriam numa viga simplesmente apoiada carregada com as cargas trapezoidais da Figura 64 66 Caso geral com substituição de variáveis O estudo apresentado nas secções anteriores e que se encontra limitado a um elemento quadrado de dimensões 2x2 pode ser facilmente estendido a elementos rectangulares de dimensões L1xL2 Toda a sua formulação seria uma extensão trivial do que foi atrás Estado Plano de Tensão Álvaro F M Azevedo 100 apresentado Nesta secção é desenvolvido um elemento finito quadrilátero de geometria arbitrária destinado à discretização de estados planos de tensão ver a Figura 65 u1 x1 x2 x1 1 2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 h x1 x2 Fig 65 Elemento finito quadrilátero de quatro nós com geometria arbitrária As coordenadas dos nós são armazenadas na matriz x cujo elemento genérico ijx corresponde à coordenada cartesiana do nó i segundo a direcção xj 42 41 32 31 22 21 12 11 x x x x x x x x x 65 De acordo com a simbologia apresentada na Secção 61 a matriz x tem dimensões nxm A espessura do elemento finito laminar representado na Figura 65 é designada por h que pode também ser uma função de x1 e de x2 A determinação da matriz de rigidez do elemento finito com a expressão 39 requer neste caso o cálculo de um integral duplo com um domínio de integração S que corresponde a um quadrilátero irregular de geometria definida pelos quatro nós do elemento Tendo em vista a sistematização deste processo de modo a facilitar a sua programação em computador revelase muito vantajoso efectuar a seguinte substituição das variáveis x1 e x2 Estado Plano de Tensão Álvaro F M Azevedo 101 2 1 2 2 2 1 1 1 s s x x x s s x 66 Na Figura 66 encontrase representado o novo domínio de integração que corresponde ao intervalo 1 1 quer para a variável s1 quer para s2 x1 1 2 3 4 x2 1 2 3 4 s1 s2 1 1 1 1 Fig 66 Substituição das variáveis x1 e x2 Os valores nodais das coordenadas s1 e s2 são os seguintes 1 1 1 1 1 1 1 1 42 41 32 31 22 21 12 11 s s s s s s s s s 67 De acordo com 66 a cada ponto s1 s2 corresponde um ponto x1 x2 A passagem do sistema de coordenadas s para o sistema de coordenadas x é efectuada com uma interpolação semelhante à que foi efectuada na Secção 62 para o campo de deslocamentos De acordo com 3 e 4 temse 41 2 1 4 31 2 1 3 21 2 1 2 11 2 1 1 2 1 1 x s s N x s s N x s s N x N s s x s s 68 42 2 1 4 32 2 1 3 22 2 1 2 12 2 1 1 2 1 2 x s s N x s s N x s s N x N s s s s x 69 No sistema de coordenadas s as funções de forma coincidem com as que foram descritas na Secção 62 bastando substituir em 6 x1 por s1 e x2 por s2 resultando Estado Plano de Tensão Álvaro F M Azevedo 102 4 1 1 4 1 1 4 1 1 4 1 1 2 1 2 1 4 2 1 2 1 3 2 1 2 1 2 2 1 2 1 1 s s s s N s s s s N s s s s N s s s s N 70 Tal como no caso do campo de deslocamentos ao atribuir a s1 s2 os valores nodais indicados em 67 obtêmse em 68 e 69 as coordenadas dos nós 65 Por exemplo para s1 s2 1 1 a função N3 vale um e as restantes são nulas obtendose em 68 31 1 11 x x e em 69 32 2 11 x x As equações 68 e 69 podem ser colocadas em forma matricial do seguinte modo 4 3 2 1 42 32 22 12 41 31 21 11 2 1 N N N N x x x x x x x x x x 71 ou 1 1 n V n m T m N x x 72 sendo 2 1 x x x 73 4 3 2 1 N N N N N V 74 Em 72 x é a matriz nxm definida em 65 Após a substituição de variáveis indicada em 66 o integral 39 passa a ser Estado Plano de Tensão Álvaro F M Azevedo 103 1 1 1 1 2 1 B D B h J d s d s K T 75 Nesta expressão J é o determinante Jacobiano que corresponde ao determinante da matriz Jacobiana J A matriz Jacobiana correspondente à transformação 66 é definida da seguinte forma 62 2 2 1 2 2 1 1 1 s x s x s x s x J 76 2 2 1 2 2 1 1 1 s x s x s x s x J J 77 Para permitir o cálculo do integral 75 todos os componentes da função integranda têm de depender de s1 e s2 Se a matriz D 31 não for constante é possível utilizar o mesmo tipo de interpolação para definir E e ν em função de s1 e s2 4 2 1 4 3 2 1 3 2 2 1 2 1 2 1 1 2 1 E s s N E s s N E s s N E N s s E s s 78 4 2 1 4 3 2 1 3 2 2 1 2 1 2 1 1 2 1 ν ν ν ν ν s s N s s N s s N N s s s s 79 Nesta expressão iE e iν são os valores no nó i do módulo de Young e do coeficiente de Poisson Na generalidade dos casos práticos E e ν são considerados constantes ao nível de cada elemento finito Quando uma estrutura apresenta mais do que um tipo de material a fronteira entre as zonas correspondentes a cada material deve coincidir com a transição entre elementos finitos Se a espessura do elemento não for constante pode ser interpolada de um modo semelhante 4 2 1 4 3 2 1 3 2 2 1 2 1 2 1 1 2 1 h s s N h s s N h s s N h N s s h s s 80 Estado Plano de Tensão Álvaro F M Azevedo 104 Nesta expressão ih é o valor da espessura no nó i Os elementos da matriz Jacobiana 76 obtêmse por derivação de 68 e 69 resultando 41 1 4 31 1 3 21 1 2 11 1 1 1 1 x s N x s N x s N x s N s x 81 41 2 4 31 2 3 21 2 2 11 2 1 2 1 x s N x s N x s N x s N s x 82 42 1 4 32 1 3 22 1 2 12 1 1 1 2 x s N x s N x s N x s N s x 83 42 2 4 32 2 3 22 2 2 12 2 1 2 2 x s N x s N x s N x s N s x 84 As equações 8184 são equivalentes à seguinte equação matricial 2 4 1 4 2 3 1 3 2 2 1 2 2 1 1 1 42 32 22 12 41 31 21 11 2 2 1 2 2 1 1 1 s N s N s N s N s N s N s N s N x x x x x x x x s x s x s x s x 85 De um modo mais compacto temse n m n m T m m s N x J 86 sendo Estado Plano de Tensão Álvaro F M Azevedo 105 2 4 1 4 2 3 1 3 2 2 1 2 2 1 1 1 s N s N s N s N s N s N s N s N s N 87 As expressões dos elementos da matriz 87 obtêmse por derivação de 70 em ordem a s1 e s2 resultando 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 1 2 1 2 1 2 1 2 s s s s s s s s s N 88 Substituindo as expressões 88 em 85 obtêmse os elementos da matriz Jacobiana em função de s1 e s2 Nota os elementos da matriz x são as coordenadas cartesianas dos nós sendo portanto constantes de valor conhecido Tendo em vista o cálculo do integral 75 encontramse já definidos em função de s1 e s2 todos os componentes da função integranda com excepção da matriz B Apresentase em seguida o procedimento para a sua obtenção No caso do elemento finito quadrilátero de quatro nós e de geometria arbitrária as equações 912 permanecem válidas ver as Secções 62 e 63 As equações 10 e 11 são em seguida reproduzidas 42 41 32 31 22 21 12 11 4 3 2 1 4 3 2 1 2 1 0 0 0 0 0 0 0 0 a a a a a a a a N N N N N N N N u u 89 Estado Plano de Tensão Álvaro F M Azevedo 106 1 1 p p m m a N u 90 No caso do elemento de geometria arbitrária as funções de forma N dependem das variáveis s1 e s2 Neste caso a interpolação dos deslocamentos 89 é efectuada de um modo coincidente com o que foi utilizado para interpolar as coordenadas cartesianas dos nós 68 e 69 Quando o método de interpolação dos deslocamentos nodais e das coordenadas cartesianas dos nós coincidem dizse que a formulação do elemento finito é isoparamétrica O campo de deformações obtémse de um modo semelhante ao que foi descrito na Secção 63 reproduzindose em seguida as equações mais significativas 1 1 m q m q u L ε 91 Por substituição de 90 em 91 chegase a 1 1 p p m q m q a N L ε 92 p m q m q p N L B 93 1 1 p p q q a B ε 94 4 3 2 1 4 3 2 1 1 2 2 1 0 0 0 0 0 0 0 0 0 0 N N N N N N N N x x x x B 95 1 4 2 4 1 3 2 3 1 2 2 2 1 1 2 1 2 4 2 3 2 2 2 1 1 4 1 3 1 2 1 1 0 0 0 0 0 0 0 0 x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N B 96 Estado Plano de Tensão Álvaro F M Azevedo 107 A matriz B depende das derivadas das funções de forma em ordem a xj j i N x De modo a ser possível calcular o integral 75 é necessário obter as expressões de j i N x em função de s1 e s2 Considerese uma das funções de forma Ni dependendo de x1 e x2 que por sua vez dependem de s1 e s2 2 1 2 2 1 1 s s x x s s N N i i 97 Pela regra da cadeia temse 2 2 2 2 1 1 2 1 2 2 1 1 1 1 s x x N s x x N s N s x x N s x x N s N i i i i i i 98 que se pode escrever da seguinte forma em notação matricial 2 2 1 2 2 1 1 1 2 1 2 1 s x s x s x s x x N x N s N s N i i i i 99 Atribuindo ao índice i os valores 1 a 4 e agrupando os quatro casos nas seguintes matrizes chegase a 2 2 1 2 2 1 1 1 2 4 1 4 2 3 1 3 2 2 1 2 2 1 1 1 2 4 1 4 2 3 1 3 2 2 1 2 2 1 1 1 s x s x s x s x x N x N x N x N x N x N x N x N s N s N s N s N s N s N s N s N 100 que de um modo mais compacto se pode escrever x J N s N 101 Estado Plano de Tensão Álvaro F M Azevedo 108 sendo J a matriz Jacobiana definida em 76 e em 86 Multiplicando ambos os membros de 101 à direita por J 1 obtémse m m n m m n J s N x N 1 102 A matriz N s foi definida em 87 e 88 sendo os seus elementos funções de s1 e s2 Em 86 pode verificarse que os elementos da matriz J são também funções de s1 e s2 Os elementos da seguinte matriz que dependem de s1 e s2 2 4 1 4 2 3 1 3 2 2 1 2 2 1 1 1 x N x N x N x N x N x N x N x N x N 103 são depois espalhados na matriz B de acordo com 96 Deste modo se alcançou o objectivo de calcular os elementos da matriz B como sendo funções de s1 e s2 Uma vez que todos os componentes da função integranda de 75 se encontram definidos em função de s1 e s2 é agora possível proceder ao cálculo da matriz de rigidez do elemento finito O facto de se tratar de um integral de difícil resolução e de os limites de integração serem 1 e 1 sugere o recurso à técnica de integração numérica que se encontra descrita no Capítulo 5 67 Algoritmo de cálculo da matriz de rigidez de um elemento isoparamétrico Um integral duplo cujos limites de integração sejam 1 e 1 para ambas as variáveis pode ser calculado pela quadratura de Gauss sendo o resultado obtido em geral um valor aproximado De acordo com o que foi exposto no Capítulo 5 a correspondente expressão é a seguinte Estado Plano de Tensão Álvaro F M Azevedo 109 1 2 1 1 2 1 1 1 1 1 2 1 GP GP n i n j j i j i f P P W W d s d s f s s 104 Nesta expressão nGP1 é o número de pontos de Gauss associado à direcção s1 e nGP2 é o número correspondente à direcção s2 Os parâmetros Wi e Wj são os pesos associados às direcções s1 e s2 A função f deve ser avaliada nos pontos de Gauss cujas coordenadas são i Pj P s s 2 1 105 De aqui em diante o segundo membro de 104 passa a ser designado por J Assim no caso de ser nGP1 2 e nGP2 2 da expansão dos somatórios em 104 resulta a seguinte expressão para J 1 1 2 2 1 1 GP n i i i i i f P P W W W W f P P J 106 2 2 2 2 1 2 1 2 2 1 2 1 1 1 1 1 f P P W W W f P P W f P P W W W W f P P J 107 De acordo com o que foi exposto no Capítulo 5 os valores dos pesos Wi e das posições Pi é neste caso 0 57735 02692 3 1 0 57735 02692 3 1 1 1 2 1 2 1 P P W W 108 passando J a ser avaliado do seguinte modo 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 f f f f J 109 O valor aproximado do integral duplo 104 depende do resultado da avaliação da função f s1 s2 em quatro pontos de Gauss cuja localização se encontra na Figura 67 Estado Plano de Tensão Álvaro F M Azevedo 110 1 2 3 4 s1 s2 1 1 1 1 Ponto de Gauss 1 2 3 4 Fig 67 Localização dos quatro pontos de Gauss no sistema de coordenadas s1 s2 De acordo com 107 para calcular o valor aproximado do integral 75 recorrendo à quadratura de Gauss com 2x2 pontos procedese do seguinte modo avaliase a sua função integranda nos quatro pontos de Gauss multiplicase o resultado correspondente a cada ponto de Gauss pelos respectivos pesos que neste caso são unitários e somamse as quatro parcelas Como se pode verificar em 75 a função integranda é um produto de matrizes de funções que é em seguida multiplicado por funções escalares Atendendo às características da quadratura de Gauss é possível avaliar todos os elementos de cada matriz em cada ponto de Gauss e só em seguida fazer o produto matricial bem como o produto pelas funções escalares avaliadas também nesse ponto de Gauss Deste modo os produtos matriciais são efectuados com valores numéricos em vez de funções facilitando assim a programação deste algoritmo em computador Apresentase em seguida a sequência de operações que têm de ser efectuadas para calcular a matriz de rigidez de um elemento finito quadrilátero recorrendo à quadratura de Gauss com 2x2 pontos Dados coordenadas cartesianas dos nós ijx espessura do elemento finito em cada nó ih Estado Plano de Tensão Álvaro F M Azevedo 113 calcular a matriz de elasticidade D com a expressão 31 calcular o vector tensão σ com a expressão 30 Depois de obtidos os valores de σ e ε nos pontos de Gauss 2x2 é possível fazer a sua interpolação ou extrapolação para outros pontos do elemento nomeadamente para os seus nós 64 Desta forma se obtêm resultados mais precisos do que aqueles que se obteriam com a avaliação directa das tensões no ponto pretendido 69 Considerações finais Neste capítulo foi apresentado o modo de obter a matriz de rigidez de um elemento finito quadrilátero de geometria arbitrária destinado à discretização de estados planos de tensão Foi apresentado com detalhe o caso do elemento de quatro nós e da quadratura de Gauss com 2x2 pontos Alguns aspectos importantes são deixados para outros capítulos tais como a assemblagem da matriz de rigidez global o desenvolvimento de elementos com mais do que quatro nós a influência do número de pontos de Gauss na qualidade dos resultados o cálculo de acções nodais equivalentes a acções concentradas distribuidas e de volume etc BIBLIOGRAFIA 61 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 62 Kreyszig E Advanced Engineering Mathematics Sixth Edition John Wiley Sons Inc 1988 63 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 64 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 Estado Plano de Tensão Álvaro F M Azevedo 114 115 CAPÍTULO 7 FUNÇÕES INTERPOLADORAS Neste capítulo são descritos diversos modos de obtenção de funções interpoladoras também designadas funções de forma São apresentados exemplos relativos a meios unidimensionais bidimensionais e tridimensionais As funções de forma obtidas por procedimentos genéricos podem depois ser utilizadas em distintas formulações do método dos elementos finitos 71 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada no âmbito da determinação de funções interpoladoras Tabela 71 Simbologia relativa à determinação de funções interpoladoras x Coordenada cartesiana x Coordenada cartesiana de um nó de um elemento finito u Campo de deslocamentos a Deslocamento nodal N Função interpoladora ou função de forma n Número de nós do elemento finito L Dimensão do elemento finito s Coordenada local curvilínea h Espessura do elemento finito laminar h Espessura do elemento finito num nó s Coordenada local de um nó de um elemento finito NV Vector das funções interpoladoras ou funções de forma Funções Interpoladoras Álvaro F M Azevedo 116 V Vector contendo os factores não constantes de um polinómio c Coeficiente de um termo de um polinómio Q Matriz cujas colunas contêm o vector V avaliado em nós do elemento finito p Número de nós de um bordo de um elemento finito Deslocamento de um nó de um elemento finito θ Rotação de um nó de um elemento finito 72 Caso unidimensional Na Figura 71 encontrase representado um elemento finito unidimensional com quatro nós colocados sobre o eixo x A posição de cada nó é definida pela respectiva coordenada cartesiana ix sendo i o número do nó u x x 1a 1 2 1x x 3 2a 3a x 2x 3x x 4 4a 4x x Fig 71 Elemento finito unidimensional de geometria arbitrária As características essenciais de uma função de forma Ni são as seguintes deve assumir o valor unitário para x ix deve anularse nos restantes nós É também desejável no caso das funções polinomiais manter o grau do polinómio tão baixo quanto possível Na Tabela 72 encontramse os valores que cada função de forma deve assumir nos nós do elemento finito Funções Interpoladoras Álvaro F M Azevedo 117 Tabela 72 Características das funções N1x N2x N3x e N4x x 1x 2x 3x 4x N1x 1 0 0 0 N2x 0 1 0 0 N3x 0 0 1 0 N4x 0 0 0 1 É fácil verificar que as seguintes funções de forma são polinómios que respeitam as condições definidas na Tabela 72 4 1 3 1 2 1 4 3 2 1 x x x x x x x x x x x x x N 1 4 2 3 2 1 2 4 3 1 2 x x x x x x x x x x x x x N 2 4 3 2 3 1 3 4 2 1 3 x x x x x x x x x x x x x N 3 3 4 2 4 1 4 3 2 1 4 x x x x x x x x x x x x x N 4 A expressão genérica para o caso de um elemento finito unidimensional com n nós é n i k k k i k i x x x x x N 1 5 A expressão 5 é designada fórmula de interpolação de Lagrange 71 sendo as expressões 14 o caso particular de 5 quando n 4 Se em 5 se considerar n 2 1 1 x e x2 1 obtêmse as funções de forma que foram determinadas no Capítulo 4 para o caso da barra de dois nós e comprimento L 2 De um modo semelhante seria possível verificar a coincidência Funções Interpoladoras Álvaro F M Azevedo 118 entre as restantes funções de forma determinadas no Capítulo 4 e as que se obtêm com 5 73 Caso bidimensional Considerese agora o elemento finito bidimensional com 16 nós representado na Figura 72 s1 1 2 23 3 4 s2 5 6 7 8 9 10 11 12 13 14 15 16 23 23 23 23 23 Fig 72 Elemento finito bidimensional com 16 nós Relativamente ao elemento finito de 16 nós pretendese obter a função de forma N7 s1 s2 Esta função deve ser unitária no nó 7 e deve anularse nos restantes nós As coordenadas do nó 7 são s1 s2 13 13 Na direcção s1 o nó 7 é o terceiro nó Por isso devese utilizar a função N3 indicada em 3 e considerar x 1s 1 1 x 1 3 2 x x3 1 3 e x4 1 Esta função é designada N31 e tem a seguinte expressão Funções Interpoladoras Álvaro F M Azevedo 119 1 3 1 1 3 1 3 3 1 1 1 1 3 1 1 1 1 1 31 s s s s N 6 1 1 3 1 16 27 1 1 1 1 31 s s s s N 7 Os índices em N31 têm o significado de função de forma unidimensional correspondente ao nó 3 e com x substituido por s1 Na direcção s2 o nó 7 é o segundo nó Por isso devese utilizar a função N2 indicada em 2 considerar x 2s e de igual forma 1 1 x 1 3 2 x x3 1 3 e x4 1 Esta função é designada N22 e tem a seguinte expressão 1 3 1 1 3 1 3 3 1 1 1 1 3 1 2 2 2 2 22 s s s s N 8 1 1 3 1 16 27 2 2 2 2 22 s s s s N 9 A função N7 s1 s2 é o produto de 7 por 9 2 22 1 31 2 1 7 s N s N s s N 10 1 1 3 1 1 1 3 1 256 729 2 2 2 1 1 1 2 1 7 s s s s s s s s N 11 Como se pode facilmente verificar esta função de forma assume o valor unitário no nó 7 e anulase nos restantes nós As funções de forma correspondentes aos restantes 15 nós poderiam ser obtidas de um modo idêntico ao que foi aqui apresentado Na Figura 73 encontrase em perspectiva o gráfico da função N7 s1 s2 Funções Interpoladoras Álvaro F M Azevedo 120 s1 s2 N7 s1 s2 Fig 73 Gráfico da função de forma N7 s1 s2 A expressão 11 é equivalente à seguinte 3 2 3 1 3 2 2 1 2 2 3 1 3 2 1 2 2 2 1 2 3 1 3 2 2 2 1 2 2 1 3 1 2 2 2 1 2 1 2 1 2 1 7 256 729 256 243 256 243 256 729 256 81 256 729 256 243 256 243 256 243 256 243 256 81 256 729 256 81 256 243 256 243 256 81 s s s s s s s s s s s s s s s s s s s s s s s s s s N 12 Funções Interpoladoras Álvaro F M Azevedo 121 O triângulo de Pascal correspondente a uma função de duas variáveis é o seguinte 6 2 5 2 1 4 2 2 1 3 2 3 1 2 2 4 1 2 5 1 6 1 5 2 4 2 1 3 2 2 1 2 2 3 1 2 4 1 5 1 4 2 3 2 1 2 2 2 1 2 3 1 4 1 3 2 2 2 1 2 2 1 3 1 2 2 2 1 2 1 2 1 1 s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s 13 Comparando 12 com o triângulo de Pascal representado em 13 pode observarse que a função de forma N7 s1 s2 é um polinómio de sexto grau incompleto em que foram utilizados apenas os 16 termos que figuram em 12 74 Procedimento genérico para determinar as funções de forma Apresentase em seguida um procedimento que permite determinar as funções de forma de um elemento finito com n nós arbitrariamente distribuídos 72 A exposição que se segue baseiase num exemplo que consiste num elemento finito de cinco nós posicionados de acordo com a Figura 74 s1 1 2 1 3 4 s2 1h h s1 s2 1 1 1 5 2h 3h 4h 5h Fig 74 Elemento finito com cinco nós As coordenadas dos cinco nós do elemento finito são no sistema de eixos s1 s2 Funções Interpoladoras Álvaro F M Azevedo 122 1 1 1 0 1 1 1 1 1 1 52 51 42 41 32 31 22 21 12 11 s s s s s s s s s s s 14 Pretendese fazer a interpolação do campo de espessuras h s1 s2 sendo utilizada a seguinte expressão em que ih representa a espessura do elemento finito no nó i 5 2 1 5 4 2 1 4 3 2 1 3 2 2 1 2 1 2 1 1 2 1 h s s N h s s N h s s N h s s N h N s s s s h 15 Recorrendo à notação matricial a equação 15 passa a 5 4 3 2 1 5 4 3 2 1 N N N N N h h h h h h 16 ou h h T NV 17 sendo 5 4 3 2 1 h h h h h h 18 5 4 3 2 1 N N N N N NV 19 Funções Interpoladoras Álvaro F M Azevedo 123 Tendo em vista a determinação das cinco funções de forma polinomiais Ni é necessário seleccionar no triângulo de Pascal um número de termos igual ao número de nós do elemento finito Por este motivo o exemplo da Figura 74 requer a escolha de cinco termos que devem ser de grau tão baixo quanto possível No triângulo de Pascal atrás apresentado 13 são assim seleccionados os seguintes termos que se agrupam num vector designado por V 2 1 2 1 2 1 1 s s s s s V 20 Na selecção efectuada foi dada preferência a termos de grau mais elevado em s1 do que em s2 devido ao facto de o elemento finito apresentar mais nós segundo a direcção s1 De acordo com a selecção de termos efectuada a função h s1 s2 vai ser aproximada com o seguinte polinómio 2 1 5 2 1 4 2 3 1 2 1 1 2 c s s c s c s c s c h s s 21 que em notação matricial se escreve 2 1 2 1 2 1 5 4 3 2 1 1 s s s s s c c c c c h 22 ou c V h T 23 sendo Funções Interpoladoras Álvaro F M Azevedo 124 5 4 3 2 1 c c c c c c 24 Ao efectuar em 22 a substituição das variáveis s1 e s2 pelas coordenadas do nó 1 pretendese obter o valor da espessura h no nó 1 1 h 12 11 2 11 12 11 5 4 3 2 1 1 1 s s s s s c c c c c h 25 Procedendo de igual forma com os restantes nós e agrupando as cinco expressões do tipo 25 numa única expressão matricial temse 52 51 42 41 32 31 22 21 12 11 2 51 2 41 2 31 2 21 2 11 52 42 32 22 12 51 41 31 21 11 5 4 3 2 1 5 4 3 2 1 1 1 1 1 1 s s s s s s s s s s s s s s s s s s s s s s s s s c c c c c h h h h h 26 ou Q c h T T 27 sendo 52 51 42 41 32 31 22 21 12 11 2 51 2 41 2 31 2 21 2 11 52 42 32 22 12 51 41 31 21 11 1 1 1 1 1 s s s s s s s s s s s s s s s s s s s s s s s s s Q 28 Funções Interpoladoras Álvaro F M Azevedo 125 No caso do exemplo da Figura 74 e de acordo com 14 os elementos de Q são 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 Q 29 Uma vez que a matriz Q é quadrada e se supõe não singular podese multiplicar à direita ambos os membros de 27 por Q1 resultando 1 Q h c T T 30 Substituindo o segundo membro de 30 em 23 resulta V Q h h T 1 31 Uma vez que são iguais os segundos membros de 17 e 31 e uma vez que o vector de espessuras h é arbitrário concluise que V Q NV 1 32 No caso do exemplo da Figura 74 a inversa da matriz Q 29 é 1 4 1 2 1 4 1 4 1 4 0 1 0 0 1 1 4 1 2 1 4 1 4 1 4 1 4 0 1 4 1 4 1 4 1 4 0 1 4 1 4 4 1 Q 1 33 resultando Funções Interpoladoras Álvaro F M Azevedo 126 2 1 2 1 2 1 1 1 4 1 2 1 4 1 4 1 4 0 1 0 0 1 1 4 1 2 1 4 1 4 1 4 1 4 0 1 4 1 4 1 4 1 4 0 1 4 1 4 4 1 s s s s s NV 34 As funções de forma são 4 1 2 1 2 1 2 1 1 s s s s N s s 35 4 1 2 1 2 1 2 1 2 s s s s s s N 36 4 2 1 2 1 2 1 2 1 2 1 3 s s s s s s s N 37 2 1 2 1 4 1 s s s N 38 4 2 1 2 1 2 1 2 1 2 1 5 s s s s s s s N 39 Existem alguns casos em que devido à localização dos nós ou devido à incorrecta selecção de termos no triângulo de Pascal a matriz Q resulta singular Nestes casos o procedimento aqui descrito não pode ser utilizado 75 Elementos bidimensionais famílias Lagrangeana e serendipity O procedimento descrito na Secção 74 encontrase bem definido com excepção do facto de ser necessário seleccionar em cada caso um adequado conjunto de termos no triângulo de Pascal Nos casos em que não existe um critério óbvio é conveniente ensaiar várias alternativas De cada conjunto de termos do triângulo de Pascal vai resultar uma distinta formulação do elemento finito sendo conveniente averiguar qual é a que conduz a resultados mais precisos Para as situações mais comuns existem já formulações que conduzem a bons resultados sendo em seguida apresentados dois desses casos que são designados de elementos da família Lagrangeana e elementos da família serendipity 72 Funções Interpoladoras Álvaro F M Azevedo 127 Os elementos bidimensionais da família Lagrangeana são quadriláteros com p2 nós sendo p o número de nós de um bordo ver a Figura 75 s1 L1 2 s2 L2 2 s1 s1 s2 s2 L1 2 L1 2 p 2 p 3 p 4 Fig 75 Elementos finitos bidimensionais da família Lagrangeana As funções de forma do elemento Lagrangeano com p 4 foram já apresentadas na Secção 73 Quando se determinam as funções de forma com o procedimento genérico descrito na Secção 74 devese seleccionar os termos do triângulo de Pascal com o critério definido na Figura 76 1 p 2 1s 2 1s 3 1s 4 1s 5 1s 6 1s 2s 2 2s 3 2s 4 2s 5 2s 6 2s 1 s2 s 2 2 s1 s 2 1 s2 s 2 3 s1 s 2 2 2 s1 s 3 1 s2 s 2 4 s1 s 2 2 3 s1 s 3 2 2 s1 s 4 1 s2 s 2 5 s1 s 2 2 4 s1 s 3 2 3 s1 s 4 2 2 s1 s 5 1 s2 s p 3 p 4 Fig 76 Selecção de termos no triângulo de Pascal para elementos finitos bidimensionais da família Lagrangeana Como se pode observar na Figura 76 o critério de selecção de termos no triângulo de Pascal é facilmente extensível a valores superiores de p Funções Interpoladoras Álvaro F M Azevedo 128 Apresentamse na Figura 77 alguns exemplos de elementos finitos da família serendipity s1 L1 2 s2 L2 2 s1 s1 s2 s2 L1 2 L1 2 p 2 p 3 p 4 Fig 77 Elementos finitos bidimensionais da família serendipity O número de nós de cada elemento da família serendipity é 4 p 1 sendo p o número de nós de um bordo Na Figura 78 encontrase o critério de selecção de termos no triângulo de Pascal para o caso de elementos da família serendipity 1 p 2 1s 2 1s 3 1s 4 1s 5 1s 6 1s 2s 2 2s 3 2s 4 2s 5 2s 6 2s 1 s2 s 2 2 s1 s 2 1 s2 s 2 3 s1 s 2 2 2 s1 s 3 1 s2 s 2 4 s1 s 2 2 3 s1 s 3 2 2 s1 s 4 1 s2 s 2 5 s1 s 2 2 4 s1 s 3 2 3 s1 s 4 2 2 s1 s 5 1 s2 s p 3 p 4 p 2 p 3 p 4 p 5 p 6 p 5 p 6 Fig 78 Selecção de termos no triângulo de Pascal para elementos finitos bidimensionais da família serendipity Funções Interpoladoras Álvaro F M Azevedo 129 De cada vez que p é incrementado uma unidade são acrescentados quatro nós ao elemento finito um em cada bordo e são seleccionados mais quatro termos no triângulo de Pascal Este critério é extensível a qualquer valor de p Na prática os elementos finitos que apresentam um bom compromisso entre o número de nós e a qualidade dos resultados obtidos são os da família serendipity com oito nós p 3 Apresentase na Figura 79 um exemplo de um destes elementos finitos no referencial x1 x2 x1 x2 p 3 Fig 79 Elemento finito de 8 nós da família serendipity Quando comparado com o quadrilátero de quatro nós o elemento finito representado na Figura 79 tem a vantagem de ser mais preciso e de se adaptar bem a fronteiras curvilíneas Apresentase em seguida um exemplo de um elemento finito que apresenta mais nós na direcção s2 do que na direcção s1 ver a Figura 710 Funções Interpoladoras Álvaro F M Azevedo 130 s1 L1 2 s2 L2 2 Fig 710 Elemento finito bidimensional com oito nós Tendo em vista a determinação das funções de forma do elemento finito representado na Figura 710 devem ser seleccionados os termos do triângulo de Pascal que se encontram assinalados na Figura 711 1 1s 2 1s 3 1s 4 1s 2s 2 2s 3 2s 4 2s 1 s2 s 2 2 s1 s 2 1 s2 s 2 3 s1 s 2 2 2 s1 s 3 1 s2 s Fig 711 Selecção de termos no triângulo de Pascal para o elemento finito bidimensional representado na Figura 710 São preferidos termos de grau mais elevado em s2 porque o elemento possui mais nós na direcção s2 do que na direcção s1 76 Propriedades das funções interpoladoras Considerese o elemento finito de três nós representado na Figura 712 Funções Interpoladoras Álvaro F M Azevedo 131 u x x 1a 1 2 1x x 3 2a 3a x 2x 3x x Fig 712 Elemento finito unidimensional de geometria arbitrária Supondo que não é efectuada qualquer substituição de variável a interpolação do campo de deslocamentos é efectuada da seguinte forma 3 3 2 2 1 1 x a N x a N N x a u x 40 Admitase agora que em todos os nós é imposto o mesmo deslocamento 3 2 1 a a a 41 Neste caso pretendese que a função interpolada u x seja uma função constante u x 42 em todos os pontos do elemento finito Substituindo 41 e 42 em 40 resulta x N x N x N 3 2 1 43 1 3 2 1 x N x N N x 44 1 1 n i Ni x 45 sendo n o número de nós do elemento finito A equação 45 constitui uma propriedade que as funções de forma devem possuir Só assim se garante que uma translação do elemento finito é correctamente interpolada com a equação 40 Funções Interpoladoras Álvaro F M Azevedo 132 É fácil constatar que todos os conjuntos de funções de forma apresentados nos Capítulos 4 e 6 possuem a propriedade 45 Outra questão que se coloca é a de definir um procedimento que garanta que as funções interpoladoras que se pretende determinar possuam a propriedade 45 Com este objectivo considerese a expressão que define as funções interpoladoras 32 V Q NV 1 46 Multiplicando ambos os membros de 46 por Q obtémse V N Q V 47 que no exemplo da Figura 74 corresponde a ver a Secção 74 2 1 2 1 2 1 5 4 3 2 1 52 51 42 41 32 31 22 21 12 11 2 51 2 41 2 31 2 21 2 11 52 42 32 22 12 51 41 31 21 11 1 1 1 1 1 1 s s s s s N N N N N s s s s s s s s s s s s s s s s s s s s s s s s s 48 Como se pode verificar em 48 se no triângulo de Pascal o elemento unitário do seu vértice for o primeiro dos termos seleccionados então o primeiro elemento do vector V é sempre unitário e a primeira linha da matriz Q tem todos os elementos também unitários A primeira das cinco equações a que 48 corresponde é 1 5 4 3 2 1 N N N N N 49 As funções de forma determinadas com 46 respeitam as condições 48 e 49 Assim fica provado que sempre que o termo unitário do triângulo de Pascal é seleccionado então as funções de forma obtidas possuem a propriedade 45 77 Interpolação Hermitiana Em todas as interpolações que foram efectuadas nas secções anteriores apenas se atendeu aos valores nodais das funções Na interpolação Hermitiana que é descrita Funções Interpoladoras Álvaro F M Azevedo 133 nesta secção são também consideradas as derivadas das funções nos nós Este tipo de interpolação tem interesse para a formulação de elementos finitos em que são consideradas as rotações eg vigas lajes Na Figura 713 encontrase um elemento finito com dois nós e comprimento L A função ux corresponde ao deslocamento vertical cujos valores nodais são 1 e 2 Nos nós 1 e 2 a rotação é θ1 e θ2 respectivamente u x x a1 1 1 2 2 1 L x x 2 2 L x x L 2 L 2 a3 2 a2 θ1 a4 θ2 Fig 713 Interpolação Hermitiana num elemento unidimensional com dois nós Os deslocamentos generalizados dos nós do elemento finito representado na Figura 713 são os seguintes 2 2 1 1 4 3 2 1 θ θ a a a a a 50 De acordo com a Figura 713 e designando dudx por u x temse 2 2 2 2 2 2 4 2 2 3 1 1 2 1 1 1 u L x u a u L u x a L u x u a L u u x a θ θ 51 Uma vez que as rotações são muito pequenas supõese Funções Interpoladoras Álvaro F M Azevedo 134 tan θ θ 52 Pretendese determinar a função ux que respeita as condições 51 Com esse objectivo admitese que a função ux é o seguinte polinómio de grau 3 3 4 2 3 2 1 c x c x c x c u x 53 que em notação matricial corresponde a 3 2 4 3 2 1 1 x x x c c c c u x 54 ou c V x u T 55 sendo 4 3 2 1 c c c c c 56 e 3 2 1 x x x V 57 Derivando ambos os membros de 53 obtémse 2 4 3 2 3 2 c x c x c u x 58 que em notação matricial corresponde a Funções Interpoladoras Álvaro F M Azevedo 135 2 4 3 2 1 3 2 1 0 x x c c c c u x 59 ou c V x u T 60 sendo 3 2 2 1 0 x x V 61 Para que 51 se verifique quando as funções u e u são 54 e 59 é necessário que 3 1 2 1 1 4 3 2 1 1 1 1 x x x c c c c u x a 62 2 1 1 4 3 2 1 1 2 3 2 1 0 x x c c c c x u a 63 3 2 2 2 2 4 3 2 1 2 3 1 x x x c c c c u x a 64 2 2 2 4 3 2 1 2 4 3 2 1 0 x x c c c c x u a 65 Funções Interpoladoras Álvaro F M Azevedo 136 Agrupando 6265 numa única expressão matricial resulta 2 2 3 2 2 1 3 1 2 2 2 1 2 1 2 1 4 3 2 1 4 3 2 1 3 3 2 2 1 1 0 1 0 1 x x x x x x x x x x c c c c a a a a 66 ou Q c a T T 67 sendo a definido por 50 c definido por 56 e Q definido por 2 2 3 2 2 1 3 1 2 2 2 1 2 1 2 1 3 3 2 2 1 1 0 1 0 1 x x x x x x x x x x Q 68 Verificase assim que as colunas da matriz Q são constituídas pelos vectores V 57 e V 61 avaliados nos pontos nodais 1x e 2x No caso da Figura 713 temse 4 3 8 4 3 8 4 4 1 2 1 2 0 1 0 1 2 3 2 3 2 2 L L L L L L L L L L Q 69 Multiplicando à direita ambos os membros de 67 por Q 1 resulta 1 Q a c T T 70 Funções Interpoladoras Álvaro F M Azevedo 137 A matriz inversa de 69 é 2 3 2 3 1 1 1 2 1 4 8 2 0 3 2 1 2 1 1 2 1 4 8 2 0 3 2 2 1 L L L L L L L L L L Q 71 Substituindo 70 em 55 obtémse V Q a x u T 1 72 A interpolação que se pretende definir deve ter as seguintes características 4 4 3 3 2 2 1 1 x a N N x a x a N N x a u x 73 que em notação matricial corresponde a 4 3 2 1 4 3 2 1 N N N N a a a a u x 74 ou aT NV u x 75 sendo a definido por 50 e 4 3 2 1 N N N N N V 76 Uma vez que são iguais os segundos membros de 72 e 75 e uma vez que o vector dos deslocamentos nodais a é arbitrário concluise que V Q NV 1 77 Funções Interpoladoras Álvaro F M Azevedo 138 No caso do exemplo da Figura 713 as funções de forma obtêmse fazendo o produto de 71 por 57 resultando 3 3 1 2 2 3 2 1 L x L x x N 78 3 2 2 2 1 2 1 4 1 8 L x L x x L x N 79 3 3 3 2 2 3 2 1 L x L x x N 80 3 2 2 4 1 2 1 4 1 8 L x L x x L x N 81 No caso particular de ser L 2 as funções de forma são as seguintes 3 1 4 1 4 3 2 1 x x x N 82 3 2 2 4 1 4 1 4 1 4 1 x x x x N 83 3 3 4 1 4 3 2 1 x x x N 84 3 2 4 4 1 4 1 4 1 4 1 x x x x N 85 Os gráficos das funções 8285 encontramse representados na Figura 714 Funções Interpoladoras Álvaro F M Azevedo 139 1 1 1 N1 x x 1 1 1 N2 x x 1 1 1 N4 x x 1 1 1 N3 x x Fig 714 Gráficos das funções Nix correspondentes ao elemento de dois nós com comprimento L 2 Apresentase em seguida o caso da interpolação Hermitiana de um elemento de três nós O elemento considerado tem comprimento L 2 e o nó intermédio centrado ver a Figura 715 u x x a1 1 1 3 1 1 x x x3 1 x L 2 a5 3 a2 θ1 a6 θ3 2 a4 θ2 a3 2 x2 0 x Fig 715 Interpolação Hermitiana num elemento unidimensional com três nós Funções Interpoladoras Álvaro F M Azevedo 140 O vector dos deslocamentos generalizados é 3 3 2 2 1 1 6 5 4 3 2 1 θ θ θ a a a a a a a 86 A função ux que respeita as condições indicadas na Figura 715 é 5 6 4 5 3 4 2 3 2 1 c x c x c x c x c x c u x 87 5 4 3 2 6 5 4 3 2 1 1 x x x x x c c c c c c u x 88 c V x u T 89 A derivada da função ux é 4 6 3 5 2 4 3 2 5 4 3 2 c x c x c x c x c u x 90 4 3 2 6 5 4 3 2 1 5 4 3 2 1 0 x x x x c c c c c c u x 91 c V x u T 92 Funções Interpoladoras Álvaro F M Azevedo 141 A matriz Q é neste caso a seguinte 4 3 5 3 4 2 5 2 4 1 5 1 3 3 4 3 3 2 4 2 3 1 4 1 2 3 3 3 2 2 3 2 2 1 3 1 3 2 3 2 2 2 1 2 1 3 2 1 5 5 5 4 4 4 3 3 3 2 2 2 1 1 1 0 1 0 1 0 1 x x x x x x x x x x x x x x x x x x x x x x x x x x x Q 93 De acordo com as coordenadas indicadas na Figura 715 temse 5 1 0 0 5 1 4 1 0 0 4 1 3 1 0 0 3 1 2 1 0 0 2 1 1 1 1 0 1 1 0 1 0 1 0 1 Q 94 1 4 1 4 1 4 1 4 0 0 3 4 1 2 5 4 1 0 0 1 0 2 0 1 0 0 1 0 2 0 1 1 4 1 4 1 4 1 4 0 0 3 4 1 2 5 4 1 0 0 Q 1 95 Atendendo a 77 temse V Q NV 1 96 5 4 3 2 6 5 4 3 2 1 1 1 4 1 4 1 4 1 4 0 0 3 4 1 2 5 4 1 0 0 1 0 2 0 1 0 0 1 0 2 0 1 1 4 1 4 1 4 1 4 0 0 3 4 1 2 5 4 1 0 0 x x x x x N N N N N N 97 Funções Interpoladoras Álvaro F M Azevedo 142 5 4 3 2 1 4 3 2 1 4 5 x x x x x N 98 5 4 3 2 2 4 1 4 1 4 1 4 1 x x x x x N 99 4 2 3 2 1 x x x N 100 5 3 4 2 x x x x N 101 5 4 3 2 5 4 3 2 1 4 5 x x x x x N 102 5 4 3 2 6 4 1 4 1 4 1 4 1 x x x x x N 103 Os gráficos das funções 98103 encontramse representados na Figura 716 78 Considerações finais Neste capítulo foram apresentados alguns procedimentos destinados à determinação de funções de forma Sempre que os procedimentos mais simples não sejam aplicáveis devese utilizar um dos métodos genéricos descritos nas Secções 74 e 77 A metodologia descrita na Secção 74 pode ser facilmente adaptada aos casos tridimensionais Neste caso no lugar do triângulo de Pascal temse uma pirâmide em cujo vértice figura o elemento unitário seguido de um segundo nível em que figuram as variáveis s1 s2 e s3 etc Funções Interpoladoras Álvaro F M Azevedo 143 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 N1 x N2 x N4 x N3 x N5 x N6 x x x x x x x Fig 716 Gráficos das funções Nix correspondentes ao elemento de três nós com comprimento L 2 BIBLIOGRAFIA 71 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 72 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 Funções Interpoladoras Álvaro F M Azevedo 144 145 CAPÍTULO 8 ASSEMBLAGEM DE ELEMENTOS FINITOS No Capítulo 3 foi apresentado com detalhe o caso da assemblagem de barras em problemas unidimensionais Neste capítulo apresentase de um modo sucinto a adaptação da técnica já descrita ao caso dos elementos finitos com mais do que dois nós e mais do que um grau de liberdade por nó 81 81 Simbologia Apresentase em primeiro lugar a simbologia adoptada na descrição da assemblagem de elementos finitos Tabela 81 Simbologia relativa à assemblagem de elementos finitos x Coordenada cartesiana a Deslocamentos nodais nos graus de liberdade da estrutura no referencial geral ag Deslocamentos nodais nos graus de liberdade do elemento finito no referencial geral K Matriz de rigidez da estrutura no referencial geral Kg Matriz de rigidez do elemento finito no referencial geral F Forças nodais equivalentes à acção exterior nos graus de liberdade da estrutura no referencial geral Fg Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial geral Assemblagem de Elementos Finitos Álvaro F M Azevedo 146 82 Assemblagem da matriz de rigidez global e do vector solicitação Depois de calculadas as matrizes de rigidez de todos os elementos finitos no referencial geral Kg é necessário proceder ao cálculo da matriz de rigidez global da estrutura K Uma operação semelhante tem de ser efectuada com os vectores solicitação dos diversos elementos finitos A assemblagem das matrizes de rigidez dos diversos elementos finitos na matriz de rigidez global é em seguida apresentada com base no exemplo da Figura 81 x1 4 5 2 1 x2 a1 a2 a3 a4 a9 a10 a7 a8 3 a5 a6 6 a11 a12 A B C Fig 81 Estrutura constituída por um elemento de 4 nós A um elemento de 2 nós B e um elemento de 3 nós C A estrutura representada na Figura 81 tem seis nós 1 a 6 e três elementos finitos A B e C O elemento A tem quatro nós o elemento B tem dois nós e o elemento C tem três nós Em cada nó existem dois graus de liberdade Em correspondência com os doze graus de liberdade da estrutura existem doze deslocamentos nodais a e doze forças nodais equivalentes à acção exterior F Assemblagem de Elementos Finitos Álvaro F M Azevedo 147 62 61 52 51 42 41 32 31 22 21 12 11 12 11 10 9 8 7 6 5 4 3 2 1 62 61 52 51 42 41 32 31 22 21 12 11 12 11 10 9 8 7 6 5 4 3 2 1 F F F F F F F F F F F F F F F F F F F F F F F F F a a a a a a a a a a a a a a a a a a a a a a a a a 1 De acordo com 1 nas considerações que se seguem é adoptada a numeração dos graus de liberdade de 1 a 12 Na relação de rigidez correspondente à estrutura F K a 2 a matriz de rigidez global K é uma matriz 12x12 Nas Figuras 82 83 e 84 encontramse representados os elementos finitos que vão ser assemblados e a respectiva numeração local nós e graus de liberdade 1 2 3 4 a7 a8 a5 a6 a3 a4 a1 a2 A Fig 82 Numerações locais do elemento finito de 4 nós A Assemblagem de Elementos Finitos Álvaro F M Azevedo 148 1 a1 a2 2 a3 a4 B Fig 83 Numerações locais do elemento finito de 2 nós B 1 a1 a2 3 a5 a6 2 a3 a4 C Fig 84 Numerações locais do elemento finito de 3 nós C São as seguintes as matrizes de rigidez dos três elementos finitos no referencial geral 88 87 86 85 84 83 82 81 78 77 76 75 74 73 72 71 68 67 66 65 64 63 62 61 58 57 56 55 54 53 52 51 48 47 46 45 44 43 42 41 38 37 36 35 34 33 32 31 28 27 26 25 24 23 22 21 18 17 16 15 14 13 12 11 A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A K A Elemento A g 3 44 43 42 41 34 33 32 31 24 23 22 21 14 13 12 11 B B B B B B B B B B B B B B B B K B Elemento B g 4 Assemblagem de Elementos Finitos Álvaro F M Azevedo 149 66 65 64 63 62 61 56 55 54 53 52 51 46 45 44 43 42 41 36 35 34 33 32 31 26 25 24 23 22 21 16 15 14 13 12 11 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C K C Elemento C g 5 Atendendo à numeração global dos graus de liberdade indicada na Figura 81 1 a 12 as matrizes de rigidez dos elementos finitos passam a ser 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 44 43 42 41 46 45 48 47 34 33 32 31 36 35 38 37 24 23 22 21 26 25 28 27 14 13 12 11 16 15 18 17 64 63 62 61 66 65 68 67 54 53 52 51 56 55 58 57 84 83 82 81 86 85 88 87 74 73 72 71 76 75 78 77 A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A K A 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 44 43 42 41 34 33 32 31 24 23 22 21 14 13 12 11 B B B B B B B B B B B B B B B B K B 7 Assemblagem de Elementos Finitos Álvaro F M Azevedo 150 44 43 42 41 46 45 34 33 32 31 36 35 24 23 22 21 26 25 14 13 12 11 16 15 64 63 62 61 66 65 54 53 52 51 56 55 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C K C 8 De acordo com o que foi exposto no Capítulo 3 a matriz de rigidez global é a soma de 6 7 e 8 resultando 44 43 42 41 46 45 34 33 32 31 36 35 24 23 22 44 21 43 42 41 26 25 46 45 48 47 14 13 12 34 11 33 32 31 16 15 36 35 38 37 24 23 22 21 26 25 28 27 14 13 12 11 16 15 18 17 64 63 62 61 66 44 65 43 42 41 54 53 52 51 56 34 55 33 32 31 64 63 62 61 24 23 22 66 21 65 68 67 54 53 52 51 14 13 12 56 11 55 58 57 84 83 82 81 86 85 88 87 74 73 72 71 76 75 78 77 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C A C A A A C C A A A A C C C A C A A A C C A A A A A A A A A A A A A A A A A A A A C C C C C B C B B B C C C C C B C B B B A A A A B B B A B A A A A A A A B B B A B A A A A A A A A A A A A A A A A A A A K K K K C B A 9 Em correspondência com os graus de liberdade indicados nas Figuras 81 a 84 têmse as forças nodais equivalentes às acções exteriores sobre a estrutura Assim e de acordo com o que foi exposto no Capítulo 3 são os seguintes os vectores solicitação correspondentes a cada elemento finito atendendo à numeração global da estrutura Assemblagem de Elementos Finitos Álvaro F M Azevedo 151 C C C C C C C B B B B B A A A A A A A A A F F F F F F F F F F F F F F F F F F F F F 4 3 2 1 6 5 4 3 2 1 4 3 2 1 6 5 8 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10 O vector F é a soma destes três vectores C C C A C A A A C B C B B A B A A A C B A F F F F F F F F F F F F F F F F F F F F F F 4 3 2 4 1 3 2 1 6 4 5 3 2 6 1 5 8 7 11 A relação de rigidez correspondente à totalidade dos graus de liberdade no referencial geral é a seguinte ver o Capítulo 3 C B A C B A F F F a K K K 12 F K a 13 Assemblagem de Elementos Finitos Álvaro F M Azevedo 152 Depois de acrescentar a 13 as condições de apoio ver o Capítulo 3 é possível resolver o sistema de equações lineares que daí resulta e obter os deslocamentos segundo todos os graus de liberdade da estrutura 83 Considerações finais Neste capítulo foi apresentada a assemblagem da matriz de rigidez global com base no armazenamento de todos os seus termos A matriz de rigidez global apresenta uma distribuição de termos particular que quando devidamente explorada conduz a significativas economias de recursos informáticos nomeadamente a redução do número de operações de cálculo e a diminuição da quantidade de memória consumida A característica mais simples de explorar é o facto de a matriz de rigidez global ser simétrica evitandose assim o cálculo e o armazenamento dos termos do seu triângulo inferior bem como todas as operações de cálculo que sobre eles teriam de ser efectuadas Considerando apenas os termos do triângulo superior é ainda vantajoso atender ao facto de muitos desses termos serem nulos O critério de selecção da técnica de armazenamento dos termos da matriz depende do método que vai ser usado para resolver o sistema de equações As técnicas de armazenamento mais comuns são as seguintes armazenamento em semibanda de largura constante armazenamento em semibanda de largura variável armazenamento em skyline e armazenamento esparso 82 BIBLIOGRAFIA 81 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 82 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 153 CAPÍTULO 9 FORÇAS NODAIS EQUIVALENTES Quando um elemento finito se encontra sujeito a acções exteriores genéricas é necessário proceder ao cálculo das forças nodais equivalentes à solicitação exterior Exemplos destas solicitações são as cargas concentradas num ponto do interior do elemento as cargas distribuídas em bordos as cargas distribuídas em faces e as forças de volume Começase por apresentar a formulação genérica do cálculo das forças nodais equivalentes seguindose um conjunto de exemplos ilustrativos dos procedimentos que em cada caso se devem adoptar 91 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada no estudo das forças nodais equivalentes a acções exteriores Tabela 91 Simbologia relativa às forças nodais equivalentes a acções exteriores x Coordenada cartesiana P Ponto onde actua uma carga concentrada L Arco onde actua uma carga distribuída por unidade de comprimento S Superfície onde actua uma carga distribuída por unidade de superfície V Volume onde actua uma carga distribuída por unidade de volume Q Carga concentrada p Carga distribuída por unidade de comprimento q Carga distribuída por unidade de superfície b Carga distribuída por unidade de volume ε Extensão Forças Nodais Equivalentes Álvaro F M Azevedo 154 γ Distorção σ Tensão normal τ Tensão tangencial u Campo de deslocamentos a Deslocamento nodal B Matriz de deformação D Matriz de elasticidade ε σ D C Elemento da matriz de elasticidade D E Módulo de elasticidade ou módulo de Young ν Coeficiente de Poisson N Função interpoladora ou função de forma K Matriz de rigidez do elemento finito no referencial geral F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial geral s Coordenada local curvilínea x Coordenada cartesiana de um nó de um elemento finito s Coordenada local de um nó de um elemento finito p Valor nodal da carga distribuída por unidade de comprimento NV Vector das funções interpoladoras ou funções de forma T Matriz de transformação nˆ Versor J Jacobiano da transformação ρ Massa específica do material g Aceleração da gravidade Forças Nodais Equivalentes Álvaro F M Azevedo 155 h Espessura do elemento finito laminar h Espessura do elemento finito num nó 92 Expressões genéricas das forças nodais equivalentes Na Figura 91 encontrase representado um corpo tridimensional sujeito a diversos tipos de acções exteriores x1 x2 x3 P x Q p x d L d L p L x q d S d S q S x b V dV b d V Fig 91 Corpo sujeito a diversos tipos de acções exteriores Os tipos de acções indicados na Figura 91 são os seguintes Força generalizada Q x concentrada no ponto P As componentes de Q x são três forças e três momentos Acção distribuída por unidade de comprimento p x Esta carga actua ao longo da linha L que se encontra definida no espaço e três dimensões As componentes de p x são três forças por unidade de comprimento e três momentos por unidade de comprimento Acção distribuída por unidade de superfície q x Esta carga actua na superfície S que se encontra definida no espaço e três dimensões As componentes de q x são três forças por unidade de superfície e três momentos por unidade de superfície Forças Nodais Equivalentes Álvaro F M Azevedo 156 Força de volume b x Esta carga actua num volume V que pode ser apenas uma parte do volume total do corpo As componentes de b x são três forças por unidade de volume Em problemas estáticos não são consideradas as componentes de momento por unidade de volume Todos os tipos de acções atrás referidos são definidos como funções das coordenadas cartesianas 3 2 1 x x x x 1 Na Figura 91 apenas foi indicado um exemplo de cada tipo de carga Nas aplicações do MEF é habitual existirem diversos exemplares de cada tipo de carga eg várias cargas concentradas em diferentes pontos do corpo várias cargas distribuídas em distintas zonas etc De acordo com o que foi exposto no Capítulo 4 o princípio dos trabalhos virtuais PTV estabelece que Trabalho Interno Trabalho Externo 2 Considerando todos os tipos de acções indicados na Figura 91 temse b V T q S T p L T Q T V T u b dV u q d S p d L u Q u V d δ δ δ δ δε σ 3 Na exposição que se segue não são consideradas as rotações nem os momentos Assim as componentes das diversas grandezas vectoriais que figuram em 3 são 12 31 23 3 2 1 12 31 23 3 2 1 τ τ τ σ σ σ σ γ γ γ ε ε ε ε 4 Forças Nodais Equivalentes Álvaro F M Azevedo 157 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 b b b b q q q q p p p p Q Q Q Q u u u u 5 Na formulação do MEF ver o Capítulo 6 o campo de deformações é interpolado a partir dos deslocamentos nodais com a seguinte expressão ε B a 6 Quando esta equação se refere aos deslocamentos virtuais e correspondentes deformações também virtuais temse B δ a δ ε 7 que é equivalente a T T T B δ a δ ε 8 No caso geral tridimensional e em materiais isotrópicos a relação entre tensões e deformações é a seguinte 91 12 31 23 3 2 1 3 3 3 1 2 2 2 1 2 2 2 1 12 31 23 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 γ γ γ ε ε ε τ τ τ σ σ σ C C C C C C C C C C C C 9 sendo ν ν ν 1 2 1 1 1 E C ν ν ν 1 2 1 2 E C ν 2 1 3 E C 10 Forças Nodais Equivalentes Álvaro F M Azevedo 158 ou de um modo mais compacto ε σ D 11 A matriz de elasticidade D depende do módulo de Young E e do coeficiente de Poisson ν Substituindo 6 em 11 resulta σ D B a 12 Na formulação do MEF ver o Capítulo 6 considerase que a interpolação do campo de deslocamentos a partir dos deslocamentos nodais é efectuada com a seguinte expressão u N a 13 A equação 13 referida à deformação virtual é a seguinte a N u δ δ 14 que é equivalente a T T T N a u δ δ 15 Substituindo todas estas equações em 3 passa a terse o PTV expresso por b V T T q S T T p L T T Q T T V T T a N b dV N q d S a a N p d L N Q a B D B a dV a δ δ δ δ δ 16 Uma vez que dV dx1 dx2 dx3 e os deslocamentos nodais não dependem das variáveis x1 x2 e x3 os vectores δ aT e a podem passar para fora dos integrais Forças Nodais Equivalentes Álvaro F M Azevedo 159 b V T T q S T T p L T T Q T T V T T N b dV a N q d S a N p d L a N Q a a B D B dV a δ δ δ δ δ 17 De acordo com o PTV a equação 17 é verdadeira para qualquer conjunto de deslocamentos virtuais concluindose assim que b V T q S T p L T Q T V T N b dV N q d S N p d L Q N a B D B dV 18 Comparando esta equação com a relação de rigidez que é utilizada no método dos deslocamentos F K a 19 temse para o caso geral indicado na Figura 91 V T B D B dV K 20 b b q q p p Q Q F F F F F 21 sendo as forças nodais equivalentes a cada carga as seguintes N Q F T Q 22 L T p p d L N F 23 S T q N q d S F 24 V T b N b dV F 25 Forças Nodais Equivalentes Álvaro F M Azevedo 160 Exceptuando casos particulares não se consegue uma precisão aceitável quando se discretiza um corpo com um único elemento finito Por este motivo devese considerar que as expressões 1825 se referem a um elemento finito e que depois se procede à habitual assemblagem da matriz de rigidez global e do vector solicitação global ver o Capítulo 8 93 Força concentrada num ponto interior O cálculo das forças nodais equivalentes a uma acção concentrada num ponto interior ao elemento finito é exemplificado com um elemento de quatro nós para estados planos de tensão ver a Figura 92 u1 x1 x2 x1 1 2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 Q1 Q2 P 2 1 Ponto sP sP P s1 s2 2 1 Ponto xP xP P Fig 92 Elemento finito de quatro nós com uma força concentrado num ponto interior De acordo com os graus de liberdade indicados na Figura 92 são os seguintes os vectores dos deslocamentos nodais e das correspondentes forças nodais Forças Nodais Equivalentes Álvaro F M Azevedo 161 42 41 32 31 22 21 12 11 42 41 32 31 22 21 12 11 F F F F F F F F F a a a a a a a a a 26 No ponto P encontrase aplicada uma força exterior com as seguintes componentes 2 1 Q Q Q 27 As coordenadas locais do ponto P são 2 1 2 1 P P P s s s s 28 As funções de forma do elemento finito são as seguintes ver o Capítulo 6 4 1 1 4 1 1 4 1 1 4 1 1 2 1 2 1 4 2 1 2 1 3 2 1 2 1 2 2 1 2 1 1 s s s s N s s s s N s s s s N s s s s N 29 As forças nodais equivalentes à carga concentrada Q são calculadas com a expressão 22 sendo a matriz N constituída pelas funções de forma 29 avaliadas no ponto 28 4 1 1 4 1 1 4 1 1 4 1 1 2 1 4 2 1 3 2 1 2 2 1 1 P P P P P P P P s s N s s N s s N s s N 30 De todas estas considerações resulta a seguinte expressão para o cálculo das forças nodais equivalentes à força Q Forças Nodais Equivalentes Álvaro F M Azevedo 162 2 1 4 4 3 3 2 2 1 1 42 41 32 31 22 21 12 11 0 0 0 0 0 0 0 0 Q Q N N N N N N N N F F F F F F F F F Q Q Q Q Q Q Q Q Q 31 A expressão 31 é facilmente avaliada desde que se conheçam as coordenadas locais s1 s2 do ponto P Contudo na generalidade dos casos práticos o ponto P é definido pelas suas coordenadas cartesianas x1 x2 Esta questão requer uma operação preliminar que consiste em calcular as coordenadas locais do ponto P a partir das suas coordenadas cartesianas Este cálculo é efectuado com base na interpolação das coordenadas cartesianas que foi apresentada no Capítulo 6 e que em seguida se reproduz 42 2 1 4 32 2 1 3 22 2 1 2 12 2 1 1 2 1 2 41 2 1 4 31 2 1 3 21 2 1 2 11 2 1 1 2 1 1 x s s N x s s N x s s N x N s s s s x x s s N x s s N x s s N x N s s s s x 32 Em 32 ijx representa a coordenada cartesiana do nó i segundo a direcção xj Substituindo em 32 x1 e x2 pelas coordenadas cartesianas do ponto P e Ni pelas funções de forma 29 resulta um sistema de duas equações não lineares com duas incógnitas s1 e s2 0 1 4 1 1 1 4 1 1 0 1 4 1 1 1 4 1 1 2 42 2 1 12 2 1 1 41 2 1 11 2 1 P P x x s s x s s x x s s x s s L L 33 que de um modo mais compacto se pode escrever da seguinte forma 0 0 2 1 2 2 1 1 s s f s s f 34 Forças Nodais Equivalentes Álvaro F M Azevedo 163 Este sistema de equações não lineares deve em geral ser resolvido por um método iterativo eg método de Newton A sua solução corresponde às coordenadas locais do ponto P 2 1 sP sP 94 Carga distribuída por unidade de comprimento Na Figura 93 encontrase representado o elemento finito de oito nós da família serendipity que neste caso se destina à discretização de estados planos de tensão Num dos bordos existe uma carga distribuída por unidade de comprimento p x x p x1 1 2 3 4 x2 s1 s2 5 6 7 8 Fig 93 Elemento finito de oito nós com uma carga distribuída por unidade de comprimento As funções de forma do elemento de oito nós são as seguintes 2 1 1 4 1 1 1 2 1 1 4 1 1 1 2 1 1 4 1 1 1 2 1 1 4 1 1 1 2 2 1 2 1 8 2 1 2 1 2 1 7 2 2 1 2 1 6 2 1 2 1 2 1 5 2 2 1 2 1 4 2 1 2 1 2 1 3 2 2 1 2 1 2 2 1 2 1 2 1 1 s s s s N s s s s s s N s s s s N s s s s s s N s s s s N s s s s s s N s s s s N s s s s s s N 35 Forças Nodais Equivalentes Álvaro F M Azevedo 164 As interpolações das grandezas correspondentes ao bordo carregado são efectuadas com as seguintes funções de forma que se obtêm substituindo s2 por 1 em 35 Notese ainda que em todos os pontos do domínio de integração do integral 23 a variável s2 assume o valor 1 0 0 0 0 0 2 1 2 1 8 1 7 1 6 1 5 1 4 1 2 1 1 3 2 1 1 2 1 2 1 1 1 s N s N s N s N s N s s s N s s N s s s N 36 Estas funções de forma coincidem com as que foram obtidas no Capítulo 4 para o elemento unidimensional de três nós Na Figura 94 está representado o eixo tangente ao bordo 1x bem como o eixo normal ao bordo 2 x O eixo tangente ao bordo segue a numeração local dos nós O eixo 2 x forma com 1x um referencial directo x1 1 2 3 x2 s1 2p 1x 2 x 1p s1 1 s1 0 s1 1 L d L Fig 94 Bordo de três nós com uma carga distribuída por unidade de comprimento Forças Nodais Equivalentes Álvaro F M Azevedo 165 A carga distribuída p x1 x2 é decomposta nas suas componentes tangencial 1 p e normal 2 p Numa análise por elementos finitos são habitualmente conhecidos os valores nodais das componentes tangencial e normal da carga distribuída que se designam por ij p ie valor da carga distribuída no nó i segundo a direcção jx Todas as grandezas relativas às cargas distribuídas são forças por unidade de comprimento de arco A interpolação das componentes tangencial e normal da carga distribuída a partir dos correspondentes valores nodais é efectuada da forma habitual recorrendo às funções de forma 36 32 1 3 22 1 2 12 1 1 1 2 31 1 3 21 1 2 11 1 1 1 1 p s N p s N p N s s p p s N p s N p N s s p 37 ou 3 2 1 32 22 12 31 21 11 2 1 N N N p p p p p p p p 38 p T NV p 39 Designando por T a matriz de transformação do referencial x1 x2 para o referencial x1 x2 temse p T p 40 e a relação inversa p T p T 41 Substituindo 39 em 41 chegase a V T T N p p T 42 Forças Nodais Equivalentes Álvaro F M Azevedo 166 A primeira linha da matriz de transformação T utilizada em 40 é constituída pelo versor 1ˆn e a segunda pelo versor 2ˆn ver a Figura 95 x1 x2 s1 1x 2 x d L 1ˆn 2ˆn d L d x1 d x2 Fig 95 Referencial tangente ao bordo O versor 1ˆn obtémse com a seguinte expressão 1 2 1 1 1 1 ˆ d s d x d s d x J n 43 sendo J a norma do vector 1 2 1 1 d s d x d s x d 2 1 2 2 1 1 d s d x d s d x J 44 O versor 2ˆn é ortogonal a 1ˆn e forma com 1ˆn um referencial directo sendo a seguinte a sua expressão 1 1 1 2 2 1 ˆ d s d x d s d x J n 45 Os elementos da matriz de transformação T são os seguintes Forças Nodais Equivalentes Álvaro F M Azevedo 167 1 1 1 2 1 2 1 1 1 d s d x d s d x d s d x d s x d J T 46 Os elementos da matriz T são calculados com base na seguinte interpolação das coordenadas de um ponto genérico do arco L 32 1 3 22 1 2 12 1 1 1 2 31 1 3 21 1 2 11 1 1 1 1 x s N x s N x N s s x x s N x s N x N s s x 47 Nesta expressão Ni são as funções de forma associadas aos nós do arco 36 e ijx representa a coordenada cartesiana do nó i segundo a direcção xj Derivando ambos os membros em ordem a s1 chegase a 32 1 3 22 1 2 12 1 1 1 2 31 1 3 21 1 2 11 1 1 1 1 x d s d N x d s d N x d s d N d s x d x d s d N x d s d N x d s d N d s x d 48 As derivadas em ordem a s1 das funções de forma 36 são 2 1 2 2 1 1 1 3 1 1 2 1 1 1 s d s N d s d s N d s d s N d 49 Para calcular as forças nodais equivalentes à carga distribuída no bordo devese utilizar a expressão 23 que em seguida se reproduz L T p p d L N F 50 Forças Nodais Equivalentes Álvaro F M Azevedo 168 Para facilitar o recurso à quadratura de Gauss ver o Capítulo 5 deve ser efectuada a seguinte mudança de variável 1 1 1 1 d s d s p d L N F T p 51 De acordo com a Figura 95 verificase que 2 2 1 2 d x d x d L 52 Atendendo a 47 temse 1 1 2 2 1 1 1 1 d s d s d x x d d s d s d x x d 53 Substituindo 53 em 52 chegase a 2 1 2 2 1 1 1 d s d x d s d x d s d L 54 Comparando 54 com 44 concluise que J d s d L 1 55 e a expressão 51 passa a ser 1 1 p J d 1s N F T p 56 Substituindo 42 em 56 obtémse 1 1 J d 1s N p T N F V T T T p 57 Uma vez que em todo o domínio de integração se verifica ser s2 1 na matriz N devem ser utilizadas as funções de forma 36 Forças Nodais Equivalentes Álvaro F M Azevedo 169 Considerando todas as expressões já deduzidas o vector das forças nodais equivalentes à carga distribuída indicada nas Figuras 93 e 94 é o seguinte 1 1 1 3 2 1 32 22 12 31 21 11 1 1 1 2 1 2 1 1 3 3 2 2 1 1 1 16 1 0 0 0 0 0 0 0 0 0 0 J d s N N N p p p p p p d s d x d s d x d s d x d s x d J N N N N N N F p M M 58 que se simplifica para 1 1 1 3 2 1 32 22 12 31 21 11 1 1 1 2 1 2 1 1 3 3 2 2 1 1 32 31 22 21 12 11 0 0 0 0 0 0 0 0 0 0 0 0 d s N N N p p p p p p d s d x d s d x d s d x d s x d N N N N N N F F F F F F p p p p p p M M M 59 No vector F p apenas as seis primeiras componentes são não nulas ie nos nós 1 2 e 3 ver a Figura 93 existem forças nodais equivalentes enquanto que nos restantes cinco nós a contribuição da carga distribuída é nula O integral 59 pode ser calculado recorrendo à quadratura de Gauss ver o Capítulo 5 Todos os componentes da função integranda de 59 são funções de s1 ou são constantes de valor conhecido como é o caso da matriz p Forças Nodais Equivalentes Álvaro F M Azevedo 170 95 Carga distribuída por unidade de superfície O cálculo das forças nodais equivalentes a uma carga distribuída por unidade de superfície é efectuado com a expressão 24 Este tipo de cargas só tem interesse prático em elementos de laje elementos de casca e faces de elementos sólidos bricks O processo de cálculo de F q é semelhante ao apresentado na Secção 94 sendo necessário adaptálo às características dos referidos elementos O domínio de integração passa a ser uma superfície 96 Carga distribuída por unidade de volume Este tipo de acção é devido à presença de forças de volume b x Estas forças estão presentes sempre que o corpo se encontra sujeito a uma aceleração O caso mais comum é o da aceleração da gravidade que se define do seguinte modo g b ρ 60 Nesta expressão ρ é a massa específica do material e g é a aceleração da gravidade No caso mais comum ie supondo que o eixo x3 é vertical e orientado para cima que a aceleração da gravidade actua segundo x3 e é negativa e que se utilizam as unidades do Sistema Internacional SI temse 9 81 0 0 3 2 1 ρ b b b 61 Em 61 a aceleração da gravidade foi considerada igual a 2 9 81 m s As unidades de b e de ρ devem ser as seguintes bi em N m3 e ρ em kg m3 ou bi em kN m3 e ρ em t m3 ou bi em MN m3 e ρ em kt m3 Ao definir o peso próprio deste modo é facilitada a sua combinação com outras componentes da aceleração g Forças Nodais Equivalentes Álvaro F M Azevedo 171 Se a única força de volume for a devida ao peso próprio então podese atribuir a ρ o valor do peso específico do material e considerar 00 1 g Deste modo fica facilitada a preparação dos dados de uma análise por elementos finitos em que é utilizado um sistema de unidades distinto do SI As forças nodais equivalentes às forças de volume são calculadas com a expressão 25 que em seguida se reproduz V T b N b dV F 62 Na Figura 96 encontrase representado um elemento finito de quatro nós destinado à discretização de estados planos de tensão u1 x1 x2 x1 1 2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 s1 s2 s1s2 h dV b d V d S Fig 96 Elemento finito de quatro nós sujeito a forças de volume No elemento representado na Figura 96 actuam forças de volume b x cujas componentes são 2 1 b b b 63 No caso do estado plano de tensão o integral 62 passa a S T b N b h d S F 64 Forças Nodais Equivalentes Álvaro F M Azevedo 172 Nesta expresão h corresponde à espessura do elemento finito que pode eventualmente ser não constante A sua interpolação a partir das espessuras nos nós ih é efectuada com a seguinte expressão ver o Capítulo 6 4 2 1 4 3 2 1 3 2 2 1 2 1 2 1 1 2 1 h s s N h s s N h s s N h N s s h s s 65 De um modo semelhante ao que foi efectuado no Capítulo 6 para a matriz de rigidez deve ser efectuada em 64 a seguinte mudança de variável 1 1 1 1 2 1 N b h J d s d s F T b 66 Nesta expressão J é o determinante Jacobiano definido no Capítulo 6 Uma vez que N é a matriz que relaciona os deslocamentos nodais com o campo de deslocamentos u N a ver o Capítulo 6 chegase à seguinte expressão final 1 1 1 1 2 1 2 1 4 4 3 3 2 2 1 1 42 41 32 31 22 21 12 11 0 0 0 0 0 0 0 0 h J d s d s b b N N N N N N N N F F F F F F F F b b b b b b b b 67 Nesta expressão os componentes da função integranda são funções de s1 e s2 ou são constantes O integral 67 pode ser calculado recorrendo à quadratura de Gauss ver o Capítulo 5 97 Considerações finais As deduções relativas a casos particulares que foram apresentadas neste capítulo podem ser facilmente adaptadas a outros casos tais como elementos finitos com mais nós outros tipos de elementos finitos meios com rotações e momentos etc Forças Nodais Equivalentes Álvaro F M Azevedo 173 BIBLIOGRAFIA 91 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 92 Hinton E Owen D R Finite Element Programming Academic Press 1980 93 Kreyszig E Advanced Engineering Mathematics Sixth Edition John Wiley Sons Inc 1988 94 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 Forças Nodais Equivalentes Álvaro F M Azevedo 174 175 CAPÍTULO 10 SÓLIDOS ESTADO PLANO DE DEFORMAÇÃO E AXISSIMETRIA Neste capítulo são descritas algumas particularidades dos elementos sólidos tridimensionais do estado plano de deformação e do estado axissimétrico Pressupõese que já é conhecida com detalhe a formulação descrita no Capítulo 6 101 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada neste capítulo Tabela 101 Simbologia relativa ao método dos elementos finitos m Número de direcções consideradas no caso tridimensional m 3 n Número de nós do elemento finito p Número de graus de liberdade do elemento finito p n x m x Coordenada cartesiana u Campo de deslocamentos a Deslocamento nodal F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial geral K Matriz de rigidez do elemento finito no referencial geral x Coordenada cartesiana de um nó de um elemento finito s Coordenada local curvilínea s Coordenada local de um nó de um elemento finito N Função interpoladora ou função de forma V Vector contendo os factores não constantes de um polinómio Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 176 NV Vector das funções interpoladoras ou funções de forma B Matriz de deformação D Matriz de elasticidade ε σ D V Volume J Jacobiano da transformação ε Extensão γ Distorção σ Tensão normal τ Tensão tangencial C Elemento da matriz de elasticidade D E Módulo de elasticidade ou módulo de Young ν Coeficiente de Poisson q Número de componentes do vector ε e do vector σ L Operador diferencial h Espessura do elemento finito laminar Q Carga concentrada p Carga distribuída por unidade de comprimento θ Ângulo direcção circunferencial P Perímetro S Superfície 102 Elementos sólidos tridimensionais bricks No desenvolvimento de elementos sólidos do tipo brick é considerada uma formulação genérica com três graus de liberdade do tipo deslocamento A exposição aqui apresentada baseiase num elemento finito sólido com oito nós ver a Figura 101 O número de graus de liberdade deste elemento é p 8 x 3 24 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 177 x1 4 5 2 1 x2 3 a52 a51 6 x3 7 8 a53 u1 x1 x2 x3 u2 x1 x2 x3 u3 x1 x2 x3 Fig 101 Elemento finito sólido de oito nós com geometria arbitrária Os vectores dos deslocamentos nodais e das forças nodais equivalentes às acções exteriores são os seguintes 83 82 81 23 22 21 13 12 11 1 83 82 81 23 22 21 13 12 11 1 F F F F F F F F F F a a a a a a a a a a p p M M 1 A matriz de rigidez do elemento K é uma matriz p x p 24 x 24 No referencial geral a matriz das coordenadas cartesianas dos nós do elemento é a seguinte 83 82 81 23 22 21 13 12 11 x x x x x x x x x x n m M M M 2 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 178 Pelos motivos referidos no Capítulo 6 é conveniente fazer a seguinte substituição de variáveis 3 2 1 3 3 3 2 1 2 2 3 2 1 1 1 s x s s x s s s x x s x s s x 3 Na Figura 102 encontrase indicado o sistema de coordenadas locais bem como o novo domínio de integração s1 4 5 2 1 s2 3 6 s3 7 8 1 1 1 1 1 1 3 2 1 s s s Fig 102 Sistema de coordenadas locais Os valores nodais das coordenadas s1 s2 e s3 são os seguintes 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 83 82 81 73 72 71 63 62 61 53 52 51 43 42 41 33 32 31 23 22 21 13 12 11 s s s s s s s s s s s s s s s s s s s s s s s s s n m 4 Uma vez que o elemento é da família Lagrangeana as funções de forma são Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 179 8 1 1 1 8 1 1 1 8 1 1 1 8 1 1 1 8 1 1 1 8 1 1 1 8 1 1 1 8 1 1 1 3 2 1 3 2 1 8 3 2 1 3 2 1 7 3 2 1 3 2 1 6 3 2 1 3 2 1 5 3 2 1 3 2 1 4 3 2 1 3 2 1 3 3 2 1 3 2 1 2 3 2 1 3 2 1 1 s s s s s s N s s s s s s N s s s s s s N s s s s s s N s s s s s s N s s s s s s N s s s s s s N s s s s s s N 5 No caso da determinação das funções de forma com o procedimento genérico descrito no Capítulo 7 é necessário seleccionar oito termos na pirâmide de Pascal que se encontra representada na Figura 103 s3 s1 s2 1 2 1s s1 s2 2 2s 2 3s s2 s3 s1 s3 Fig 103 Pirâmide de Pascal Os termos seleccionados são agrupados no vector V que é neste caso o seguinte 3 2 1 3 1 3 2 2 1 3 2 1 1 s s s s s s s s s s s s V 6 Nota o termo de terceiro grau não se encontra representado na Figura 103 O procedimento para a determinação das funções de forma é semelhante ao que foi descrito no Capítulo 7 A interpolação das coordenadas cartesianas é efectuada com a seguinte expressão Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 180 8 2 1 83 23 13 82 22 12 81 21 11 3 2 1 N N N x x x x x x x x x x x x M L L L 7 que corresponde a 1 1 n V n m T m N x x 8 De acordo com o que foi apresentado no Capítulo 9 a matriz de rigidez do elemento finito é calculada com a expressão genérica V T B D B dV K 9 Após a substituição de variáveis 3 passa a terse 1 1 1 1 1 1 3 2 1 B D B J d s d s d s K T 10 No caso tridimensional a matriz Jacobiana J é a seguinte 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 s x s x s x s x s x s x s x s x s x J 11 e o determinante Jacobiano é J J 12 No caso dos materiais isotrópicos é a seguinte a relação entre tensões e deformações 101 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 181 12 31 23 3 2 1 3 3 3 1 2 2 2 1 2 2 2 1 12 31 23 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 γ γ γ ε ε ε τ τ τ σ σ σ C C C C C C C C C C C C 13 sendo ν ν ν 1 2 1 1 1 E C ν ν ν 1 2 1 2 E C ν 2 1 3 E C 14 Em 13 o número de componentes de σ e ε q é 6 podendo escreverse de um modo mais compacto 1 1 q q q q D ε σ 15 A matriz de elasticidade D depende do módulo de Young E e do coeficiente de Poisson ν A matriz Jacobiana obtémse com a seguinte expressão que resulta da derivação de 7 em ordem a s1 s2 e s3 3 8 2 8 1 8 3 2 2 2 1 2 3 1 2 1 1 1 83 23 13 82 22 12 81 21 11 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 s N s N s N s N s N s N s N s N s N x x x x x x x x x s x s x s x s x s x s x s x s x s x M M M L L L 16 De um modo mais compacto temse Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 182 n m n m T m m s N x J 17 Os elementos da matriz N s são as derivadas de 5 em ordem a s1 s2 e s3 No caso geral tridimensional a relação entre o campo dos deslocamentos e o campo das deformações é 101 3 2 1 1 2 1 3 2 3 3 2 1 12 31 23 3 2 1 0 0 0 0 0 0 0 0 0 u u u x x x x x x x x x γ γ γ ε ε ε 18 ou de modo mais compacto 1 1 m q m q u L ε 19 A interpolação do campo de deslocamentos é efectuada com 83 82 81 23 22 21 13 12 11 8 2 1 8 2 1 8 2 1 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a a a a a a a a a N N N N N N N N N u u u M L 20 ou Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 183 1 1 p p m m a N u 21 Substituindo 21 em 19 chegase a 1 1 p p m q m q a N L ε 22 sendo p m q m q p N L B 23 8 2 1 8 2 1 8 2 1 1 2 1 3 2 3 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 N N N N N N N N N x x x x x x x x x B L L L 24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 8 2 8 1 2 2 2 1 1 2 1 1 8 3 8 1 2 3 2 1 1 3 1 2 8 3 8 2 2 3 2 2 1 3 1 3 8 3 2 3 1 2 8 2 2 2 1 1 8 1 2 1 1 x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N x N B L L L L L L 25 Os elementos da matriz B obtêmse com a seguinte expressão que é uma generalização do que foi exposto no Capítulo 6 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 184 1 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 3 8 2 8 1 8 3 2 2 2 1 2 3 1 2 1 1 1 3 8 2 8 1 8 3 2 2 2 1 2 3 1 2 1 1 1 s x s x s x s x s x s x s x s x s x s N s N s N s N s N s N s N s N s N x N x N x N x N x N x N x N x N x N M M M M M M 26 Em notação matricial temse m m n m m n J s N x N 1 27 Uma vez que todos os componentes da função integranda de 10 se encontram definidos em função de s1 s2 e s3 já é possível calcular o integral triplo recorrendo à quadratura de Gauss descrita no Capítulo 5 O algoritmo de cálculo da matriz de rigidez do elemento sólido tridimensional é semelhante ao que foi apresentado no Capítulo 6 sendo apenas necessário remover todas as operações relativas à espessura do elemento e adaptar as dimensões das matrizes envolvidas no cálculo O cálculo das tensões e deformações no elemento finito são também efectuadas de acordo com o que foi exposto no Capítulo 6 desde que se efectuem as necessárias adaptações ao caso tridimensional 103 Estado plano de deformação As características de um estado plano de deformação encontramse descritas em 101 Na Figura 104 está representado um corpo prismático de secção constante e eixo segundo x3 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 185 x1 x2 x3 Plano médio 12 12 h 1 O Fig 104 Estado plano de deformação De acordo com as hipóteses consideradas para o estado plano de deformação são efectuadas as seguintes simplificações 12 2 1 12 31 23 3 2 1 1 0 0 0 γ ε ε γ γ γ ε ε ε ε q 28 Atendendo a 13 temse 12 2 1 3 3 3 1 2 2 2 1 2 2 2 1 12 31 23 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 γ ε ε τ τ τ σ σ σ C C C C C C C C C C C C 29 Os parâmetros C1 C2 e C3 encontramse definidos em 14 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 186 A equação 29 simplificase da seguinte forma 12 2 1 3 1 2 2 1 12 2 1 0 0 0 0 γ ε ε τ σ σ C C C C C 30 sendo a tensão normal σ 3 dada por 2 1 2 3 ε ε σ C 31 As tensões tangenciais τ23 e τ31 são nulas Em 101 encontrase a relação inversa de 13 cuja terceira equação é 3 2 1 3 1 σ ν σ ν σ ε E E E 32 Como no estado plano de deformação ε3 0 resulta 2 1 3 σ ν σ σ 33 Uma vez que se supõe que os campos de tensões e deformações não variam com x3 é suficiente estudar o comportamento de um troço de comprimento unitário cujo plano médio passa pela origem e é perpendicular a x3 ver a Figura 104 Nestas circunstâncias apenas é necessário discretizar com elementos finitos a secção transversal que é perpendicular a x3 e que passa pela origem ver a Figura 105 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 187 u1 x1 x2 x1 1 2 3 4 x2 u2 x1 x2 a41 a42 a31 a32 a21 a22 a11 a12 h 1 Q p Fig 105 Elemento finito para estados planos de deformação A formulação de um elemento finito para estados planos de deformação coincide com a que foi apresentada no Capítulo 6 com excepção do seguinte em todas as expressões deve ser considerado h 1 a relação constitutiva σ D ε passa a ser a equação 30 sendo os elementos da matriz D definidos em 14 a força Q que actua no corpo ver a Figura 105 passa a ser uma força por unidade de comprimento segundo x3 a força por unidade de comprimento p que actua no corpo ver a Figura 105 passa a ser uma força por unidade de superfície 104 Estado axissimétrico Na Figura 106 encontrase representada uma secção transversal que ao ser rodada em torno de x2 gera um sólido de revolução que corresponde a um reservatório axissimétrico A secção transversal está discretizada em elementos finitos quadriláteros contidos no plano x1 x2 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 188 x2 Eixo de axissimetria x1 u1 u2 uθ θ A Fig 106 Reservatório axissimétrico Se além da geometria for também axissimétrica a distribuição do material que constitui o reservatório bem como as características das acções exteriores verificase que o comportamento é também axissimétrico Devido ao facto de as acções axissimétricas serem na generalidade dos casos de natureza gravítica admitese que o eixo de axissimetria x2 é sempre vertical O eixo x1 bem como os componentes segundo x1 de todas as grandezas recebem a designação de radiais O componente uθ do campo de deslocamentos designase deslocamento circunferencial A deformação circunferencial εθ corresponde ao quociente entre a variação do perímetro de uma fibra e o perímetro original ver a Figura 107 x2 x1 A A xA 1 uA 1 Fig 107 Definição da deformação circunferencial Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 189 O perímetro inicial da circunferência que passa pelo ponto A é 1 0 2 xA P π 34 O perímetro da circunferência após a deformação é 1 1 2 A A u x P π 35 A deformação circunferencial é 0 0 P P P θε 36 Substituindo 34 e 35 em 36 resulta 1 1 1 1 1 1 2 2 2 2 A A A A A A x u x x u x π π π π εθ 37 No caso de um ponto genérico temse 1 1 x u θε 38 Por não corresponderem a uma deformação axissimétrica são nulas as distorções γ 1θ e γ 2θ ver a Figura 106 Nestas circunstâncias o vector ε é o seguinte 1 2 2 1 1 1 2 2 1 1 12 2 1 x u x u x u x u x u γ ε ε ε θ 39 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 190 De um modo análogo ao que foi exposto no Capítulo 6 temse 2 1 1 2 1 2 1 12 2 1 0 1 0 0 u u x x x x x γ ε ε ε θ 40 que de um modo mais compacto corresponde a ε L u 41 Nestas circunstâncias a relação constitutiva é 12 2 1 3 1 2 2 2 1 2 2 2 1 12 2 1 0 0 0 0 0 0 γ ε ε ε τ σ σ σ θ θ C C C C C C C C C C 42 Os parâmetros C1 C2 e C3 encontramse definidos em 14 Na Figura 108 estão representados os componentes não nulos do tensor das tensões em problemas axissimétricos No cálculo da matriz de rigidez K deve ser considerado o seguinte ver o Capítulo 9 V T B D B dV K 43 S T d S x B D B K 2π 1 44 Nesta expressão 2 π x1 é o perímetro da fibra correspondente aos pontos de abcissa x1 Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 191 x2 x1 σ 1 σ 2 σ θ τ 12 τ 12 γ 1θ γ 2θ 0 τ 1θ τ 2θ 0 Fig 108 Componentes não nulos do tensor das tensões em problemas axissimétricos Após a mudança de variável temse 1 1 1 1 2 1 1 2 J d s d s x B D B K T π 45 Uma vez que se pretende que todos os componentes da função integranda de 45 sejam funções de s1 e s2 devese calcular x1 com a seguinte expressão 1 2 1 11 2 1 1 2 1 1 n n x s s N x N s s x s s L 46 em que n é o número de nós do elemento finito A matriz B e o determinante Jacobiano J são calculados tal como foi exposto no Capítulo 6 supondo que os elementos finitos do problema axissimétrico se situam no plano x1 x2 ver a Figura 106 O integral duplo 45 pode ser calculado recorrendo à quadratura de Gauss ver o Capítulo 5 Os integrais correspondentes ao cálculo das forças nodais equivalentes às acções exteriores têm de ser adaptados ao caso axissimétrico de um modo semelhante ao que foi apresentado para a matriz de rigidez K Sólidos Estado Plano de Deformação e Axissimetria Álvaro F M Azevedo 192 Em problemas axissimétricos os apoios assentamentos de apoio e as cargas concentradas em nós estão presentes numa linha que se obtém por rotação do correspondente ponto em torno de x2 ver a Figura 106 105 Considerações finais Neste capítulo a descrição das formulações foi efectuada de um modo mais sucinto uma vez que se supunha um bom conhecimento da formulação do estado plano de tensão exposta no Capítulo 6 Os elementos sólidos tridimensionais bricks são os que aparentemente permitem modelar qualquer geometria No entanto apresentam os inconvenientes de necessitarem de uma preparação dos dados mais trabalhosa requererem um maior esforço computacional e apresentarem maior dificuldade na interpretação dos resultados Os estados planos de deformação surgem tipicamente em problemas geotécnicos eg muros de suporte barragens gravidade túneis Os estados axissimétricos apresentam o inconveniente de necessitarem de uma acção axissimétrica que na prática nem sempre ocorre Para ultrapassar esta limitação em 102 é desenvolvida uma formulação correspondente a sólidos de revolução com acções não axissimétricas Contudo o seu elevado grau de dificuldade constitui um entrave a uma aplicação generalizada BIBLIOGRAFIA 101 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 102 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 103 Zienkiewicz O C Taylor R L The Finite Element Method Fourth Edition McGrawHill 1988 193 CAPÍTULO 11 FLEXÃO DE VIGAS Antecedendo a apresentação da formulação de diversos tipos de elementos de viga efectuase em seguida uma revisão dos fundamentos da flexão de vigas Apenas são consideradas as deformações devidas às tensões normais 111 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação da flexão de vigas Tabela 111 Simbologia relativa à flexão de vigas g Referencial geral l Referencial local i Primeiro nó de uma barra j Segundo nó de uma barra x Coordenada cartesiana u Campo de deslocamentos θ Ângulo G Centro de gravidade ε Extensão σ Tensão normal E Módulo de elasticidade ou módulo de Young S Superfície N Esforço axial Flexão de Vigas Álvaro F M Azevedo 194 M Momento flector k Curvatura A Área da secção transversal da barra prismática I Momento de inércia da secção transversal da barra prismática 112 Flexão composta plana Considerese uma barra prismática de eixo rectilíneo e secção variável representada na Figura 111 conjuntamente com os referenciais geral g e local l g 1 g 2 g 3 l 1 l 2 l 3 i j i j Fig 111 Barra i j referencial geral g e referencial local l O eixo l1 do referencial local coincide com o eixo da barra e está orientado do nó i para o nó j correspondendo i e j à numeração global dos nós da malha Os eixos l2 e l3 são os eixos principais centrais de inércia da secção transversal Apesar de a secção ser variável considerase que a localização destes eixos é constante ao longo da barra A transformação dos deslocamentos generalizados e das forças generalizadas entre os referenciais l e g encontrase descrita nos Capítulos 2 e 3 No estudo que se segue apenas é considerado o referencial l que passa a ser designado por x Supõese também que todas as acções estão contidas no plano x1 x3 sendo os Flexão de Vigas Álvaro F M Azevedo 195 deslocamentos segundo x2 considerados nulos Na Figura 112 está representado o eixo da barra indeformado e deformado para o caso da flexão plana Encontrase também representada a secção transversal da barra cujos eixos principais centrais de inércia são x2 e x3 x1 x3 u3 x1 u3 u1 θ G x2 x3 Fig 112 Barra deformada e secção transversal Uma vez que os deslocamentos são muito pequenos considerase tan θ θ 1 sendo θ o ângulo de rotação do eixo da barra A função u3 ver a Figura 112 corresponde ao deslocamento do eixo da barra que apenas depende de x1 De acordo com 1 e com a equação da deformada que se encontra representada na Figura 112 temse 1 3 1 d x d u x θ 2 Na Figura 113 encontramse representados os deslocamentos de três pontos de uma secção transversal A O e B Flexão de Vigas Álvaro F M Azevedo 196 O x3 u3 u1 x1 u3 x1 A B A B θ x3 θ u1 x1 x3 A A B B O O O uO1 x1 Fig 113 Barra deformada e deslocamentos da secção transversal Admitese que uma secção plana se mantém plana após a deformação Admitese também que uma secção perpendicular ao eixo da barra mantém esta característica após a deformação O ponto O apresenta coordenada nula segundo x3 O deslocamento do ponto O segundo x1 é designado uO1 e depende apenas de x1 O deslocamento de um ponto genérico da secção transversal segundo x1 depende de x1 e de x3 e é definido pela seguinte expressão 1 3 1 1 3 1 1 x x x u x x u O θ 3 ou 3 θ 1 1 x u u O 4 A extensão segundo x1 é dada por 1 3 1 1 3 1 1 1 1 1 d x d x d x d u x u x x u O O θ θ ε 5 A extensão ε1 é positiva quando existe um alongamento Substituindo 2 em 5 obtémse 2 1 3 2 3 1 1 1 d x x d u d x d u O ε 6 Flexão de Vigas Álvaro F M Azevedo 197 Uma vez que se consideram muito pequenas as dimensões da secção transversal em comparação com o comprimento da barra podese desprezar o efeito das tensões normais σ2 e σ3 ficando a lei de Hooke reduzida a ε1 σ1 E ou 1 3 1 1 1 1 d x d E x d x E d u E O θ ε σ 7 sendo E o módulo de Young que é sempre positivo 111 A um valor positivo da tensão σ1 corresponde uma tracção Substituindo 2 em 7 resulta 2 1 3 2 3 1 1 1 d x E x d u d x d u E O σ 8 A resultante das tensões normais na secção transversal é ver a Figura 114 S O S d S d x d E x d x E d u d S N 1 3 1 1 1 θ σ 9 sendo S a área da secção transversal ver a Figura 112 De um modo semelhante se define o momento flector como sendo S O S x d S d x d E x d x E d u x d S M 3 1 3 1 1 3 1 θ σ 10 Considerase que um momento flector é positivo quando provoca tracções nas fibras que têm coordenada x3 positiva ver a Figura 114 Flexão de Vigas Álvaro F M Azevedo 198 x3 x1 N M Curvatura k positiva Momento flector positivo Esforço axial positivo Fig 114 Definição do esforço axial e do momento flector Supondo que o módulo de Young E é constante em todos os pontos da barra e passando para fora do integral tudo o que não depende de x2 nem de x3 resulta de 9 e 10 S S O x d S d x E d d S d x E d u N 3 1 1 1 θ 11 e S S O x d S d x E d x d S d x E d u M 2 3 1 3 1 1 θ 12 Uma vez que os eixos x2 e x3 são principais centrais de inércia o seguinte momento estático é nulo 0 3 S x d S 13 A área e o momento de inércia em relação a x2 são definidos com as seguintes expressões ver a Figura 112 S d S A 14 S x d S I 2 3 15 Flexão de Vigas Álvaro F M Azevedo 199 É assim possível simplificar 11 e 12 para 1 1 d x E A d u N O 16 d 1x E I d M θ 17 Designando a extensão correspondente ao eixo da barra por εo temse 1 1 d x d u O εo 18 E A o N ε 19 que corresponde à expressão clássica relativa à tracção de barras 112 Substituindo 2 em 17 temse 2 1 3 2 d x E I d u M 20 Designando por k a curvatura da barra 2 1 3 2 d x k d u 21 resulta E I k M 22 ou E I M k 23 que corresponde a uma das expressões clássicas da flexão de vigas 112 Na Figura 112 todos os pontos da linha que representa o eixo da barra deformada apresentam ordenada positiva u3 0 primeira derivada positiva 0 1 3 d x d u e Flexão de Vigas Álvaro F M Azevedo 200 segunda derivada também positiva 0 2 1 2 3 d x d u Atendendo à definição do sentido positivo do momento flector ver a Figura 114 verificase que um momento flector positivo provoca uma curvatura negativa Esta questão encontrase expressa na equação 23 em que E e I são sempre positivos A expressões 16 e 17 são equivalentes às seguintes E A N d x d u o 1 24 E I M d x d 1 θ 25 Substituindo 24 e 25 em 7 temse E I M E x E A N E 3 1 σ 26 resultando 3 1 x I M A σ N 27 que corresponde à expressão clássica da flexão composta 112 113 Considerações finais A formulação apresentada neste capítulo apresenta a vantagem de ser facilmente estendida a problemas em que os eixos x2 e x3 não são principais centrais de inércia Apresenta também a vantagem de recorrer a um conjunto de convenções coerente com o que é habitual considerar no método dos elementos finitos Fica assim facilitada a formulação de elementos finitos de viga bem como a sua combinação com outros tipos de elementos Flexão de Vigas Álvaro F M Azevedo 201 BIBLIOGRAFIA 111 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 112 Massonnet C Résistance des Matériaux Dunod Paris 1968 Flexão de Vigas Álvaro F M Azevedo 202 203 CAPÍTULO 12 VIGA DE EULERBERNOULLI Designase por EulerBernoulli a formulação do elemento finito de viga em que se considera que as secções se mantêm planas e normais ao eixo da barra após a deformação Deste modo não é considerada a deformação devida ao corte 121 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação do elemento de viga de EulerBernoulli Tabela 121 Simbologia relativa ao elemento de viga de EulerBernoulli L Comprimento da barra prismática x Coordenada cartesiana u Campo de deslocamentos a Deslocamento generalizado nodal Deslocamento nodal θ Rotação nodal x Coordenada cartesiana de um nó de um elemento finito N Função interpoladora ou função de forma ε Deformação B Matriz de deformação σ Tensão normal E Módulo de elasticidade ou módulo de Young p Acção exterior distribuída por unidade de comprimento Viga de EulerBernoulli Álvaro F M Azevedo 204 F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial local V Volume S Superfície I Momento de inércia da secção transversal da barra prismática K Matriz de rigidez do elemento finito no referencial local M Momento flector s Coordenada local s Coordenada local de um nó de um elemento finito J Jacobiano da transformação J d x1 d s 122 Viga de dois nós sem substituição de variável Na Figura 121 encontrase representado um elemento de viga com dois nós e com comprimento L ver os Capítulos 7 e 11 u3 x1 x1 a1 1 1 2 2 11 1 L x x 2 21 1 L x x L 2 L 2 a3 2 a2 θ1 a4 θ2 x3 Fig 121 Elemento de viga com dois nós Os deslocamentos generalizados dos nós do elemento finito representado na Figura 121 são os seguintes Viga de EulerBernoulli Álvaro F M Azevedo 205 2 2 1 1 4 3 2 1 θ θ a a a a a 1 De acordo com o que foi exposto no Capítulo 11 o deslocamento lateral é u3x1 A coordenada cartesiana ijx corresponde ao nó i e referese ao eixo xj A interpolação do campo de deslocamentos é efectuada com a seguinte expressão ver o Capítulo 7 4 1 4 3 1 3 2 1 2 1 1 1 1 3 a x N a N x a x N a N x x u 2 que em notação matricial se escreve 4 3 2 1 4 3 2 1 1 3 a a a a N N N N x u 3 ou N a u 3 4 sendo 4 3 2 1 N N N N N 5 As funções de forma N1 a N4 são as correspondentes à interpolação Hermitiana e têm as seguintes expressões ver o Capítulo 7 3 1 3 1 1 1 2 2 3 2 1 L x L x x N 6 3 1 2 2 1 1 1 2 1 2 1 4 1 8 L x L x x L x N 7 Viga de EulerBernoulli Álvaro F M Azevedo 206 3 1 3 1 1 3 2 2 3 2 1 L x L x x N 8 3 1 2 2 1 1 1 4 1 2 1 4 1 8 L x L x x L x N 9 Considerando apenas os deslocamentos laterais u3x1 ie considerando constante a componente uO1 do campo de deslocamentos temse de acordo com o que foi exposto no Capítulo 11 2 1 3 2 3 1 d x ε x d u 10 Designando por 1ε a seguinte componente da expressão 10 2 1 3 2 1 d x ε d u 11 passa a terse 1 3 1 ε ε x 12 Substituindo 2 em 11 chegase a 4 3 2 1 2 1 4 2 2 1 3 2 2 1 2 2 2 1 1 2 1 a a a a d x d N d x d N d x d N d x d N ε 13 Definindo a matriz B da seguinte forma 2 1 4 2 2 1 3 2 2 1 2 2 2 1 1 2 d x d N d x d N d x d N d x d N B 14 passa a escreverse B a 1ε 15 Viga de EulerBernoulli Álvaro F M Azevedo 207 Substituindo 15 em 12 obtémse ε1 x3 B a 16 Atendendo às funções de forma 6 a 9 são os seguintes os componentes da matriz B 1 2 1 3 1 2 3 1 6 1 12 6 1 12 L x L L x L x L L x B 17 Também de acordo com o que foi exposto no Capítulo 11 temse 1 1 ε σ E 18 e atendendo a 16 B a x E 3 σ1 19 Considerese que na viga da Figura 121 actua a carga uniformemente distribuída representada na Figura 122 x3 x1 F1 1 2 L 2 L 2 F2 p F3 F4 Fig 122 Carga uniformemente distribuída e respectivas forças nodais equivalentes As forças nodais equivalentes à acção exterior encontramse também representadas na Figura 122 e apresentam os mesmos sentidos positivos que foram considerados para os deslocamentos generalizados ai O princípio dos trabalhos virtuais PTV que foi apresentado no Capítulo 4 corresponde ao seguinte Viga de EulerBernoulli Álvaro F M Azevedo 208 L T V T u p d L dV δ δε σ 20 No caso da viga representada nas Figuras 121 e 122 a equação 20 passa a 2 2 1 3 2 2 1 1 1 L L L L S u p d x d S d x δ δ ε σ 21 Nesta equação S é a superfície correspondente à secção transversal da barra ver o Capítulo 11 e 3 2 d x d x d S 22 A equação 16 referida à deformação virtual é a seguinte a x B δ δ ε 3 1 23 que é equivalente a 3 1 x B a T T δ δ ε 24 A equação 4 referida à deformação virtual é a seguinte a N u δ δ 3 25 que é equivalente a aT NT u δ δ 3 26 Substituindo todas estas equações em 21 passa a terse o PTV expresso por 2 2 1 2 2 1 2 3 L L T T L L S T T p d x N a a B E x B a d S d x δ δ 27 Passando para fora de cada integral tudo o que não depende da respectiva variável chegase a Viga de EulerBernoulli Álvaro F M Azevedo 209 2 2 1 2 2 1 2 3 L L T T L L S T T d x N p a x d S d x a B B E a δ δ 28 Nesta expressão considerase que o módulo de Young E é constante dentro da secção transversal e variável ao longo do eixo da barra O momento de inércia em relação ao eixo x2 é definido da seguinte forma sendo designado por I2 S x d S I 2 3 2 29 De acordo com o PTV a equação 28 é verdadeira para qualquer conjunto de deslocamentos virtuais concluindose assim que 2 2 1 2 2 1 2 L L T L L T d x N p B B E I d x a 30 A matriz de rigidez do elemento de viga é 2 2 1 2 L L T B B E I d x K 31 e o vector de solicitação é 2 2 1 L L T d x N p F 32 Supondo o módulo de Young e o momento de inércia constantes em toda a barra a expressão da matriz de rigidez passa a 2 2 1 2 L L T B B d x E I K 33 Substituindo em 33 a expressão 17 temse Viga de EulerBernoulli Álvaro F M Azevedo 210 2 2 1 1 2 1 3 1 2 1 3 1 2 1 3 1 2 1 3 2 6 1 12 6 1 12 6 1 12 6 1 12 L L d x L x L L x L x L x L L x L x L L x L x L E I K 34 Depois de efectuar o cálculo dos integrais presentes em 34 chegase a L L L L L L L L L L L L L L L L E I K 4 6 2 6 6 12 6 12 2 6 4 6 6 12 6 12 2 2 2 3 2 3 2 2 2 3 2 3 2 35 Esta matriz coincide com a que se obtém pelos métodos clássicos da teoria das estruturas reticuladas 121 Substituindo 5 em 32 e considerando as funções de forma 69 obtémse 2 2 1 3 1 2 2 1 1 3 1 3 1 3 1 2 2 1 1 3 1 3 1 1 2 1 4 1 8 2 2 3 2 1 1 2 1 4 1 8 2 2 3 2 1 L L x d L x L x x L L x x L L x L x x L L x x L p F 36 Depois de efectuar o cálculo dos integrais presentes em 36 chegase a 12 2 12 2 2 2 L L L L p F 37 que também corresponde ao que se obtém por métodos clássicos 121 Viga de EulerBernoulli Álvaro F M Azevedo 211 No Capítulo 11 encontrase deduzida a seguinte expressão para o cálculo do momento flector na viga quando o módulo de Young é constante 2 1 3 2 2 d x d u E I M 38 Atendendo a 11 15 e 38 concluise que o momento flector pode ser obtido com E I B a M 2 39 A matriz B é avaliada no ponto em que se pretende calcular o momento flector Devese notar que em geral esta expressão não fornece valores para os momentos flectores coincidentes com os da teoria clássica porque quando os deslocamentos a são nulos o momento flector calculado com 39 é nulo em toda a barra sendo assim ignorada a contribuição das cargas que actuam no seu interior eg carga distribuída p Esta questão obriga a que seja efectuada uma discretização de cada barra de um pórtico em vários elementos finitos ver a Figura 123 Fig 123 Exemplo discretização das barras de um pórtico em 22 elementos finitos O procedimento aqui apresentado para o cálculo da matriz de rigidez e do vector solicitação apresenta a vantagem de ser mais facilmente estendido a outras situações mais elaboradas eg elementos finitos com mais do que dois nós barras de secção variável barras não rectilíneas Viga de EulerBernoulli Álvaro F M Azevedo 212 123 Viga de três nós sem substituição de variável A formulação da matriz de rigidez e vector solicitação da viga de três nós é efectuada de um modo semelhante ao que foi exposto na Secção 122 As únicas diferenças são o aumento da dimensão de todos os vectores e matrizes envolvidos e o recurso às expressões da interpolação Hermitiana com três nós ver o Capítulo 7 124 Viga de dois nós com substituição de variável Quando é utilizada a interpolação Hermitiana e se faz uma substituição de variável surgem algumas questões que são apresentadas com base no exemplo da Figura 124 u3 x1 x1 a1 1 1 2 2 11 1 L x x 2 21 1 L x x L 2 L 2 a3 2 a2 θ1 a4 θ2 x3 u3 s s 1 1 a 1 2 1 1 s s 1 s s2 1 1 2 θ1 a 2 3 a 4 θ2 a Fig 124 Substituição de variável num elemento de viga com dois nós A transformação entre a coordenada x1 e a coordenada s é neste caso simples efectuada com a seguinte expressão Viga de EulerBernoulli Álvaro F M Azevedo 213 L s x 2 1 40 O vector dos deslocamentos generalizados na viga real é 2 2 1 1 4 3 2 1 θ θ a a a a a 41 sendo 1 3 d x θ d u 42 Após a substituição de variável e atendendo à coordenada local s temse 2 2 1 1 4 3 2 1 θ θ a a a a a 43 sendo d s d u 3 θ 44 A derivada em ordem a s da função u3x1s é pela regra da cadeia d s d x d x d u d s d u 1 1 3 3 45 De acordo com 40 temse neste caso 2 1 L d s d x 46 Designando por J a seguinte derivada Viga de EulerBernoulli Álvaro F M Azevedo 214 2 1 L d s d x J 47 temse d x J d u d s u d 1 3 3 48 e atendendo a 42 e 44 θ θ J 49 Nos nós temse J J 2 2 1 1 θ θ θ θ 50 Substituindo 50 em 43 obtémse J J a a a a a 2 2 1 1 4 3 2 1 θ θ 51 Atendendo a 41 chegase a a J a a J a a a a a a 4 3 2 1 4 3 2 1 52 A interpolação do campo de deslocamentos pode ser efectuada com base na coordenada s sendo utilizada a seguinte expressão ver a Figura 124 4 4 3 3 2 2 1 1 3 N s a N s a N s a N s a u s 53 Viga de EulerBernoulli Álvaro F M Azevedo 215 As funções de forma i N são definidas com as seguintes expressões que correspondem à interpolação Hermitiana numa viga com comprimento L 2 ver o Capítulo 7 3 1 4 1 4 3 2 1 s s N s 54 3 2 2 4 1 4 1 4 1 4 1 s s s s N 55 3 3 4 1 4 3 2 1 s s s N 56 3 2 4 4 1 4 1 4 1 4 1 s s s s N 57 Substituindo 52 em 53 chegase a 4 4 3 3 2 2 1 1 3 N s J a N s a N s J a N s a u s 58 Uma vez que se pretende que a interpolação de u3 seja efectuada da seguinte forma 4 3 2 1 4 3 2 1 3 a a a a N s N s N s N s u s 59 ou N a u s 3 60 concluise que N J N N J N N N N N N 4 3 2 1 4 3 2 1 61 sendo 3 1 1 4 1 4 3 2 1 s s s N s N 62 Viga de EulerBernoulli Álvaro F M Azevedo 216 J s s s s J N s N 3 2 2 2 4 1 4 1 4 1 4 1 63 3 3 3 4 1 4 3 2 1 s s s N s N 64 J s s s s J N s N 3 2 4 4 4 1 4 1 4 1 4 1 65 Atendendo às equações 1014 existe a necessidade de calcular a seguinte matriz 2 1 4 2 2 1 3 2 2 1 2 2 2 1 1 2 d x d N d x d N d x d N d x d N B 66 Para calcular as derivadas de Ni em ordem a x1 quando apenas se conhecem as funções Nis 6265 devese recorrer à regra da cadeia d s d x d x d N d s d N i i 1 1 67 Atendendo a 47 fica J d x d N d s d N i i 1 68 Derivando outra vez em ordem a s e considerando de novo a regra da cadeia temse d s d x d s d N d x d d s d N d s d i i 1 1 69 Considerando 47 e 68 chegase a J J d x d N d x d d s N d i i 1 1 2 2 70 Uma vez que J é constante temse Viga de EulerBernoulli Álvaro F M Azevedo 217 2 2 1 2 2 2 J d x d N d s N d i i 71 que é equivalente a 2 2 2 2 1 2 1 J d s d N d x N d i i 72 Derivando duas vezes as funções de forma 6265 em ordem a s temse s d s d N 2 3 2 2 1 73 J s d s d N 2 3 2 1 2 2 2 74 s d s d N 2 3 2 3 2 75 J s d s d N 2 3 2 1 2 4 2 76 Estas expressões são substituídas em 72 obtendose assim as segundas derivadas de Ni em ordem a x1 2 2 1 1 2 2 3 s J d x d N 77 J s d x d N 2 3 2 1 2 1 2 2 78 2 2 1 3 2 2 3 s J d x d N 79 J s d x d N 2 3 2 1 2 1 4 2 80 Atendendo a 47 temse Viga de EulerBernoulli Álvaro F M Azevedo 218 L s d x N d 2 2 1 2 1 6 81 L s L d x d N 3 1 2 1 2 2 82 L s d x N d 2 2 1 2 3 6 83 L s L d x d N 3 1 2 1 4 2 84 De acordo com 66 são os seguintes os elementos da matriz B em função da variável s L s L L s L s L L s B 3 1 6 3 1 6 2 2 85 De acordo com 33 ie supondo o módulo de Young e a secção constantes temse a seguinte expressão para a matriz de rigidez do elemento finito de viga no referencial local 2 2 1 2 L L T B B d x E I K 86 Após a substituição da variável x1 pela variável s 86 passa a 1 1 1 2 d s d s B B d x E I K T 87 Atendendo a 47 e a 85 temse 1 1 2 2 2 2 2 2 3 1 6 3 1 6 3 1 6 3 1 6 L d s L s L L s L s L s L L s L s L L s L s L E I K 88 Viga de EulerBernoulli Álvaro F M Azevedo 219 Notese que todos os elementos da matriz que constitui a função integranda são funções de s O comprimento da barra L é um parâmetro fixo Após o cálculo dos integrais obtémse a seguinte matriz L L L L L L L L L L L L L L L L E I K 4 6 2 6 6 12 6 12 2 6 4 6 6 12 6 12 2 2 2 3 2 3 2 2 2 3 2 3 2 89 Considerando a carga uniformemente distribuída representada na Figura 122 temse o seguinte vector solicitação que é calculado com a expressão 32 2 2 1 L L T d x N p F 90 Após a substituição da variável x1 pela variável s 90 passa a 1 1 1 d s d s d x N p F T 91 Atendendo a 47 e às funções de forma 6265 temse 1 1 3 2 3 3 2 3 2 2 4 1 4 1 4 1 4 1 4 1 4 3 2 1 2 4 1 4 1 4 1 4 1 4 1 4 3 2 1 d s L L s s s s s L s s s s s p F 92 Do cálculo destes integrais resulta Viga de EulerBernoulli Álvaro F M Azevedo 220 12 2 12 2 2 2 L L L L p F 93 As expressões 89 e 93 coincidem com as que se obtêm recorrendo à teoria clássica da flexão de vigas 121 A formulação aqui apresentada possui contudo a vantagem de ser extensível a casos mais genéricos tais como vigas curvas e vigas de secção variável em que pode ser vantajosa a utilização de funções de interpolação de grau mais elevado e consequentemente o recurso a elementos finitos com mais do que dois nós O cálculo do momento flector num ponto definido pela coordenada s é efectuado com a expressão 39 sendo a matriz B calculada com 85 Nas condições do elemento atrás descrito é possível demonstrar que os valores mais correctos do campo de momentos flectores se encontram nos pontos cuja coordenada s é 3 s 1 94 Se se pretender conhecer os valores do campo de momentos noutros pontos é em geral mais vantajoso efectuar uma extrapolação ou interpolação simples a partir dos pontos 94 O campo de esforços transversos pode ser obtido por derivação do campo de momentos em ordem a x 125 Considerações finais A formulação da viga de EulerBernoulli aqui apresentada não é mais desenvolvida porque na prática é preferível utilizar uma formulação que entre em linha de conta com a deformação por esforço transverso Esta formulação é apresentada no Capítulo 13 Viga de EulerBernoulli Álvaro F M Azevedo 221 BIBLIOGRAFIA 121 Correia de Araújo F Cálculo Matricial das Estruturas Contínuas pelo Método dos Deslocamentos Revista Engenharia Publicação dos Alunos da FEUP Ano XIX Número 43 NovembroDezembro 196566 122 Hinton E Owen D R J An Introduction to Finite Element Computations Pineridge Press Swansea UK 1979 123 Hughes T J R The Finite Element Method Linear Static and Dynamic Finite Element Analysis PrenticeHall Inc 1987 Viga de EulerBernoulli Álvaro F M Azevedo 222 223 CAPÍTULO 13 VIGA DE TIMOSHENKO Na formulação do elemento de viga de Timoshenko 131 é considerado que as secções planas se mantêm planas Contudo supõese que uma secção normal ao eixo da viga não mantém essa característica após a deformação Deste modo é possível considerar a deformação devida ao corte 131 Simbologia Apresentase em primeiro lugar um resumo da simbologia adoptada na formulação do elemento de viga de Timoshenko Tabela 131 Simbologia relativa ao elemento de viga de Timoshenko L Comprimento da barra prismática x Coordenada cartesiana u Campo de deslocamentos a Deslocamento generalizado nodal Deslocamento nodal θ Rotação nodal x Coordenada cartesiana de um nó de um elemento finito s Coordenada local s Coordenada local de um nó de um elemento finito J Jacobiano da transformação J d x1 d s N Função interpoladora ou função de forma G Centro de gravidade Viga de Timoshenko Álvaro F M Azevedo 224 φ Rotação correspondente à deformação por esforço transverso ε Extensão Bb Matriz de deformação relativa ao termo de flexão bending σ Tensão normal E Módulo de elasticidade ou módulo de Young γ Distorção Bs Matriz de deformação relativa ao termo de corte shear τ Tensão tangencial G Módulo de distorção W Trabalho V Volume S Superfície I Momento de inércia da secção transversal da barra prismática A Área da secção transversal da barra prismática A Área efectiva de corte relativa à secção transversal da barra prismática α Coeficiente de redução da área da secção transversal para atender ao corte F Forças nodais equivalentes à acção exterior nos graus de liberdade do elemento finito no referencial local K Matriz de rigidez do elemento finito no referencial local M Momento flector V Esforço transverso 132 Viga de dois nós com substituição de variável Na Figura 131 encontrase representado um elemento de viga com dois nós e com comprimento L ver o Capítulo 11 Supõese que nos nós não há deslocamentos segundo x1 Deste modo apenas se considera o comportamento à flexão da viga Viga de Timoshenko Álvaro F M Azevedo 225 u3 x1 x1 a1 1 1 2 2 11 1 L x x 2 21 1 L x x L 2 L 2 a3 2 a2 θ1 a4 θ2 x3 u3 s s 1 2 1 1 s s 1 s s2 1 1 θ s Fig 131 Substituição de variável num elemento de viga com dois nós Os deslocamentos generalizados dos nós do elemento finito representado na Figura 131 são os seguintes 2 2 1 1 4 3 2 1 θ θ a a a a a 1 A transformação entre a coordenada x1 e a coordenada s é neste caso simples efectuada com a seguinte expressão L s x 2 1 2 sendo a derivada em ordem a s a seguinte 2 1 L d s d x J 3 Viga de Timoshenko Álvaro F M Azevedo 226 A interpolação do deslocamento lateral u3 e da rotação θ é efectuado separadamente para cada uma destas variáveis Assim e uma vez que u3 e θ apresentam dois valores nodais cada é utilizada a seguinte interpolação unidimensional com dois nós 3 2 1 1 3 s a N s a N s u 4 4 2 2 1 s a N s a N s θ 5 Neste exemplo com dois nós as funções de forma são as seguintes ver o Capítulo 4 2 1 1 s s N 6 2 1 2 s s N 7 Na Figura 132 está representado o eixo da viga na sua posição inicial sobre x1 e a correspondente deformada ver o Capítulo 11 Está também representada a secção transversal cujos eixos são x2 e x3 Uma vez que se consideram pequenas deformações supõese que o declive da recta tangente ao eixo coincide com o ângulo de rotação do eixo da barra x1 x3 u3 x1 u3 u1 1 3 d x d u G x2 x3 Fig 132 Barra deformada e secção transversal Na Figura 133 estão indicados os seguintes ângulos rotação do eixo da barra 1 3 d x d u rotação da secção transversal θ e rotação correspondente à deformação por esforço transverso φ Encontrase também representado o campo de deslocamentos u1 na secção transversal Viga de Timoshenko Álvaro F M Azevedo 227 φ x3 u3 u1 x1 u3 x1 A B A B θ x3 θ u1 x1 x3 A A B B O O O O 1 3 d x u d Fig 133 Barra deformada e deslocamentos da secção transversal Na formulação da viga de EulerBernoulli considerase que o ângulo φ é nulo sendo os ângulos 1 3 d x d u e θ coincidentes Na formulação da viga de Timoshenko o ângulo φ é considerado não nulo sendo φ θ 1 3 d x d u 8 Estes três ângulos dependem de x1 De acordo com a Figura 133 temse 1 3 3 1 1 x x u x x θ 9 A extensão ε1 é definida por 132 θ ε 3 1 1 1 1 x x x u 10 sendo Viga de Timoshenko Álvaro F M Azevedo 228 1 3 1 d x x d θ ε 11 Designando por 1ε a seguinte componente da expressão 11 1 1 d x dθ ε 12 passa a terse 1 3 1 ε ε x 13 Derivando 5 em ordem a x1 chegase a 4 1 2 2 1 1 1 a d x d N d x a d N d x d θ 14 Substituindo 14 em 12 obtémse 4 3 2 1 1 2 1 1 1 0 0 a a a a d x d N d x d N ε 15 Considerese agora uma matriz de deformação que é designada B b pelo facto de estar associada à flexão bending A sua definição é a seguinte 1 2 1 1 0 0 d x d N d x d N B b 16 Atendendo a 1 15 passa a escreverse B b a 1ε 17 Substituindo 17 em 13 obtémse a B x b 3 ε1 18 Viga de Timoshenko Álvaro F M Azevedo 229 Tal como no Capítulo 12 considerase a lei de Hooke referida apenas à tensão normal σ1 e à extensão ε1 1 1 ε σ E 19 Substituindo 18 em 19 temse a E x B b 3 σ1 20 A distorção γ13 é definida por 132 1 3 3 1 13 d x d u d x d u γ 21 Atendendo a 9 temse a seguinte distorção média 133 1 3 3 3 13 d x d u x d x d θ γ 22 1 3 13 d x d u θ γ 23 Substituindo 8 em 23 obtémse φ γ 13 24 Derivando 4 em ordem a x1 chegase a 3 1 2 1 1 1 1 3 a d x d N d x a d N d x d u 25 Substituindo 5 e 25 em 23 obtémse 3 1 2 1 1 1 4 2 2 1 13 a d x d N d x a d N N a N a γ 26 Em notação matricial temse Viga de Timoshenko Álvaro F M Azevedo 230 4 3 2 1 2 1 2 1 1 1 13 a a a a N d x d N N d x d N γ 27 Considerese agora uma matriz de deformação que é designada B s pelo facto de estar associada ao corte shear A sua definição é a seguinte 2 1 2 1 1 1 N d x d N N d x d N B s 28 Atendendo a 1 27 passa a escreverse B s a γ13 29 Uma vez que de acordo com a lei de Hooke para materiais isotrópicos 132 13 13 γ τ G 30 temse depois de substituir 29 em 30 G B a s τ13 31 De acordo com o Princípio dos Trabalhos Virtuais PTV ver o Capítulo 4 admitese que Trabalho Interno Trabalho Externo 32 Considerando δ Wbi Trabalho interno associado à flexão bending 33 δ Wsi Trabalho interno associado ao corte shear 34 δ W e Trabalho externo 35 De acordo com 32 temse Viga de Timoshenko Álvaro F M Azevedo 231 e i s bi W W W δ δ δ 36 Considerando que a contribuição da flexão para o trabalho apenas depende da tensão normal σ1 temse ver o Capítulo 4 V bi dV W δ ε1 σ1 δ 37 A equação 18 referida à deformação virtual é a seguinte a B x b δ δ ε 3 1 38 sendo equivalente a 3 1 x B a T b T δ δ ε 39 Substituindo 39 e 20 em 37 obtémse 2 2 1 3 3 L L S b T b T bi a d S d x x E x B B a W δ δ 40 Nesta equação S é a superfície correspondente à secção transversal da barra ver o Capítulo 11 De acordo com a Figura 132 temse 3 2 d x d x d S 41 Supondo o módulo de Young E constante em todos os pontos do elemento de viga e passando para fora dos integrais tudo o que não depende da respectiva variável de integração resulta x d S d x a B B E a W L L S b T b T bi 2 2 1 2 3 δ δ 42 O momento de inércia em relação ao eixo x2 é definido da seguinte forma sendo designado por I2 ver a Figura 132 S x d S I 2 3 2 43 Viga de Timoshenko Álvaro F M Azevedo 232 Substituindo 43 em 42 e supondo que a barra é de secção constante passa a terse d x a B B E I a W L L b T b T bi 2 2 1 2 δ δ 44 Considerando que a contribuição do corte para o trabalho apenas depende da tensão tangencial τ13 temse V si dV W δ γ13 τ13 δ 45 A equação 29 referida à deformação virtual é a seguinte B s δ a δ γ 13 46 sendo equivalente a T δ aT Bs δ γ 13 47 Substituindo 47 e 31 em 45 obtémse 2 2 1 L L S s T s T si G B a d S d x B a W δ δ 48 Supondo o módulo de distorção G constante em todos os pontos do elemento de viga e passando para fora dos integrais tudo o que não depende da respectiva variável de integração resulta d S d x a B B G a W L L S s T s T si 2 2 1 δ δ 49 A área da secção transversal da barra é S d S A 50 Na expressão 49 é necessário introduzir o factor correctivo de corte α sendo a área reduzida de corte definida por 133 Viga de Timoshenko Álvaro F M Azevedo 233 A A α 51 Substituindo 50 em 49 e considerando a área reduzida de corte passa a terse no caso de uma barra de secção constante B d x a B G A a W L L s T s T si 2 2 1 δ δ 52 Por uma questão de simplificação desta exposição considerase que o trabalho externo associado às forças exteriores δ W e inclui apenas a contribuição das forças generalizadas concentradas nos nós da barra Nestas condições temse F a W T e δ δ 53 As componentes do vector F são forças generalizadas forças e momentos em correspondência com os quatro graus de liberdade dos nós da barra ver a Figura 131 Substituindo 44 52 e 53 em 36 obtémse F a d x a B B G A a d x a B B E I a T L L s T s T L L b T b T δ δ δ 2 2 1 2 2 1 2 54 Uma vez que 54 tem de se verificar para qualquer deformação virtual δ a chegase à habitual equação F K a 55 sendo a matriz de rigidez K calculada com a seguinte expressão 2 2 1 2 2 1 2 L L s T s L L b T b d x B B G A d x B B E I K 56 Depois de efectuar em 56 a substituição de variável definida em 2 temse 1 1 1 1 1 1 2 d s d s d x B B G A d s d s d x B B E I K s T s b T b 57 Viga de Timoshenko Álvaro F M Azevedo 234 Substituindo 3 em 57 chegase a 1 1 1 1 2 2 2 d s B B G A L d s B B L E I K s T s b T b 58 Para se obter os elementos das matrizes B b 16 e B s 28 em função da variável s é necessário calcular a derivada das funções de forma em ordem a x1 Para isso é suficiente recorrer à regra da cadeia ficando d s d x d x d N d s d N i i 1 1 59 Atendendo a 3 passa a terse 2 1 L d x d N d s d N i i 60 que é equivalente a d s d N L d x d N i i 2 1 61 Da derivação de 6 e 7 em ordem a s resulta 2 1 1 d s d N 62 2 2 1 d s d N 63 Atendendo a 61 temse L d x d N 1 1 1 64 L d x d N 1 1 2 65 Substituindo 64 e 65 em 16 obtémse Viga de Timoshenko Álvaro F M Azevedo 235 L L B b 1 0 1 0 66 Substituindo 6 7 64 e 65 em 28 obtémse s L s L B s 2 1 1 1 2 1 1 1 67 Os elementos das matrizes Bb e Bs apenas dependem de L que se considera um parâmetro fixo e da variável s Depois de substituir estas expressões em 58 e de calcular os integrais em ordem a s resulta a seguinte expressão para a matriz de rigidez do elemento finito 3 2 1 6 2 3 2 1 2 1 1 0 0 1 0 1 0 0 0 0 2 2 2 2 L SIM L L L L L L L A G SIM L E I K 68 O facto de a expressão 68 ser aproximada obriga a que na análise de um pórtico cada uma das suas barras tenha de ser discretizada em vários elementos finitos Esta questão foi já referida no Capítulo 12 No Capítulo 11 encontrase deduzida a seguinte expressão para o cálculo do momento flector na viga quando o módulo de Young é constante 1 2 d x d E I M θ 69 Atendendo a 12 17 e 69 concluise que o momento flector pode ser obtido com E I B a M b 2 70 A matriz Bb é avaliada no ponto em que se pretende calcular o momento flector As expressões 29 e 31 referemse à distorção média e à tensão tangencial média O esforço transverso V é calculado com a seguinte expressão Viga de Timoshenko Álvaro F M Azevedo 236 S d S x V 3 τ13 71 Nesta expressão τ13x3 representa a tensão tangencial real que depende da coordenada x3 ver a Figura 132 Uma vez que na presente formulação apenas se dispõe da tensão tangencial média τ13 é necessário calcular o esforço transverso V com base na área efectiva de corte A A sua expressão é a seguinte 134 1 3 d x d u G A V θ 72 Atendendo a 23 temse V G A γ13 73 Substituindo 24 em 73 obtémse G A φ V 74 O ângulo φ está representado na Figura 133 Substituindo 29 em 73 chegase a G A B a V s 75 A matriz Bs é avaliada no ponto em que se pretende calcular o esforço transverso A expressão que fornece a área efectiva de corte A depende da forma de secção transversal 133 Devese ter em consideração que quer o momento flector quer o esforço transverso apenas apresentam valores com precisão aceitável em determinados pontos do elemento finito 135 Se se pretender conhecer os valores dos esforços noutros pontos é em geral preferível efectuar uma extrapolação ou interpolação simples a partir dos pontos em que os resultados são mais correctos Viga de Timoshenko Álvaro F M Azevedo 237 133 Considerações finais A formulação da viga de Timoshenko aqui apresentada pode ser estendida aos seguintes casos barras com mais do que dois nós barras curvilíneas barras de secção variável barras tridimensionais sujeitas a flexão desviada inclusão da torção consideração do centro de corte distinto do centro de gravidade barras em que as propriedades do material variam ao longo do eixo da barra ou dentro da secção transversal etc 134 BIBLIOGRAFIA 131 Oñate E Cálculo de Estructuras por el Método de los Elementos Finitos Análisis Estático Lineal Segunda Edición CIMNE Barcelona 1995 132 Azevedo A F M Mecânica dos Sólidos Faculdade de Engenharia da Universidade do Porto 1996 133 Massonnet C Résistance des Matériaux Dunod Paris 1968 134 Barros J A O Método dos Elementos Finitos Aplicado a Estruturas Reticuladas Relatório 01DECE99 Universidade do Minho 2001 135 Cook R D Malkus D S Plesha M E Witt R J Concepts and Applications of Finite Element Analysis Fourth Edition John Wiley Sons Inc 2002 Viga de Timoshenko Álvaro F M Azevedo 238 239 ANEXO A UTILIZAÇÃO DO PROGRAMA FEMIX 31 Neste capítulo é efectuada uma breve descrição das tarefas que é necessário empreender para analisar uma estrutura com o programa FEMIX Versão 31 A documentação completa bem como as instruções para o download do programa FEMIX 31 encontramse no seguinte URL httpcivilfeupptSoftwareFemix31Femix31Manualhtm Referese em primeiro lugar o modo de instalação seguindose um exemplo de aplicação Nesta publicação não se pretende repetir o conteúdo do manual do programa A1 devendo o leitor recorrer à documentação completa sempre que surgirem dúvidas A1 Instalação Descarregar do URL acima referido o seguinte ficheiro femixV310031zip Fazer a extracção de todo o conteúdo deste ficheiro para um directório qualquer Sugerese a instalação em C sendo automaticamente criado neste local um directório chamado Cfemix Dentro deste directório surgem outros subdirectórios É aconselhável acrescentar ao PATH o directório Cfemixbin Em Windows 2000 ou XP esta operação pode ser efectuada clicando em My Computer com o botão da direita e seleccionando Properties Advanced Environment Variables Seleccionar em seguida o PATH do utilizador corrente a carregar no botão Edit Em seguida devese acrescentar no fim da lista de directórios o seguinte texto Cfemixbin Aconselhase também a criação no Desktop de Shortcuts para os seguintes programas Utilização do Programa Femix 31 Álvaro F M Azevedo 240 Cfemixbins3dcadexe Cfemixbindrawmeshexe Habitualmente a invocação dos diversos módulos é feita a partir da linha de comandos Para obter uma janela que suporta a invocação de comandos devese seleccionar Start Run e em seguida escrever na janela de texto cmd Para aumentar o número de linhas de texto devese clicar no canto superior esquerdo da janela de comandos e seleccionar Properties Layout Em seguida aumentar o parâmetro Window Size Height Para testar a instalação devese fazer o seguinte Abrir uma janela de comandos cmd cd emp md femix cd femix s3dcad Se a instalação tiver sido feita correctamente deve ser possível arrancar o programa s3dcad a partir do directório corrente Tendo em vista uma familiarização com os diversos ficheiros que fazem parte da instalação aconselhase uma inspecção ao conteúdo de todos os directórios que se encontram dentro de Cfemix A2 Preparação dos dados Apresentase em seguida uma descrição dos principais passos a dar para se chegar aos resultados de uma análise com o programa FEMIX 31 Todas as fases são exemplificadas com base na estrutura representada na Figura A1 Utilização do Programa Femix 31 Álvaro F M Azevedo 241 x1 x2 50 kN m 100 m 10 m E 20 GPa ν 015 h 03 m Fig A1 Viga que se pretende analisar pelo MEF usando o programa FEMIX Para criar um ficheiro de dados contendo a quase totalidade da informação que descreve o problema devese escrever o seguinte a partir da linha de comandos A2 Nota o caracter e todos os que se encontram à sua direita são comentários não sendo necessário digitálos s3dcad csm Create a simple mesh 2 Rectangle 10 Size in x 1 1 Size in x 2 gen Generate a refined mesh 2 Surfaces 4 N of nodes of the generated elements 4 N of divisions for all the elements in s1 1 N of divisions for all the elements in s2 ren Renumber elements nodes and special nodes 1 Default answer 2 Default answer 3 Default answer y Default answer y Default answer y Default answer 10e5 Default answer wri Write a s3d file viga44 Job name elementos de 4 nós malha com 4 elementos gld Write a gldat file femix viga44 Job name Utilização do Programa Femix 31 Álvaro F M Azevedo 242 1 Plane stress 1 From the coordinates x1x2 end End s3dcad Da execução destes comandos resultam os seguintes ficheiros viga44s3d Ficheiro com a geometria tendo em vista a sua visualização gráfica viga44gldat Ficheiro com os dados para a análise com o FEMIX Uma vez que se tratam de ficheiros do tipo texto o seu conteúdo pode ser inspeccionado por exemplo com os programas Notepad ou Word O ficheiro viga44s3d destinase ao programa drawmesh ver a Figura A2 Os principais comandos deste programa são os seguintes File Import Importar o ficheiro de extensão s3d View Set View Angles XY Visualizar o plano XY Options Markers Colocar tudo Visible Options Numbers Colocar tudo Visible Options Lines Alterar o Shrink factor para 90 View Shading Fazer a coloração dos elementos Nota para muitos dos comandos existem botões nas barras de ferramentas bem como teclas de atalho fazer Help Keyboard Commands Fig A2 Visualização da malha com o programa drawmesh O ficheiro viga44gldat que foi gerado com o programa s3dcad ainda não se encontra completo Referemse em seguida as alterações que devem ser efectuadas Utilização do Programa Femix 31 Álvaro F M Azevedo 243 Substituir o título Rectangular mesh por Viga discretizada com 4 elementos de 4 nos kNm Substituir o bloco de parâmetros pelo seguinte 4 nelem n of elements in the mesh 10 npoin n of points in the mesh 2 nvfix n of points with fixed degrees of freedom 1 ncase n of load cases 1 nmats n of sets of material properties 1 nspen n of sets of element nodal properties 1 ntype problem type 4 nnode n of nodes per element 2 ngaus n of Gauss points in the integration rule element stiffness 2 ngstr n of Gauss points in the integration rule stresses 2 ndime n of geometric dimensions 2 ndofn n of degrees of freedom per node 0 nnscs n of points with specified coordinate system 0 nsscs n of specified coordinate systems 0 npspr n of springs 0 nsspv n of spring vectors 4 nprop n of material properties used in the formulation 1 npren n of element nodal properties used in the formulation 0 nwink n of element faces with Winkler coefficients Acrescentar as definições das características dos apoios ao seguinte bloco de dados Points with fixed degrees of freedom and fixity codes 1fixed0free ivfix nofix ifpre 1 1 1 1 2 9 0 1 Remover os seguintes blocos de dados Points with specified coordinate system Specified coordinate system index Spring index point number type of spring vector spring constant value and Spring vector index Modificar as propriedades do material para o seguinte Utilização do Programa Femix 31 Álvaro F M Azevedo 244 Sets of material properties Young modulus Poisson ratio mass per unit volume and thermic coeff imats young poiss dense alpha 1 200e6 015 00 00 kPa Modificar o bloco das espessuras nodais para o seguinte Sets of element nodal properties ispen 1 inode thickness 1 03 2 03 3 03 4 03 Substituir os casos de carga que surgem por defeito pelas seguintes linhas Title of the first load case Carga distribuida de 50 kNm Load parameters 0 nplod n of point loads in nodal points 0 ngrav gravity load flag 1yes0no 4 nedge n of edge loads FEM only 0 nface n of face loads FEM only 0 nteme n of elements with temperature variation FEM only 0 nudis n of uniformly distributed loads 3d frames and trusses only 0 ntral n of trapezoidal distributed loads 3d frames and trusses only 0 nepoi n of bar point loads 3d frames and trusses only 0 ntemb n of bars with temper variation 3d frames and trusses only 0 nprva n of prescribed and non zero degrees of freedom Edge load loaded element loaded points and load value local coordinate system iedge loele 1 1 lopoe fe1 fe2 2 00 500 Utilização do Programa Femix 31 Álvaro F M Azevedo 245 4 00 500 iedge loele 2 2 lopoe fe1 fe2 4 00 500 6 00 500 iedge loele 3 3 lopoe fe1 fe2 6 00 500 8 00 500 iedge loele 4 4 lopoe fe1 fe2 8 00 500 10 00 500 ENDOFFILE A3 Execução do programa Depois de ter o ficheiro viga44gldat completamente definido escrever na linha de comandos o seguinte prefemix viga44 Verificar a coerência dos dados femix viga44 d Calcular a matriz de rigidez global calcular o vector solicitação global e resolver o sistema de equações lineares posfemix viga44 Gravar diversos tipos de ficheiros de resultados Depois de executar as diversas opções do programa posfemix podese inspeccionar os ficheiros que foram criados dos quais se destacam os seguintes viga44gllpt dados formatados viga44rslpt resultados formatados viga44mes3d malha indeformada Utilização do Programa Femix 31 Álvaro F M Azevedo 246 viga44dms3d malha deformada viga44pss3d tensões principais viga44dipva campo de deslocamentos viga44ds3d malha indeformada desconectada viga44dstpva campo de tensões relativo à malha desconectada A4 Visualização gráfica Para visualizar os ficheiros de extensão s3d devese fazer no drawmesh File Import Para visualizar os campos escalares contidos em ficheiros de extensão pva devese fazer no drawmesh PVA Import Este opção deve ser seleccionada depois de se ter lido a correspondente malha num ficheiro de extensão s3d Para capturar o conteúdo de uma janela do drawmesh podese fazer File Export View Image Desta forma é criado um ficheiro com extensão bmp que pode em seguida ser inserido num documento Word ou em qualquer outra aplicação Windows Para combinar a malha indeformada com a malha deformada devese escrever na linha de comandos s3djoin o viga44medm viga44me viga44dm Em seguida importar o ficheiro viga44medms3d com o drawmesh ver a Figura A3 Fig A3 Visualização da malha deformada com o programa drawmesh Utilização do Programa Femix 31 Álvaro F M Azevedo 247 Para combinar as tensões principais com a malha indeformada devese escrever na linha de comandos s3djoin o viga44meps viga44me viga44ps Em seguida importar o ficheiro viga44mepss3d com o drawmesh ver a Figura A4 Fig A4 Visualização das tensões principais com o programa drawmesh Na Figura A5 encontrase representado o campo escalar correspondente aos deslocamentos segundo x1 Fig A5 Visualização do campo de deslocamentos horizontais com o programa drawmesh Utilização do Programa Femix 31 Álvaro F M Azevedo 248 Na Figura A6 encontrase representado o campo escalar correspondente às tensões normais segundo x1 Fig A6 Visualização do campo de tensões normais σx1 com o programa drawmesh A5 Considerações finais Neste capítulo foi apresentado um exemplo muito simples de aplicação do programa FEMIX à análise de uma estrutura pelo MEF Para fazer aplicações a outros tipos de estruturas aconselhase a leitura da correspondente documentação A1 A2 BIBLIOGRAFIA A1 Azevedo A F M Barros J A O Manual de Utilização do Programa FEMIX Versão 31 Porto 2000 httpcivilfeupptSoftwareFemix31Femix31Manualhtm A2 Azevedo A F M Barros J A O Manual de Utilização do Programa S3DCAD Versão 30 Porto 1998 httpcivilfeupptSoftwareFemix31pdfS3dcadpdf