Métodos Numéricos para Engenharia: Introdução ao Método dos Elementos Finitos

CEFETPR CENTRO FEDERAL DE EDUCAÇÃO TECNOLÓGICA DO PARANÁ DAMEC DEPARTAMENTO ACADÊMICO DE MECÂNICA MÉTODOS NUMÉRICOS PARA A ENGENHARIA INTRODUÇÃO AO MÉTODO DOS ELEMENTOS FINITOS Fundamentos Teóricos MARCO ANTÔNIO LUERSEN Março 2000 PARANÁ 2 Introdução ao MEF M A Luersen ÍNDICE CAPÍTULO 1 Conceitos Iniciais e Método dos Resíduos Ponderados 3 11 Introdução 3 12 Conceitos Iniciais 5 13 Método dos Resíduos Ponderados 7 131 Método da Colocação 11 132 Subdomínio 12 133 Mínimos Quadrados 12 134 Colocação dos Mínimos Quadrados 13 135 Galerkin 13 CAPÍTULO 2 Elementos Finitos Unidimensionais 21 21 Método dos Elementos Finitos de Galerkin 21 22 Obtenção da Matriz de Rigidez e Vetor de Carga Globais 30 23 Imposição das Condições de Contorno 33 24 Resolução do Sistema Algébrico Linear e Determinação dos Deslocamentos Nodais 35 25 Cálculo de Deformações e Tensões Quantidades Derivadas 37 26 Cálculo das Reações 41 27 Condução de Calor Unidimensional em Regime Permanente 42 28 Formulação de Elementos Finitos Utilizando o Princípio da Mínima Energia Potencial 47 29 Método Direto 51 210 Obtenção de Funções de Interpolação Unidimensionais C0 54 211 Transformação de Coordenadas para o Elemento de Barra 60 REFERÊNCIAS BIBLIOGRÁFICAS 65 APÊNDICE 67 Cap 1 Conceitos Introdutórios 3 CAPÍTULO 1 Conceitos Iniciais e Método dos Resíduos Ponderados 11 INTRODUÇÃO Para a representação de fenômenos físicos os engenheiros e físicos geralmente estabelecem um sistema de equações diferenciais válidas em certa região domínio e impõem nesse sistema condições de contorno e condições iniciais Até esse estágio o modelo matemático está completo e para aplicações práticas necessitase somente a solução para um conjunto particular de dados numéricos Aqui surge a maior dificuldade pois apenas formas muito simples de equações dentro de contornos triviais são possíveis de serem resolvidas exatamente através de métodos matemáticos Equações diferenciais ordinárias com coeficientes constantes são alguns exemplos para os quais soluções analíticas são possíveis Para solucionar tais dificuldades há a necessidade da utilização de métodos numéricos obtendose assim soluções aproximadas da equação diferencial que se deseja resolver sujeita a certas condições de contorno e condições iniciais A idéia básica dos métodos numéricos para solucionar tal classe de problemas é discretizar o problema contínuo Ou seja o conjunto infinito de números que representam a função ou funções desconhecidas é substituído por um conjunto finito de parâmetros desconhecidos sendo que esse processo requer alguma forma de aproximação Os parâmetros desconhecidos são encontrados através de solução algébrica ou seja obtémse um sistema algébrico de equações do tipo b A x 11 4 Introdução ao MEF M A Luersen onde A é uma matriz quadrada e nãosingular b é o vetor independente e x é o vetor que contém os parâmetros desconhecidos da solução aproximada ou seja é o vetor solução Por vezes a matriz A pode depender do vetor x caracterizando um problema nãolinear Ou ainda A x e b podem ser todos dependentes do tempo desejandose a solução não somente para um instante mas para um determinado intervalo de tempo necessitando solucionar a equação 11 várias ou inúmeras vezes Dentre os métodos numéricos mais conhecidos e utilizados podese destacar o Método das Diferenças Finitas o Método dos Volumes Finitos o Método dos Elementos Finitos e o Método dos Elementos de Contorno Os três primeiros métodos são métodos de domínio isto é a discretização é ao longo de toda a região Já o Método dos Elementos de Contorno como o próprio nome diz é um método de contorno necessitando apenas fazer a discretização ao longo do contorno da região de estudo Outra classe de métodos numéricos que está se destacando são os chamados métodos sem malha meshless methods onde contrário ao Método dos Elementos Finitos não necessita da definição de uma malha mas somente de pontos ao longo do domínio Maiores detalhes sobre métodos sem malha podem ser encontrados em Duarte 1995 e Duarte e Oden 1995 Na engenharia outras formas de abordar a obtenção de soluções aproximadas são utilizadas tais como os princípios de equilíbrio conservação e métodos energéticos Mas a qualquer uma dessas formas de abordar o problema podese sempre também associar uma equação diferencial Assim de uma forma geral o Método dos Elementos Finitos MEF pode ser encarado como um método numérico aproximado para solucionar equações diferencias com precisão aceitável para engenheiros Seu desenvolvimento inicial foi em aplicações de análise estrutural em aeronaves utilizando o princípio de equilíbrio de forças e análise matricial Posteriormente foi aplicado a outros fenômenos sendo atualmente utilizado na solução de muitos problemas dentre os quais podese citar análise de tensões vibrações acústica escoamento de fluidos distribuição de temperatura campo eletromagnético injeção de plástico etc tanto para modelar fenômenos lineares como nãolineares Existem inúmeros softwares comerciais no mercado Cap 1 Conceitos Introdutórios 5 Ansys análise estrutural linear e nãolinear acústica eletromagnetismo escoamento de fluidos distribuição de temperatura ProMechanica análise estrutural linear e distribuição de temperatura Abaqus análise estrutural linear e nãolinear Nastran análise estrutural Cosmos análise estrutural Algor análise estrutural linear e nãolinear análise térmica escoamento de fluido LSDyna conformação mecânica de metais crash IDeas análise estrutural linear e nãolinear análise térmica escoamento de fluidos CMold simulação do processo de injeção de plástico Moldflow simulação do processo de injeção de plástico 12 CONCEITOS INICIAIS No estudo de soluções aproximadas de equações diferenciais é necessário introduzir alguns conceitos e definições Funcional expressão integral escalar que contém implicitamente as equações diferenciais que descrevem o problema Na mecânica estrutural o funcional mais largamente utilizado é a expressão da energia potencial total Funcionais também são avaliados para problemas de condução de calor escoamento de fluidos acústica etc De uma forma genérica um funcional no espaço ℜ2 pode ser escrito da seguinte forma I F x y u v u x u y v x v x dxdy 2 2 12 sendo x e y as variáveis independentes u e v funções de x e y O termo funcional indica que I depende não somente de u v e suas derivadas em um ponto mas também do seu efeito integrado sobre a região de interesse Na mecânica 6 Introdução ao MEF M A Luersen dos sólidos um exemplo de funcional é a expressão da energia potencial total π de um corpo elástico que escrita em notação indicial é dada por Γ Ω Ω σ ε π Γ Ω Ω t u d b u d d 2 1 i i i i ij ij 13 onde σ ij são as componentes de tensão ε ij as componentes de deformação bi as forças de corpo it as forças de superfície ui os deslocamentos Ω o domínio de interesse e Γ o contorno do domínio de interesse Cálculo Variacional parte da matemática que se preocupa em encontrar valores estacionários extremos de funcionais Ou seja para o funcional I da expressão 12 procurase as funções uxy e vxy de modo que I seja estacionário As funções u e v que satisfizerem esta condição são também a solução das equações de EulerLagrange associadas ao funcional As equações de EulerLagrange representam em formato de equações diferenciais o mesmo problema representado de forma integral pelo funcional a elas associado Noções sobre cálculo variacional podem ser encontrados em Dym Shames 1978 Tauchert 1974 Bassanezi Ferreira Jr 1988 entre outros trabalhos Princípio Variacional usa uma expressão integral chamada funcional que possui implicitamente as equações diferenciais e as condições de contorno não essenciais de um problema e nele são aplicados os princípios do cálculo variacional procurando uma posição estacionária geralmente um mínimo para o funcional Freqüentemente o funcional é uma forma quadrática O Princípio da Mínima Energia Potencial é um dentre os vários princípios utilizados na mecânica Assim no sentido clássico quando se fala em formulação variacional pressupõese a minimização de uma forma quadrática Utilizando um princípio variacional o Método dos Elementos Finitos pode ser apresentado como uma generalização do Método de RayleighRitz Detalhes sobre o Método de RayleighRitz podem ser encontrados no livro Energy Principles in Structural Mechanics Tauchert 1974 Em outras áreas que não seja a mecânica estrutural um princípio variacional pode não ser obtido Isto acontece se a equação diferencial do problema contém derivadas de ordem ímpar ocorrendo na mecânica dos fluidos para alguns tipos de escoamento Cap 1 Conceitos Introdutórios 7 Mesmo assim o Método dos Elementos Finitos pode ser aplicado utilizando o Método dos Resíduos Ponderados Tanto o Método de RayleighRitz como o Método dos Resíduos Ponderados usam expressões integrais que contém as equações diferenciais do problema físico Ambas as formulações funcional e residual são conhecidas como forma fraca das equações que governam o problema A equação diferencial é conhecia como forma forte A forma fraca força as condições de continuidade na média integral sense ao passo que a forma forte força as condições de continuidade em cada ponto Note que a solução da forma forte é sempre também solução da forma fraca ao passo que a afirmação contrária nem sempre é verdadeira A solução da forma fraca é conhecida também como solução generalizada 13 MÉTODO DOS RESÍDUOS PONDERADOS Para se aplicar o Método dos Resíduos Ponderados necessitase primeiramente saber a equação diferencial que rege o problema físico em questão Seja a equação diferencial que governa um determinado problema físico com condições especificadas no contorno problema de valor no contorno Du f 0 no domínio Ω Bu g 0 no contorno Γ condições de contorno 14 onde u variável dependente Por exemplo deslocamentos em um ponto temperatura em um ponto potencial elétrico x variável independente Por exemplo coordenadas de um ponto f g funções de x constantes ou zero dependendo do problema D B operadores diferenciais Se o operador diferencial D é de ordem 2m definese dois tipos de condições de contorno de acordo com a ordem do operador diferencial B Esta classificação é importante no estudo de soluções de equações diferenciais pois cada tipo de condição de contorno é abordada de uma forma 8 Introdução ao MEF M A Luersen Condições de contorno essenciais ou condições de contorno de Dirichlet condições de contorno que envolvem derivadas de ordem 0 a m1 Em análise estrutural as condições de contorno essenciais também são conhecidas como condições de contorno geométricas ou cinemáticas Condições de contorno naturais ou não essenciais ou condições de contorno de Neumann condições de contorno que envolvem derivadas de ordem m a 2m1 Uma equação que combine uma condição de contorno essencial e uma natural é chamada de condição de contorno mista Para um mesmo ponto apenas um tipo de condição de contorno pode ser especificada isto é se é especificada uma condição de contorno essencial a natural nesse mesmo ponto ou região não pode ser especificada pois depende da condição de contorno já fornecida e viceversa ao se especificar uma condição de contorno natural O que é possível é terse uma condição de contorno mista isto é para um mesmo ponto ou região uma equação que envolve os dois tipos de condições de contorno Por exemplo para uma equação diferencial de segunda ordem condições de contorno essenciais são as derivadas de ordem zero ou seja a própria função e condições de contorno naturais são as derivadas primeira da função A seguir são mostrados alguns exemplos de equações diferenciais associadas ao problema físico que representam Flexão de viga material elásticolinear e teoria de pequenas deformações EI d w dx p x 4 4 para este caso D EI d dx 4 4 u w deflexão f p carregamento transversal e as condições de contorno naturais Bu g 0 EI d w dx M B 2 2 0 momento fletor EI d w dx VB 3 3 0 esforço cortante Cap 1 Conceitos Introdutórios 9 Como condições de contorno essenciais para este problema temse a deflexão w e sua derivada dw dx esta última representando a rotação Condução de calor unidimensional em regime permanente d dx k x A x dT x dx Q x A x 0 onde T temperatura k condutividade térmica A área Q fonte de calor k dT dx fluxo de calor Como condição de contorno natural temse a especificação do fluxo q k dT dx condição de contorno essencial a temperatura T e para este problema há a ocorrência comum de uma condição de contorno mista representada pela troca de calor por convecção k dT dx h T T 0 onde h é o coeficiente de troca de calor por convecção e T a temperatura do meio envolvendo a superfície do contorno Condução de calor bidimensional regime permanente sem fonte de calor Q0 k e A constantes 2 0 T Equação de Laplace 10 Introdução ao MEF M A Luersen onde 2 2 2 x y operador laplaciano bidimensional Escoamento potencial bidimensional 2 0 φ onde φ x u velocidade na direção x φ y v velocidade na direção y Eletrostática d dx x A x d x dx x A x ε φ ρ 0 onde φ potencial elétrico ε permissividade A área ρ densidade de carga elétrica d dx φ campo elétrico Considerando o caso unidimensional em geral desconhecese a solução ux do problema em questão e procurase uma solução aproximada u x Tipicamente u x é um polinômio que satisfaz as condições de contorno essenciais e contém coeficientes a determinar a1 a 2 a n Assim para se obter a solução aproximada u x devese determinar os coeficientes ai tal que u e u sejam suficientemente próximas segundo Cap 1 Conceitos Introdutórios 11 um determinado critério estabelecido Ou seja a x a x a u 2 2 1 0 sendo os coeficientes ai determinados segundo critérios que serão vistos Substituindo u no lugar de u nas equações diferenciais 14 temse dois tipos de erros ou resíduos R R a x Du f D D i resíduo interior no domínio R R a x Bu g C C i resíduo no contorno 15 Em alguns casos existem somente condições de contorno essenciais Assim ao se aplicar o Método dos Resíduos Ponderados RC não necessita ser considerado pois na escolha da solução aproximada u esta deve a priori satisfazer as condições de contorno essenciais tal como no método de RayleighRitz Os resíduos podem se anular para alguns valores de x mas só serão nulos para todos os valores de x se a solução aproximada u for a solução exata isto é se u x u x Presumese que u é uma boa aproximação de u e os resíduos sejam pequenos Resíduos pequenos podem ser alcançados de várias maneiras cada uma delas resultando num sistema de equações algébricas de ordem n a ser resolvido onde as incógnitas são os coeficientes ai Algumas dessas maneiras são apresentadas a seguir 131 Método da Colocação Para n diferentes valores de x o resíduo é imposto como sendo nulo R a x D i i 0 para i1 2 j1 R a x C i i 0 para ij j1 n 16 12 Introdução ao MEF M A Luersen xi xi1 ux ux Figura 11 Representação do método da colocação onde a solução aproximada é imposta igual à solução exata nas coordenadas xi e xi 1 132 Subdomínio Sobre n diferentes regiões do domínio Ω e do contorno Γ a integral do resíduo é imposta nula R a x d D i i Ω Ω 0 para i1 2 j1 R a x d C i i Γ Γ 0 para ij j1 n 17 133 Mínimos Quadrados Os coeficientes ai são escolhidos de forma a minimizar a função I I ai 0 para i1 2 n 18 A função I é formada integrando os quadrados dos resíduos Cap 1 Conceitos Introdutórios 13 I R a x d R a x d D i C i Ω Γ Ω Γ 2 2 α 19 onde α é um escalar arbitrário que pondera a importância de RC em relação a RD 134 Colocação dos Mínimos Quadrados ou Mínimos Quadrados pontuais Similar ao anterior sendo que a função I é dada de forma discreta isto é n j i 2 i i C 1 j 1 i 2 i i D R a x R a x I α 110 135 Galerkin o mais importante Selecionamse funções peso W W x i i e impõese que a média ponderada do resíduo RD com relação às funções peso é igual a zero Em termos matemáticos RD é feito ortogonal às funções peso o produto interno entre i W e D R é nulo ver figura 12 R W x R a x d i i D i Ω Ω 0 para i1 2 n 111 RD u espaço gerado pelas funções peso u Figura 12 O resíduo RD é ortogonal ao espaço gerado pelas funções peso Wi 14 Introdução ao MEF M A Luersen No método de BubnovGalerkin ou comumentemente chamado simplesmente de Galerkin as funções peso são os coeficientes das coordenadas generalizadas ai Assim W u a i i 112 Ou seja a base de funções para aproximar u e para aproximar Wi são as mesmas No método chamado de PetrovGalerkin outras formas de Wi são utilizadas ou seja o conjunto de funções peso é diferente do conjunto de funções utilizadas para a aproximação No método de Galerkin o resíduo no contorno RC é usado em combinação com integração por partes para a imposição das condições de contorno naturais Se existir um princípio variacional associado à equação diferencial Galerkin e RayleighRitz darão soluções idênticas quando utilizase a mesma função aproximada u Exemplo Método de Galerkin Seja a equação diferencial 1 0 c x u xx para 0 x L 113a com condições de contorno 2 u 0 em x 0 Dirichlet 0 b u x em x L Neumann 113b 1 Na equação diferencial que segue e frequentemente no decorrer do texto será utilizada a notação comma u seja x representa x xy representa 2 x y etc Cap 1 Conceitos Introdutórios 15 associada com o problema de uma barra submetida a um carregamento distribuído ao longo do comprimento e um carregamento de tração na extremidade x L veja Apêndice A A barra possui seção transversal constante igual a A e módulo de elasticidade também constante igual a E A figura 13 esquematiza o problema sendo u a função deslocamento ao longo do eixo x L EA qxq ox σ0 x Figura 13 Barra submetida a tração e carregamento distribuído Os valores das constantes c e b são dados por c q EA 0 b E σ 0 114 Escolhese como função aproximada u por exemplo um polinômio de segundo grau da forma 2 2 1 a x a x u 115 Note que u satisfaz independentemente dos coeficientes ai a serem determinados a condição de contorno de Dirichlet essencial que é o deslocamento nulo em x 0 Ou seja dizse que u é um campo de deslocamentos cinematicamente admissível Aplicandose o método de Galerkin onde 0 W x R a x d R i D i i Ω Ω i1 2 n 116 16 Introdução ao MEF M A Luersen para o problema em questão temse 0 c x dx W x u R L 0 xx i i i1 2 117 O Método de Galerkin inicia com uma integração por partes reduzindo a ordem da diferenciação dentro da integral relaxando assim a ordem da continuidade requerida para a função aproximada A integração por partes também serve para introduzir as condições de contorno naturais ou de Neumann conforme será visto a seguir Integrando por partes apenas a parcela W x u i xx e lembrando a expressão da integração por partes v du u v u dv 118 fazse u Wi e dx u dv xx tendose 0 W c x dx u W W x u R L 0 i x xi L 0 x i i i 1 2 119 Como x a u W 1 1 e 2 2 2 x a u W e o termo fora da integral 0 0 u W x x i para i1 e i2 pois W x W x 1 2 0 0 0 e L b W x Lu W x i x i ver condições de contorno eq 113b temse 0 W b W cx dx W u R x L i L 0 i x xi i i 1 2 120 Cap 1 Conceitos Introdutórios 17 e 1 W x1 2x W 2 x 2a x a u 2 1 x obtendose as duas componentes do produto interno iR i 1 L 0 2 2 1 1 0 bL dx cx 2a x a R 0 bL 3 x c a x x a L 0 3 2 2 1 bL 3 L c L a La 3 2 2 1 121 i 2 L 0 2 3 2 2 1 2 0 bL dx cx 4a x 2a x R 0 bL 4 x c 3 a x 4 a x 2 L 0 4 3 2 1 2 2 4 2 3 2 1 bL 4 L c 3 L a 4 L a 122 As equações 121 e 122 formam um sistema algébrico de equações lineares cujas incógnitas são a1 e a 2 18 Introdução ao MEF M A Luersen 2 4 4 1 3 3 1 2 1 3 3 4 2 2 bL cL bL cL a a L L L L 123 sendo a solução deste sistema 12 7cL b a 2 1 e 4 cL a2 124 Assim a solução aproximada u fica 2 2 x 4 cL x 12 7cL b u 125 A solução exata deste problema é facilmente obtida sendo dada por 3 2 x 6 c x 2 cL b u 126 Cap 1 Conceitos Introdutórios 19 Comparação de Resultados software Mathcad Solução Analítica X Solução via Galerkin utilizando como função aproximada um polinômio de ordem 2 σo 45 E 21 105 A 150 qo 65 L 900 b σo E c qo E A Solução analítica u x b c L2 2 x c 6 x3 Solução via Galerkin uap x b 7 c L2 12 x c L 4 x2 x 0 05 900 0 225 450 675 900 0 15 30 45 60 u x uap x x Cálculo de Tensões σ E ε E dudx Solução Exata tensão σ x E d dx b c L2 2 x c 6 x3 σ x E b 1 2 c L2 1 2 c x2 Solução via Galerkin tensão σap x E d dx b 7 c L2 12 x c L 4 x2 σap x E b 7 12 c L2 1 2 c L x 20 Introdução ao MEF M A Luersen x 0 05 900 0 225 450 675 900 0 6250 125 104 1875 104 25 104 σ x σap x x Fim da listagem do arquivo do Mathcad Exercício Resolver o problema anterior através do Método de RayleighRitz a utilizando como função de aproximação 2 2 1 a x a x u b utilizando como função de aproximação 3 3 2 2 1 a x a x a x u Comparar as duas soluções com a solução analítica já dada anteriormente pela expressão 126 e com a solução via Galerkin utilizando 2 2 1 a x a x u expressão 125 Cap 2 Elementos Finitos Unidimensionais 21 CAPÍTULO 2 Elementos Finitos Unidimensionais 21 MÉTODO DOS ELEMENTOS FINITOS DE GALERKIN Seja o mesmo problema da barra resolvido no capítulo anterior representado pela equação diferencial 0 q x AE u xx para 0 x LT 21a e com condições de contorno u 0 em x 0 e 0 P AEu x em x LT 21b LT EA qx x P Figura 21 Barra submetida a tração Para aplicação do Método dos Elementos Finitos ao problema acima subdivide se o corpo em análise no presente caso a barra em uma série de pedaços chamados elementos finitos e aplicase o método de Galerkin em um pedaço elemento finito característico da estrutura A função aproximada no Método dos Elementos Finitos tem como incógnitas os deslocamentos em determinados pontos do elemento sendo esses 22 Introdução ao MEF M A Luersen pontos denominados nós O número de nós depende do tipo de função de aproximação que se está arbitrando Se for uma equação de primeiro grau polinômio de grau 1 isto é o deslocamento é linear ao longo do comprimento do elemento terseá dois nós localizados nas duas extremidades do elemento Se se arbitrar uma função quadrática polinômio de grau 2 terseá três nós um em cada extremidade e um terceiro intermediário e assim por diante Deste modo obtémse um sistema algébrico de equações característico de um elemento finito sendo as incógnitas os deslocamentos nodais do elemento Montase esse sistema para cada um dos elementos que formam a estrutura sendo feita depois a superposição isto é juntamse os pedaços do corpo de modo a preservar a continuidade ou seja deslocamentos em nós coincidentes serão iguais obtendose um sistema algébrico de equações para toda a estrutura No estudo de condução de calor ao invés de deslocamentos nodais como incógnitas temse temperaturas nodais No caso da barra temse uma incógnita associada com cada nó portanto temse um grau de liberdade por nó Posteriormente será visto que podese ter mais incógnitas graus de liberdade associadas a apenas um nó como é o caso da elasticidade plana onde temse como incógnitas dois deslocamentos um em cada direção e portanto o elemento finito terá dois graus de liberdade por nó Para elementos finitos unidimensionais sua representação é apenas uma linha unindo os seus nós figura 24 Na análise do problema da figura 21 dividese a barra em três elementos finitos de comprimentos L1 L2 e L3 conforme figura 22 O processo de divisão do corpo em elementos finitos chamase discretização x LT EA qx P L1 L2 L3 Figura 22 Discretização em elementos finitos Cap 2 Elementos Finitos Unidimensionais 23 Analisase um elemento característico de comprimento L conforme mostra a figura 23 L x e Figura 23 Elemento finito característico Fazse agora a aproximação u para o elemento Esta função aproximada terá como incógnitas os deslocamento nos pontos nodais Como exemplo arbitrarseá uma função aproximada de primeiro grau isto é o deslocamento ao longo do elemento será uma interpolação linear entre os deslocamentos das extremidades Portanto o elemento e terá dois nós localizados nas extremidades sendo seus graus de liberdade incógnitas os deslocamentos nessas posições ou seja u1 e u2 conforme mostra a figura 24 nó 1 nó 2 u u e L 1 2 Figura 24 Nós e graus de liberdade do elemento finito Sem perda de generalidade considerando o sistema de coordenadas com origem no nó 1 u pode ser escrito da seguinte forma u x x L u x L u 1 1 2 22 ou de uma forma mais geral 2 i 1 i i N x u u x 23 24 Introdução ao MEF M A Luersen onde Ni são as chamadas funções de interpolação sendo que i varia de 1 ao número de nós do elemento que no presente caso é igual a dois As funções de interpolação Ni portanto são dadas por 2 N x L 1 1 e N x L 2 24 Verifique que para as coordenadas nodais u é igual aos deslocamentos nodais ou seja para x 0 u u 1 e para x L u u 2 A expressão 23 pode ser expressa em forma matricial u x N u 25 onde N é denominada de matriz das funções de interpolação do elemento e u o vetor de deslocamentos nodais do elemento A expressão 25 é uma expressão geral onde N contém funções de interpolação não necessariamente lineares e u representa todos os graus de liberdade do elemento Esta forma será extensivamente utilizada no decorrer do presente texto Para uma interpolação linear N e u são dados por N 1 1 2 x L x L X e u u u X 1 2 2 1 26 Em uma forma gráfica a equação 22 pode ser representada pela figura 25 onde os deslocamentos ao longo do elemento são interpolados de forma linear pelos deslocamentos das extremidade u1 e u2 2 A obtenção das funções de interpolação será vista no item 28 Cap 2 Elementos Finitos Unidimensionais 25 x2 x3 ux x 1 u 1 u 2 u 3 Figura 25 Interpolação linear de deslocamentos Definida a função de aproximação u pelo método de Galerkin as funções peso Wi são dadas por W u u N i i i ou W N x L W N x L 1 1 2 2 1 27 e as equações de Galerkin são dadas por R W x AEu q dx i i xx L 0 0 28 Assim 0 q dx AE u N L 0 xx i 29 Integrando por partes a primeira parcela da expressão 29 26 Introdução ao MEF M A Luersen N AEu dx N AEu N AEu dx i xx L i x L i x x L 0 0 0 210 Note que no contorno AEu AE A P x ε σ carga concentrada aplicada e portanto são as condições de contorno naturais Assim temse R N A Eu N q dx N P i i x x i L i L 0 0 0 211 Para i1 N AE N u N u N q dx N x L P N x P x x x L 1 1 1 2 2 1 0 1 1 0 0 N AE N u N u dx N qdx N x L P N x P x x x L L 1 1 1 2 2 0 1 0 1 1 0 0 N AEN u N AEN u dx N qdx P x x x x L L 1 1 1 1 2 2 0 1 0 1 212 como N L 1 x 1 e N L 2 x 1 213 e multiplicando toda a expressão por 1 temse 1 L 0 L 0 2 1 P L qdx x 1 dx u L AE 1 L 1 u L AE 1 L 1 214 Integrando obtémse Cap 2 Elementos Finitos Unidimensionais 27 EA L u EA L u P x L qdx L 1 2 1 0 1 215 Procedendo de forma similar para i2 obtémse a segunda equação EA L u EA L u P x L qdx L 1 2 2 0 216 O segundo termo da direita nas expressões 215 e 216 são chamados de forças nodais consistentes pois são forças aplicadas nos nós provenientes de uma força distribuída ao longo do elemento São ditas consistentes pois são ponderadas pelas funções de interpolação Em forma matricial as duas últimas expressões podem ser escritas como EA L u u F F 1 1 1 1 1 2 1 2 217 onde F P x L qdx L 1 1 0 1 F P x L qdx L 2 2 0 218 ou de forma compacta K u F E E E 219 onde K E EA L 1 1 1 1 é a matriz de rigidez do elemento u E u u 1 2 é o vetor deslocamentos nodais do elemento F E F F 1 2 é o vetor de carregamentos nodais do elemento 28 Introdução ao MEF M A Luersen A expressão 217 ou 219 é a equação de elementos finitos para o elemento de barra considerando a utilização de funções de interpolação lineares Agora necessitase montar a equação para toda a estrutura Isto é feito superpondose todos os elementos que compõem a estrutura ou seja K u F E E NE E E NE 1 1 ou K u F 220 onde K é a matriz de rigidez global da estrutura isto é uma matriz que engloba todas as matrizes elementares u é o vetor deslocamentos nodais da estrutura F é o vetor de carregamentos nodais da estrutura A seguir é mostrado como são obtidas as matrizes e vetores globais que representam toda a estrutura Voltando agora à barra que foi divida em três elementos numerase todos os elementos e nós da estrutura A estrutura portanto terá quatro nós figura 26 x 1 2 3 4 1 2 3 Figura 26 Barra discretizada em três elementos finitos Observando a figura 26 verificase que o elemento 1 é composto pelos nós 1 e 2 o elemento 2 pelos nós 2 e 3 e o elemento 3 pelos nós 3 e 4 Esta relação entre número do elemento e numeração dos nós globais dáse o nome de conectividade ou incidência ver tabela 21 Cap 2 Elementos Finitos Unidimensionais 29 Tabela 21 Tabela da conectividade Elemento Nó 1 Nó 2 1 1 2 2 2 3 3 3 4 Montando a matriz de rigidez e vetor de carga para cada um dos elementos e para generalizar o problema considerase que cada um dos elementos tenha diferentes propriedades isto é módulo de elasticidade E1 E 2 E 3 área da seção transversal A1 A2 A3 e comprimento L1 L2 L3 Para efeitos de demonstração da superposição do vetor carga será considerado que a força distribuída é constante isto é q x q 0 Elemento 1 composto pelos nós 1 e 2 Matriz de rigidez e vetor de carga do elemento associados aos graus de liberdade u 1 e u 2 u u E A L u u 1 2 1 1 1 1 1 2 1 1 1 1 K F1 1 1 2 1 0 1 1 0 1 2 2 F F q L R q L 221 Onde R1 é a reação na direção x devido à restrição no nó 1 e também desconhecia a priori Elemento 2 composto pelos nós 2 e 3 Matriz de rigidez e vetor de carga do elemento associados aos graus de liberdade u 2 e u 3 30 Introdução ao MEF M A Luersen u u E A L u u 2 3 2 2 2 2 2 3 1 1 1 1 K F 2 1 2 2 2 0 2 0 2 2 2 F F q L q L 222 Elemento 3 composto pelos nós 3 e 4 Matriz de rigidez e vetor de carga do elemento associados aos graus de liberdade u 3 e u 4 u u E A L u u 3 4 3 3 3 3 3 4 1 1 1 1 K F 3 1 3 2 3 0 3 0 3 2 2 F F q L q L P 223 22 OBTENÇÃO DA MATRIZ DE RIGIDEZ E VETOR DE CARGA GLOBAIS A matriz de rigidez global K e do vetor de carga global F da estrutura são obtidos a partir da superposição das matrizes e vetores elementares criandose uma matriz ou vetor que englobe todas as matrizes elementares Os efeitos dos graus de liberdade coincidentes são somados da seguinte forma 3 3 3 3 3 3 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 L E A L E A 0 0 L E A L E A L E A L E A 0 0 L E A L E A L E A L A E 0 0 L E A L A E K 224 o vetor de deslocamentos globais é dado por Cap 2 Elementos Finitos Unidimensionais 31 u u u u u 1 2 3 4 225 e o vetor de carga global F F F F F F F F F F F q L R q L q L q L q L q L P 1 2 3 4 1 1 2 1 1 2 22 13 2 3 0 1 1 0 1 0 2 0 2 0 3 0 3 2 2 2 2 2 2 226 Note que a matriz de rigidez global 224 é simétrica e apresenta valores não nulos somente na diagonal principal e nas diagonais adjacentes caracterizando a chamada matriz banda Em programas de elementos finitos dependendo do tipo de solução estas propriedades diminuem o espaço de armazenamento da matriz Por exemplo podese armazenar somente os elementos da matriz que estiverem na semibanda superior ou inferior conforme representa a figura 27 onde na primeira coluna é armazenada a diagonal principal 0 d g c f b e a d g 0 0 g c f 0 0 f b e 0 0 e a Figura 27 Representação do armazenamento em banda de uma matriz simétrica Considerando os três elementos finitos do mesmo tamanho L L L L 1 2 3 mesma área de seção transversal A A A A 1 2 3 e mesmo módulo de elasticidade E E E E 1 2 3 temse o seguinte sistema linear a ser resolvido 32 Introdução ao MEF M A Luersen EA L u u u u q L R q L q L q L P 1 1 0 0 1 2 1 0 0 1 2 1 0 0 1 1 2 2 1 2 3 4 0 1 0 0 0 227 O sistema de equações acima ainda não pode ser solucionado pois a matriz de rigidez é singular Para se poder obter uma solução é necessário impor as condições de contorno isto é impor ao modelo sua vinculação suportes Caso o número de vinculações seja insuficiente a matriz continuará sendo singular Fisicamente isto indica que algum movimento de corpo rígido ainda é possível Sem vinculações suficientes a estrutura pode assumir infinitas configurações não sendo possível obter uma única solução Cap 2 Elementos Finitos Unidimensionais 33 23 IMPOSIÇÃO DAS CONDIÇÕES DE CONTORNO Para o exemplo da figura 26 temse como única vinculação o deslocamento nulo do nó 1 ou seja u1 0 228 e o sistema linear da expressão 227 pode ser reescrito como EA L u u u u q L q L q L P 1 0 0 0 0 2 1 0 0 1 2 1 0 0 1 1 0 2 1 2 3 4 0 0 0 229 onde a primeira equação de 229 foi eliminada sendo substituída pela equação u1 0 e como sabese que o deslocamento no nó 1 é nulo eliminouo das outras equações Ao invés de se resolver o sistema de 4 equações a 4 incógnitas 229 podese agora resolver o seguinte sistema de 3 equações a 3 incógnitas EA L u u u q L q L q L P 2 1 0 1 2 1 0 1 1 2 2 3 4 0 0 0 230 No caso de se ter um valor de variável prescrita deslocamento 1 u diferente de zero ou seja u u 1 1 comum em problema de distribuição de temperatura onde temse temperaturas prescritas diferentes de zero a imposição desta restrição é feita da seguinte forma 34 Introdução ao MEF M A Luersen P 2 q L L q u L EA L q L u EA u u u u 1 1 0 0 1 2 1 0 0 1 2 0 0 0 0 1 L EA 0 0 1 0 1 4 3 2 1 231 Assim de uma forma geral quando se tem o sistema de equações k k k k k k k k k k k k k k k k u u u u F F F F i n i n i i ii in n n ni nn i n i n 11 12 1 1 21 22 2 2 1 2 1 2 1 2 1 2 L L L L M M L M L M L L M M L M L M L L M M M M 232 com vinculação u u i i o sistema acima pode ser escrito como i ni n i i i2 2 i i1 1 n i 2 1 nn n2 1 n 2n 22 21 1n 12 11 k u F u k u F k u F u u u u k 0 k k 0 1 0 0 k 0 k k k 0 k k M M M M L L M L M L M M L L M L M L M M L L L L 233 A forma de imposição das condições de contorno mostrada acima é uma dentre as várias existentes Cap 2 Elementos Finitos Unidimensionais 35 24 RESOLUÇÃO DO SISTEMA ALGÉBRICO LINEAR E DETERMINAÇÃO DOS DESLOCAMENTOS NODAIS Após feita a imposição das condições de contorno resolvese por algum método numérico o sistema linear de equações obtendose assim os deslocamentos nodais Dentre os métodos de solução de sistema de equações lineares podemos destacar os métodos diretos da eliminação de Gauss decomposição LU decomposição de Cholesky método de solução frontal wavefront e os métodos iterativos de Gauss Seidel Jacobi IOR Gradiente Otimizado Steepest Descent Gradientes Conjugados entre outros Cada um desses métodos tem vantagens e desvantagens dependendo do tipo de matriz obtida simétrica nãosimétrica positiva definida matriz banda matriz esparsa condicionamento etc da disponibilidade de memória do computador precisão de resultados e tempo necessário para obtenção da solução Como para o problema estudado o sistema de equações possui pequena dimensão representado pelas expressões 229 ou 230 podese resolvêlo facilmente obtendose os seguintes resultados u u q L EA PL EA u L EA PL EA u q L EA PL EA 1 2 0 2 3 0 2 4 0 2 0 5 2 4q 2 9 2 3 234 A solução analítica para este problema é dada por u x P x EA x dx q EA L x x Px EA q EA Lx x Px EA x T 0 0 2 0 2 2 3 2 235 que nas posições correspondentes aos nós temse os mesmos valores obtidos pela solução numérica 234 36 Introdução ao MEF M A Luersen A coincidência de resultados entre a solução numérica e a solução exata nos pontos nodais ocorre apenas para problemas simples e como já mencionado na maioria dos casos não se tem uma solução analítica para comparação de resultados Nesses casos a qualidade da solução numérica obtida é analisada de outras formas utilizandose estimadores de erros Szabó Babuska 1991 Zienkiewicz Zhu 1987 Substituindo em 234 os seguintes valores numéricos 150mm A 12 x10 MPa E 6750 N P N mm 56 q 900mm L 2 5 o T 236 obtémse os valores de deslocamento nos pontos nodais u u mm u mm u mm 1 2 3 4 0 0 11071 0 20286 0 27643 237 Graficamente as duas soluções podem ser visualizadas na figura 28 0 3 00 60 0 9 00 Posição x mm 000 040 080 120 160 Deslocamento u mm Solução exata Solução numérica aproximada Figura 28 Solução exata e solução através do MEF Cap 2 Elementos Finitos Unidimensionais 37 Note que a solução analítica é um polinômio de segundo grau parábola enquanto a solução numérica MEF é linear dentro de cada elemento devido à escolha do tipo das funções de interpolação Se fossem utilizadas funções de interpolação de ordem maior um polinômio quadrático ou cúbico a solução numérica e a solução exata seriam coincidentes ao longo de toda a barra Note também que se o carregamento distribuído constante q0 for nulo a solução numérica e a solução analítica são coincidentes ao longo de toda a barra e não somente nos pontos nodais Isto era de se esperar pois a solução analítica é linear e as funções de interpolação utilizadas também são lineares 25 CÁLCULO DE DEFORMAÇÕES E TENSÕES QUANTIDADES DERIVADAS Analisarseá agora as derivadas da função aproximada u x e da solução analítica u x Essas derivadas estão associadas às deformações e consequentemente às tensões Para uma barra submetida a um esforço normal a deformação normal εx é dada por ε x du dx 238 e a tensão normal utilizando a Lei de Hooke σ ε x E x 239 Para a solução via elementos finitos uma vez calculados os deslocamentos nodais basta calcular para cada elemento finito a deformação x ε dada pela derivada dx du Como os deslocamentos são dados por 2 1 2 2 1 1 u L x u L x 1 N u N u u x 239 38 Introdução ao MEF M A Luersen A deformação no elemento é dada por L u u Lu x L u x 1 dx d dx du 1 2 2 1 x ε 240 E a tensão no elemento fica L u E u E 1 2 x x ε σ 241 onde 1 u e 2 u são os deslocamentos nodais do elemento A derivada da solução deformação e tensão via MEF é contínua apenas ao longo do elemento finito possuindo descontinuidades nos pontos nodais A derivada da solução analítica equação 235 é dada por EA P x EA 3L q EA P x EA L q dx du 0 T 0 x ε 242 e a tensão A P x A 3L q E 0 x x ε σ 243 A figura 29 mostra graficamente para os valores numéricos utilizados na expressão 236 a derivada da solução analítica e da solução numérica sendo esta última constante ao longo de um elemento finito Cap 2 Elementos Finitos Unidimensionais 39 0 30 0 6 00 90 0 Posição x mm 000 02 0 000 02 5 000 03 0 000 03 5 000 04 0 Deformação dudx Deformação derivada dudx Solução exata Solução numérica Figura 29 Derivada da solução exata e da solução via MEF Note que na interface entre os elementos finitos as deformações e consequentemente as tensões são descontínuas Para obtenção das tensões nesses pontos várias técnicas podem ser utilizadas Cook Malkus Plesha 1989 Hinton Campbell 1974 sendo a mais simples a média ponderada entre os valores na descontinuidade Essa descontinuidade das derivadas acontece devido à escolha das funções de interpolação que são da classe C 0 C zero isto é apenas a derivada de ordem zero ou seja a própria função de interesse é contínua ao longo de todo o corpo Se fossem utilizadas funções de interpolação da classe C 1 C um como por exemplo os polinômios de Hermite Ziekiewicz Taylor 1991 não só a função seria contínua como também sua derivada primeira mas em contrapartida teriase mais graus de liberdade e o sistema linear a ser resolvido seria de maior ordem Certos tipos de problemas como por exemplo na modelagem de vigas utilizando a teoria clássica viga de EulerBernoulli exigese para se ter uma adequada convergência a escolha de funções da classe C 1 Esta exigência do grau de continuidade está associada com a ordem do operador diferencial da equação que rege o problema Se o operador for de ordem 2m exigese pelo menos funções de continuidade Cm 1 Cook Malkus Plesha 1988 As afirmações acima sobre a continuidade das variáveis de interesse suas derivadas e das funções de interpolação requeridas são válidas para os chamados elementos finitos baseados em campos de deslocamentos isto é onde somente as 40 Introdução ao MEF M A Luersen variáveis primais são aproximadas deslocamentos na mecânica dos sólidos temperatura em problemas de condução de calor campo eletromagnético em problemas de eletromagnetismo etc Para outros tipos de formulações tais como as formulações híbridas e mistas as afirmações não são válidas Formulações de elementos híbridos e mistos fogem ao escopo deste texto Detalhes sobre tais elementos podem ser encontrados em Pian Tong 1987 Fazendo uma análise para o caso de q0 0 temse que a derivada da solução analítica e da solução numérica são coincidentes sendo constante ao longo de todo o comprimento da barra figura 210 Obviamente para a variável primal deslocamento u as soluções também são coincidentes ao longo de todo o corpo 0 30 0 6 00 90 0 Posição x mm 000 01 9 000 02 0 000 02 1 000 02 2 000 02 3 000 02 4 Deformação dudx Figura 210 Coincidência entre a derivada da solução exata e a derivada da solução numérica para o caso do carregamento distribuído q0 ser nulo Cap 2 Elementos Finitos Unidimensionais 41 26 CÁLCULO DAS REAÇÕES Voltando ao sistema de equações antes da imposição das condições de contorno expressão 227 tinhase EA L u u u u q L R q L q L q L P 1 1 0 0 1 2 1 0 0 1 2 1 0 0 1 1 2 2 1 2 3 4 0 1 0 0 0 244 Como já foram determinados os deslocamentos nodais podese utilizar a primeira equação para a determinação da reação R1 EA L u EA L u u u q L R 1 2 3 4 0 1 0 0 2 245 De 234 temse que u1 0 e u q L EA PL EA 2 0 2 5 2 246 Assim R q L P 1 3 0 que é o mesmo resultado obtido por equilíbrio verifique Portanto para se obter as reações primeiramente determinamse os deslocamentos nodais e em posse deles retornase ao sistema de equações original antes da imposição das condições de contorno 42 Introdução ao MEF M A Luersen Exercício Seja a barra escalona mostrada abaixo com seções transversais circulares de 10 mm e 6 mm de diâmetro e módulo de elasticidade E 7x10 MPa 4 Calcule as tensões nas barras e as reações utilizando elementos finitos 200 mm 150 mm 8 kN 27 CONDUÇÃO DE CALOR UNIDIMENSIONAL EM REGIME PERMANENTE O problema de condução de calor unidimensional com área A e condutividade térmica k constantes em regime permanente é regido pela seguinte equação diferencial k d T x dx Q x 2 2 0 247 sendo Q x a fonte interna de calor por unidade de volume Note que a equação diferencial 247 é da mesma forma que a equação da barra 21a e portanto para um elemento utilizando funções de interpolação lineares a equação de elementos finitos é dada por k L T T x L x L Q dx q q L R R 1 1 1 1 1 1 2 0 1 2 248 onde T1 e T2 são as temperaturas nos pontos nodais do elemento e qr está associado com o fluxo de calor no contorno O exemplo a seguir ilustra um problema de condução unidimensional solucionado pelo MEF Cap 2 Elementos Finitos Unidimensionais 43 Exemplo Usando cinco elementos finitos resolva o problema de distribuição de temperaturas no sólido apresentado na figura abaixo Sólido Ambiente 1 Ambiente 2 h 2000 Wm K h 5000 Wm K 2 T 100 C 1 2 ο T 50 C 2 ο 01 m oo oo 1 2 k 250 Wm K Figura 211 Problema de distribuição de temperaturas Dividese o corpo em cinco elementos finitos de mesmo comprimento conforme figura 212 01 m x x 01 m 4 2 1 3 5 6 nós Figura 212 Discretização do problema em cinco elementos finitos 44 Introdução ao MEF M A Luersen O modelo de elementos finitos possui 5 elementos e 6 nós sendo sua conectividade dada na tabela 22 Tabela 22 Conectividade para o modelo de elementos finitos da figura 212 Elemento Nó 1 Nó 2 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 A equação de elementos finitos para cada elemento pode ser obtida à partir de 248 sendo Q x 0 pois não há fonte interna de calor Como o tamanho dos elementos é igual a matriz de rigidez elementar para os cinco elementos toma a mesma forma K E k L 1 1 1 1 12500 1 1 1 1 249 Fazendose a superposição para todos os elementos finitos da mesma forma que foi feito para o problema da barra temse o seguinte sistema global antes de serem impostas as condições de contorno 6 R 5 R 4 R 3 R 2 R 1 R 6 5 4 3 2 1 q q q q q q T T T T T T 1 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 0 1 1 12500 250 Na posição x 0 nó 1 e na posição x 0 1 o fluxo de calor é dado através da condição de convecção Cap 2 Elementos Finitos Unidimensionais 45 q h T T T R1 1 1 1 1 200000 2000 q h T T T R6 2 2 6 6 250000 5000 251 Os outros valores de q R são todos nulos pois o calor que sai de um elemento é o mesmo que entra no elemento adjacente Introduzindo 248 em 247 obtémse 12500 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 12500 200000 2000 0 0 0 0 250000 5000 1 2 3 4 5 6 1 6 T T T T T T T T 252 Rearranjando de modo que o vetor independente não tenha variáveis temperaturas temse 14500 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 25000 12500 0 0 0 0 12500 17500 200000 0 0 0 0 250000 1 2 3 4 5 6 T T T T T T 253 A solução analítica para este problema é dada por T x T T T x L h h h L k 1 1 2 1 1 2 1 1 1 254 Para a posição dos pontos nodais a solução analítica 254 e a solução numérica do sistema linear 253 coincidem dando os seguintes resultados 46 Introdução ao MEF M A Luersen 59091 T 62727 T 66364 T 70000 T 73636 T 77273 T 6 5 4 3 2 1 255 Para este problema a coincidência de resultados não se dá somente na variável primal T mas também nas suas derivadas Esta coincidência entre as soluções era de se esperar pois como a fonte interna de calor é nula Q 0 a solução analítica é uma equação linear e as funções de interpolação utilizadas na aproximação via elementos finitos também são lineares similar ao caso da barra quando o carregamento distribuído q0 é nulo As derivadas de T em relação a x agora estão associadas ao fluxo de calor q k dT dx Exercício Dado o sólido da figura abaixo formado por dois materiais A e B com geração interna de calor Calcule as temperaturas no centro dos materiais e o fluxo total de calor que sai pelo lado direito da parede Despreze a resistência de contato 2 cm 01 m 2 5 cm Parede isolada T 50 C ο P k 175 W m C ο A k 35 W m C ο B A B q 500 kW m3 Cap 2 Elementos Finitos Unidimensionais 47 28 FORMULAÇÃO DE ELEMENTOS FINITOS UTILIZANDO O PRINCÍPIO DA MÍNIMA ENERGIA POTENCIAL generalização do Método de Rayleigh Ritz O Método dos Elementos Finitos pode ser definido como a aplicação do método de RayleighRitz em subregiões do domínio os elementos finitos onde os coeficientes a determinar são os graus de liberdade deslocamentos temperatura etc Na mecânica dos sólidos a energia potencial total π para um sólido elástico linear em notação indicial é expressa por Γ Ω Ω ε σ π Γ Ω Ω t u d b u d d 2 1 i i i i ij ij 256 onde ij σ são as componentes de tensão ijε as componentes de deformação ib as forças de corpo it as forças de superfície iu os deslocamentos Ω o domínio de interesse e Γ o contorno do domínio de interesse Para o problema da barra temse um estado uniaxial de tensões isto é apenas a tensão normal ao longo do eixo da barra é diferente de zero E a expressão da energia potencial total fica Γ Ω Ω σ ε π Γ Ω Ω t u d bud d 2 1 x x 257 Utilizando a Lei de Hooke σ ε x E x 258 temse Γ Ω Ω ε π Γ Ω Ω t ud bud d E 2 1 x 2 259 48 Introdução ao MEF M A Luersen Mas como ε x du dx elasticidade linear e q sendo uma força de corpo por unidade de comprimento então Γ Ω π Γ Ω t ud qudx d dx E du 2 1 L 0 2 260 A partir da expressão 260 aplicase o Método de RayleighRitz à um elemento do domínio onde as variáveis a serem determinadas são os graus de liberdade nodais Como no método de RayleighRitz a função aproximada deve satisfazer as condições de contorno cinemáticas para um elemento esta condição é representada pelo fato de que nos pontos nodais o valor da função aproximada deverá ser igual ao valor do próprio grau de liberdade nodal de modo a terse continuidade entre os elementos Assim para um elemento finito de comprimento L e de área A temse t u dx q u dx dx dx EA du 2 1 L 0 L 0 2 Γ π 261 Aproximandose o deslocamento u através de funções lineares com dois nós por elemento como já feito anteriormente na utilização do Método de Galerkin ou seja u x x L u x L u 1 1 2 ou u x N x u i i i 1 2 262 onde as funções de interpolação lineares são dadas por N x L 1 1 e N x L 2 263 a derivada da função aproximada u x em relação a x fica sendo Cap 2 Elementos Finitos Unidimensionais 49 du dx u L u L 1 2 264 A energia potencial total agora aproximada π a toma a forma tu dx qu dx dx dx EA du 2 1 L 0 L 0 2 a Γ π 265 Substituindo 264 em 265 e verificando que as forças de superfície estão associadas à forças no contorno do elemento ou seja forças aplicadas nos nós temse π a L L EA L u u dx q x x L u x L u dx P u P u 1 2 1 2 2 1 2 0 1 2 0 1 1 2 2 266 Em notação matricial 266 toma a forma 2 1 2 1 L 0 2 1 2 1 L 0 2 2 1 a P P u u xdx q L x L x 1 u u u dx u 1 1 1 1 L EA u u 2 1 π 267 Aplicase agora o Princípio da Mínima Energia Potencial PMEP que diz Dentre todas as equações de deslocamentos que satisfazem a compatibilidade interna e as condições de contorno em um sistema estável aquelas que também satisfazem as equações de equilíbrio fazem da energia potencial um mínimo Em termos matemáticos o PMEP é expresso por δπ a 0 ou π a iu 0 pois as incógnitas coeficientes no Método de RayleighRitz são os graus de liberdade nodais Assim 50 Introdução ao MEF M A Luersen π a L L L u EA L dxu EA L dxu x L qdx P 1 2 0 1 2 0 2 0 1 1 0 268 π a L L L u EA L dxu EA L dxu x L qdx P 2 2 0 1 2 0 2 0 2 0 e a equação de elementos finitos em forma matricial fica 2 1 L 0 2 1 P P dx q L x L x 1 u u 1 1 1 1 L EA 269 que é a mesma obtida pelo Método dos Elementos Finitos de Galerkin Exercício Deduzir a equação de elementos finitos do elemento de barra através do Princípio dos Trabalhos Virtuais PTV utilizando funções de interpolação lineares PTV diz que O trabalho virtual realizado pelas forças internas é igual ao trabalho virtual realizado pelas forças externas ou seja Γ δ Ω δ Ω δε σ Γ Ω Ω t u d b u d d i i i i ij ij E para a obtenção da equação de elementos finitos necessitase aproximar tanto os deslocamentos reais ui como os deslocamentos virtuais δui Cap 2 Elementos Finitos Unidimensionais 51 29 MÉTODO DIRETO O método direto consiste em obter a equação de elementos finitos à partir dos conceitos de equilíbrio Para elementos simples este método é uma alternativa em relação à utilização dos métodos energéticos ou resíduos ponderados Mas para elementos mais complexos sua utilização é praticamente inviável Para exemplificar a utilização do método direto ele será aplicado para a obtenção da equação do elemento finito para análise estrutural mais simples o elemento de barra Analisandose um elemento finito de barra de comprimento L com dois nós atuando sobre eles as forças F1 e F2 conforme esquematiza a figura 213 fazendo o equilíbrio de forças temse F F 1 2 0 ou F F 1 2 270 L F1 F2 Figura 213 Forças atuantes sobre o elemento de barra Como o elemento de barra está submetido apenas a um carregamento axial tem se somente a tensão normal σ x atuando ao longo da barra e σ x Px A 271 onde Px é o esforço normal e A a área da seção transversal do elemento Assim pode se escrever F 1 σ x A e F 2 σ x A 272 52 Introdução ao MEF M A Luersen Utilizando a lei de Hooke para o estado uniaxial de tensões σ ε x E x temse F EA x 1 ε F EA x 2 ε 273 A deformação axial ε x é dada em função do deslocamento axial ux através de sua derivada em relação à coordenada x ε x du x dx 274 Assumindo que a deformação ε x é constante ao longo do elemento podese escrever ε x L L 275 onde L é a variação do comprimento do elemento e L o comprimento original do elemento L pode ser dado em função dos deslocamentos nodais u1 e u2 L u u 2 1 276 e a deformação ε x toma a forma ε x u u L 2 1 277 Substituindose 277 em 273 temse F EA L u u 1 2 1 F EA L u u 2 2 1 278 Cap 2 Elementos Finitos Unidimensionais 53 Em notação matricial 278 toma a forma EA L u u F F 1 1 1 1 1 2 1 2 279 que é a mesma equação já deduzida anteriormente através da utilização de outros métodos para o elemento de barra Exercício Deduza utilizando o método direto a equação de elementos finitos de uma barra circular sujeita apenas à momentos torçores nas extremidades sendo os graus de liberdade as rotações ϕ 1 e ϕ 2 ver figura A matriz de rigidez será em termos de J G e L sendo J o momento polar de inércia e G o módulo de elasticidade transversal L T1 ϕ1 T2 ϕ2 Viuse até o presente momento três maneiras de se formular um problema através da utilização do Método dos Elementos Finitos sendo que para todas elas seguiuse uma metodologia similar Assim o Método dos Elementos Finitos pode ser subdividido em seis pontos básicos 1 Especificação do tipo de aproximação linear quadrática etc 2 Discretização do domínio em elementos finitos localização dos pontos nodais especificando as coordenadas bem como a conectividade 3 Desenvolvimento do sistema de equações Utilizando o Método dos Resíduos Ponderados Princípio da Energia Potencial Estacionária equilíbrio ou outro princípio obtémse um sistema de equações para cada elemento obtenção das matrizes de rigidez e vetores de carga para cada elemento A partir disso fazse a adequada superposição das matrizes e vetores obtendose a matriz de rigidez global e o vetor de cargas global 4 Imposição das condições de contorno na matriz de rigidez e vetor de carga global 54 Introdução ao MEF M A Luersen 5 Solução do sistema de equações 6 Cálculo de outras quantidades de interesse Essas quantidades estão geralmente relacionadas com as derivadas das incógnitas e são as componentes de tensão fluxo de calor velocidade de fluidos campo elétrico dependendo do problema que se está resolvendo 210 OBTENÇÃO DE FUNÇÕES DE INTERPOLAÇÃO UNIDIMENSIONAIS C 0 Neste item será mostrado como a partir de uma aproximação linear para a função de interesse φ a a x 1 2 280 obtémse as funções de interpolação lineares N i podendose assim escrever φ em função dos valores nodais φ N iφ i i 281 como feito nas expressões 22 para o problema da barra No caso de utilizarse uma aproximação linear na forma da expressão 282 o somatório na expressão 281 vai de 1 a 2 sendo as funções de interpolação N 1 e N 2 dadas como já visto anteriormente por N x L 1 1 e N x L 2 282 O método aqui apresentado apesar de ser aplicado para a obtenção de funções de interpolação lineares pode ser aplicado a funções de ordem mais alta e até mesmo à elementos bidimensionais e tridimensionais Cap 2 Elementos Finitos Unidimensionais 55 Para as coordenadas nodais φ deve apresentar os valores nodais ou seja para x 0 φ φ 1 e para x L φ φ 2 Assim temse a a 1 2 1 0 φ a a L 1 2 2 φ 283 Em forma matricial 283 pode ser escrito como 1 0 1 1 2 1 2 L a a φ φ ou A a φφ 284 onde A 1 0 1 L a a a 1 2 e φφ φ φ 1 2 Pré multiplicando 284 pela inversa da matriz A A 1 obtémse A A a A 1 1 φφ ou a A 1 φφ 285 A 1 é facilmente obtida resultando em A 1 1 0 1 1 L L 286 e os coeficientes ai são expressos em função dos valores nodais que a partir de agora são as incógnitas Assim 56 Introdução ao MEF M A Luersen a a L L 1 2 1 2 1 0 1 1 φ φ ou a1 1 φ e a L L 2 1 2 1 1 φ φ 287 Substituindo os valores de a1 e a2 na expressão 280 temse φ φ φ φ 1 1 2 1 1 L L x ou φ φ φ 1 1 2 x L x L 288 Note que 288 está escrita na forma da expressão 281 onde as funções de interpolação são exatamente aquelas expressas em 282 A representação gráfica das funções de interpolação lineares unidimensionais pode ser visualizada na figura 214 Note que Ni é unitária no nó i e nula no nó j e viceversa para a função de interpolação N j Para polinômios interpoladores de maior ordem acontece o mesmo isto é a função Ni é unitária na posição nodal correspondente ao nó i e nula nos demais nós do elemento A figura 215 representa a interpolação linear da variável φ através de funções de interpolação lineares equação 288 0 1 xi xj Ni L 0 1 xi xj Nj L Figura 214 Funções de interpolação lineares unidimensionais Cap 2 Elementos Finitos Unidimensionais 57 xi xj φi L φj x φ φ a1 a2 x Figura 215 Variação linear da variável φ Da mesma forma como foram obtidas as funções de interpolação lineares pode se obter funções de ordem maior Por exemplo para funções de interpolação quadráticas fazse φ a a x a x 1 2 3 2 289 e para as posições nodais x1 x 2 e x 3 o elemento quadrático possui 3 nós φ assume os valores nodais isto é φ φ x1 1 φ φ x2 2 290 φ φ x3 3 e a partir disto procedese da mesma forma como feito para as funções de interpolação lineares obtendose os valores dos coeficientes a1 a 2 e a 3 em função dos valores nodais φ 1 φ 2 e φ 3 Uma forma geral de obtenção de funções de interpolação polinomiais de continuidade C 0 é através da fórmula de interpolação de Lagrange Para um elemento de n nós ou seja polinômio interpolador de grau n1 a função de interpolação associada ao nó i é dada por N x x x x i j j i j j i n 1 291 58 Introdução ao MEF M A Luersen Assim por exemplo um elemento quadrático grau 2 terá 3 nós e consequentemente 3 funções de interpolação dadas por N x x x x x x x x 1 2 3 2 1 3 1 N x x x x x x x x 2 1 3 1 2 3 2 292 N x x x x x x x x 3 1 2 1 3 2 3 Graficamente as funções de interpolação quadráticas podem ser visualizadas na figura 216 A figura 217 representa a interpolação quadrática da variável φ 0 1 x1 x2 N1 x3 0 1 x1 x2 N2 x3 0 1 x1 x2 N3 x3 Figura 216 Funções de interpolação quadráticas unidimensionais Cap 2 Elementos Finitos Unidimensionais 59 x1 x2 φ1 φ2 x φ φ a1 a2 x x3 a3 x2 φ3 Figura 217 Variação quadrática da variável φ Como a expressão 288 é uma expressão geral as funções lineares também podem ser obtidas a partir dela Para este caso n2 assim N x x x x 1 2 2 1 e N x x x x x x x x 2 1 1 2 1 2 1 293 Sem perda de generalidade podese fazer x1 0 e x 2 L e assim 293 toma a forma da expressão 282 Algumas características podem ser notadas nas funções de interpolação unidimensionais de continuidade C 0 1 A função de interpolação Ni é igual a 1 no nó i e zero nos demais nós do elemento 2 a soma de todas as funções de interpolação do elemento é igual a 1 Explique por quê 3 a soma das derivadas de todas as funções de interpolação do elemento em relação a x é igual a zero consequência óbvia a partir de 2 Exercício Mostre que a matriz de rigidez para o elemento de barra utilizando funções de interpolação quadráticas expressões 292 e para os três nós igualmente espaçados é dada por K E EA L 3 7 1 8 1 7 8 8 8 16 60 Introdução ao MEF M A Luersen 211 TRANSFORMAÇÃO DE COORDENADAS PARA O ELEMENTO DE BARRA Nas deduções da matriz de rigidez e vetor de carga do elemento de barra ele estava referenciado em relação a um sistema com eixo coordenado x alinhado com o seu eixo axial figura 22 No caso de não se ter esta situação necessitase fazer uma transformação de coordenadas Com esta transformação temse a possibilidade de modelar treliças planas pois estas tem a característica de seus membros sofrerem apenas tração ou compressão A figura 218 apresenta o elemento finito em um sistema cujos eixos coordenados não estão alinhados com seu eixo axial Na figura 218a são apresentados os graus de liberdade do elemento Note que o elemento não alinhado apesar de possuir um grau de liberdade por nó ui com direção coincidente com seu eixo centroidal possuirá dois graus de liberdade por nó em relação ao sistema de coordenadas global figura 218b provenientes da decomposição daquele grau de liberdade nos eixos globais x e y L θ 1 2 u1 u2 x y θ u1 u2 x y v1 v2 a b Figura 218 Elemento de barra rotacionado em relação ao sistema de coordenadas global xy Sendo u o deslocamento ao longo do elemento u o deslocamento ao longo do eixo x e v o deslocamento ao longo do eixo y u u v cos sen θ θ 294 Cap 2 Elementos Finitos Unidimensionais 61 A partir dessa decomposição podese escrever as relações entre os deslocamentos nodais e suas componentes nas direções x e y como u u u v u v 1 2 1 1 2 2 0 0 0 0 cos sen cos sen θ θ θ θ ou u T u 295 onde u u u 1 2 é o vetor deslocamentos no sistema local u u v u v 1 1 2 2 é o vetor deslocamentos no sistema global xy e T é a matriz de transformação dada por T cos sen cos sen θ θ θ θ 0 0 0 0 Uma força P ao longo do elemento mesma direção de u pode ser decomposta em componentes nas direções x e y P x P cos θ e P y P senθ 296 e assim a relação entre as forças nodais e suas componentes nas direções x e y atuantes nos dois nós do elemento é escrita como 2 1 2 y 2 x 1 y 1 x P P sen 0 cos 0 0 sen 0 cos P P P P θ θ θ θ ou P T P T 297 onde 62 Introdução ao MEF M A Luersen P P P 1 2 é o vetor carregamento nodal no sistema local P P P P P x y x y 1 1 2 2 é o vetor carregamento no sistema global xy e T T é a transposta da matriz de transformação T dada por T T cos sen cos sen θ θ θ θ 0 0 0 0 Inserindo as relações 295 e 297 na equação de elementos finitos do elemento de barra agora expressa em termos de u1 e u2 ou seja E A L u u P P 1 1 1 1 1 2 1 2 298 e pré multiplicando por T T obtémse cos sen cos sen cos sen cos sen θ θ θ θ θ θ θ θ 0 0 0 0 1 1 1 1 0 0 0 0 1 1 2 2 1 1 2 2 E A L u v u v P P P P x y x y 299 ou K u P 2100 onde K é a matriz de rigidez do elemento de barra expressa no sistema de coordenadas global xy e dada por Cap 2 Elementos Finitos Unidimensionais 63 K T K T T EA L simé trica cos cos sen cos cos sen sen sen cos sen cos cos sen sen 2 2 2 2 2 2 θ θ θ θ θ θ θ θ θ θ θ θ θ θ 2101 sendo K a matriz de rigidez do elemento de barra no sistema local A equação 2100 ou 299 é a equação para o elemento de barra em relação a um sistema de coordenadas xy para o caso geral em que seu eixo centroidal está inclinado em relação ao eixo x de um ângulo θ A força axial ao longo de cada membro pode ser obtida a partir dos deslocamentos no sistema de coordenadas local equação 298 Inserindo 295 em 298 podese obter as forças axiais nos membros em função dos deslocamentos nodais no sistema xy que foram obtidos diretamente da resolução do sistema linear 2100 P T u AE L 1 1 1 1 ou P P AE L u v u v 1 2 1 1 2 2 1 1 1 1 0 0 0 0 cos sen cos sen θ θ θ θ 2102 Desenvolvendo as multiplicações matriciais da expressão 2100 temse P AE L u u v v P AE L u u v v 1 1 2 1 2 2 1 2 1 2 cos sen cos sen θ θ θ θ 2103 Note que P P 1 2 e assim apenas uma das equações 2103 necessita ser calculada Geralmente é escolhida a força associada ao nó 2 P2 pois pode ser interpretado que o membro está sob tração se for positiva e sob compressão se for negativa A partir de P2 podese também calcular a tensão normal de cada elemento apenas dividindose P2 pela área correspondente σ P A 2 2104 64 Introdução ao MEF M A Luersen Para montar a matriz de rigidez necessitase das propriedades E e A de cada elemento do comprimento L e do ângulo θ sendo que estes dois últimos dados geralmente não são fornecidos diretamente Podese calcular o comprimento L o coseno e o seno do ângulo θ para cada elemento a partir de suas coordenadas nodais ou seja L x x y y 2 1 2 2 1 2 cosθ x x L 2 1 e senθ y y L 2 1 2105 onde x y 1 1 são as coordenadas do primeiro nó do elemento e x y 2 2 do segundo nó do elemento Exercício Dada a estrutura apresentada na figura 219 calcular o deslocamento no ponto A as forças axiais e as tensões normais em cada membro utilizando o MEF A área de cada membro está mostrada na figura e os dois membros são de aço cujo módulo de elasticidade E é igual a 2 0 10 5 X MPa Utilizar unidades compatíveis 60 cm 80 cm A F 4000 kgf 30 o área 2 cm área 6 cm 2 2 Figura 219 Treliça plana a ser resolvida pelo MEF Referências Bibliográficas 1 Bassanezi Ferreira Jr 1988 2 Bathe KJ 1982 Finite Element Procedures in Enginnering Analysis Prentice Hall 3 Buchanan G R 1995 Theory and Problems of Finite Element Analysis Shawms Outlne Series McGraw Hill Inc 4 Carey GF Oden JT 1981 Finite Elements vol I PrenticeHall New Jersey 5 Cook RD Malkus DS Plesha ME 1988 Concepts and Applications of Finite Element Analysis 3nd Edition John Wiley Sons New York 6 Davis JP e Rabinowitz P 1984 Methods of Numerical Integration 2nd Edition Academic Press Inc Ltd London 7 Duarte CAM Oden JT 1995 Hp Clouds A Meshless Method to Solve BoundaryValue Problems TICAM Report 9505 8 Duarte CAM 1995 A Review of Some Meshless Methods to Solve Partial Differential Equations TICAM Report 9506 9 Dym Shames 1978 10 Fagan MJ 1992 Finite Element Analysis Theory and Practice Logman Scientific and Technical 11 Hinton E Campbell J S 1974 Local and Global Smoothing of Discontinuous Finite Element Functions Using a Least Square Method Int J Num Meth Eng vol 8 pp 461480 12 Hughes TJR 1987 The Finite Element Method Linear Static and Dynamic Finite Element Analysis PrenticeHall Englewood Cliffs New Jersey 13 Pian THH Tong P 1987 Mixed and Hybrid FiniteElement Methods em Finite Element Handbook Part 2 FEM Fundamentals H Kardestuncer Cap 5 McGrawHill Book Company 14 Szabó B Babuska I 1991 Finite Element Analysis John Wiley Sons Co New York USA 15 Tauchert T 1974 Energy Principles in Structural Mechanics 16 Timoshenko SP Goodier JN 1980 Teoria da Elasticidade 3a Ed Guanabara Dois Rio de Janeiro 17 Zienkiewicz Zhu 1987 A Simple Error Estimator and Adaptive Procedure for Pratical Engineering Analysis Int J Num Meth Eng vol 24 pp 335357 Apêndice A Dedução da equação diferencial para o problema da barra L EA qx σ0 x Analisandose um comprimento x da barra qx x F x x A σ F x x x x A σ x Fazendo o equilíbrio de forças na direção x temse Fx 0 0 q x F x x F x Dividindo por x 0 q x F x x F x Fazendo o limite x 0 temse 0 q dx dF x 0 q dx d x A σ Utilizando a Lei de Hooke dx E du E ε σ 0 q dx dx d EA du Se a área A e o módulo de elasticidade E forem constantes temse 0 EA q dx d u 2 2 para 0 x L Condição de contorno em xL o L x σ σ o dx E du σ 0 E dx du σ o