UNIVERSIDADE FEDERAL DO PARÁ CENTRO TECNOLÓGICO

UNIVERSIDADE FEDERAL DO PARÁ CENTRO TECNOLÓGICO DEPARTAMENTO DE ENGENHARIA CIVIL NÚCLEO DE INSTRUMENTAÇÃO E COMPUTAÇÃO APLICADA À ENGENHARIA O Método dos Elementos Finitos Aplicado ao Problema de Condução de Calor Autor: Prof. Remo Magalhães de Souza, M.Sc., Ph.D Belém 05/2003 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 1. Introdução O Método dos Elementos Finitos (MEF) consiste em um método numérico aproximado para análise de diversos fenômenos físicos que ocorrem em meios contínuos, e que são descritos através de equações diferenciais parciais, com determinadas condições de contorno (Problemas de Valor de Contorno), e possivelmente com condições iniciais (para problemas variáveis no tempo). O MEF é bastante genérico, e pode ser aplicado na solução de inúmeros problemas da engenharia. 1.1. Idéia básica do Método dos Elementos Finitos A idéia principal do Método dos Elementos Finitos consiste em se dividir o domínio (meio contínuo) do problema em sub-regiões de geometria simples (formato triangular, quadrilaeral, cúbico, etc.), conforme ilustra esquematicamente a Figura 1.1. Esta idéia é bastante utilizada na engenharia, onde usalmente tenta-se resolver um problema complexo, subdividindo-o em uma série de problemas mais simples. Logo, trata-se de um procedimento intuitivo para os engenheiros. pontos nodais elementos finitos contorno original Figura 1.1 – Malha de Elementos Finitos (para problema plano) 1 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Devido ao fato das sub-regiões apresentarem dimensões finitas, estas sub-regiões são chamadas “elementos finitos”, em contraste com os elementos infinitesimais utilizados no cálculo diferencial e integral. Advém daí, o nome “Método dos Elementos Finitos”, estabelecido por Ray Clough, na década de 50. Os elementos finitos utilizados na discretização (subdivisão) do domínio do problema são conectados entre si através de determinados pontos, denominados nós ou pontos nodais, conforme indica a Figura 1.1. Ao conjunto de elementos finitos e pontos nodais, dá-se, usualmente o nome de malha de elementos finitos. Diversos tipos de elementos finitos já foram desenvolvidos. Estes apresentam formas geométricas diversas (por exemplo, triangular, quadrilateral, cúbico, etc) em função do tipo e da dimensão do problema (se uni, bi, ou tridimensional). A Figura 1.2 apresenta a geometria de vários tipos de elementos finitos. Elemento de barra com dois nós Elemento de barra com três nós Elemento triangular com três nós Elemento quadrilateral com quatro nós Elemento triangular com seis nós Elemento tetraédrico com quatro nós Elemento quadrilateral com nove nós Elemento hexaédrico com oito nós Figura 1.2 – Diferentes tipos de elementos finitos 2 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza A precisão do método depende da quantidade de nós e elementos, e do tamanho e tipo dos elementos presentes na malha. Um dos aspectos mais importantes do MEF diz respeito a sua convergência. Embora trata-se de um método aproximado, pode-se demonstrar que em uma malha consistente, a medida que o tamanho dos elementos finitos tende a zero, e conseqüentemente, a quantidade de nós tende a infinito, a solução obtida converge para a solução exata do problema. Ou seja, quanto menor for o tamanho e maior for o número de elementos em uma determinada malha, mais precisos serão os resultados da análise. 1.2. Campos de aplicação O número de áreas de aplicação para o MEF tem crescido de forma considerável recentemente. Dentre os inúmeros campos de aplicação possíveis, podem se citar: Indústria da Construção Civil; Indústria automobilística, naval, aeronáutica e aeroespacial; Metalurgia; Mineração; Exploração de petróleo; Setor energético; Telecomunicações; Forças Armadas; Meio ambiente; Recursos Hídricos; Saúde. As primeiras aplicações do MEF foram em problemas de engenharia estrutural, mais especificametne, sobre análise de tensões. Neste tipo de problema, busca-se determinar as tensões, deformações e deslocamentos em um corpo sólido sujeito a determinadas ações tais como cargas (forças aplicadas) e recalques (deslocamentos impostos). Exemplos de tais aplicações compreendem o estudo do comportamento de estruturas civis, tais como edifícios, pontes, barragens, e túneis, onde os elementos finitos são utilizados na discretização de vigas, lajes, treliças, paredes, fundações, etc. O estudo de análise de tensões também é importante em outras áreas da engenharia, tais como engenharia mecânica, naval, aeronáutica, aeroespacial, onde são necessários análises das estruturas e peças mecânicas de máquinas, automóveis, caminhões, navios, aviões, espaçonaves, etc. Dentro da área de mecânica dos sólidos, podem ser realizadas: análise estática, análise modal (problemas de auto valor e auto-vetor, para estudo de vibrações e instabilidade estrutural), e análise dinâmica. Além da aplicação clássica do MEF na solução de problemas da mecânica dos sólidos, várias outras áreas da engenharia empregam atualmente o MEF como uma poderosa ferramenta na análise de diversos fenômenos físicos, e no projeto e análise de diversos equipamentos, dispositivos, processos industriais, etc. 3 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza A quantidade de problemas físicos que podem ser analisados com o MEF é bastante grande. A título de ilustração podem-se citar as seguintes áreas: • Transferência de calor; • Elastostática; • Elastodinâmica; • Eletroestática; • Eletromagnetismo; • Acústica; • Fadiga; • Mecânica da fratura; • Hidráulica; • Hidrodinâmica; • Aerodinâmica; • Biomecânica; • Lubrificação; • Problemas de interação fluído-estrutura; • Problemas de propagação de ondas; • Dispersão de contaminantes; Vários dos fenômenos listados acima podem ser agrupados em uma categoria especial de problema físico, denominado problema de campo (ou, mais particularmente, problema de potencial). Exemplos comums de problemas de campo são: • Condução de calor, • Condução elétrica; • Campos gravitacionais; • Campos eletroestáticos; • Campos magnetoestáticos; • Fluxo irrotacional de fluidos ideais; • Percolação através de um meio poroso; • Torsão de barras prismáticas; 4 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Os fenômenos de campo descritos acima têm em comum o mesmo tipo de equação diferencial governante, qual seja a equação quasi-harmônica. Casos particulares da equação quasi-harmônica são as conhecidas equações de Poisson, e de Laplace. No capítulo 2, apresenta-se o desenvolvimento da equação de Poisson, com aplicação do problema de condução de calor. Entretanto, o mesmo desenvolvimento pode ser aplicado a outros problemas de campos com poucas alterações. 1.3. O conceito de grau de liberdade no MEF Além dos conceitos de “elementos finitos” e “nós” no MEF, um outro conceito muito importante refere-se ao conceito de “grau de liberdade” (degree of freedom) ou, “gdl” (dof). A idéia de grau de liberdade tem sua origem na idéia do movimento de partículas em problemas da Mecânica, onde se considera que, conforme ilustra a Figura 1.3: • Um ponto apresenta, no espaço tridimensional, três graus de liberdade, quais sejam três possíveis movimentos de translação. • Mais genericamente, um corpo rígido apresenta, no espaço tridimensional, seis graus de liberdade, quais sejam, três possíveis movimentos de translação e três possíveis movimentos de rotação. (a) (b) Figura 1.3 – Graus de liberdade. a) graus de liberdade de um ponto; b) graus de liberdade de um corpo rígido. O comportamento de um elemento é praticamente definido pelo número e posicionamento dos nós, e pelo número de graus de liberdade (gdl) por nó. O mesmo elemento finito (com a mesma forma e mesmo número de nós), como por exemplo, o elemento triangular de três nós pode ser utilizado com diferentes graus de liberdade, dependendo da dimensão e tipo do problema em questão. 5 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Em problemas de mecânica dos sólidos (análise de tensões), os graus de liberdade dos nós correspondem aos possíveis movimentos que estes podem sofrer. Por exemplo, o problema de análise de tensões em um meio tridimensional apresenta três graus de liberdade por nó (três translações). No caso plano, existem dois graus de liberdade por nó (duas translações). Estes movimentos ou deslocamentos dos nós são as incógnitas principais da análise pelo método tradicional de Elementos Finitos do problema geral da Mecânica dos sólidos. Por um outro lado, no problema de condução de calor, por exemplo, embora não se estude o movimento de partículas, utiliza-se comumente o termo “grau de liberdade” para fazer referência à incógnita principal do problema, qual seja o valor do campo de temperatur nos nós da malha. 6 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 2. Problema de condução de calor O problema de condução de calor é estudado como motivação inicial para aplicação do MEF. Escolheu-se este problema pela fácil interpretação física das equações, e da sua relevância prática para diversos setores da engenharia. A apresentação é feita inicialmente, utilizando-se a forma tradicional das equações diferencias que governam o problema (“forma forte”), por tratar-se de um procedimento mais simples. Posteriormente, para possibilitar a aplicação direta do MEF na solução do problema de condução de calor, apresenta-se também a “forma fraca” da equação governante. 2.1. Forma forte das equações governantes do problema Equation Section (Next) Considere um corpo bidimensional (de espessura constante) com domínio Ω e contorno Γ , com referência a um sistema de coordenadas cartesianas (x, y) conforme ilustra a Figura 2.1. Γ Ω y x Figura 2.1 - corpo bidimensional com domínio Ω e contorno Γ , com referência a um sistema de coordenadas cartesianas (x, y). Seja Q( x, y ) a taxa de geração de calor interna ou fonte1 (calor por unidade de volume e tempo) e qx ( x, y ) e q y ( x, y ) as componentes do vetor fluxo de calor (calor por unidade de área e tempo) em um ponto (x,y) do corpo Ω 1 Fontes de calor Q são proporcionadas, por exemplo, por resistência à corrente elétrica e reações químicas. 7 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza  qx ( x, y )  q = q ( x, y ) =  .  q y ( x, y )  (2.1) A equação que governa o problema de condução de calor em um meio bidimensional em equilíbrio (regime estacionário, sem variação no tempo) pode ser facilmente deduzida considerando-se um elemento diferencial de lados dx e dy , e com fluxo de calor atravessando o contorno do elemento, conforme ilustra a Figura 2.2 qy + qx ∂q y ∂y dy qx + dy Q ∂q x dx ∂x dx qy Figura 2.2 – elemento diferencial com fluxo de calor atravessandoo contorno do elemento Considerando, sem perda de generalidade, que a espessura do corpo é unitária, a taxa de calor gerado no corpo é igual a Qd x dy . Se as faces anterior e posterior indicadas na figura forem isoladas termicamente, então a seguinte condição deve ser satisfeita Qdxdy + q x dy + q y dx = (q x + ∂q y ∂q x dx)dy + (q y + dy )dx . ∂x ∂y (2.2) Cancelando os termos repetidos, e dividindo a equação resultante por d x dy chega-se a equação que governa o problema estacionário de condução de calor − ou, de forma mais compacta ∂q x ∂q y − + Q = 0 em Ω , ∂x ∂y −divq + Q = 0 em Ω , (2.3) (2.4) ou, ainda 8 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza −∇T q + Q = 0 em Ω , onde ∂  ∂x  ∇=  ∂  ∂y  (2.5) (2.6) denota o operador diferencial nabla (ou del), tal que ∇Tq = ∂ ∂x ∂  qx  ∂qx ∂q y + = divq .  = ∂y  q y  ∂x ∂y (2.7) No caso de fluxo unidimensional, observa-se fisicamente que o fluxo de calor em uma direção é proporcional à taxa de variação da temperatura T naquela direção (Lei de Fourier). Assim, qx = −κ x ∂T , ∂x (2.8) onde κ x é o coeficiente de condutividade térmica (calor por unidade de comprimento, tempo e temperatura). Para o caso mais geral (bi ou tridimensional), observa-se que o vetor fluxo de calor é função do gradiente de temperatura T onde, para o caso bidimensional, q = −κ∇T , κ xx ( x, y ) κ xy ( x, y )  κ = κ ( x, y ) =   κ xy ( x, y ) κ yy ( x, y )  (2.9) (2.10) é a matriz de condutividade térmica, e ∂  ∂T   ∂x   ∂x  ∇T =   T =   = gradT . ∂  ∂T   ∂y   ∂y  (2.11) Assim, a eq. (2.9) pode-ser escrita como 9 O MEF aplicado ao Problema de Condução de Calor ou seja, Remo M. de Souza  ∂T  κ xx κ xy   ∂x   q x   , q  = −  ∂T  y  κ xy κ yy     ∂y  (2.12) ∂T ∂T   + κ xy q x = −  κ xx ∂x ∂y   ∂T ∂T   + κ yy q y = −  κ xy ∂x ∂y   (2.13) Substitutindo as eqs. (2.13) em (2.3), chega-se a ∂T ∂T ∂  κ xx + κ xy  ∂x ∂y ∂x  ∂T ∂T  ∂   + ∂y  κ xy ∂x + κ yy ∂y     +Q = 0.  (2.14) Se as direções cartesianas (x,y) coincidirem com as direções principais do material, então κ xy = 0 . Além disso, no caso particular de um meio isotrópico (com mesma condutividae térmica em todas as direções), tem-se que κ xx = κ yy = κ . Neste caso, a matriz de condutividade térmica, pode ser escrita como 0  κ ( x, y ) 1 0  κ = κ ( x, y ) =  = κ ( x, y )    = κ ( x, y ) I , κ ( x, y )   0 0 1  (2.15) com I sendo a matriz identidade de ordem 2. No caso de um meio homogêneo, a condutividade térmica não depende das coordenadas (x ,y), ou seja, κ xx e κ xx são constantes. Para um meio isotrópico e homogêneo, tem-se κ xy = 0 e κ xx = κ yy = κ = cte . Neste caso, a eq. (2.14) fica  ∂ 2T ∂ 2T +  ∂x 2 ∂y 2  κ  +Q = 0   (2.16) a qual é conhecida como Equação de Poisson. Esta equação governa vários dos problemas de campo importantes na engenharia. Pode-se também obter a eq. (2.16) de forma mais direta e compacta, substituindo-se a eq. (2.9) na eq. (2.5), o que resulta em 10 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza ∇ T κ ∇T + Q = 0 (2.17) κ∇ 2T + Q = 0 (2.18) Considerando-se novamente um meio isotrópico e homogêneo (com κ = κ I constante), esta equação fica onde ∇ 2 é o operador Laplaciano, tal que ∂ ∂x ∇ 2T = ∇T ∇T = ∂ ∂  ∂x  ∂ 2T ∂ 2T = + T   ∂y  ∂  ∂x 2 ∂y 2  ∂y  (2.19) Para o caso particular em que Q = 0 , ou seja, sem nenhuma fonte de calor interna, a eq. (2.18), fica ∇ 2T = 0 ou ∂ 2T ∂x 2 + ∂ 2T ∂y 2 =0 (2.20) a qual é conhecida como Equação de Laplace. 2.1.1. Condições de contorno Em geral, três diferentes tipos de condições de contorno podem ser considerados para o problema de condução de calor, quais sejam: a) Imposição de temperatura; b) Imposição de fluxo de calor; c) Imposição da relação entre temperatura e o fluxo de calor (ocorrendo na parte do contorno sujeita a convecção). Por simplicidade, serão consideradas na discussão a seguir, apenas os tipos de condições de contorno (a) e (b). Para isto, considera-se que o contorno Γ é subdivido em duas subregiões, ΓT e Γ q , conforme indica a Figura 2.3, tal que ΓT ∪ Γ q = Γ ΓT ∩ Γ q = ∅ (2.21) 11 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Γq Ω y x ΓT Figura 2.3 – subdivisão do contorno do corpo As regiões ΓT e Γ q são definidas de acordo com o tipo de condição de contorno considerada, quais sejam: a) Imposição de temperatura. Este caso corresponde ao tipo mais simples de condição de contorno, e consiste basicamente em se especificar o valor da temperatura na região ΓT do contorno, ou seja T =T onde T é a temperatura conhecida no contorno ΓT . em ΓT (2.22) b) Imposição de fluxo de calor. Neste caso, considera-se o “equilíbrio” de fluxo de calor em um elemento infinitesimal na região Γ q do contorno, conforme indica a Figura 2.4.a n̂ Ω y x Γq qx ∆s cos α ∆s α Q ∆s senα ΓT (a) qn n̂ α qy (b) Figura 2.4 – Equilíbrio de fluxo no contorno. a) Corpo com detalhe do elemento infinitesimal no contorno; b) fluxos de calor no elemento infinitesimal; 12 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza A Figura 2.4.b mostra um detalhe do elemento diferencial com as componentes de fluxo de calor que atuam no contorno do elemento. Nesta figura,  nˆ x  cos α  nˆ =   =   nˆ y  senα  (2.23) é o vetor normal unitário a superfície do contorno, com α sendo o ângulo que este vetor normal forma com o eixo das abscissas; ∆s é o comprimento da face do elemento triangular referente ao contorno do corpo; e qn é o valor conhecido do fluxo normal à superfície no contorno Γ q . De acordo com a Figura 2.4.b , para que haja equilíbrio de fluxo de calor no contorno, a seguinte equação deve ser satisfeita 1 Q ∆s 2 cos α senα + qx ∆s cos α + q y ∆s senα + qn ∆s = 0 . 2 (2.24) 1 Q ∆s cos α senα + qx cos α + q y senα + qn = 0 . 2 (2.25) Dividindo a equação por ∆s , tem-se Levando esta expressão ao limite quando o lado ∆s tende a zero, observa-se que o primeiro termo desaparece. Assim, a equação fica − q x nˆ x − q y nˆ y = qn , (2.26) ou de forma mais compacta, −qT nˆ = qn em ΓT (2.27) 2.1.2. Resumo das equações Por conveniência, as equações que governam o problema de condução de calor, na forma forte, são resumidamente apresentadas no quadro abaixo: Forma forte da equação que governa o problema −∇T q(T ) + Q = 0 em Ω (2.28) 13 O MEF aplicado ao Problema de Condução de Calor Relação constitutiva do meio (Lei de Fourier) q(T ) = −κ∇T Remo M. de Souza em Ω (2.29) Condições de contorno T =T −qT nˆ = qn em ΓT   em Γ q  (2.30) O problema de condução de calor consiste em se resolver a equação diferencial parcial (PDE, Partial Differential Equation) (2.28), considerando a relação constitutiva (2.29) do material, e satisfazendo as codições de contorno (2.30). Este tipo de problema é comumente denominado problema de valor de contorno (BVP, Boundary Value Problem) As equações apresentadas no quadro acima são expressas na chamada forma forte, o que significa que estas equações devem ser satisfeitas pontualmente, ou seja, a solução do problema consiste em satisfazes estas equações, para qualquer ponto (x,y) do meio. 2.1.3. Forma fraca das equações governantes do problema A obtenção da forma fraca das equações que governam o problema consiste no estabelecimento de equações integrais sobre o domínio Ω e o contorno Γ do corpo, referentes à satisfação destas equações em um sentido “médio” (ao contrário do sentido restrito pontual da forma forte). Nos desenvolvimentos seguintes, utiliza-se a notação compacta de diferenciação em relação as variáveis x e y, tal que ϕ, x = ∂ϕ ∂ϕ ∂ 2ϕ ∂ 2ϕ , ϕ, y = , ϕ, xx = , ϕ, xy = ∂x ∂y ∂y∂x ∂y 2 (2.31) A forma fraca das equações governantes pode ser obtida seguindo-se os seguintes passos: 1o Passo) Multiplica-se a eq. (2.28) por uma função arbitrária w( x, y ) , denominada, função teste, tal que ( ) w( x, y ) −∇T q(T ) + Q = 0 (2.32) 14 O MEF aplicado ao Problema de Condução de Calor ou seja, w( x, y ) ( qx, x + q y, y − Q ) = 0 Remo M. de Souza (2.33) 2o Passo) Integra-se a equação acima sobre o domínio Ω ∫ w( x, y) ( qx, x + q y, y − Q ) dΩ = 0 Ω (2.34) Observa-se que caso se determine o campo de temperatura T ( x, y ) que resolva a eq. (2.28), ou seja, caso a eq. (2.28) seja satisfeita, então a eq. (2.34) também é automaticamente satisfeita para qualquer que seja a função w( x, y ) . Por outro lado, pode-se demonstrar que caso se determine um campo de temperatura T ( x, y ) que satisfaça a eq. (2.34) para qualquer que seja a função w( x, y ) , então este campo é a solução da equação (2.28). 3o Passo) Faz-se integração por partes da equação acima, utilizando-se o teorema integral de Gauss. Para isso, inicialmente transfere-se as derivadas da função qx e q y para a função teste w , utilizandose a regra de derivada do produto wqx, x = ( wqx ), x − w, x qx wq y , y = ( wq y ), y − w, y q y (2.35) ∫ w ( qx, x + q y, y − Q ) dΩ = tal que Ω ∫ ( (wqx ), x + (wq y ), y − w, x qx − w, y q y − wQ ) dΩ (2.36) Ω ∫ ( (wqx ), x + (wq y ), y ) dΩ = ∫ ( (wqx )nˆ x + (wq y )nˆ y ) dΓ Utiliza-se, então, o teorema integral de Gauss, tal que Ω Γ = ∫ w ( qx nˆ x + q y nˆ y ) dΓ Γ = ∫ wqT nˆ dΓ (2.37) Γ 15 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza ∫ w ( qx, x + q y, y − Q ) dΩ Assim, substitutindo-se a eq. (2.37) em (2.36), chega-se a 0= = Ω T ∫ ( − w, x qx − w, y q y − wQ ) dΩ + ∫ wq nˆ dΓ Ω (2.38) Γ 4o Passo) Considera-se as condições de contorno (2.30), tal que ∫ wq Γ n dΓ = Tˆ = ∫ ΓT ∫ ΓT wqT nˆ dΓ + wqT nˆ dΓ − ∫ wq Γq ∫ Γq n dΓ Tˆ wqn dΓ (2.39) Observa-se que o termo qT nˆ é igual a qn em Γ q , sendo portanto perfeitamente conhecido nesta região do contorno. Porém, o termo qT nˆ é desconhecido em ΓT . Esta dificuldade pode ser resolvida eliminando-se a incógnita q , através da seguinte restrição na função teste w( x, y ) = 0 tal que o termo ∫ ΓT em ΓT wqT nˆ dS = 0 (2.40) (2.41) As funções teste w( x, y ) que satisfazem a condição (2.40) são denominadas funções admissíveis. Assim, substituindo-se a eq. (2.41) em (2.39), e o resultado em (2.38), esta pode ser reescrita como ∫ ( − w, x qx − w, y q y − wQ ) dΩ − ∫ wqn dΓ = 0 Ω ou seja, de forma mais compacta Γq − ∫ ( ∇w )T q dΩ = Ω ∫ wQ dΩ + ∫ wqn dΓ Ω Γq (2.42) (2.43) Esta equação consiste na forma fraca do problema de condução de calor. Observa-se que esta equação independe das propriedades do material (relações constitutivas) do meio. 16 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Considerando a relação constitutiva do meio dada pela eqs. (2.29) (Lei de Fourrier), e levando em conta ainda que w é uma função escalar ( w = wT ), a eq. (2.43) pode ser escrita como ∫ ( ∇w ) Ω T κ (∇T )dΩ = ∫w Ω T Q dΩ + ∫w Γq T qn dΓ (2.44) Esta equação consiste na forma fraca do problema de condução de calor para um material obedecendo as relações constituivas do material referentes às leis de Fourrier. A solução do problema nesta forma fraca consiste em se determinar o campo de temperaturas T ( x, y ) satisfazendo a eq. (2.44), para toda função teste w( x, y ) admissível. Uma solução aproximada para este problema pode ser obtida através do método dos Elementos Finitos, descrito na próxima seção. 17 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 3. Aplicação do MEF ao problema de Condução de Calor Conforme discutido anteriormente, a idéia básica do MEF consiste em discretizar (subdivir) o domínio do problema utilizando-se uma malha de elementos finitos. Na malha, os elementos são interligados através dos nós, conforme indica a Figura 1.1. O passo inicial para utilização do MEF consiste na etapa de criação da malha de elementos finitos pelo usuário. Para isso, o usuário especifica a localização dos nós, utilizando-se um sistema de coordenadas cartesianas, em um posicionamento arbitrário, conforme ilustra a Figura 3.1. 7 3 10 ( x3 , y3 ) 4 ( x4 , y4 ) 2 ( x2 , y2 ) 5 9 6 1 11 ( x1, y1) y 8 12 x Figura 3.1 – Especificação da posição dos nós da malha. Em geral, o número de graus de liberdade por nó da malha está relacionado com o tipo e a dimensão do problema em questão. No caso de problema de potencial, o objetivo inicial é a determinação de um campo escalar correspondente à solução do problema. Por exemplo, no problema de condução do calor, objetiva-se determinar o campo de temperaturas, o qual consiste em um campo escalar. Neste caso, os elementos empregados na análise devem possui um grau de liberdade (gdl) por nó, independentemente da dimensão do problema (se uni, bi ou tridimensional). No problema de condução de calor, quando se utiliza o MEF, as incógnitas principais do problema são as temperaturas nodais, ou seja, são os valores do campo de temperaturas avaliados nos nós da malha. Essas temperaturas nodais podem ser armazenadas em um arranjo unidimensional (vetor) da seguinte maneira 18 O MEF aplicado ao Problema de Condução de Calor  T1  T   2    T =  T3  ,  #    TN g  Remo M. de Souza (2.45) onde T1 é a temperatura correspondente ao gdl 1, T2 é a temperatura correspondente ao gdl 2, e assim por diante, até o número de graus de liberdade N g da malha. Através do MEF, a equação diferencial que governa o problema é transformada em um sistema de equações algébricas do tipo KT = F (2.46) Onde K é uma matriz de condutividade do problema, (em geral denominada matriz de rigidez), de ordem N g × N g , e F é um vetor de coeficientes (em geral denominado vetor de forças), de ordem N g ×1 , e T é o vetor de incógnitas. No caso do problema de condução de calor, F tem o sentido de fontes concentradas de calor (calor por unidade de tempo) nos nós da malha  F1   F   2    F =  F3   #     FN g  (2.47) onde F1 é a fonte de calor correspondente ao gdl 1, F2 é a fonte correspondente ao gdl 2, e assim por diante, até o número de graus de liberdade N g da malha. A idéia de fonte de calor concentrada pode ser inserida no contexto do problema contínuo desenvolvido anteriormente considerando-se uma função Q ( x, y ) taxa de geração de calor, pontual, ou seja, uma função nula em todo o domínio, exceto em um determinado ponto P (função singular, ou delta de Dirac). Utilizando-se esta idéia, tem-se que (ver eq. (2.44)) ∫ w( x, y)Q ( x, y) dΩ = w( xP , yP ) Fc P = wP Fc P Ω (2.48) 19 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza onde FcP corresponde a fonte de calor concentrada no ponto P, e wP corresponde à função teste avaliada neste ponto. Como podem haver várias fontes de calor concentras na malha de elementos finitos, é conveniente reescrever a eq. (2.44) como, T T ∫ ( ∇w) κ (∇T )dΩ = ∫ w Q dΩ + Ω onde w é um vetor tal que = Ω ∫w Ω T Q dΩ + ∫ Γq wT qn dΓ + ∫w Γq T ∑ wP Fc P NP P =1 qn dΓ + w Fc (2.49) T  w( x1, y1 )   w1   w( x , y )   w  2 2   2   w( x , y )   w  w= 3 3 = 3    #  #      w( xg , y g )   wg  (2.50)  Fc1   F   c2    Fc =  Fc3   #     Fc N g  (2.51) e é o vetor contendo as fontes de calor concentradas nos nós da malha. Caso não exista fonte de calor concentrada em um determinado nó, então a respectiva componente no vetor Fc é nula. 3.1. Formulação do elemento finito triângular linear para o problema de condução de calor bidimensional A apresentação a seguir utiliza a idéia intuitiva, comum na engenharia, de se resolver um problema complexo, sudvidindo-o em problemas menores mais simples. Assim, o desenvolvimento a seguir mostra inicialmente a formulação de um elemento finito simples, e posteriormente, apresenta 20 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza como este elemento é considerado na solução do problema global, considerando toda a malha de elementos finitos. A partir dos valores das temperaturas nos nós de um elemento pode-se determinar o valor do campo de temperatura em um ponto qualquer no interior do elemento, realizando-se uma interpolação dos valores nodais. Esta interpolação pode ser linear, quadrática, ou referente a qualquer outra função polinomial, dependendo do número de nós do elemento. Na verdade, pode-se também utilizar outras funções de interpolação além das funções polinomiais, tais como funções trigonométricas, exponenciais, etc. Um dos elementos finitos mais simples já desenvolvidos é o elemento finito triangular com interpolação linear. Este elemento apresenta uma forma triângular, com três nós I, J, e K posicionados nos vértices do triângulo, conforme indica a Figura 3.2. K ( xK , y K ) I ( xI , y I ) y J ( xJ , y J ) x Figura 3.2 – Elemento finito triangular linear, com referência ao sistema de eixos cartesianos. Na Figura 3.2 estão indicadas as coordenadas ( xI , yI ) , ( xJ , y J ) e ( xK , y K ) , dos nós I, J, e K, respectivamente, do elemento triangular. Estas coordenadas são fornecidas como dados de entrada do problema. O elemento triangular linear, quando utilizado em problemas de condução de calor, possui um grau de liberdade por nó, totalizando três graus de liberdade, quais sejam os valores TI , TJ , e TK . Estes graus de liberdade correspondem ao valor do campo de temperatura avaliado nos nós I, J, e K do elemento. Estes graus de liberdade são armazenados no vetor de temperaturas nodais Te do elemento  TI    T =  TJ    TK  e (2.52) 21 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza A seguir determinam-se as funções de intepolação do elemento, as quais permitem calcular o valor do campo de temperatura T em um ponto ( x, y ) qualquer no interior deste elemento. 3.1.1. Funções de forma do elemento A formulação do elemento triangular linear baseia-se na hipótese de que, no interior do elemento, o campo de temperatura seja uma função linear das coordenadas ( x, y ) . Assim, assume-se o seguinte campo de temperatura T ( x, y ) = a1 + a2 x + a3 y (2.53) onde a1 , a2 e a3 são constantes a serem determinadas. A eq. (2.53) pode ser escrita de forma mais compacta como T ( x, y ) = x( x, y )a onde x= 1 x e y  a1    a = a2     a3  (2.54) (2.55) (2.56) O vetor a, contendo as constantes, pode ser determinado através da imposição do valor da temperatura em cada nó, ou seja T ( xI , yI ) = a1 + a2 xI + a3 yI = TI T ( xJ , y J ) = a1 + a2 xJ + a3 y J = TJ T ( xK , yK ) = a1 + a2 xK + a3 yK = TK (2.57) As eqs. (2.57) podem ser reescritas na forma matricial como 1 xI 1 x J  1 xK yI   a1   TI      y J   a2  =  TJ  yK   a3  TK  (2.58) ou, em forma, mais compacta Ga = Te (2.59) onde 22 O MEF aplicado ao Problema de Condução de Calor 1 xI G = 1 xJ  1 xK Remo M. de Souza yI  y J  y K  (2.60) é uma matriz contendo as coordenadas dos nós do elemento. Pode-se determinar o vetor de constantes a, invertendo-se a eq. (2.59) a = G −1Te onde  ( x J y K − xK y J ) ( x K y I − xI y K ) ( xI y J − x J y I )  1  (y − y ) ( yK − yI ) ( yI − y J )  G −1 = J K  det(G )  ( xK − xJ ) ( x I − xK ) ( xJ − xI )  com det(G ) = ( xJ yK + xI y J + xK yI ) − ( xJ yI + xK y J + xI yK ) (2.61) (2.62) (2.63) Conforme apresentado no Apêndice A, o determinante da matriz G corresponde à duas vezes a área At do elemento, ou seja 2 At = det(G ) (2.64) Substituindo a eq. (2.61) na eq. (2.54), chega-se a ou ainda, onde T ( x, y ) = x( x, y )G −1Te (2.65) T ( x, y ) = N ( x, y )Te (2.66) N( x, y ) = x( x, y )G −1 (2.67) Considerando a eq. (2.64), e desenvolvendo o produto dado na equação acima, chega-se a N( x, y ) = N1 ( x, y ) N 2 ( x, y ) N3 ( x, y ) (2.68) onde 23 O MEF aplicado ao Problema de Condução de Calor N1 ( x, y ) = N 2 ( x, y ) = N 3 ( x, y ) = Remo M. de Souza 1 ( ( x J y K − xK y J ) + ( y J − y K ) x + ( xK − x J ) y ) 2 At 1 ( ( x K y I − xI y K ) + ( y K − y I ) x + ( xI − xK ) y ) 2 At (2.69) 1 ( ( xI y J − x J y I ) + ( y I − y J ) x + ( x J − xI ) y ) 2 At A matriz N( x, y ) é uma matriz contendo funções de interpolação dos graus de liberdade nodais, neste caso, das temperaturas nodais. Esta matriz é usualmente denominada matriz de funções de forma. Observando a eq. (2.66), conclui-se que a matriz de funções de forma permite determinar a temperatura T em um ponto ( x, y ) qualquer do elemento, a partir dos valores das temperaturas nodais Te . Uma interpretação geométrica dos termos da matriz de funções de forma é apresentada no Apêndice A. 3.1.2. Derivadas das funções de forma do elemento Na solução do problema de condução de calor, torna-se necessário o cálculo de derivadas do campo de temperatura T ( x, y ) , conforme se observa na eq. (2.44) onde considera-se o gradiente deste campo. As derivadas do campo de temperatura, na formulação do elemento finito, pode ser facilmente calculada a partir da eq. (2.66) onde ∇T ( x , y ) = ∇N ( x , y ) T e = B ( x , y ) T e (2.70) ∂  ∂x  B( x, y ) = ∇N( x, y ) =   N1 ( x, y ) N 2 ( x, y ) N3 ( x, y ) ∂  ∂y  (2.71)  ∂N1 ( x, y )  ∂x =  ∂N1 ( x, y )  ∂y ∂N 2 ( x, y ) ∂x ∂N 2 ( x, y ) ∂y ∂N3 ( x, y )   ∂x  ∂N3 ( x, y )   ∂y Calculando as derivadas das funções de forma em relação às coordenadas cartesianas (ver eqs. (2.68) e (2.69)), chega-se a 24 O MEF aplicado ao Problema de Condução de Calor B= 1 2 At Remo M. de Souza ( y J − yK ) ( yK − yI ) ( yI − y J )  (x − x ) (x − x ) (x − x )  J I K J I   K (2.72) 3.1.3. Interpolação da função teste Analogamente ao campo de temperatura T ( x, y ) , a função teste w( x, y ) no interior de cada elemento também pode ser obtida por interpolação de valores nodais  wI    w =  wJ  w   K e (2.73) utilizando-se a mesma matriz de funções de forma N ( x, y ) (ver eq. (2.66)). Ou seja, w( x, y ) = N( x, y )w e (2.74) Da mesma forma, o gradiente de w( x, y ) , necessário na eq. (2.44), também pode ser obtido de forma similar ao gradiente da temperatura T ( x, y ) (ver eq. (2.70)), resultando em ∇w( x, y ) = B( x, y )w e (2.75) 3.1.4. Matriz de condutividade do elemento A eq. (2.49) é válida para um domínio Ω de formato qualquer, e também, para qualquer subdomíno dentro do domínio original. Assim, pode-se particularizar esta equação para um subdomínio referente a um elemento finito, ou seja, ∫e ( ∇w) Ω T κ (∇T )dΩ = ∫e w T Ω Q dΩ + ∫e w T Γq qn dΓ + w e Fce T (2.76) onde Ωe , e Γe correspondem, respectivamente, ao domínio e contorno de um elemento, e  Fce   I e  e  Fc =  Fc J   e   Fc K  (2.77) 25 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza é o vetor de fluxos nodais do elemento. A substituição das eqs. (2.70), (2.75), e (2.74) em (2.76), leva a we T ∫e B T Ω κB dΩ Te = w e ∫e N T T Ω Q dΩ + w e T ∫e N T Γq qn dΓ + w e Fce T (2.78) T Colocando-se o termo w e em evidência, chega-se a w e  ∫ BT κB dΩ Te − ∫ N T Q dΩ − ∫ N T qn dΓ − Fce  = 0  e  e e Ω  Ω Γ q   T (2.79) O vetor w e deve ser completamente arbitrário; portanto, conclui-se que ∫e B Ω T ∫e N κB dΩ Te = Ω T Q dΩ + ∫e N Γq T qn dΓ + Fce (2.80) ou, de forma mais compacta, K e Te = F e onde Ke = ∫e B Ω T (2.81) κB dΩ (2.82) é matriz de condutividade (ou de rigidez) do elemento, e Fe = ∫e N Ω T Q dΩ + ∫e N T Γq qn dΓ + Fce (2.83) é o vetor de fontes nodais total (ou vetor de forças) do elemento. Os dois primeiros termos do lado direito da eq. (2.83), correspondem às parcelas de fluxo nodal equivalentes aos fluxos de calor “distribuído” no domínio (por unidade de volume) e no contorno do elemento (por unidade de superfície). Estas parcelas podem ser agrupadas em um único vetor, Fqe = ∫e N Ω T Q dΩ + ∫e N Γq T qn dΓ (2.84) tal que F e = Fqe + Fce (2.85) 26 O MEF aplicado ao Problema de Condução de Calor 3.2. Remo M. de Souza Montagem da matriz de condutividade e do vetor de fontes nodais do modelo A partir das matrizes de condutividade e vetores de fontes nodais dos elementos que formam a malha de elementos finitos, pode-se obter a matriz de condutividade e a matriz de fontes nodais do modelo. Para isso, será utilizado como exemplo uma malha simples ilustrada na Figura 3.3. 4 5 d 1 b c a 2 3 Figura 3.3 – Malha simples de Elementos Finitos Na discussão seguinte, considera-se a conectividade dos elementos para a malha da Figura 3.3, conforme especificada na Tabela 3.1. Elemento Nó I Nó J Nó K a 2 3 1 b 2 1 4 c 3 5 1 d 1 5 4 Tabela 3.1 – Tabela de incidência nodais dos elementos 27 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 3.2.1. Relação entre os vetores de temperaturas nodais do elemento e do modelo Seja T o vetor de temperaturas nodais da malha representada na Figura 3.3  T1  T   2  T = T3  T   4 T5  (2.86) O vetor de temperaturas nodais do elemento a é TIa    T2  a  a   T = TJ  = T3   a  T  TK   1  (2.87) Observa-se que a equação acima pode ser escrita como o seguinte produto matricial  T1  TIa      0 1 0 0 0  T2      Ta = TJa  = 0 0 1 0 0  T3   a  1 0 0 0 0  T    4 T  K  T5  ou ainda, onde Ta = H a T 0 1 0 0 0 H =  0 0 1 0 0  1 0 0 0 0  a (2.88) (2.89) (2.90) é denominada matriz de incidência do elemento a. Procedendo-se de forma análoga, chega-se as seguintes relações para os outros elementos Tb = Hb T Tc = H c T Td = H d T (2.91) onde 28 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 0 1 0 0 0  H = 1 0 0 0 0   0 0 0 1 0  0 0 1 0 0 c  H =  0 0 0 0 1  1 0 0 0 0  b (2.92) 1 0 0 0 0  H =  0 0 0 0 1   0 0 0 1 0  d são as matrizes de incidência dos elementos b, c e d, respectivamente. Deve-se notar que o número de linhas da matriz de incidência de um elemento é igual ao número de graus de liberdade do elemento (neste caso, três, para o elemento triangular linear), e o número de colunas é igual ao número de graus de liberdade do modelo (neste caso, cinco, para a malha mostrada na Figura 3.3). A matriz de incidência de cada elemento é facilmente determinada pelas seguintes regras simples: 1) Para a primeira linha, tem-se o valor um na coluna correspondente ao nó I do elemento, com as demais colunas iguais a zero; 2) Para a segunda linha, tem-se o valor um na coluna correspondente ao nó J do elemento, com as demais colunas iguais a zero; 3) Para a terceira linha, tem-se o valor um na coluna correspondente ao nó K do elemento, com as demais colunas iguais a zero; ou seja, de forma mais geral, tem-se para cada linha n, o valor 1 na coluna correspondente ao grau de liberdade do nó n, com as demais colunas iguais a zero. 3.2.2. Relação entre os vetores de fontes nodais do elemento e do modelo Seja F o vetor de fontes nodais da malha mostrada na Figura 3.3 29 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza  F1  F   2  F3  F=   F4   F5     F6  (2.93) O vetor de fontes nodais do elemento a é  FIa    a  a F =  FJ   a F  K  Analogamente, para os outros elementos, têm-se  FIb    b  b F =  FJ   b  FK   FIc    c  c F =  FJ   c  FK  (2.94)  FId    d  d F =  FJ   d  FK  (2.95) Uma condição de “equilíbrio” a ser satisfeita, é que a soma das fontes nodais, referentes a um nó comum, de cada elemento que estão conectados a este nó comum, deve ser igual à fonte nodal total aplicada deste nó. Ou seja, para o nó 1, tem-se F1 = FKa + FJb + FKc + FId (2.96) Analogamente, têm-se para os demais nós F2 = FIa + FIb F3 = FJa + FIc F4 = FKb + FKd (2.97) F5 = FJc + FJd As eqs. (2.96) e (2.97) podem ser escritas na seguinte forma matricial  F1  0  F  1  2    F3  = 0  F  0  4   F5  0 0 1 0  FIa    0 0    1   1 0   FJa  +  0   0 0  F a  0  K  0 0 0  1 0 0  FIb    0 0   0   0 0   FJb  + 1   0 1   F b  0  K  0 0 0  0 1 1  FIc    0 0   0   0 0   FJc  + 0   0 0  F c  0  K  0 1 0  0 0  d 0 0   FI    0 0   FJd   0 1  F d   K  1 0  (2.98) 30 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Observa-se que as matrizes mostradas na equação acima correspondem às matrizes transpostas das matrizes de incidência de cada elemento. Assim, a equação acima pode ser escrita, de forma compacta, como F = H a F a + H b Fb + H c F c + H d F d (2.99) F = H1 F1 + H 2 F 2 + H3 F3 + H 4 F 4 (2.100) T T T T Para escrever a equação acima na forma de somatório, pode-se definir a ≡ 1 , b ≡ 2 , c ≡ 3 e d ≡ 4 . Assim, T T ou seja, T T F = ∑ He Fe ne T (2.101) e=1 onde e representa um elemento genérico e ne é o número de elementos da malha (neste caso, ne = 4 ). Através da eq. (2.101), pode-se, portanto, determinar o vetor de fontes nodais da malha. 3.2.3. Obtenção da matriz de condutividade do modelo A eq. (2.81) apresenta a relação entre os vetores de temperaturas e fontes nodais de um elemento genérico e. Assim, para a malha da Figura 3.3, pode-se escrever F a = K a Ta Fb = K b Tb F c = K c Tc F d = K d Td (2.102) Substituindo as equações acima na eq. (2.99), tem-se F = H a K a Ta + Hb K b Tb + H c K c Tc + H d K d Td T T T T (2.103) Substituindo, agora, as eqs. (2.89) e (2.91) na equação acima, tem-se F = H a K a H a T + Hb K b HbT + H c K c H c T + H d K d H d T T T T T (2.104) Colocando o vetor de temperaturas nodais do modelo T em evidência, chega-se a F = ( H a K a H a + H b K b H b + H c K c H c + H d K d H d )T (2.105) F = KT (2.106) T ou ainda, T T T onde 31 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza K = H a K a H a + Hb K b Hb + H c K c H c + H d K d H d T T T T (2.107) é a matriz de condutividade (ou de rigidez) do modelo. Esta equação também pode ser escrita na forma de somatório, seguindo a idéia utilizada na eq. (2.101). Assim, em geral, pode-se obter a matriz de rigidez do modelo, através do seguinte somatório K = ∑ He K e He ne T (2.108) e=1 Através da eq.(2.108), pode-se, portanto, determinar a matriz de condutividade da malha. Embora as eqs. (2.101) e (2.108) tenham sido deduzidas para o malha mostrada na Figura 3.3, estas equações são completamente genéricas, podendo ser utilizadas para qualquer outra malha de elementos finitos. Para isso, deve-se apenas determinar a matriz de incidência associada a cada elemento para a malha em questão. Deve-se ressaltar que embora as eqs. (2.101) e (2.108) representem teoricamente a forma de obtenção do vetor de fontes nodais, e da matriz de condutividade do modelo, na prática, este processo torna-se ineficiente para malhas refinadas (com muitos elementos), em função da grande quantidade de zeros presente nas matrizes de incidência. Assim, na prática, utiliza-se um algoritmo computacional para montagem do vetor de fontes nodais e matriz de condutividade da malha, o qual evita a multiplicação desnecessária dos números zero presentes nas matrizes de incidência dos elementos. 3.3. Imposição das condições de contorno e solução do sistema de equações Caso as temperaturas nodais de todos os nós da malha fossem conhecidas, as fontes nodais poderiam ser facilmente determinadas através da eq. (2.106). Entretanto, em situações práticas, se conhece a temperatura nodal de alguns nós, e a fonte nodal dos demais nós. Para a determinação dos valores desconhecidos das fontes e temperaturas nodais, deve-se numerar os graus de liberdade da malha, de tal maneira que os nós com fonte nodal prescrita (conhecida) sejam numerados primeiro, e os nós com temperatura prescrita (conhecida) sejam numerados por último. Por exemplo, considera-se a malha da Figura 3.3. Neste caso, a eq. (2.106) pode ser escrita na forma expandida como 32 O MEF aplicado ao Problema de Condução de Calor  K11 K  21  K31   K 41  K51 K12 K13 K14 K 22 K 23 K 24 K32 K33 K34 K 42 K 43 K 44 K52 K53 K54 Remo M. de Souza K15   T1   F1  K 25  T2   F2      K35  T3  =  F3   K 45  T4   F4      K55  T5   F5  (2.109) Considera-se agora, que as fontes nodais sejam prescritas para os nós 1, 2 e 3, e que as temperaturas sejam prescritas para os nós 4 e 5. Com isso, pode-se particionar o sistema acima, da seguinte maneira,   K11    K 21   K31    K 41    K51 K12 K 22 K32 K 42 K52 ou de forma mais compacta onde  K11 K 00 =  K 21  K31 K K10 =  41  K51 K13  K 23  K33  K 43   K 44 K53   K54  K 00 K  10 K12 K 22 K32 K 42 K52  K14 K  24  K34 K15    T1   F1        K 25    T2   F2   =  F  K35   T3   3     K 45   T4   F4      K55   T5   F5  K 01  T0  F0   =  K11   T1   F1  K13   K14  K 23  K 01 =  K 24 K33   K34 K 43  K K11 =  44  K53   K54 K15   T1   F1       K 25  T0 = T2  F0 =  F2      K35  T3   F3  K 45  T  F  T1 =  4  F1 =  4   K55  T5   F5  (2.110) (2.111) (2.112) Desta forma, deve-se notar que os vetores T1 e F0 são conhecidos, ao passo que os vetores T0 e F1 são desconhecidos. Desenvolvendo a eq. (2.111), tem-se K 00T0 + K 01T1 = F0 K10T0 + K11T1 = F1 (2.113) Com isso, o vetor de temperaturas nodais T0 pode ser calculado, a partir da primeira das equações acima, T0 = K 00−1 ( F0 − K 01T1 ) (2.114) 33 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Após a determinação de T0 , o vetor de fontes nodais F1 , pode ser calculado diretamente utilizando-se a segunda das eqs. (2.113) 3.4. F1 = K10T0 + K11T1 (2.115) Resumo das etapas de análise pelo MEF As etapas de análise pelo Método dos elementos finitos são descritas resumidamente abaixo: 1) Montagem da matriz de condutividade do material (eq. (2.10)) para cada elemento κ xx ( x, y ) κ xy ( x, y )  κ = κ ( x, y ) =   κ xy ( x, y ) κ yy ( x, y )  (2.116) 2) Montagem da matriz com as derivadas das funções de forma (eq. (2.72)) para cada elemento B= 1 2 At ( y J − yK ) ( yK − yI ) ( yI − y J )     ( x K − x J ) ( xI − xK ) ( x J − x I )  (2.117) 3) Determinação da matriz de condutividade para cada elemento, através da eq. (2.82) Ke = ∫e B Ω T κB dΩ (2.118) Para o caso particular do elemento triangular linear, com material homogêneo, as matrizes κ e B são constantes (independentes de x e y). Assim, a matriz de condutividade do elemento pode ser obtida como K e = BT κBAt t (2.119) 4) Determinação do vetor de fontes ou fluxos nodais para cada elemento, através da eq. (2.83) Fe = ∫e N Ω T Q dΩ + ∫e N Γq T qn dΓ + Fce (2.120) Para o caso particular do elemento triangular linear, com fonte Q constante, e fluxo normal prescrito no contorno do elemento qn nulo, o vetor F e pode ser obtido como 34 O MEF aplicado ao Problema de Condução de Calor 1  F = 1 QAt t + Fce 1  e Remo M. de Souza (2.121) 5) Determinação da matriz de incidência H e para cada elemento, conforme explicado na seção 3.2.1. 6) Montagem da matriz de condutividade do modelo de acordo com a eq. (2.108) K = ∑ He K e He ne T e =1 (2.122) 7) Montagem do vetor de fontes nodais do modelo de acordo com a eq. (2.101) F = ∑ He Fe ne T e =1 (2.123) Na verdade, apenas se conhece uma parte deste vetor, denominada F0 , em função das condições de contorno. Após a montagem do vetor F total como mostrado acima, extrai-se a parte F0 deste vetor, e ignora-se a parte F1 , a qual será recalculada posteriormente. 8) Montagem da parte conhecida T1 do vetor de temperaturas nodais 9) Partição do sistema de equações KT = F , considerando as condições de contorno, de acordo com a seção 3.3.  K 00 K  10 K 01  T0  F0   =  K11   T1   F1  (2.124) 10) Solução do sistema de equações, de acordo com as eqs. (2.114) e (2.115) T0 = K 00−1 ( F0 − K 01T1 ) F1 = K10T0 + K11T1 (2.125) (2.126) 11) Montagem do vetor de temperaturas nodais do modelo T  T =  0  T1  (2.127) 12) Determinação do vetor de temperaturas nodais de cada elemento, utilizando-se a matriz de incidência. 35 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Te = H e T (2.128) 13) Determinação do gradiente de temperatura no interior de cada elemento (eq. (2.70) ∇T = BTe (2.129) 14) Determinação do fluxo de calor no interior de cada elemento (eq. (2.9)). q = −κ∇T (2.130) Com isto, tem-se a solução do problema de condução de calor por elementos finitos, utilizandose o elemento triangular linear. 36 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza A. Interpretação geométrica das funções de forma do elemento Equation Section 1 Para possibilitar uma intepretação geométrica dos termos da matriz de funções de forma N( x, y ) , deve-se inicialmente, calcular a área do elemento triangular. Isto pode ser feito facilmente calculando-se a norma do produto vetorial entre dois vetores r e s definidos arbitrariamente por duas arestas do elemento, conforme ilustra a Figura A.1. K ( xK , y K ) s I ( xI , y I ) y r J ( xJ , y J ) x Figura A.1 – Vetores definidos pelas arestas do elemento, para determinação da área do triângulo.  rx   xJ − xI      r = ry  =  y J − yI     0    rz   e  s x   xK − xI      s =  s y  =  yK − yI     0    sz   (A.1) A área At do elemento pode ser calculada como At = 1 r×s 2 (A.2) Assim, ˆi r × s = rx ˆj ry sx sy ˆi ˆj kˆ kˆ rz = ( xJ − xI ) ( y J − yI ) 0 s z ( xK − xI ) ( y K − y I ) 0 = 0ˆi + 0ˆj + ( ( xJ − xI )( yK − yI ) − ( xK − xI )( y J − yI ) ) kˆ (A.3) onde î , ĵ e k̂ são os vetores base unitários (versores) do sistema de coordenadas cartesianas ( x, y, z ) conforme indica a Figura A.1. 37 O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza Substituindo-se a eq. (A.3) na eq. (A.2), chega-se a 1 1 r × s = ( ( xJ − xI )( y K − yI ) − ( xK − xI )( y J − yI ) ) 2 2 1 = ( ( x J y K + x I y J + xK y I ) − ( x J y I + xK y J + xI y K ) ) 2 At = (A.4) Comparando as equações (2.63) e (A.4), conclui-se que det(G ) = 2 At (A.5) ou seja, o determinante da matriz G corresponde a duas vezes a área do elemento. Para interpretação geométrica dos outros termos da matriz N( x, y ) considera-se um ponto P de coordenadas ( x, y ) no interior do elemento, e divide-se o elemento em três triângulos com vértice em P, conforme mostra a Figura A.2. K ( xK , y K ) A2 A1 I ( xI , y I ) y P ( x, y ) A3 J ( xJ , y J ) x Figura A.2 – Sub-áreas no interior do elemento, definidas por um ponto P de coordenadas ( x, y ) . Desta forma, o cálculo da área A1 do triângulo definido pelos pontos P, J e K pode ser calculada utilizando-se a eq. (A.4), com as variáveis x e y, no lugar de xI e yI , respectivamente 1 ( ( xJ yK + xy J + xK y ) − ( xJ y + xK yJ + xyK ) ) 2 1 = ( ( x J y K − xK y J ) + ( y J − y K ) x + ( xK − x J ) y ) 2 A1 ( x, y ) = (A.6) De forma similar, tem-se, para as áreas A2 e A3 A2 ( x, y ) = 1 ( ( xK y I − x I y K ) + ( y K − y I ) x + ( xI − xK ) y ) 2 (A.7) 38 O MEF aplicado ao Problema de Condução de Calor A3 ( x, y ) = Remo M. de Souza 1 ( ( xI y J − x J y I ) + ( y I − y J ) x + ( x J − xI ) y ) 2 (A.8) Comparando as eqs. (A.6), (A.7) e (A.8) com os elementos da matriz N( x, y ) na eq. (2.67), e considerando ainda a eq. (A.5), conclui-se que a matriz de funções de forma pode ser expressa como N ( x, y ) = ou ainda como, A1 ( x, y ) At A2 ( x, y ) At A3 ( x, y ) At N( x, y ) = N1 ( x, y ) N 2 ( x, y ) N3 ( x, y ) (A.9) (A.10) com A ( x, y ) N1 ( x, y ) = 1 At A ( x, y ) N 2 ( x, y ) = 2 At A ( x, y ) N 3 ( x, y ) = 3 At (A.11) sendo denominadas coordenadas naturais do triângulo. É interessante observar que a medida que o ponto P “se aproxima” do nó I por exemplo, N1 → 1 , N 2 → 0 e N3 → 0 . Quando o ponto P se situa exatamente sobre o nó I, isto é, quando x = xI e y = yI , tem-se que N1 ( xI , yI ) = 1 , N 2 ( xI , yI ) = 0 , e N3 ( xI , yI ) = 0 . Em geral, esta importante propriedade das funções de forma pode ser expressa como 1 se L = M N L ( xM , y M ) =  0 se L ≠ M (A.12) onde considera-se que a numeração dos nós do elemento é tal que, I ≡ 1 , J ≡ 2 e K ≡ 3 . 39