O Método dos Elementos Finitos Aplicado ao Problema …luciac/fem/livros-fem... · O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 2 Devido ao fato das sub-regiões

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

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.

contorno original

elementos finitospontos nodais

Figura 1.1 – Malha de Elementos Finitos (para problema plano)


2

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 triangularcom três nós

Elemento hexaédricocom oito nós

Elemento quadrilateralcom nove nós

Elemento quadrilateralcom quatro nós

Elemento tetraédricocom quatro nós

Elemento triangularcom seis nós

Elemento de barracom dois nós

Elemento de barracom três nós

Figura 1.2 – Diferentes tipos de elementos finitos


3

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.


4

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;


5

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.


6

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.


7

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.

x

y

Γ

Ω

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 ( , )xq x y e ( , )yq 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.


8

( , )( , )

( , )x

y

q x yx y

q x y = =

q q . (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

xq

yq

dx

dydx

xqq xx

∂+∂

Q

dyyy

qq

y∂

+∂

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 d dQ x y . Se as faces anterior e posterior indicadas na figura forem isoladas

termicamente, então a seguinte condição deve ser satisfeita

d d d d ( + d )d ( + d )dyxx y x y

qqQ x y q y q x q x y q y xx y

∂∂+ + = +

∂ ∂. (2.2)

Cancelando os termos repetidos, e dividindo a equação resultante por d dx y chega-se a equação

que governa o problema estacionário de condução de calor

0yx qq Qx y

∂∂− − + =∂ ∂

em Ω , (2.3)

ou, de forma mais compacta

div 0Q− + =q em Ω , (2.4)

ou, ainda


9

T 0Q−∇ + =q em Ω , (2.5)

onde

x

y

∂ ∂ ∇ = ∂ ∂

(2.6)

denota o operador diferencial nabla (ou del), tal que

T divx yxy

q qqqx y x y

∂ ∂∂ ∂ ∇ = = + = ∂ ∂ ∂ ∂ q q . (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,

x xTqx

κ ∂= −

∂, (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

T= − ∇q κ , (2.9)

onde, para o caso bidimensional,

( , ) ( , )( , )

( , ) ( , )xx xy

xy yy

x y x yx y

x y x y

κ κ

κ κ

= =

κ κ (2.10)

é a matriz de condutividade térmica, e

grad

Tx xT T T

Ty y

∂ ∂ ∂ ∂ ∇ = = = ∂ ∂ ∂ ∂

. (2.11)

Assim, a eq. (2.9) pode-ser escrita como


10

xx xyx

y xy yy

Tq x

Tqy

κ κ

κ κ

∂ ∂ = − ∂ ∂

, (2.12)

ou seja,

xx xyx

xy yyy

T Tq

x yT T

qx y

κ κ

κ κ

∂ ∂ += − ∂ ∂ ∂ ∂ += − ∂ ∂

(2.13)

Substitutindo as eqs. (2.13) em (2.3), chega-se a

0xx xy xy yyT T T T

Qx y x yx y

κ κ κ κ∂ ∂ ∂ ∂∂ ∂ + ++ + = ∂ ∂ ∂ ∂∂ ∂ . (2.14)

Se as direções cartesianas (x,y) coincidirem com as direções principais do material, então

0xyκ = . 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 1 0( , ) ( , ) ( , )

0 ( , ) 0 1x y

x y x y x yx y

κκ κ

κ

= = = =

κ κ I , (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 0xyκ = e ctexx yyκ κ κ= = = . Neste caso, a

eq. (2.14) fica

2 2

2 2 0T T Qx y

κ ∂ ∂ + =+ ∂ ∂

(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


11

T 0T Q∇ ∇ + =κ (2.17)

Considerando-se novamente um meio isotrópico e homogêneo (com κ=κ I constante), esta

equação fica

2 0T Qκ∇ + = (2.18)

onde 2∇ é o operador Laplaciano, tal que

2 22 T

2 2T TxT T T

x y x yy

∂ ∂ ∂ ∂ ∂∂ ∇ = ∇ ∇ = = + ∂∂ ∂ ∂ ∂ ∂

(2.19)

Para o caso particular em que 0Q = , ou seja, sem nenhuma fonte de calor interna, a eq. (2.18),

fica

2 0T∇ = ou 2 2

2 2 0T Tx y

∂ ∂+ =

∂ ∂(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)


12

x

y

Ω

ΓT

Γq

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

x

y

Ω

ΓT

Γq

n

xq

yq

coss α∆ Q

nq

n α

α

s∆

sens α∆

(a) (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;


13

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,

ˆ cosˆ

ˆ senx

y

nn

αα

= =

n (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 nq é 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

21 cos sen cos sen 02 x y nQ s q s q s q sα α α α∆ + ∆ + ∆ + ∆ = . (2.24)

Dividindo a equação por s∆ , tem-se

1 cos sen cos sen 02 x y nQ s q q qα α α α∆ + + + = . (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

ˆ ˆx x y y nq n q n q− − = , (2.26)

ou de forma mais compacta,

T ˆ nq− =q n 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 ( ) 0 emT Q−∇ + = Ωq (2.28)


14

Relação constitutiva do meio (Lei de Fourier)

( ) emT T= − ∇ Ωq κ (2.29)

Condições de contorno

T

em

ˆ emT

n q

T T

q

= Γ

− = Γ q n(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

2 2, , , ,2, , ,x y xx xyx y y xy

ϕ ϕ ϕ ϕϕ ϕ ϕ ϕ∂ ∂ ∂ ∂= = = =∂ ∂ ∂ ∂∂

(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

( )T( , ) 0( )w x y T Q =−∇ +q (2.32)


15

ou seja,

( ), ,( , ) 0x x y yq q Qw x y + − = (2.33)

2o Passo) Integra-se a equação acima sobre o domínio Ω

( ), ,( , ) d 0x x y yq q Qw x yΩ

+ − Ω =∫ (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 xq e yq para a função teste w , utilizando-

se a regra de derivada do produto

, , ,

, , ,

( )

( )x x x x x x

y y y y y y

wq wq w q

wq wq w q

= −

= −(2.35)

tal que

( )

( )

, ,

, , , ,

d

( ) ( ) d

x x y y

x x y y x x y y

q q Qw

wq wq w q w q wQΩ

Ω

+ − Ω =

+ − − − Ω

∫

∫(2.36)

Utiliza-se, então, o teorema integral de Gauss, tal que

( ) ( )

( )

, ,

T

ˆ ˆ( ) ( ) ( ) ( )d d

ˆ ˆ d

ˆ d

x x y y x x y y

x x y y

wq wq wq n wq n

q n q nw

w

Ω Γ

Γ

Γ

+ +Ω = Γ

+= Γ

= Γ

∫ ∫

∫

∫ q n

(2.37)


16

Assim, substitutindo-se a eq. (2.37) em (2.36), chega-se a

( )

( )

, ,

T, ,

0 d

ˆd d

x x y y

x x y y

q q Qw

w q w q wQ wΩ

Ω Γ

+ −= Ω

− − −= Ω+ Γ

∫

∫ ∫ q n(2.38)

4o Passo) Considera-se as condições de contorno (2.30), tal que

T T T

T

ˆ ˆ ˆd d d

ˆ d d

T q

T q

n

w w w

w wq

Γ Γ Γ

Γ Γ

Γ = Γ + Γ

= Γ − Γ

∫ ∫ ∫

∫ ∫

q n q n q n

q n(2.39)

Observa-se que o termo T ˆq n é igual a nq em qΓ , sendo portanto perfeitamente conhecido nesta

região do contorno. Porém, o termo T ˆq 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

( , ) 0w x y = em TΓ (2.40)

tal que o termo

T ˆ d 0T

w SΓ

=∫ q n (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

( ), , d d 0q

x x y y nw q w q wQ wqΩ Γ

− − − Ω− Γ =∫ ∫ (2.42)

ou seja, de forma mais compacta

( )T d d dq

nwQ wqwΩ Ω Γ

− Ω = Ω+ Γ∇∫ ∫ ∫q (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.


17

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 ( Tw w= ), a eq. (2.43) pode ser escrita como

( ) T TT ( )d d dq

nT w Q w qwΩ Ω Γ

∇ Ω = Ω+ Γ∇∫ ∫ ∫κ (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.


18

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.

29

12

10

11

8

7

4

3

65

1

x

y1 1( , )x y

2 2( , )x y

3 3( , )x y

4 4( , )x y

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


19

1

2

3

gN

TTT

T

=

T , (2.45)

onde 1T é a temperatura correspondente ao gdl 1, 2T é a temperatura correspondente ao gdl 2, e assim

por diante, até o número de graus de liberdade gN 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 g gN N× , e F é um vetor de coeficientes (em geral denominado vetor de forças), de ordem

1gN × , 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

1

2

3

gN

FFF

F

=

F (2.47)

onde 1F é a fonte de calor correspondente ao gdl 1, 2F é a fonte correspondente ao gdl 2, e assim por

diante, até o número de graus de liberdade gN 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))

( , ) ( , )d ( , )P P c P cP Pw x y Q x y w x y F w FΩ

Ω = =∫ (2.48)


20

onde PcF corresponde a fonte de calor concentrada no ponto P, e Pw 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 TT

1

T T T

( )d d d

d d

P

q

q

N

n P cPP

n c

T w Q w q w Fw

w Q w q

=Ω Ω Γ

Ω Γ

∇ Ω = Ω + Γ +∇

= Ω + Γ +

∑∫ ∫ ∫

∫ ∫

κ

w F(2.49)

onde w é um vetor tal que

1 1 1

2 2 2

3 3 3

( , )( , )( , )

( , )g g g

w x y ww x y ww x y w

w x y w

= =

w (2.50)

e

1

2

3

g

c

c

cc

cN

FFF

F

=

F (2.51)

é 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 cF é 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


21

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.

x

y( , )I II x y

( , )K KK x y

( , )J JJ x y

Figura 3.2 – Elemento finito triangular linear, com referência ao sistema de eixos cartesianos.

Na Figura 3.2 estão indicadas as coordenadas ( , )I Ix y , ( , )J Jx y e ( , )K Kx y , 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 IT , JT , e KT .

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 eT do elemento

Ie

J

K

TTT

=

T (2.52)


22

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

1 2 3( , )T x y a a x a y= + + (2.53)

onde 1a , 2a e 3a são constantes a serem determinadas. A eq. (2.53) pode ser escrita de forma mais

compacta como

( , ) ( , )T x y x y= x a (2.54)

onde1 x y=x (2.55)

e

1

2

3

aaa

=

a (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

1 2 3

1 2 3

1 2 3

( , )( , )( , )

I I I I I

J J J J J

K K K K K

T x y a a x a y TT x y a a x a y TT x y a a x a y T

= + + == + + == + + =

(2.57)

As eqs. (2.57) podem ser reescritas na forma matricial como

1

2

3

111

II I

JJ J

KK K

a Tx ya Tx ya Tx y

=

(2.58)

ou, em forma, mais compacta

e=Ga T (2.59)

onde


23

111

I I

J J

K K

x yx yx y

=

G (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)

1 e−=a G T (2.61)

onde

1( ) ( ) ( )

1 ( ) ( ) ( )det( )

( ) ( ) ( )

J K K J K I I K I J J I

J K K I I J

K J I K J I

x y x y x y x y x y x yy y y y y yx x x x x x

−− − −

= − − − − − −

GG

(2.62)

com

det( ) ( ) ( )J K I J K I J I K J I Kx y x y x y x y x y x y= + + − + +G (2.63)

Conforme apresentado no Apêndice A, o determinante da matriz G corresponde à duas vezes a

área tA do elemento, ou seja

2 det( )tA = G (2.64)

Substituindo a eq. (2.61) na eq. (2.54), chega-se a

1( , ) ( , ) eT x y x y −= x G T (2.65)

ou ainda,

( , ) ( , ) eT x y x y= N T (2.66)

onde1( , ) ( , )x y x y −=N x G (2.67)

Considerando a eq. (2.64), e desenvolvendo o produto dado na equação acima, chega-se a

1 2 3( , ) ( , ) ( , ) ( , )x y N x y N x y N x y=N (2.68)

onde


24

( )

( )

( )

1

2

3

1( , ) ( ) ( ) ( )2

1( , ) ( ) ( ) ( )21( , ) ( ) ( ) ( )

2

J K K J J K K Jt

K I I K K I I Kt

I J J I I J J It

N x y x y x y y y x x x yA



= − + − + −

= − + − + −

= − + − + −

(2.69)

A matriz ( , )x yN é 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

eT .

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)

( , ) ( , ) ( , )e eT x y x y x y∇ =∇ =N T B T (2.70)

onde

1 2 3

31 2

31 2

( , ) ( , ) ( , ) ( , ) ( , )

( , )( , ) ( , )

( , )( , ) ( , )

xx y x y N x y N x y N x y

yN x yN x y N x y

x x xN x yN x y N x y

y y y

∂ ∂ = ∇ = ∂ ∂

∂∂ ∂ ∂ ∂ ∂ =

∂∂ ∂ ∂ ∂ ∂

B N

(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


25

( ) ( ) ( )1( ) ( ) ( )2

J K K I I J

K J I K J It

y y y y y yx x x x x xA

− − − = − − −

B (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

Ie

J

K

www

=

w (2.73)

utilizando-se a mesma matriz de funções de forma ( , )x yN (ver eq. (2.66)). Ou seja,

( , ) ( , ) ew x y x y= N w (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

( , ) ( , ) ew x y x y∇ = B w (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,

( )TT TT ( )d d d

e e eq

e en cT w Q w qw

Ω Ω Γ

∇ Ω = Ω+ Γ +∇∫ ∫ ∫κ w F (2.76)

onde eΩ , e eΓ correspondem, respectivamente, ao domínio e contorno de um elemento, e

ec I

e ec c J

ec K

F

F

F

=

F (2.77)


26

é o vetor de fluxos nodais do elemento.

A substituição das eqs. (2.70), (2.75), e (2.74) em (2.76), leva a

T T T TT T Td d de e e

q

e e e e e en cQ q

Ω Ω Γ

Ω = Ω+ Γ +∫ ∫ ∫w B κB T w N w N w F (2.78)

Colocando-se o termo Tew em evidência, chega-se a

T T T Td d d 0e e e

q

e e en cQ q

Ω Ω Γ

Ω − Ω− Γ − = ∫ ∫ ∫w B κB T N N F (2.79)

O vetor ew deve ser completamente arbitrário; portanto, conclui-se que

T T Td d de e e

q

e en cQ q

Ω Ω Γ

Ω = Ω+ Γ +∫ ∫ ∫B κB T N N F (2.80)

ou, de forma mais compacta,

e e e=K T F (2.81)

onde

T de

e

Ω

= Ω∫K B κB (2.82)

é matriz de condutividade (ou de rigidez) do elemento, e

T Td de e

q

e en cQ q

Ω Γ

= Ω+ Γ +∫ ∫F N N F (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,

T Td de e

q

eq nQ q

Ω Γ

= Ω+ Γ∫ ∫F N N (2.84)

tal que

e e eq c= +F F F (2.85)


27

3.2. 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.

5

2 3

4

1

a

b c

d

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


28

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

1

2

3

4

5

TTTTT

=

T (2.86)

O vetor de temperaturas nodais do elemento a é

2

3

1

aI

a aJaK

T TT T

TT

= =

T (2.87)

Observa-se que a equação acima pode ser escrita como o seguinte produto matricial

1

2

3

4

5

0 1 0 0 00 0 1 0 01 0 0 0 0

aI

a aJaK

TT TT T

TTT

= =

T (2.88)

ou ainda,a a=T H T (2.89)

onde

0 1 0 0 00 0 1 0 01 0 0 0 0

a =

H (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

b b

c c

d d

=

=

=

T H T

T H T

T H T

(2.91)

onde


29

0 1 0 0 01 0 0 0 00 0 0 1 0

0 0 1 0 00 0 0 0 11 0 0 0 0

1 0 0 0 00 0 0 0 10 0 0 1 0

b

c

d

= = =

H

H

H

(2.92)

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,


3) Para a terceira linha, tem-se o valor um na coluna correspondente ao nó K do elemento,


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


30

1

2

3

4

5

6

FFFFFF

=

F (2.93)

O vetor de fontes nodais do elemento a é

aI

a aJaK

F

F

F

=

F (2.94)

Analogamente, para os outros elementos, têm-se

b c dI I I

b b c c d dJ J Jb c dK K K

F F F

F F F

F F F

= = =

F F F (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

1a b c dK J K IF F F F F= + + + (2.96)

Analogamente, têm-se para os demais nós

2

3

4

5

a bI Ia cJ Ib dK Kc dJ J

F F F

F F F

F F F

F F F

= +

= +

= +

= +

(2.97)

As eqs. (2.96) e (2.97) podem ser escritas na seguinte forma matricial

1

2

3

4

5

0 0 1 0 1 0 0 0 1 1 0 01 0 0 1 0 0 0 0 0 0 0 00 1 0 0 0 0 1 0 00 0 0 0 0 1 0 0 00 0 0 0 0 0 0 1 0

a b cI I Ia b cJ J Ja b cK K K

FF F FFF F FF

F F F FF

= + + +

0 0 00 0 10 1 0

dIdJdK

F

F

F

(2.98)


31

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 Ta a b b c c d d= + + +F H F H F H F H F (2.99)

Para escrever a equação acima na forma de somatório, pode-se definir 1a ≡ , 2b ≡ , 3c ≡ e

4d ≡ . Assim,

T T T T1 1 2 2 3 3 4 4= + + +F H F H F H F H F (2.100)

ou seja,

T

1

ene e

e== ∑F H F (2.101)

onde e representa um elemento genérico e en é o número de elementos da malha (neste caso, 4en = ).

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

a a a b b b c c c d d d= = = =F K T F K T F K T F K T (2.102)

Substituindo as equações acima na eq. (2.99), tem-se

T T T Ta a a b b b c c c d d d= + + +F H K T H K T H K T H K T (2.103)

Substituindo, agora, as eqs. (2.89) e (2.91) na equação acima, tem-se

T T T Ta a a b b b c c c d d d= + + +F H K H T H K H T H K H T H K H T (2.104)

Colocando o vetor de temperaturas nodais do modelo T em evidência, chega-se a

T T T T( )a a a b b b c c c d d d= + + +F H K H H K H H K H H K H T (2.105)

ou ainda,

=F KT (2.106)

onde


32

T T T Ta a a b b b c c c d d d= + + +K H K H H K H H K H H K H (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

T

1

ene e e

e== ∑K H K H (2.108)

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


33

11 12 13 14 15 1 1

21 22 23 24 25 2 2

31 32 33 34 35 3 3

41 42 43 44 45 4 4

51 52 53 54 55 5 5

K K K K K T FK K K K K T FK K K K K T FK K K K K T FK K K K K T F

=

(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,

11 12 13 14 15 1 1

21 22 23 24 25 2 2

31 32 33 34 35 3 3

41 42 43 44 45 4 4

51 52 53 54 55 5 5

K K K K K T FK K K K K T FK K K K K T F

K K K K K T FK K K K K T F

=

(2.110)

ou de forma mais compacta

00 01 0 0

10 11 1 1

=

K K T FK K T F

(2.111)

onde

11 12 13 14 15 1 1

00 21 22 23 01 24 25 0 2 0 2

31 32 33 34 35 3 3

41 42 43 44 45 4 410 11 1 1

51 52 53 54 55 5 5

K K K K K T FK K K K K T FK K K K K T F

K K K K K T FK K K K K T F

= = = =

= = = =

K K T F

K K T F

(2.112)

Desta forma, deve-se notar que os vetores 1T e 0F são conhecidos, ao passo que os vetores 0T

e 1F são desconhecidos.

Desenvolvendo a eq. (2.111), tem-se

00 0 01 1 0

10 0 11 1 1

+ =+ =

K T K T FK T K T F

(2.113)

Com isso, o vetor de temperaturas nodais 0T pode ser calculado, a partir da primeira das

equações acima,

( )10 00 0 01 1

−= −T K F K T (2.114)


34

Após a determinação de 0T , o vetor de fontes nodais 1F , pode ser calculado diretamente

utilizando-se a segunda das eqs. (2.113)

1 10 0 11 1= +F K T K T (2.115)

3.4. 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 xy

xy yy

x y x yx y

x y x y

κ κ

κ κ

= =

κ κ (2.116)

2) Montagem da matriz com as derivadas das funções de forma (eq. (2.72)) para cada elemento

( ) ( ) ( )1( ) ( ) ( )2

J K K I I J

K J I K J It

y y y y y yx x x x x xA

− − − = − − −

B (2.117)

3) Determinação da matriz de condutividade para cada elemento, através da eq. (2.82)

T de

e

Ω

= Ω∫K B κB (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

TetA t=K B κB (2.119)

4) Determinação do vetor de fontes ou fluxos nodais para cada elemento, através da eq. (2.83)

T Td de e

q

e en cQ q

Ω Γ

= Ω+ Γ +∫ ∫F N N F (2.120)

Para o caso particular do elemento triangular linear, com fonte Q constante, e fluxo normal

prescrito no contorno do elemento nq nulo, o vetor eF pode ser obtido como


35

111

e et cQA t

= +

F F (2.121)

5) Determinação da matriz de incidência eH 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)

T

1

ene e e

e== ∑K H K H (2.122)

7) Montagem do vetor de fontes nodais do modelo de acordo com a eq. (2.101)

T

1

ene e

e== ∑F H F (2.123)

Na verdade, apenas se conhece uma parte deste vetor, denominada 0F , em função das condições

de contorno. Após a montagem do vetor F total como mostrado acima, extrai-se a parte 0F deste vetor,

e ignora-se a parte 1F , a qual será recalculada posteriormente.

8) Montagem da parte conhecida 1T 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.

00 01 0 0

10 11 1 1

=

K K T FK K T F

(2.124)

10) Solução do sistema de equações, de acordo com as eqs. (2.114) e (2.115)

( )10 00 0 01 1

−= −T K F K T (2.125)

1 10 0 11 1= +F K T K T (2.126)

11) Montagem do vetor de temperaturas nodais do modelo

0

1

=

TT

T(2.127)

12) Determinação do vetor de temperaturas nodais de cada elemento, utilizando-se a matriz de

incidência.


36

e e=T H T (2.128)

13) Determinação do gradiente de temperatura no interior de cada elemento (eq. (2.70)

eT∇ = BT (2.129)

14) Determinação do fluxo de calor no interior de cada elemento (eq. (2.9)).

T= − ∇q κ (2.130)

Com isto, tem-se a solução do problema de condução de calor por elementos finitos, utilizando-

se o elemento triangular linear.


37

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

( , )x yN , 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.

x

y( , )I II x y

( , )K KK x y

( , )J JJ x yr

s

Figura A.1 – Vetores definidos pelas arestas do elemento, para determinação da área do triângulo.

e0 0

x J I x K I

y J I y K I

z z

r x x s x xr y y s y y

r s

− − = = − = = −

r s (A.1)

A área tA do elemento pode ser calculada como

12tA = ×r s (A.2)

Assim,

( )

ˆ ˆ ˆ ˆ ˆ ˆ

( ) ( ) 0( ) ( ) 0

ˆ ˆ ˆ0 0 ( )( ) ( )( )

x y z J I J I

K I K Ix y z

J I K I K I J I

r r r x x y yx x y ys s s

x x y y x x y y

× = = − −− −

= + + − − − − −

i j k i j kr s

i j k

(A.3)

onde i , j e k são os vetores base unitários (versores) do sistema de coordenadas cartesianas ( , , )x y z

conforme indica a Figura A.1.


38

Substituindo-se a eq. (A.3) na eq. (A.2), chega-se a

( )

( )

1 1 ( )( ) ( )( )2 21 ( ) ( )2

t J I K I K I J I

J K I J K I J I K J I K

A x x y y x x y y

x y x y x y x y x y x y

= = − − − − −×

= + + − + +

r s(A.4)

Comparando as equações (2.63) e (A.4), conclui-se que

det( ) 2 tA=G (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 ( , )x yN 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.

x

y( , )I II x y

( , )K KK x y

( , )J JJ x y

( , )P x y2A

3A

1A

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 1A 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 Ix e Iy , respectivamente

( )

( )

11( , ) ( ) ( )21 ( ) ( ) ( )2

J K J K J K J K

J K K J J K K J

A x y x y xy x y x y x y xy

x y x y y y x x x y

= + + − + +

= − + − + −(A.6)

De forma similar, tem-se, para as áreas 2A e 3A

( )21( , ) ( ) ( ) ( )2 K I I K K I I KA x y x y x y y y x x x y= − + − + − (A.7)


39

( )31( , ) ( ) ( ) ( )2 I J J I I J J IA x y x y x y y y x x x y= − + − + − (A.8)

Comparando as eqs. (A.6), (A.7) e (A.8) com os elementos da matriz ( , )x yN 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

31 2 ( , )( , ) ( , )( , )t t t

A x yA x y A x yx yA A A

=N (A.9)

ou ainda como,

1 2 3( , ) ( , ) ( , ) ( , )x y N x y N x y N x y=N (A.10)

com

31 21 2 3

( , )( , ) ( , )( , ) ( , ) ( , )t t t

A x yA x y A x yN x y N x y N x yA A A

= = = (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, 1 1N → , 2 0N → e 3 0N → . Quando o ponto P se situa

exatamente sobre o nó I, isto é, quando Ix x= e Iy y= , tem-se que 1( , ) 1I IN x y = , 2 ( , ) 0I IN x y = , e

3( , ) 0I IN x y = . Em geral, esta importante propriedade das funções de forma pode ser expressa como

1se ( , )

0 se L M ML M

N x yL M=

= ≠(A.12)

onde considera-se que a numeração dos nós do elemento é tal que, 1I ≡ , 2J ≡ e 3K ≡ .

Documents

O Método dos Elementos Finitos Aplicado ao Problema …luciac/fem/livros-fem... · O MEF aplicado ao Problema de Condução de Calor Remo M. de Souza 2 Devido ao fato das sub-regiões