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
dy
∂q y
∂y
dy
qx +
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
−
∂q x ∂q y
−
+ Q = 0 em Ω ,
∂x
∂y
(2.3)
ou, de forma mais compacta
−divq + Q = 0 em Ω ,
(2.4)
ou, ainda
8
O MEF aplicado ao Problema de Condução de Calor
Remo M. de Souza
−∇T q + Q = 0 em Ω ,
(2.5)
∂
 ∂x 
∇= 
∂
 ∂y 
(2.6)
onde
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
q = −κ∇T ,
(2.9)
κ xx ( x, y ) κ xy ( x, y ) 
κ = κ ( x, y ) = 

κ xy ( x, y ) κ yy ( x, y ) 
(2.10)
onde, para o caso bidimensional,
é 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
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)
ou seja,
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
∇ T κ ∇T + Q = 0
Remo M. de Souza
(2.17)
Considerando-se novamente um meio isotrópico e homogêneo (com κ = κ I constante), esta
equação fica
κ∇ 2T + Q = 0
(2.18)
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
em ΓT
(2.22)
onde T é a temperatura conhecida no contorno ΓT .
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̂
∆s
Γq
Ω
y
qx
∆s cos α
qn
α
n̂ α
Q
∆s senα
x
ΓT
(a)
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)
Dividindo a equação por ∆s , tem-se
1
Q ∆s cos α senα + qx cos α + q y senα + qn = 0 .
2
(2.25)
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
Remo M. de Souza
Relação constitutiva do meio (Lei de Fourier)
q(T ) = −κ∇T
em Ω
(2.29)
em ΓT 

em Γ q 
(2.30)
Condições de contorno
T =T
−qT nˆ = qn
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
Remo M. de Souza
ou seja,
w( x, y ) ( qx, x + q y, y − Q ) = 0
(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)
tal que
∫ w ( qx, x + q y, y − Q ) dΩ =
Ω
∫ ( (wqx ), x + (wq y ), y − w, x qx − w, y q y − wQ ) dΩ
(2.36)
Ω
Utiliza-se, então, o teorema integral de Gauss, tal que
∫ ( (wqx ), x + (wq y ), y ) dΩ = ∫ ( (wqx )nˆ x + (wq y )nˆ y ) dΓ
Ω
Γ
= ∫ w ( qx nˆ x + q y nˆ y ) dΓ
(2.37)
Γ
= ∫ wqT nˆ dΓ
Γ
15
O MEF aplicado ao Problema de Condução de Calor
Remo M. de Souza
Assim, substitutindo-se a eq. (2.37) em (2.36), chega-se a
0=
∫ w ( qx, x + q y, y − Q ) dΩ
Ω
=
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
Tˆ
n dΓ =
Γ
∫
wqT nˆ dΓ +
∫
wqT nˆ dΓ −
ΓT
=
∫ wq
Tˆ
n dΓ
Γq
ΓT
∫
wqn dΓ
(2.39)
Γq
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
em ΓT
(2.40)
tal que o termo
∫
wqT nˆ dS = 0
(2.41)
ΓT
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
Ω
(2.42)
Γq
ou seja, de forma mais compacta
− ∫ ( ∇w )T q dΩ =
Ω
∫ wQ dΩ + ∫ wqn dΓ
Ω
(2.43)
Γq
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Ω =
T
∫w
Ω
Q dΩ +
T
∫w
qn dΓ
(2.44)
Γq
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
( x3 , y3 )
3
10
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Ω +
Ω
Ω
=
T
∫w
Ω
Q dΩ +
∫
wT qn dΓ +
Γq
T
∫w
NP
∑ wP Fc P
P =1
(2.49)
T
qn dΓ + w Fc
Γq
onde w é um vetor tal que
 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 )
y
I ( xI , y I )
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
(2.54)
onde
x= 1 x
y
(2.55)
e
 a1 
 
a = a2 
 
 a3 
(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)
Ga = Te
(2.59)
ou, em forma, mais compacta
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
(2.61)
 ( 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 ) 
(2.62)
det(G ) = ( xJ yK + xI y J + xK yI ) − ( xJ yI + xK y J + xI yK )
(2.63)
onde
com
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
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)
ou ainda,
onde
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 ) =
1
( ( x J y K − xK y J ) + ( y J − y K ) x + ( xK − x J ) y )
2 At
N 2 ( x, y ) =
1
( ( x K y I − xI y K ) + ( y K − y I ) x + ( xI − xK ) y )
2 At
N 3 ( x, y ) =
1
( ( xI y J − x J y I ) + ( y I − y J ) x + ( x J − xI ) y )
2 At
Remo M. de Souza
(2.69)
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)
∇T ( x , y ) = ∇N ( x , y )T e = B ( x, y ) T e
(2.70)
onde
∂
 ∂x 
B( x, y ) = ∇N( x, y ) =   N1 ( x, y ) N 2 ( x, y ) N3 ( x, y )
∂
 ∂y 
 ∂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
(2.71)
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Ω =
T
∫e w
Q dΩ +
Ω
T
∫e w
T
qn dΓ + w e Fce
(2.76)
Γq
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
T
Ω
∫e N
T
Q dΩ + w e
T
Ω
∫e N
T
T
qn dΓ + w e Fce
(2.78)
Γq
T
Colocando-se o termo w e em evidência, chega-se a
T
w e  ∫ BT κB dΩ Te − ∫ N T Q dΩ − ∫ N T qn dΓ − Fce  = 0
 e

e
e
Ω

Ω
Γ
q


(2.79)
O vetor w e deve ser completamente arbitrário; portanto, conclui-se que
∫e B
T
κB dΩ Te =
∫e N
Ω
T
∫e N
Q dΩ +
Ω
T
qn dΓ + Fce
(2.80)
Γq
ou, de forma mais compacta,
K e Te = F e
(2.81)
onde
Ke =
∫e B
T
κB dΩ
(2.82)
Ω
é matriz de condutividade (ou de rigidez) do elemento, e
Fe =
∫e N
T
Q dΩ +
Ω
∫e N
T
qn dΓ + Fce
(2.83)
Γq
é 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
T
qn dΓ
(2.84)
Γq
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 
(2.88)
Ta = H a T
(2.89)
0 1 0 0 0
H =  0 0 1 0 0 
1 0 0 0 0 
(2.90)
ou ainda,
onde
a
é 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
(2.91)
Td = H d T
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
 F 
F =  3
 F4 
 F5 
 
 F6 
(2.93)
O vetor de fontes nodais do elemento a é
 FIa 
 
a  a
F =  FJ 
 a
F
 K 
(2.94)
Analogamente, para os outros elementos, têm-se
 FIb 
 
b  b
F =  FJ 
 b
 FK 
 FIc 
 
c  c
F =  FJ 
 c
 FK 
 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
(2.97)
F4 = FKb + FKd
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
T
T
T
T
F = H a F a + H b Fb + H c F c + H d F d
(2.99)
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
T
T
F = H1 F1 + H 2 F 2 + H3 F3 + H 4 F 4
(2.100)
ou seja,
ne
T
F = ∑ He Fe
(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
T
T
T
T
F = H a K a Ta + Hb K b Tb + H c K c Tc + H d K d Td
(2.103)
Substituindo, agora, as eqs. (2.89) e (2.91) na equação acima, tem-se
T
T
T
T
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
(2.104)
Colocando o vetor de temperaturas nodais do modelo T em evidência, chega-se a
T
T
T
T
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)
ou ainda,
onde
31
O MEF aplicado ao Problema de Condução de Calor
T
T
T
Remo M. de Souza
T
K = H a K a H a + Hb K b Hb + H c K c H c + H d K d H d
(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
ne
T
K = ∑ He K e He
(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
 K
  51
K12
K 22
K32
K13 
K 23 
K33 
 K14
K
 24
 K34
K 42
K52
K 43   K 44
K53   K54
K15    T1   F1 
     
K 25    T2   F2 
 =  F 
K35   T3 
 3  


K 45   T4   F4 
 
 
K55   T5   F5 
(2.110)
ou de forma mais compacta
 K 00
K
 10
K 01  T0  F0 
 = 
K11   T1   F1 
(2.111)
onde
 K11

K 00 =  K 21
 K31
K12
K 22
K32
K
K10 =  41
 K51
K 42
K52
K13 
 K14


K 23  K 01 =  K 24
 K34
K33 
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.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)
F1 = K10T0 + K11T1
3.4.
(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 − x ) (x − x ) (x − x ) 
J
I
K
J
I 
 K
(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
T
qn dΓ + Fce
(2.120)
Γq
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)
ne
T
K = ∑ He K e He
(2.122)
e =1
7) Montagem do vetor de fontes nodais do modelo de acordo com a eq. (2.101)
ne
T
F = ∑ He Fe
(2.123)
e =1
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 )
(2.125)
F1 = K10T0 + K11T1
(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
(A.3)
= 0ˆi + 0ˆj + ( ( xJ − xI )( yK − yI ) − ( xK − xI )( y J − yI ) ) kˆ
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
A3
P ( x, y )
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
A3 ( x, y )
At
(A.9)
N( x, y ) = N1 ( x, y ) N 2 ( x, y ) N3 ( x, y )
(A.10)
N ( x, y ) =
A1 ( x, y )
At
A2 ( x, y )
At
ou ainda como,
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