Apostila de Elementos Finitos

Prévia do material em texto

1
Capítulo III 
 
 
 
3.2. Método das Diferenças Finitas 
 
 
Este item 3.2 foi extraído da referência [55] Carnahan 1969, sofrendo modificações. 
Em todo o desenvolvimento da análise numérica, utiliza-se as diferenças finitas, com um 
roteiro a ser seguido numa seqüência lógica. Partindo da definição de operadores 
numéricos de diferenças finitas (∆, ∇, µδ e δ), introduz-se o conceito de interpolação 
através das fórmulas de Gregory-Newton e Stirling, que utilizam estes operadores. Em 
seguida, introduz-se a derivação numérica e a integração numérica (Quadratura) por meio 
da derivação e integração da fórmula de Gregory-Newton, chegando às fórmulas de 
Newton-Cotes. 
 
Seguindo essa técnica, praticamente todos os tópicos da análise numérica podem 
ser introduzidos por meio das diferenças finitas, e depois, eles são desenvolvidos para 
além das diferenças finitas. No estudo de análise numérica de equações diferenciais não 
é diferente. Sugere-se uma introdução por meio dos métodos das diferenças finitas (do 
inglês: Finite Difference Methods ou FDM), e posteriormente, o desenvolvimento do 
assunto para além das diferenças finitas, como por exemplo, a introdução do Método dos 
Elementos Finitos. Com isso, segue-se um desenvolvimento didático muito próximo do 
desenvolvimento histórico, uma vez que, conforme foi citado no item 1.2 da Introdução, o 
método das diferenças finitas surgiu antes do Método dos Elementos Finitos. 
 
O método das diferenças finitas pode ser utilizado para resolver problemas de valor 
de contorno ou valor inicial, envolvendo equações diferenciais ordinárias ou parciais. 
Assim, este método pode ser usado para solucionar as equações de modelos a 
parâmetros concentrados ou distribuídos. A técnica consiste em substituir cada derivada 
ou diferencial das equações diferenciais por aproximação de diferenças finitas ou 
acréscimo finitos das variáveis, como mostra as equação 3.1 abaixo: 
 
�
�
�
�
�
�
�
�
�
�
�
∆∆
∆
≈
∂∂
∂
∆
∆
≈
∂
∂
∆
∆
≈
∂
∂
∆
∆
≈
∂
∂
∆
∆
≈
∆
∆
≈
∆
∆
≈
∆≈∆≈
yx
u
yx
u
x
u
x
u
x
u
x
u
x
u
x
u
x
y
dx
yd
x
y
dx
yd
x
y
dx
dy
ydyxdx
.
,,,
,,
,
22
3
3
3
3
2
2
2
2
3
3
3
3
2
2
2
2
 (3.1) 
 
O Método de Elementos Finitos é bem mais recente que o anterior, sendo mais 
genérico, e podendo ser aplicado a complexas estruturas geométricas e a ambientes com 
várias mudanças de meio. Ele possui uma formulação matemática mais trabalhada, 
sendo, portanto, um conjunto de técnicas e métodos que se baseia na discretização do 
problema em elementos pequenos e na aproximação de cada elemento por um conjunto 
de polinômios. 
 
Considere, primeiramente, o problema formado por equações diferenciais 
ordinárias (EDO’s). Existem dois tipos. Um deles, é o problema de valor inicial, que 
 2
assume a forma geral abaixo: 
 
F[ t, y(t), y´(t) ] =0, t>0 
t=t0, y=yo, y´=y´0 (3.2) 
 
Onde, o t é a variável independente, usualmente o tempo; y é um vetor de variáveis 
dependentes; y´ é a sua derivada em relação a t; F é um vetor de funções de t, y, e y´; e, 
finalmente, yo e y´o são vetores que representam as condições iniciais do problema. Note-
se, que o domínio da variável t é semi-infinito, e que a solução deste problema deverá ser 
obtida marchando-se no tempo, a partir da condição inicial. Caso exista, pelo menos, uma 
função dentro do vetor F que não dependa de nenhum elemento do vetor y´, a equação 
representa um sistema de equações algébrico-diferenciais (sistema de EAD). [68] 
 
O outro tipo de problema é o do valor de contorno, que assume a seguinte forma 
geral para sistemas de segunda ordem: 
 
F[ x, y(x), y´(x), y´´(x) ] = 0, xo<x<xf 
x=xo, go(x,y,y´)=0 (3.3) 
x=xf, gf(x,y,y´)=0 
 
Onde, o x é a variável independente, usualmente, uma coordenada espacial; y é o 
vetor de variáveis dependentes; y´ e y´´ são as suas derivadas: primeira e segunda, 
respectivamente, em relação a x; F é um vetor de funções; e go e gf são vetores de 
funções que representam as condições de contorno nos limites do domínio do sistema de 
equações. 
 
O objetivo do Método das Diferenças Finitas é transformar um problema composto 
por equações diferenciais em um problema formado por equações algébricas. O primeiro 
passo, nesta direção, é a chamada discretização do domínio da variável independente. A 
discretização consiste em dividir o domínio de cálculo em um certo número de 
subdomínios. Para um domínio semi-infinito, existem infinitos subdomínios. Quando o 
domínio é finito, o número de subdomínios também o é, e digamos que seja J. Em 
qualquer caso, estipulam-se os pontos que delimitam os subdomínios, que, no caso de 
um domínio finito, são iguais a (J+1), em número. 
 
Note-se, que os subdomínios podem ter o mesmo tamanho, gerando uma malha 
uniforme, ou então, formando uma malha não-uniforme. Embora as discretizações 
baseadas no primeiro tipo de malha sejam mais simples, existem vantagens numéricas, 
em muitos casos, no uso de malhas não-uniformes. 
 
O segundo passo é gerar aproximações para as derivadas das variáveis 
dependentes que aparecem nas diferenciais, nos pontos discretos, xj ou tj, isto é, obter y´j 
e y´´j, utilizando apenas os valores de y nestes pontos discretos: yj. Finalmente, aplicam-
se as equações diferenciais ordinárias aos pontos discretos xj, substituindo as 
aproximações obtidas para y´j e y´´j. Isto gera sistemas de equações algébricas na forma: 
 
f( yj ) =0, (3.4) 
 
Onde, o f é um vetor de equações algébricas que depende dos valores 
desconhecidos yj, sendo que esta dependência varia conforme o tipo de problema, de 
contorno ou inicial. Este sistema de equações, quer seja ele linear ou não linear, pode ter 
 3
a sua solução obtida. Note-se, que a solução assim obtida para o problema consistirá em 
uma seqüência de pontos, xj ou tj, onde se conhecem os valores de y, yj. 
 
Ficam claras, agora, duas características do Método de Diferenças Finitas: a 
aplicação das equações diferenciais é local, isto é, em cada ponto, xj ou tj, e a solução 
obtida é composta por um conjunto enumerável de pontos onde os valores da solução são 
conhecidos. 
 
Um dos passo necessários na solução de equações diferenciais por diferenças 
finitas é a aproximação das derivadas presentes nestas equações, aplicadas a um dado 
ponto arbitrário, xj ou tj. Uma maneira simples de se obter estas aproximações, é por meio 
do uso da expansão de uma função em série de Taylor, em torno de um dado ponto. Seja 
xj este ponto base, podemos escrever o valor de y( xj+1 )= yj+1, pela seguinte série infinita: 
 
...
!4
)(
!3
)(
!2
)(
)(
4
1
3
1
2
1
11 +
−
′′′′
+
−
′′′
+
−
′′
+−′+= +++++
jjjjjjjjj
jjjjj
xxyxxyxxy
xxyyy (3.5) 
 
Enquanto que o valor de y(xj-1)=yj-1 é dado por: 
 
...
!4
)(
!3
)(
!2
)(
)(
4
1
3
1
2
1
11 +
−′′′′
+
−′′′
−
−′′
+−′−=
−−−
−−
jjjjjjjjj
jjjjj
xxyxxyxxy
xxyyy (3.6) 
 
Considere, agora, a necessidade de se aproximar o valor de y´j , o que será feito, 
utilizando-se as expansões acima. Estas equações podem ser escritas de forma mais 
compacta por meio da definição do comprimento do domínio j : hj=xj-xj-1. 
 
Dessa forma, multiplicando a segunda expansão (3.6) por hj+1
2, e diminuindo da 
primeira expansão (3.5) multiplicada por hj
2, obtemos a seguinteexpressão, na qual y´´j foi 
eliminado: 
 
)(
)(
])([
2
2
11
2
1
2
1
22
11
2
hO
hhhh
yhyhhyh
y
jjjj
jjjjjjj
j +
+
−−+
=′
++
−+++ (3.7) 
 
Nela, o O(z) indica que a aproximação tem ordem de grandeza de z, isto é, o valor 
exato da derivada da função, no ponto considerado, é obtido, a partir da expressão 
aproximada, no limite, quando z→0. Esta ordem de grandeza é oriunda do termo de 
menor ordem (ou primeiro termo) entre aqueles que envolvem as derivadas de maior 
ordem. O conjunto deste termos, ou a sua forma simplificada de representação por ordem 
de grandeza, é denominado de erro de truncamento. 
 
Para uma malha uniforme hj=h, qualquer que seja j, a aproximação dada fica com a 
simplificação: 
 
)(
2
211 hO
h
yy
y jjj +
−
=′
−+ (3.8) 
 
Ela é chamada aproximação por diferença central da derivada primeira de y. 
 
Podemos, ainda, usar as expansões, para obter mais duas aproximações para a 
 4
derivada primeira de y, que para uma malha uniforme são dadas por: 
 
)(1 hO
h
yy
y jjj +
−
=′
− (3.9) 
 
Ela é obtida, a partir da segunda expansão (3.6), sendo chamada de aproximação 
por diferença descendente (backward differentiation): 
 
)(1 hO
h
yy
y jjj +
−
=′
+ (3.10) 
 
Ela é obtida, a partir da primeira expansão (3.5), sendo chamada de aproximação 
por diferença ascendente (forward differentiation). 
 
Para uma aproximação da derivada segunda, pode-se somar as duas expansões 
(3.5) e (3.6), obtendo-se: 
 
)(
2
2
2
11 hO
h
yyy
y jjjj +
+−
=′′
−+ (3.11) 
 
Ela é chamada aproximação por diferenças centrais da derivada segunda de y. 
 
A equação anterior envolve três valores funcionais, para descrever uma 
aproximação da derivada segunda da função, o que representa o mínimo necessário para 
isto, já que a derivada primeira tem que ser eliminada da forma final, e, portanto, pelo 
menos duas expansões em série de Taylor têm que ser consideradas. 
 
Nada impede, que seja utilizada uma outra expansão em série de Taylor, para 
melhorar a ordem de aproximação das equações acima. Por exemplo, poder-se-ia utilizar 
a expansão para o valor funcional yj-2 ( ou yj+2 ), para eliminar o primeiro termo do erro de 
truncamento da equação, obtendo-se, assim, uma aproximação de ordem h3. Entretanto, 
aproximações envolvendo mais de três valores funcionais, em pontos adjacentes, 
apresentam uma maior dificuldade de solução das equações algébricas obtidas pelo 
processo de discretização. 
 
Como exemplo, considere o problema de valor de contorno abaixo: 
 
y´´+y´-2y=0 com x=0, y(0)=0 e para x=1, y(1)=1 (3.12) 
 
Seja o domínio discretizado por uma malha uniforme com J, e subdomínios de 
comprimento h e com xo=0 e xf=1. Aplicando a equação diferencial acima, nos pontos 
onde não se conhecem os valores funcionais de y, temos: 
 
y´´j+y´j-2yj=0, j=1,2,3,4....J-1 (3.13) 
 
Utilizando as aproximações das derivadas primeiras e segunda por diferenças 
centrais e arrajandos os termos, tem-se: 
 
2(yj+1 –2 yj + yj-1) + h(yj+1-yj-1) – 4 h
2 yj =0, j=1,2,3,4....,J-1 (3.14) 
 
 5
ou 
 
(2+h) yj+1 - 4(h
2+1) yj + (2-h) yj-1 =0, j=1,2,3,4,.....,J-1 (3.15) 
e yo=0 e yJ=1 
 
Logo, caímos em um sistema linear de equações algébricas. As condições de 
contorno do problema dado acima são chamadas de primeiro tipo, isto é, definem o valor 
da variável no contorno, sendo facilmente incorporadas ao sistema algébrico das 
equações discretizadas. Diversos outros tipos de condições de contorno são possíveis, 
sendo que, a sua utilização no sistema de equações algébricas discretizadas, torna-se um 
pouco mais elaborada. 
 
Em geral, a condição de contorno pode ser não-linear, na qual go e gt são funções 
arbitrárias de x, y e y´. Entretanto, são três, os tipos existentes de condições de contorno 
lineares. 
 
Diz-se, que a condição de contorno é do primeiro tipo, quando o valor da variável 
dependente é dado no contorno, sendo facilmente utilizada nas equações discretizadas. 
O problema acima serve de exemplo, e a forma geral é dada por: 
 
x=xc, y=yc (3.16) 
 
Quando a condição de contorno é de segundo tipo, o valor da derivada da variável 
dependente é dado no contorno, isto é: 
 
x=xc, y´=y´c (3.17) 
 
Esta condição de contorno tem que ser discretizada, para ser combinada com o 
sistema algébrico discretizado, fazendo: 
 
h
yy
y jjj 2
11 −+ −
=′ (3.18) 
 
A condição de contorno é dita de terceiro tipo, quando tem a seguinte forma geral: 
 
x=xc, ay´+by=c (3.19) 
 
Nela, a, b e c são constantes conhecidas. O seu tratamento é similar ao dado às 
condições de contorno de segundo tipo. 
 
Considere, agora, o seguinte problema de valor inicial que envolve apenas uma 
diferencial ordinária na sua forma normal: 
 
y´=f(t,y) com t=0, y(0)=yo (3.20) 
 
Sendo o intervalo genérico entre tj e tj+1, considere diferentes formas de aproximar 
a derivada primeira, utilizando, como informação conhecida, apenas o ponto j, isto é 
y(tj)=yj . 
 
Utilizando a aproximação de diferenças finitas para frente, para y´j, obtém-se : 
 
 6
Yj+1=yj+hf(t,yj) + O(h
2), h=tj+1-tj (3.21) 
 
Fórmula que permite calcular yj+1 a partir de yj , com erro da ordem de h
2. Esta 
equação é explícita no valor desconhecido de yj+1, sendo, pois, o método denominado de 
explícito. Especificamente, representa o método explícito de Euler. Não caímos em um 
sistema linear, tendo a solução direta. 
 
Caso, por outro lado, resolvermos utilizar a aproximação de diferenças finitas para 
trás de y´j+1 , dada com j+1 no lugar de j, podemos aplicar no ponto tj+1 e escrever: 
 
yj+1=yj + h f(tj+1, yj+1) + O(h
2), h=tj+1-tj (3.22) 
 
Fórmula que calcula yj+1 a partir de yj, com erro da ordem de h
2. Note-se que, no 
caso geral, a equação é não linear no valor desconhecido de yj+1, sendo, pois, necessário, 
utilizar-se um método adequado à solução de problemas não lineares, para se obter o 
valor de yj+1.Assim, como yj+1 não pode ser explicitado a partir da equação, o método é 
denominado de implícito. Mais especificamente, este é chamado método implícito de 
Euler. 
 
Além dos dois métodos acima, podemos, ainda, obter um terceiro, a partir da 
aproximação por diferença central de y´j+1/2, no intervalo considerado. Aplicando ao meio 
do intervalo, temos: 
 
yj+1=yj + h f(tj+h/2, yj+1/2) + O(h
3), h=tj+1-tj (3.23) 
 
Fórmula que ainda não pode ser usada para obter yj+1, porque o valor de f no ponto 
considerado, não é conhecido. Entretanto, expandindo fj e fj+1 em série de Taylor, em 
torno do ponto tj+1/2=tj+h/2, tem-se que: 
 
yj+1=yj+h/2 . [ f(tj,yj) + f(tj+1,yj+1) ] + O(h
3), h= tj+1-tj(3.24) 
 
Fórmula que permite calcular yj+1, ainda que de forma implícita. Este método é 
ainda implícito, sendo denominado de método trapezoidal (ou de Crank-Nicholson). 
 
Considere a solução numérica do problema de valor inicial, abaixo: 
 
y´=-y2, t>0 com t=0, y=1 (3.25) 
 
Ela utiliza os métodos de Euler, até o ponto t=1, com passos uniformes de 
integração de 0,2. A solução analítica deste problema é: 
 
y(t)=1/(1+t) (3.26) 
 
Aplicando os métodos, temos: 
 7
 
Tabela 3.1 – Comparação dos métodos de diferenças finitas aplicados à equação 3.25. 
[68]. 
T Euler 
Explícito 
Euler 
Implícito 
Trapezoidal Solução 
Analítica 
0,0 1,0000 1,0000 1,0000 1,0000 
0,2 0,8000 0,8541 0,8310 0,8333 
0,4 0,6720 0,7434 0,7113 0,7143 
0,6 0,5817 0,6572 0,6220 0,6250 
0,8 0,5140 0,5880 0,5528 0,5556 
1,0 0,4612 0,5316 0,4975 0,5000 
 
 
No caso de problema de valor inicial com equações de ordem maior do que 1, uma 
redução de ordem deve ser feita para uma ordem menor, com a inserção de uma variável 
no sistema. Exemplo: 
 
y´´=3y´+y+5x com y(0)=1 e y´(0)=2 (3.27) 
 
Introduz-se a variável z com z=y´, derivando z´=y´´. Assim, temos agora, dois P.V.I. 
(Problema de Valor Inicial) de primeira ordem: 
 
�
�
�
=
=′
�
�
�
=
++=′
1)0(
2)0(
53
y
zy
z
xyzz
 (3.28) 
 
Resolvemos, simultaneamente, os dois P.V.I., construindo uma tabela de x, y e z. 
[52] 
 
Os métodos de Runge-Kutta são de ponto simples, explícitos, mas com diversos 
estágios, de modo a se obter uma maior ordem de aproximação. A idéia básica deste tipo 
de método é definir que a variação da variável dependente, no passo em questão, é dada 
por uma média ponderada de variações desta variável, calculadas com avaliações 
diferentes da função derivada, isto é: 
 
�
�
�
�
�
�
�
>++∆=∆
∆=∆
∆=− �
=
+
1 com , ),(.
),(.1
1
1
ibyatfty
ytfty
yCyy
ijiji
jj
n
i
iijj
 (3.29) 
 
Nela, Ci, ai e bi são coeficientes a serem determinados. Usando a equação acima, 
e ou igualando-a à série de Taylor, determinamos estes coeficientes para a diversa ordem 
de Runge-Kutta. Mostra-se, que Runge-Kutta é, na verdade, uma variação do Método das 
Diferenças Finitas. Por exemplo: 
 
Para segunda ordem: C1+C2=1, aC2=1/2 e bC2=1/2. Se escolhermos C2=1/2, 
chegamos a Euler Modificado, onde C1=C2=0,5 e a=b=1. 
 8
 
Para quarta ordem: C1=C4=1/6, C2=C3=1/3, a2=a3=h/2, a4=h, b2=b3=1/2 e b4=1. 
 
Os problemas matemáticos descritos nesta seção correspondem a modelos mais 
simples, onde existe apenas uma variável independente, seja ela o tempo ou uma 
coordenada espacial. Entretanto, modelos físicos mais elaborados originam equações 
diferenciais parciais (EDP’s), com duas ou mais variáveis independentes. As equações 
diferenciais parciais, com suas condições auxiliares, formam tanto problemas de valor 
inicial quanto problemas de valor de contorno. 
 
A discretização de problemas, em mais de uma variável dependente, segue um 
procedimento similar ao visto para problemas unidimensionais. O primeiro passo, aqui, é, 
como antes, a discretização do domínio de cálculo. Serão considerados, neste estudo, 
problemas com, no máximo, duas coordenadas espaciais, já que estes apresentam todas 
as características dos problemas multidimensionais. A extensão do procedimento para 
problemas tridimensionais é facilmente obtida. 
 
Considere uma função u(t,x,y) definida em um domínio 0≤x≤1, 0≤y≤1 e t ≥0, que 
podem ser coordenadas adimensionais ou não. A discretização do domínio pode ser feita 
com malhas uniformes ou não uniformes. Como não há nenhuma característica 
fundamental do procedimento de discretização, que seja dependente do tipo da malha, a 
atenção especial é dada a malhas uniformes. 
 
 i-1, j+1 i, j+1 i+1, j+1 
 y 
 
 
 i-1, j i, j i+1, j 
 
 i-1, j-1 i, j-1 i+1, j-1 
 
 x 
Figura 3.1 – Malha de diferenças finitas em problemas com duas dimensões espaciais x e 
y. Note-se, que, em geral, as malhas são quadradas, e se adota um índice para cada 
variável i para x e j para y. 
 
Temos assim: ui,j
n=u(tn,xi,yj), onde tn, xi e yj representam cada ponto da malha de 
discretização. 
 
O segundo passo é, também, a aproximação por diferenças finitas das derivadas, 
que aparecem na equação diferencial parcial, e que podem ser obtidas das expansões da 
variável dependente, em série de Taylor, em relação a uma ou mais variáveis 
independentes. [55] 
 
Dada a expansão de Taylor para u(x,y), em um formato clássico: 
 
....),()(
!
1
...),()(
!3
1
),()(
!2
1
),()(),(),(
3
2
+
∂
∂
+
∂
∂
++
∂
∂
+
∂
∂
+
∂
∂
+
∂
∂
+
∂
∂
+
∂
∂
+=++
yxu
y
k
x
h
n
yxu
y
k
x
h
yxu
y
k
x
hyxu
y
k
x
hyxukyhxu
n
 (3.30) 
 
Ou fazendo a expansão até a derivada segunda : 
 9
 
...
),(
2
1),(),(
2
1
),(),(
),(),(
2
22
2
2
+
∂
∂
+
∂∂
∂
+
∂
∂
+
∂
∂
+
∂
∂
+=++
y
yxu
k
yx
yxu
hk
x
yxu
h
y
yxu
k
x
yxu
hyxukyhxu
 (3.31) 
 
Assumindo-se o índice i para x e j para y, tem-se: u(x,y)=ui,j , u(x+h,y)=ui+1,j, 
u(x,y+k)=ui,j+1, e assim por diante, onde cada espaçamento de h em x corresponde a 
somar ou subtrair 1 no índice i, e cada espaçamento de k em y corresponde a somar ou 
subtrair 1 no índice j. Com isso, faz-se a expansão em série, de Taylor, até a derivada 
segunda de u(x+h,y) e u(x-h,y), usando a equação 3.31, abaixo: 
 
...
2
1
....
),(
2
1),(
),(),(
,2
2
2
,,,1
2
2
2
+
∂
∂
+
∂
∂
+=
+
∂
∂
+
∂
∂
+=+
+ jijijiji u
x
hu
x
huu
x
yxu
h
x
yxu
hyxuyhxu
 (3.32) 
e 
...
2
1
....
),(
2
1),(
),(),(
,2
2
2
,,,1
2
2
2
−
∂
∂
+
∂
∂
−=
−
∂
∂
+
∂
∂
−=−
+ jijijiji u
x
hu
x
huu
x
yxu
h
x
yxu
hyxuyhxu
 (3.33) 
 
Da equação 3.32, truncando todos os termos maiores que h2, ou seja, com erro 
O(h2), e isolando a derivada primeira, chega-se a uma primeira aproximação por 
diferenças finitas da derivada parcial de u em relação a x: 
 
x
uu
h
uu
x
u jijijijiji
∆
−
=
−
=
∂
∂ ++ ,,1,,1, (3.34) 
 
Da equação 3.33, truncando todos os termos maiores que h2 , ou seja, com erro 
O(h2), e isolando a derivada primeira, chega-se a uma segunda aproximação por 
diferenças finitas da derivada parcial de u em relação a x: 
 
x
uu
h
uu
x
u jijijijiji
∆
−
=
−
=
∂
∂
−− ,1,,1,, (3.35) 
 
Somando-se a equação 3.32 com a equação 3.33, truncando todos os termos 
maiores que h3, ou seja, com erro O(h3), e isolando a derivada primeira, chega-se a uma 
terceira aproximação por diferenças finitas da derivada parcial de u em relação a x: 
 
x
uu
h
uu
x
u jijijijiji
∆
−
=
−
=
∂
∂
−+−+
22
,1,1,1,1,(3.36) 
 
Subtraindo-se a equação 3.33 da equação 3.32, truncando todos os termos 
maiores que h3, ou seja, com erro O(h3), e isolando a derivada segunda, chega-se a uma 
aproximação por diferenças finitas da derivada parcial de segunda ordem de u em relação 
a x: 
 10
 
2
,1,,1
2
,1,,1
2
,
2
)(
22
x
uuu
h
uuu
x
u jijijijijijiji
∆
+−
=
+−
=
∂
∂
−+−+ (3.37) 
 
A última etapa para a solução de uma equação diferencial parcial por diferenças 
finitas, é substituir as aproximações das derivadas na equação e nas suas condições de 
contorno, gerando um sistema algébrico, cuja solução fornece a solução aproximada do 
problema original. 
 
Dado o exemplo abaixo, considere uma barra fina metálica de 1 metro de 
comprimento, e estude a condução de calor, segundo a equação: 
 
2
2
x
u
t
u
∂
∂
=
∂
∂
 (3.38) 
 
u(x,t) é a temperatura da barra na posição x e instante t. 
 
Condição Inicial: u(x,0)=0oC 
Condições de Contorno: u(0,t)=u(1,t)=100oC 
Passos: ∆x=0,2metros e ∆t=0,01segundos 
 
Substituindo, na equação diferencial, as aproximações de diferenças finitas das 
equações 3.34 e 3.37, e trocando o índice i por j para a derivada em relação a t, tem-se: 
 
2
,1,,1,1,
)(
2
x
uuu
t
uu jijijijiji
∆
+−
=
∆
−
−++ (3.39) 
 
Isolando o termo ui,j+1 , que corresponde ao tempo seguinte, e substituindo os 
valores de ∆x e ∆t, tem-se : 
 
4
2 ,1,,1
1,
jijiji
ji
uuu
u −++
++
= (3.40) 
 
A equação acima fornece a temperatura do tempo j+1 em função do tempo anterior 
j. Quando se tem uma equação deste tipo, onde o novo valor é calculado em função de 
valores conhecidos, dá-se o nome de método explícito, pois pode-se explicitar um valor 
desconhecido em função de outros já conhecidos. Em seguida, pode-se chegar na tabela 
abaixo: 
 
 
Tabela 3.2: Solução numérica do problema de condução de calor. [55] 
 I 0 1 2 3 4 5 
J t \ x 0,0 m 0,2 m 0,4 m 0,6 m 0,8 m 1,0 m 
0 0,00 0 0 0 0 0 0 
1 0,01 100 0 0 0 0 100 
2 0,02 100 25 0 0 25 100 
3 0,03 100 37,5 6,25 6,25 37,5 100 
4 0,04 100 45,3 14,0 14,0 45,3 100 
5 0,05 100 51,1 21,8 21,8 51,1 100 
6 0,06 100 56,0 29,1 29,1 56,0 100 
 11
 
Veja-se, agora, um outro exemplo: considere uma placa quadrada fina metálica de 
1 metro de comprimento, e encontre a temperatura de equilíbrio (após um tempo muito 
grande), segundo a equação de Laplace: 
0
2
2
2
2
=
∂
∂
+
∂
∂
y
u
x
u
 (3.41) 
u(x,y) é a temperatura da barra na posição x e y como na figura abaixo: 
 
 y 
 
 1 
 
 
 x 
0 1 
 Figura 3.2 – Placa Quadrada do Exemplo 
 
Condições de Contorno: u(0,y)=u(x,0)=100oC e u(1,y)=u(x,1)=0oC 
Passos: ∆x=∆y=0,25 metros 
 
 
Substituindo, na equação diferencial, as aproximações de diferenças finitas da 
equação 3.37, e trocando o índice i por j para a derivada em relação a y, tem-se: 
 
0
)(
2
)(
2
2
1,,1,
2
,1,,1
=
∆
+−
+
∆
+−
−+−+
y
uuu
x
uuu jijijijijiji (3.42) 
 
Lembrando que ∆x = ∆y, corta os denominadores e isolando o termo ui,j , tem-se : 
 
4
1,1,,1,1
,
−+−+ +++
=
jijijiji
ji
uuuu
u (3.43) 
 
A equação acima fornece a temperatura na posição i, j , em função da média 
aritmética das temperatura de cima, de baixo, da direita e da esquerda. Quando se tem 
uma equação deste tipo, onde o novo valor é calculado em função de valores 
desconhecidos, dá-se o nome de método implícito. Em seguida, pode-se chegar a um 
sistema linear, onde a equação 3.43 é expandida para cada ponto do interior da placa, 
seguindo a numeração dada na figura abaixo: 
 y 
 100 0 0 0 0 
 100 u7 u8 u9 0 
 
 100 u4 u5 u6 0 
 
 100 u1 u2 u3 0 
 
 100 100 100 100 100 
 
 0 0,25 0,5 0,75 1 x 
 
 Figura 3.3 : Malha da discretização da placa. 
 
Expandindo a equação 3.43 para cada ponto do interior, de u1 a u9, chega-se às 
equações: 
 12
 
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
=+−−
=−+−−
=−+−
=−+−−
=−−+−−
=−−+−
=−+−
=−−+−
=−−
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
+++
=
+++
=
+++
=
+++
=
+++
=
+++
=
+++
=
+++
=
+++
=
04
04
1004
04
04
1004
1004
1004
2004
4
00
4
0
4
0100
4
0
4
4
100
4
1000
4
100
4
100100
986
9875
874
9653
86542
7541
632
5321
421
68
9
579
8
48
7
395
6
2846
5
175
4
62
3
513
2
42
1
uuu
uuuu
uuu
uuuu
uuuuu
uuuu
uuu
uuuu
uuu
uu
u
uuu
u
uu
u
uuu
u
uuuu
u
uuu
u
uu
u
uuu
u
uu
u
 (3.44) 
 
Resolvendo o sistema linear da equação 3.44, chega-se à solução abaixo: 
 
000,50428,71714,85
571,28000,50428,71
285,14571,28000,50
321
654
987
===
===
===
uuu
uuu
uuu
 (3.45) 
 
Observa-se a simetria em relação à diagonal principal da placa, e constata-se que 
os valores da temperatura, nos vértices da placa, não entraram no cálculo. [55] 
 
 
 13
3.3. Método dos Elementos Finitos 
 
 
O primeiro método numérico com o intuito de resolver equações diferenciais 
parciais (PDE - Partial Differential Equations) foi o método das diferenças finitas. Neste 
método, o domínio da solução é dividido em uma malha de pontos ou nós discretos. O 
PDE é então aplicado para cada nó e suas derivadas substituídas por diferenças finitas 
divididas. Embora tal aproximação seja, conceitualmente, de fácil compreensão, registra-
se alguns inconvenientes. Em particular, torna-se difícil sua aplicação em sistemas com 
geometria irregular, condições de contorno não usuais ou composição heterogênea. 
Figura 3.4 – (a) Uma peça industrial com geometria irregular e composição não 
homogênea. (b) Tal sistema é muito difícil de modelar com a aproximação por diferenças 
finitas. Isso se deve ao fato de complicadas aproximações serem requeridas nos 
contornos do sistema e na fronteira entre as regiões de diferentes composições. (c) Uma 
discretização de elementos finitos é muito melhor aplicada a tal sistema. [54] 
 
O Método dos Elementos Finitos fornece uma alternativa melhor a tais sistemas. 
Em contraste à técnica das diferenças finitas, o MEF divide o domínio da solução em 
formas simples de regiões ou “elementos”. Uma solução aproximada do PDE pode ser 
desenvolvida para cada um destes elementos. A solução total é então gerada, colocando-
as juntas ou montando-as. Utilizando-se as soluções individuais, toma-se o cuidado de 
assegurar a continuidade dos contornos entre os elementos. Assim, o PDE é satisfeito em 
forma de fatias. 
 
O uso de elementos, em vez de malhas retangulares, fornece melhor aproximação 
em sistemas comformatos irregulares, além de valores desconhecidos poderem ser 
gerados continuamente por meio do domínio da solução inteira, em vez de pontos 
isolados. 
 
 Em razão de uma descrição detalhada ir além do escopo deste texto, este capítulo 
se conterá a uma introdução geral do Método dos Elementos Finitos. O objetivo é mostrar, 
de forma simples e fácil, suas características, princípios e capacidades. Assim, a seção 
seguinte abordará uma visão geral dos passos envolvidos na solução de um problema, 
utilizando o MEF. Este é seguido por um simples exemplo: molas ligadas em séries. 
Embora este exemplo não envolva PDE, permite-se desenvolver e mostrar os principais 
aspectos da aproximação de elementos finitos, sem ser desencorajados por fatores 
complicados. Pode-se, então, discutir algumas características, envolvendo o emprego do 
método de elementos finitos em PDE. 
 
 14
 
 
 
 
Etapas para aplicação do Método dos Elementos Finitos 
 
- Pré-Processamento: 
- Definição do problema e do domínio. 
- Discretização ou divisão do domínio em elementos. 
 
- Processamento: 
- Obter as equações dos elementos [k]{u}={f}. 
- Escolha da função de aproximação. 
- Ajuste ótimo da função de aproximação. 
- Formulação Direta.(ou) 
- Método dos Resíduos Ponderados.(ou) 
- Método Colocacional 
- Método de Subdomínios 
- Método dos Mínimos Quadrados 
- Método de Galerkin 
- Técnica Variacional. 
- Método Rayleigh-Ritz 
- Montagem ou colocação das equações dos elementos juntas [K]{u´}={F´}. 
- Acréscimo das condições iniciais e de contorno [ ] { } { }Fuk ′=′ . 
- Solução do sistema linear (ou não linear) {u´}. 
 
- Pós-Processamento: 
- Apresentação dos resultados ou visualização gráfica. 
- Determinação de variáveis secundárias. 
 
 
3.3.1. Visão Geral 
 
 
Este desenvolvimento foi extraído da referência [54], Chapra 1997. Embora as 
particularidades irão variar, utiliza-se, usualmente, na implementação do Método dos 
Elementos Finitos, um padrão de procedimentos, passo a passo. A seguir, é apresentada 
uma breve visão geral de cada um desses passos, cuja aplicação, nos contextos de 
Engenharia, irá ser demonstrada em itens subseqüentes. 
 
 
3.3.1.1 Discretização (Pré-Processamento) 
 
Este passo envolve a divisão do domínio solução em elementos finitos, e podem 
ser em uma, duas ou três dimensões. Os pontos de interseção das linhas que descrevem 
os lados dos elementos são referenciados como nós, e os lados são chamados de linhas 
ou planos nodais. 
 15
 
 
Figura 3.5 - Exemplos de elementos empregados em (a) uma, (b) duas, e (c) três 
dimensões. 
 
 
3.3.1.2. Equações dos Elementos (Processamento) 
 
A seguir, desenvolvem-se equações, a fim de aproximar a solução de cada 
elemento. Isso envolve dois sub-passos. Primeiro, escolhe-se uma função apropriada com 
coeficientes desconhecidos, que serão usados para aproximar a solução. Por último, 
avaliam-se os coeficientes, em que as funções se aproximam da solução, de forma 
considerada ótima. [54] 
 
Escolha das Funções de Aproximação 
 
Considera-se que, por serem de fácil manipulação matemática, os polinômios são 
freqüentemente empregados para este propósito. Para o caso unidimensional, a 
 Elemento Linear
 (a) Unidimensional
 Nó
 Linha
 Nodal
 Elemento Elemento
 Quadrílateral Triangular
 (b) Bidimensional
 Elemento
 Hexaédrico
 Plano Nodal
 Elemento Tetraédrico
 (c) Tri-dimensional
 16
alternativa mais simples é um polinômio de primeira ordem, ou uma linha reta: 
 
u(x) = a0 + a1 x (3.46) 
 
 
Nesta fórmula, u(x) é a variável dependente; a0 e a1 são constantes; e x é a 
variável independente. Essa função deve passar através dos valores u(x) nos pontos 
finais dos elementos em x1 e x2. Portanto: 
 
u1 = a0 + a1 x1 
u2 = a0 + a1 x2 
 
Onde u1 = u (x1) e u2 = u (x2). Estas equações podem ser resolvidas, usando a 
regra de Cramer, onde: 
 
a0 = (u1 x2 – u2 x1) / (x2 – x1) a1 = (u2 – u1) / (x2 – x1) 
 
Este resultado pode, então, ser substituído na Eq. (3.46), a qual, depois de se 
arrumar os termos, pode ser escrita como: 
 
u = N1 u1 + N2 u2 (3.47) 
 
onde N1 = (x2 – x) / (x2 – x1) (3.48) 
 
e N2 = (x – x1) / (x2 – x1) (3.49) 
 
A equação (3.47) é chamada função “de aproximação” ou “de forma”, e N1 e N2 são 
chamados de funções de interpolação. Inspecionando melhor, percebe-se que a Eq. 
(3.47) é, de fato, o polinômio interpolador de primeira ordem de Lagrange. Ela fornece um 
significado, para predizer valores intermediários (que é interpolar) entre valores dados u1 
e u2 nos nós. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 17
 
Figura 3.6 (b) Uma função de aproximação ou forma para (a) um elemento linear. As 
correspondentes funções de interpolação são mostradas em (c) e (d). Note-se, que a 
soma das funções de interpolação (N1+N2) são iguais a 1. 
 
Em adição, lidar com equações lineares facilita operações como a diferenciação e 
integração. Tais manipulações, mais à frente, serão importantes em outros itens. A 
derivação da Eq. (3.47) é: 
 
2
2
1
1 u
dx
dN
u
dx
dN
dx
du
+= (3.50) 
 
De acordo com as Eq. (3.48) e (3.49), as derivadas de N1 e N2 podem ser 
calculadas como: 
 
12
2
12
1 1 
1
xxdx
dN
xxdx
dN
−
=
−
−= (3.51) 
 
E, portanto, a derivada de u é: 
 
)(
1
21
12
uu
xxdx
du
+−
−
= (3.52) 
 
Em outras palavras, essa é a diferença dividida, representando a inclinação da reta 
conectada nos nós. 
 
A integral pode ser expressa como: 
 
 Nó 1 Nó 2
 (a)
 u1 u
 u2
 x1 x2
 (b)
 1 N1
 x1 x2
 (c)N2 1
 x1 x2
 (d)
 18
dxuNuNudx
x
x
x
x
 2211
2
1
2
1
+=� � 
 
Cada termo, no lado direito, é somente a integral de um triângulo reto com base x2 
– x1 e altura u. Isto é: 
 
uxxNudx
x
x
)(
2
1
12
2
1
−=� 
 
Assim, a integral inteira é: 
 
)(
2 12
21
2
1
xx
uu
udx
x
x
−
+
=� (3.53) 
 
Ou seja, simplesmente, a regra dos trapézios. 
 
 
Obtenção de um Ajuste Ótimo da Função de Aproximação 
 
Após a escolha da função de interpolação, as equações que governam o 
comportamento dos elementos deve ser desenvolvida. Elas representam um ajuste da 
função de aproximação, com a finalidade de solução da subjacente equação diferencial. 
Vários métodos são disponíveis para este propósito. Entre os mais comuns, cita-se a 
aproximação direta, o método dos resíduos ponderados, e a técnica variacional. O 
resultado desses métodos é análogo para o ajuste de curvas. Contudo, em vez de ajustar 
funções para dados, eles especificam relacionamentos entre a desconhecida Eq. (3.47) 
para satisfazer as subjacentes PDE, em uma forma apropriada. 
 
Matematicamente, o resultado das equações dos elementos irá, freqüentemente, 
consistir de um conjunto de equações lineares algébricas, podendo ser expressa na forma 
matricial: 
 
[ k ] { u } = { F } (3.54) 
 
Nela, [ k ] é uma matriz propriedade ou rigidez do elemento; { u } é um vetor coluna 
de valores desconhecidos dos nós; e { F } é um vetor coluna, refletindo o efeito de 
quaisquer influências externas aplicadas nos nós. Percebe-se, que, em alguns casos, as 
equações podem ser não lineares. Contudo, nos exemplos elementares descritos aqui e 
em muitos dos problemas práticos, os sistemas são lineares. 
 
 
3.3.1.3. Montagem (Processamento) 
 
 Depois de se obter as equações dos elementos individuais, elas devem ser 
colocadas juntas ou montadas, para caracterizar o comportamento unificado do sistema 
inteiro. O processo de montagem é governado pelo conceito de continuidade. Isto é, as 
soluções de elementos contíguos são combinadas, e os valores desconhecidos (algumas 
vezes, as derivadas) de seus comuns nós são equivalentes. Assim, a solução total será 
contínua. 
 19
 
Quando todas as versões individuais da Eq. (3.54) são, finalmente, montadas, o 
sistema inteiro é expresso sob forma matricial, como: 
 
[ K ] { u´ } = { F´ } (3.55) 
 
Nela, [ K ] é a matriz propriedade montada e { u´ } e { F´ } são vetores colunas de 
valores desconhecidos dos nós e forças externas feitas com apóstrofos, para denotar uma 
montagem dos vetores { u } e { F } dos elementos individuais. 
 
 
3.3.1.4. Condições de Contorno e Iniciais (Processamento) 
 
Antes da Eq. (3.55) poder ser resolvida, deve-se modificá-la, para considerar as 
condições iniciais e de contorno do sistema. Estes ajustes resultam em: 
 [ ] { } { }Fuk ′=′ (3.56) 
 
Nela, as barras significam as condições de contorno incorporadas. 
 
3.3.1.5. Solução (Processamento) 
 
Pode-se obter a solução da Eq. (3.56), com técnicas para a resolução de sistemas 
lineares e não lineares. Em muitos casos, os elementos serão configurados, de modo que 
as equações resultantes possam ser unidas, diminuindo o tamanho do sistema. Assim, o 
esquema de eficiência mais alto disponível a cada sistema é possível de ser empregado. 
O uso de simetrias, onde uma parte do sistema é igual à outra, possibilita a redução da 
ordem do sistema linear. 
 
3.3.1.6 Apresentação dos Resultados (Pós-Processamento) 
 
Obtida a solução, esta será exibida na forma de tabelas ou gráficos. Em adição, 
variáveis secundárias serão determinadas e expressas. 
 
Apesar dos passos precedentes serem muito genéricos, eles são comuns na 
maioria das implementações do Método dos Elementos Finitos. No item seguinte, ilustra-
se como eles podem ser aplicados na obtenção do resultado numérico de dois sistemas 
físicos simples – o primeiro, molas ligadas em série, e depois, uma haste sendo aquecida. 
 
3.3.2. Exemplo da Temperatura de Equilíbrio 
 
Nos item seguinte vai-se exemplificar o Método dos Resíduos Ponderados usando 
a Técnica de Galerkin para um caso bidimensional. Usa-se um exemplo semelhante ao do 
item 3.2 do MDF. Considere uma placa quadrada fina metálica de 1 metro de 
comprimento, e encontre a temperatura de equilíbrio (após um tempo muito grande), 
segundo as equação de Laplace: 
 
0
2
2
2
2
=
∂
∂
+
∂
∂
y
u
x
u
 (3.57) 
 
 20
u(x,y) é a temperatura da barra na posição x e y como na figura abaixo: 
 
 
 
 
 y 
 
 1 
 
 
 x 
0 1 
 Figura 3.7 – Placa Quadrada do Exemplo. 
 
Condições de Contorno: u(x,0)=0oC, u(x,1)=100oC, u(0,y)=75oC e u(1,y)=50oC 
 
 
3.3.3. Problema Bidimensional 
 
Este item foi extraído da referência [54]. Embora o número de operações 
matemáticas cresça bastante, a extensão da aproximação de elementos finitos para duas 
dimensões é, conceitualmente, similar à aplicação unidimensional discutida. Assim, ela 
segue os mesmos passos, como mostrado no subitem 3.2.1. 
 
3.3.3.1. Discretização 
 
Uma variedade de simples elementos, como triângulos e quadriláteros, são 
usualmente empregados para fracionar os elementos finitos em duas dimensões. Na 
discussão presente, limitar-se-á a elementos triangulares do tipo descrito na Fig. 3.8. 
 y 
 3 
 
 
 2 1 
 
 x 
Figura 3.8 Um elemento triangular 
 
3.3.3.2. Equações dos Elementos 
 
Como no caso unidimensional, o próximo passo é desenvolver uma equação para 
aproximar a solução ao elemento. Em um elemento triangular, a aproximação mais 
simples é um polinômio linear [Compare com a Eq. (3.46)]. 
 
u(x, y) = a0 + a1,1x + a1,2y (3.58) 
 
Onde u é a variável dependente, os a’s representam os coeficientes, e x e y são 
variáveis independentes. Essa função deve passar pelos valores de u(x, y) nos nós do 
triângulo (x1, y1), (x2, y2), e (x3, y3). Portanto: 
 
u1(x, y) = a0 + a1,1x1 + a1,2 y1 
u2(x, y) = a0 + a1,1x2 + a1,2 y2 
 21
u3(x, y) = a0 + a1,1x3 + a1,2 y3 
 
Ou na forma matricial: 
 
�
�
�
�
�
�
�
�
�
�
=
�
�
�
�
�
�
�
�
�
�
�
�
�
	
�
�
�
�
3
2
1
2,1
1,1
0
33
22
11
1
1
1
u
u
u
a
a
a
yx
yx
yx
 
 
Este pode ser resolvida a: 
 
)]()()([
2
1
1221331132233210 yxyxuyxyxuyxyxuA
a
e
−+−+−= (3.59) 
)]()()([
2
1
2131323211,1 yyuyyuyyuA
a
e
−+−+−= (3.60) 
)]()()([
2
1
1233122312,1 xxuxxuxxuA
a
e
−+−+−= (3.61) 
 
Onde Ae é a área do elemento triangular: 
 
)]()()[(
2
1
122131132332 yxyxyxyxyxyxAe −+−+−= 
 
As equações (3.59) até (3.61) podemser substituídas na Eq. (3.58). Depois de 
agrupar os termos, o resultado pode ser expresso como: 
 
u = N1u1 + N2u2 + N3u3 (3.62) 
 
Onde: 
 
])()()[(
2
1
233223321 yxxxyyyxyxA
N
e
−+−+−= 
])()()[(
2
1
311331132 yxxxyyyxyxA
N
e
−+−+−= 
])()()[(
2
1
122112213 yxxxyyyxyxA
N
e
−+−+−= 
 
A equação (3.62) fornece uma maneira de predizer valores intermediários para o 
elemento, com base nos valores dos nós. A Figura 3.9 mostra a função de aproximação 
com as correspondentes funções de interpolação. Percebe-se, que a soma das funções 
de interpolação é sempre igual a um. 
 
 22
 
Figura 3.9 (a) Uma linear função de aproximação para um elemento triangular. As 
correspondentes funções de interpolação são mostradas em (b) até (d). Fonte [54] 
 
Também no caso unidimensional, vários métodos são disponíveis para desenvolver 
as equações dos elementos, baseadas na subjacente PDE e nas funções de 
aproximação. As equações resultantes são, consideravelmente, mais complicadas que as 
do caso unidimensional. Contudo, devido às funções de aproximação serem usualmente 
polinômios de mais baixa ordem, como na Eq. (3.58), os termos da matriz final do 
elemento consistirá a polinômios de baixa ordem e constantes. 
 
 23
3.3.3.3. Ajuste Ótimo da Função de Aproximação 
 
O desenvolvimento foi extraído da referência [37]. Seja em R2, espaço onde se 
situa o maior número de problemas físicos, a ortogonalidade de duas funções f e g é dada 
por: 
 
�
Ω
= 0 . dvgf 
 
Logo, na resolução de um problema, onde as derivadas parciais podem ser 
traduzidas pela pesquisa de uma função v , tal como os operadores L sobre o domínio e B 
sobre a fronteira, que verifica-se: L(v) – f = 0 e B(v) – g =0. [37] 
 
O método denominado Método dos Resíduos Ponderados consiste, na pesquisa de 
funções v , em que satisfaçam a condição de contorno, ponderadas por funções u, tais 
que, para toda função u que satisfaça condições de continuidade determinadas, se possa 
escrever: 
 
�
Ω
=− 0))(.( dwfvLu (3.63) 
 
Se o conjunto de funções u é de dimensão infinita, então, é possível obter uma 
equivalência entre o problema, as derivadas parciais e sua formulação integral. 
Entretanto, nas aplicações práticas, as funções u formam um espaço de dimensão finita, e 
a fórmula (2.63) constitui uma única aproximação caracterizada pela função dada e por 
este conjunto de funções. 
 
A vantagem do Método dos Resíduos Ponderados, em relação à formulação 
variacional, é poder se aplicar a qualquer equação, independentemente, da existência e 
do conhecimento de uma formulação variacional do problema; por outro lado, de início, 
existe um erro de método caracterizado pela escolha das funções u; no entanto, este 
último ponto é de importância secundária, pois este erro e o de aproximação se 
conjugam, para resultar, sob certas condições, os mesmos resultados nas duas 
formulações. 
 
Tem-se, como exemplo, o problema térmico (Para o exemplo do item 3.3.2 tem-se 
k=1 e Q=0): 
 
sITT
Q
y
T
k
yx
T
k
x
0
0
=
=+��
�
�
��
�
�
∂
∂
∂
∂
+�
�
�
�
�
�
∂
∂
∂
∂
 
 
A formulação, em termos do Método dos Resíduos Ponderados, consiste em 
escolher funções T, que verificam as condições de contorno, onde u∀ é: 
 
0=��
�
�
�
�
�
�
+��
�
�
��
�
�
∂
∂
∂
∂
+�
�
�
�
�
�
∂
∂
∂
∂
��٠dQy
T
k
yx
T
k
x
u 
 
Uma integração, por partes, permite transformar esta integral: 
 24
 
0=
∂
∂
+��
�
�
�
�
�
�
+��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
− ��� dSn
T
kuduQ
y
u
y
T
x
u
x
T
k
S
 
 
Se impusermos u=0 sobre o contorno, o segundo termo desaparece e resulta a: 
 
0=��
�
�
�
�
�
�
−��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
�� dQuy
u
y
T
x
u
x
T
k 
 
A seqüência de funções de projeção ui tomada e as funções de aproximação αi 
escolhida para T caracterizarão, portanto, o método. [37] 
 
O princípio geral consiste em determinar os coeficientes u1, u2,..., uNN da 
aproximação u* , pela realização de um certo número de condições a impor. Se a função u 
é substituída por sua aproximação u* sobre todo domínio Ω, obtém-se: 
 
���
===
′=′′=′=
NN
i
iyiyix
NN
i
ix
NN
i
ii NuuNuuNuu
1
*
1
*
1
* ; ; , 
 
Onde os ui são os coeficientes numéricos. Como foi visto, anteriormente, são, de 
fato, valores de u* em cada um dos nós da discretização. Logo, o funcional vem a ser uma 
função exclusiva dos coeficientes u1, u2,..., uNN. Pode-se, então, escrever a função 
pesquisada u* aproximada pela combinação linear: 
 
NNNN NANANAu +++= ...2211
* 
 
Onde os coeficientes A1, A2,..., ANN serão determinados pelo método, de forma a 
realizar a melhor aproximação possível de u sobre a base de funções N1, N2,..., NNN. 
 
No Método dos Elementos Finitos, o domínio de estudo é discretizado em 
subdomínios chamados Elementos Finitos, sobre os quais, a função procurada é 
aproximada por um polinômio. A discretização realizada é uma partição do domínio, sem 
buracos nem recobrimentos. Os polinômios próprios de cada elemento devem respeitar, 
na fronteira, as condições de continuidade compatíveis com aquelas impostas pela 
natureza do problema. Este vínculo permite determinar o conjunto de funções Ni, a partir 
das funções Ni
(e) definidas sobre cada elemento. 
 
Um exemplo em R2, é uma divisão em dois sub-domínios triangulares de primeira 
ordem, sobre os quais, a função procurada é aproximada por um polinômio de primeira 
ordem: P=a+bx+cy. Em cada elemento, existem três coeficientes a determinar, e, 
portanto, três monômios, perfazendo um total de seis coeficientes não conhecidos; mas 
as condições de continuidade sobre a aresta, que une os vértices 2 e 3 aos seis 
coeficientes, restringem a 4 o número real de coeficientes não conhecidos. Há, então, 
uma função de aproximação que afeta cada vértice, como mostrado na Figura 3.10. 
 
 
 
 
 
 
 25
 
 Figura 3.10 - Domínio a dois elementos triangulares. 
 
A escolha de cada elemento das funções de aproximação definirá o seu tipo, e 
caracterizará, pela seqüência, a natureza das funções de aproximação. Há uma ligação 
matemática rigorosa entre a escolha da natureza das funções de aproximação (lineares, 
quadráticas, cúbicas) e a forma dos elementos, sempre definidas por trechos sobre a 
discretização. 
 
Na aplicação deste método, é preciso escolher um conjunto de funções de projeção 
(ou funções de ponderação) Φ1, Φ2,..., ΦNN antes de escrever as equações de projeção 
L(u*) sobre cada uma destas funções. Na técnica de Galerkin tem-se N1=Φ1 , N2=Φ2 e 
assim por diante.[37] 
 
 
0
1
1 =Ω�
�
�
�
�
�
�
�
−�
�
�
�
�
�
�
�
Φ ���
=
dfNuL
NN
j
jj 
0
1
2 =Ω�
�
�
�
�
�
�
�
−�
�
�
�
�
�
�
�
Φ ���
=
dfNuL
NN
j
jj (3.64) 
... 
0
1
=Ω
�
�
�
�
�
�
�
�
−�
�
�
�
�
�
�
�
Φ ���
=
dfNuL
NN
j
jjNN 
 
Obtém-se, de novo, um sistema de NN equações algébricas a resolver para se 
encontrar as NN incógnitas u1, u2,..., uNN. 
 
 
Observações Importantes: 
 
A aplicação do método dos elementos finitos conduz à substituição de uma 
equação ou sistema de equações a derivadas parciais por um sistema de equações 
algébricos, contento coeficientes da função de aproximação, que são, na realidade, os 
valoresdas funções de aproximação nos nós do domínio discretizado. Como no Método 
dos Resíduos Ponderados as funções de ponderação são idênticas às funções de 
aproximação, o sistema de equações obtido é idêntico àquele obtido pela formulação 
variacional. 
 
Se o operador L é linear, então a funcional é quadrática, implicando na linearidade 
do sistema de equações algébricas obtido. Da mesma forma, no Método dos Resíduos 
Ponderados, a linearidade de L implica na linearidade das equações (3.64), pois, em cada 
uma delas, os coeficientes ui podem ser colocadas em evidência na integral. 
 3
 1
 4
 2
 26
 
( )( )� �������
=
Ω
=
ΩΦ−ΩΦ=Ω
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
−Φ
NN
i
iijj
NN
j
jji fddNLudfNuL
11
 
 
 
 
3.3.3.4. Condições de Contorno e Montagem 
 
A incorporação das condições de contorno e a montagem do sistema matricial 
também se tornam mais complicados, quando a técnica dos elementos finitos é aplicada 
em problemas de duas ou três dimensões. Contudo, como na derivação da matriz 
elemento, as dificuldades relativas ao mecanismo dos processos são maiores do que a 
complexidade conceitual. Por exemplo, o estabelecimento da topologia do sistema, o qual 
era trivial para o caso de uma dimensão, torna-se um problema de grande importância em 
duas ou três dimensões. Em particular, a escolha do esquema de numeração irá ditar a 
estrutura do sistema matricial resultante e, portanto, a eficiência com a qual ele pode ser 
resolvido. A Figura 3.11 mostra um esquema desenvolvido para uma placa plana em 
equilíbrio térmico, solucionada por meio do método dos elementos finitos. 
 
 
 
 
Figura 3.11 - Um esquema de numeração dos nós e elementos para uma aproximação 
por elementos finitos, para uma placa plana em equilíbrio térmico. [54] 
 
Logo para o exemplo tem-se: 
 
2
2
2
2
y
T
x
T
R
∂
∂
+
∂
∂
= 
 
Integrando pelo Método dos Resíduos ponderados e associado à técnica de 
Galerkin, tem-se: 
 
0
2
2
2
2
=��
�
�
��
�
�
∂
∂
+
∂
∂
��
Ω
d
y
T
x
T
Ni 
Y 
X 
100 oC 
O oC
50 oC 
75 oC 
 27
 
onde Ni são as funções de teste no caso as funções de Lagrange N1, N2 e N3 da função 
de aproximação T = T1 . N1 + T2 . N2 + T3 . N3. Após uma integração por partes (Teorema 
de Green), obtém-se: 
 
dS
n
T
Nd
y
N
y
T
x
N
x
T
S
i
ii
��� ∂
∂
−=∆��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
∆
 
 
para cada elemento triangular. O termo do lado direito representa o fluxo sobre o contorno 
do elemento triangular (Percorrendo no sentido anti-horário). Assim para um elemento 
genérico 1,2 e 3 pode-se escrever: 
 
[ ]( ) [ ]( ) [ ]( )312312 F 31 F 23 F 12 ladonoFluxoladonoFluxoladonoFluxodSnTNS i ++=∂
∂
− � 
 
Pode-se aplicar as equações acima em cada elemento da malha: 
Para o Elemento (1) tem-se: 
 
 Nó 7 (0,25;0,25) 
 
 
 
 
 
 Nó 1 (0;0) Nó 2 (0,25;0) 
 
 
Figura 3.12 – Elemento 1 da Malha com as coordenadas (x,y) de cada nó. 
 
Para i=1: 
 
dS
n
T
Ndydx
y
N
y
T
x
N
x
T
S
x
�� � ∂
∂
−=��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
1
25,0
0 0
11 
 
Substituindo a integral do lado direito pelo fluxo em cada lado do triangulo e 
lembrando que N1 é zero para o lado oposto ao vértice 1 (F27=0) tem-se: 
 
7112
25,0
0 0
11 FFdydx
y
N
y
T
x
N
x
Tx
+=��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
� � 
 
Resolvendo o lado esquerdo tem-se: 
 
711221 5,05,0 FFTT +=− (3.65) 
 
Para i=2 tem-se de maneira similar: 
 
 28
1227
25,0
0 0
22 FFdydx
y
N
y
T
x
N
x
Tx
+=��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
� � 
ou 
 
122771 5,05,0 FFTT +=+− (3.66) 
 
Para i=3 tem-se: 
 
2771
25,0
0 0
33 FFdydx
y
N
y
T
x
N
x
Tx
+=��
�
�
��
�
�
∂
∂
∂
∂
+
∂
∂
∂
∂
� � 
 
ou 
 
277172 5,05,0 FFTT +=+− (3.67) 
 
As equações de (3.65) até (3.67) formam as equações elemento para o primeiro 
triangulo. De maneira similar para o segundo elemento tem-se: 
 
�
�
�
�
�
+=+−
+=−
+=−
177676
766176
611761
5,05,0
5,05,0
5,05,0
FFTT
FFTT
FFTT
 
 
Cada equação acima foi obtida para N1, N2 e N3, respectivamente. Deve notar que 
F71=-F17 ou seja o fluxo que sai do elemento 1 e vai para o elemento 2 através do lado 
modal 17 é o mesmo que sai do elemento 2 e vai para o elemento 1 só que com sentido 
oposto. Isso é fundamental na fase de montagem onde os termos de fluxo interno na 
placa irão se anular. 
 
Vai-se procedendo assim para os 32 elementos da malha do exemplo. Para os 
triângulos inferiores (Elementos impares 1,3,5,7...31) tem-se: 
 
�
�
�
�
�
+=+−
+=+−
+=−
233132
231231
311221
5,05,0
5,05,0
5,05,0
FFTT
FFTT
FFTT
 
 
Para os triângulos superiores (Elementos pares 2,4,6,8...32) tem-se: 
 
�
�
�
�
�
+=+−
+=−
+=−
133232
322132
211321
5,05,0
5,05,0
5,05,0
FFTT
FFTT
FFTT
 
 
Montando as equações elementos tem-se 25 equações e 41 incógnitas, das quais 
25 são de T1, T2,..., T25 mais 16 correspondente ao fluxo em cada lado externo dos 
elemento (F12, F23, F34, F45, F5,10 .... F61). 
 
 
 29
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
=−−++−−
+=+−+−=−−++−−
=+++−+=+−+−
=+++−+=+−+−
+=+−+−=−−+−−−
+=+++−=−−+−−−
+=+−+−=−−+−+−
+=+++−=+−+−
=−+−+−=−−+−
=++−+−+=++−+−
+=+−+−=−−+−
+=++−−=+−+−
=−+−−+=+−+
;05,05,05,05,05,0
;5,05,05,0;05,05,05,05,15,05,0
;5,05,05,0;5,05,05,0
;5,05,05,0;5,05,05,0
;5,05,05,0;05,05,05,15,05,0
;5,05,05,0;05,05,05,15,05,0
;5,05,05,0;05,05,05,05,05,0
;5,05,05,0;5,05,0
;05,05,05,05,0;5,05,05,0
;05,05,05,0;5,05,05,0
;5,05,05,0;5,05,05,0
;5,05,05,0;5,05,05,0
;05,05,05,05,0;5,05,0
1915141387
25,2025,24252420191817131276
24,232524191816,1111,61712116
23,222423181715,1010,51514109
23,2222,212322171615159843
22,2121,162221171614138732
20,1525,202520151413128732
24,2525,20251915141612761
242318132110,591054
232217131210,545109543
21,1616,1121161211349432
20,1510,52015109238421
19181410916128721
TTTTTT
FFTTTTTTTTTT
FTTTTFFTTTT
FTTTTFFTTTT
FFTTTTTTTTTT
FFTTTTTTTTTT
FFTTTTTTTTTT
FFTTTTFTTTT
TTTTTFTTTT
TTTTTFFTTTTT
FFTTTTFTTTT
FFTTTTFTTTT
TTTTTFFTTTT
 
 
Pode-se agora colocar as condições de contorno: 
 
�
�
�
�
�
�
�
=====
===
===
=====
;100;100;100;100;100
;50;50;50
;75;75;75
;0;0;0;0;0
2524232221
201510
16116
54321
TTTTT
TTT
TTT
TTTTT
 
 
Resolvendo o sistema final de 25 equações e 25 incógnitas tem-se : 
 
�
�
�
�
�
===
===
===
93,3326,3386,42
46,5225,5617,63
64,6912,7657,78
987
141312
191817
TTT
TTT
TTT
 
 
Solução para as temperaturas internas na placas em equilíbrio térmico. As outras 
incógnitas tem um interesse menorpois são o fluxo nas bordas da placa. 
 
 
 
 
 30
 
3.3.3.5. Solução e Apresentação do Resultado (Pos-Processamento) 
 
Embora o mecanismo seja complicado, o sistema matricial é meramente um 
conjunto de n equações simultâneas, que pode ser solucionada, para se encontrar os 
valores das variáveis dependente nos n nós. A Figura 3.13 mostra uma solução possível, 
no caso de uma placa em equilíbrio térmico. 
 
 
 
Figura 3.13 - A distribuição da temperatura de uma placa em equilíbrio térmico, sendo 
calculada por meio do Método dos Elementos Finitos. [54] 
 
 
 3.3.3.6. Comparação com a Solução pelo MDF 
 
 Mesmo neste exemplo com meio homogêneo e geometria simples, nota-se uma 
vantagem no MEF, pois não se é obrigado a usar malhas retangulares de mesmo 
tamanho como no MDF. Em outros problemas com meios não homogêneos, geometrias 
complexas e parâmetros não lineares as vantagens do MEF aumentam bastante. O que 
torna o MEF um método bem mais amplo e com uma aplicabilidade prática muito maior 
que o MDF. Assim o aumento na complexidade matemática dos elementos finitos é 
compensada por sua portabilidade e eficiência na solução de equações diferenciais 
parciais encontradas na engenharia. Por esta razão o MEF é o método mais comumente 
implementado pelas ferramentas de CAE e para muitos autores é a base da Engenharia 
Assistida por Computador. 
 
 
 
 
 31
Apêndice A 
Ajuste Ótimo da Função de Aproximação em uma 
Dimensão 
 
A.1. Solução de Elementos Finitos para Molas em Séries 
 
 
 Figura A.1 (a) Uma série de molas interconectadas. Uma extremidade é fixada na 
parede, enquanto a outra é submetida a uma força constante F. (b) Representação em 
elementos finitos. Cada mola é representada por um elemento. Portanto, o sistema 
consiste de quatro elementos e cinco nós. [54] 
 
Descrição do problema: a figura A.1 mostra uma série de molas Interconectadas. 
Uma extremidade é fixada a uma parede, enquanto a outra é sujeita a uma força 
constante F. Usando, passo a passo, os procedimentos do Método dos Elementos Finitos, 
determina-se o deslocamento das molas. 
 
Solução − Discretização: o modo de particionar esse sistema é, obviamente, tratar 
cada mola como um elemento. Assim, o sistema consiste de quatro elementos e cinco 
nós (Fig. A.1b). 
 
Equações dos Elementos: como este sistema é muito simples, suas equações dos 
elementos podem ser escritas diretamente, sem o recurso da aproximação matemática. 
Este é um exemplo de aproximação direta para os elementos derivados. 
 
 
 
 32
 
 
 
 
Figura A.2 - Um diagrama de corpo livre de um sistema de molas. [54] 
 
A Figura A.2 mostra um elemento individual. O relacionamento entre a força F e o 
deslocamento x pode ser representado, matematicamente, pela lei de Hooke: 
 
F = k x 
 
Onde, k representa a constante da mola, a qual pode ser interpretada como a força 
requerida para causar uma unidade de deslocamento. Se uma força F1 é aplicada no nó 
1, o seguinte balanço de força (reação) deve segurar: 
 
F = k (x1 – x2) 
 
Onde, x1 é deslocamento do nó um da sua posição de equilíbrio; e x2 o 
deslocamento do nó dois da sua posição de equilíbrio. Assim, x2 – x1 representa o quanto 
a mola é alongada ou comprimida e relativa ao equilíbrio (Fig. A.2). 
 
Essa equação pode também ser escrita como: 
 
F1 = k x1 – k x2 
 
Para um sistema estacionário, um balanço de forças também necessita que F1=F2 
e, portanto: 
 
F2 = -k x1 + k x2 
 
Estas duas simultâneas equações especificam o comportamento do elemento em 
resposta às forças aplicadas. Podem ser escritas numa forma matricial, como: 
 
�
�
�
�
�
�
=
�
�
�
�
�
�
�
�
	
�
�
−
−
2
1
2
1
F
F
x
x
kk
kk
 
 
Ou então: 
 
[ k ] { x } = { F } , 
 
Onde a matriz [ k ] é a matriz propriedade do elemento. Neste caso, é também 
referenciada com a matriz rigidez do elemento. Note-se, que esta última equação tem sido 
 33
moldada no formato da Eq. (3.54). Assim, obteve-se sucesso na geração de uma equação 
matricial, que descreve o comportamento de um elemento típico no sistema. 
 
Antes de proceder para o próximo passo, montagem da solução total, introduzir-se-
á alguma notação. Os elementos [ k ] e { F } são convencionalmente colocados 
sobreescrito e subescrito, como em: 
 
 
�
�
�
�
�
�
=
�
�
�
�
�
�
�
�
	
�
�
−
−
)(
2
)(
1
2
1
)(
22
)(
21
)(
12
)(
11
e
e
ee
ee
F
F
x
x
kk
kk
 , 
 
 
Onde, o sobreescrito (e) designa que estas são equações elemento; os k’ s são 
também colocados subescritos, e kij denota sua localização na linha i e coluna j da matriz. 
Para o presente caso, elas são também fisicamente interpretadas como representando a 
força requerida no nó i, para induzir uma unidade de deslocamento no nó j. 
 
Montagem − Antes das equações elementos serem montadas, todos os elementos 
e nós devem ser numerados. Esse esquema global de numeração especifica a 
configuração ou topologia do sistema (o presente caso usa um esquema idêntico ao da 
tabela A.1). Ou seja, mostra-se o nó que pertence a cada elemento. Uma vez que a 
topologia é especificada, as equações, para cada elemento, podem ser escritas com 
referência às coordenadas globais. 
 
As equações elementos podem então ser adicionadas, uma de cada vez, para 
montar o sistema total. O resultado final pode ser expresso na forma matricial como 
[lembrando Eq. (3.55)]: 
 
[ K ] { x´ } = { F´ } 
 
Onde: 
 
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
−
−+−
−+−
−+−
−
=
)4(
22
)4(
21
)4(
12
)4(
11
)3(
22
)3(
21
)3(
12
)3(
11
)2(
22
)2(
21
)2(
12
)2(
11
)1(
22
)1(
21
)1(
12
)1(
11
][
kk
kkkk
kkkk
kkkk
kk
K 
 
e 
 
�
�
�
�
��
�
�
�
�
�
�
�
��
�
�
�
=′
)4(
2
)1(
1
0
0
0
}{
F
F
F 
 
E { x´ } e { F´ } são os vetores deslocamento e força expandida. Quanto às 
 34
equações que foram montadas, as forças internas se cancelaram. Assim, o resultado final 
para { F´ } é zero, em todas as linhas, exceto no primeiro e último nó. 
 
Antes de proceder o próximo passo, deve-se comentar a estrutura da matriz 
propriedade montada. Ela é tridiagonal. Isso é um resultado direto do esquema particular 
de numeração escolhido (Tabela A.1) antes da montagem. Embora não seja muito 
importante no contexto presente, com a realização de tal união, sistemas esparsos podem 
ser uma vantagem na colocação de problemas mais complicados. Isso é devido a 
esquemas eficientes disponíveis para a resolução de tal sistema. 
 
Condições de Contorno − O sistema presente é sujeito a simples condições de 
contorno e x1=0. Introduzindo essas condições, e aplicando o esquema de remuneração, 
reduz-se o sistema para (k´s=1): 
 
�
�
�
�
�
�
�
�
�
�
�
�
�
�
=
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
	
�
�
�
�
�
−
−−
−−
−
Fx
x
x
x
0
0
0
11
121
121
12
5
4
3
2
 
 
O sistema está agora na forma da Eq. (3.56), e está pronto para ser resolvido. 
 
Embora a redução das equações é, certamente, uma valiosa aproximação 
incorporada nas condições de contorno, usualmente, prefere-se deixar o número de 
equações intactas, quando a solução é realizada por computador. Neste caso, temos de 
F1 = -F (reação de F da parede) e x1=0. Assim fica o sistema: 
 
�
�
�
�
��
�
�
�
�
�
�
�
��
�
�
�−
=
�
�
�
�
��
�
�
�
�
�
�
�
��
�
�
�
�
�
�
�
�
�
	
��
�
�
�
�
�
−
−−
−−
−−
−
F
F
x
x
x
x
x
0
0
0
11
121
121
121
11
5
4
3
2
1
 
 
Uma vez incorporadas as condições de contorno, muda-se ao próximo passo: a 
solução. 
 
Gerando a Solução − Resolvendo o sistema linear com uma das técnicas 
numéricas, onde todos os k´s = 1 e F = 1, temos: 
 
X1 = 0, x2 = 1, x3 = 2, x4 = 3 e x5 = 4. 
 
Apresentação dos Resultados (Pós-processamentos) − O resultado pode agora ser 
desenhado graficamente. Na Figura A.3, os resultados estão como esperados. Cada mola 
alongada é uma unidade de deslocamento. 
 
 35
 
Figura A.3 (a) O diagrama do sistema original. (b) O sistema depois da aplicação da força 
constante. Os deslocamentos são indicados no espaço entre os dois sistemas. [54] 
 
 
A.2. Exemplo de uma Haste Sendo Aquecida 
 
A figura A.4 mostra um sistema modelado pela equação de Poisson, na forma 
unidimensional: 
 
)(
2
2
xf
dx
Td
−= (A.1) 
 
Onde, f(x) é uma função que define a fonte de calor ao longo da haste, e onde, as 
pontas da haste, são mantidas numa temperatura fixa 
 
T(0, t) = T1 e T(L, t) = T2. 
 
Essa não é uma equação diferencial parcial, mas uma equação diferencial ordinária com 
condições de contorno. Esse modelo simples é usado, porque ele permitirá introduzir a 
aproximação por elementos finitos, sem algumas das complicações, como por exemplo, 
um PDE de duas dimensões. 
 
Figura A.4 (a) Uma longa e fina haste sujeita a fixas condições de contorno e a uma 
contínua fonte de calor ao longo do seu eixo. (b) A representação de elementos finitos 
consistindo de quatro elementos de igual comprimento e cinco nós.[54] 
 
Exemplo: solução analítica para uma haste sendo aquecida. 
Descrição do Problema − Resolver a Eq. (A.1) para uma haste de 10 cm, com as 
seguintes condições de contorno: T(0, t) = 40 oC e T(10, t) = 200 oC e uma fonte de calor 
 36
uniforme f(x) = 10. 
 
Solução − A equação a ser resolvida é: 
 
10
2
2
−=
dx
Td
 
 
Assume-se uma solução na forma: 
 
T = a x2 + b x + c 
 
Ela é diferenciada duas vezes, resultando T´´ = 2 a. Substituindo este na equação 
diferencial, temos: a = -5. As condições de contorno são usadas para avaliar os 
coeficientes restantes. Da primeira condição, para x=0: 
 
40 = -5 (0)2 + b(0) + c ou c = 40. 
 
Similarmente, para a segunda condição: 
 
200 = -5(10)2 + b(10) + 40 
 
A qual é solucionada dando b = 66. Portanto, a solução final é: 
 
T = -5 x2 + 66 x + 40 
 
O resultado é mostrado na figura A.5 
 
Figura A.5 - Distribuição de temperatura ao longo de uma haste aquecida por uma fonte 
uniforme de calor e mantida fixa a temperatura nos extremos da haste. 
 
 
Temperatura da Haste
0
40
80
120
160
200
240
280
x 5
x (cm)
T
 (
G
ra
u
s 
C
el
su
s)
 37
A.3. Discretização 
 
Uma configuração simples modeladora do sistema é uma série de elementos de 
igual comprimento (Fig. A.4b). Assim, o sistema é tratado com quatro elementos de igual 
comprimento e cinco nós. 
 
 
A.4. Equações dos Elementos 
 
Um elemento individual é mostrado na Fig. A.6a. A distribuição de temperatura 
para o elemento é representada pela função de aproximação: 
 
2211
~
TNTNT += (A.2) 
 
Onde N1 e N2 são funções lineares de interpolação especificadas pelas Eq. (3.48) e 
(3.49), respectivamente. Assim, como descrito na Fig. A.6b, a função de aproximação é 
uma interpolação linear entre as duas temperaturas dos nós. 
 
 Nó 1 Nó 2 
 
 (a) 
 
 T
~
 
 T2 
 
 T1 
 
 
 
 
 
 x1 x2 
 (b) 
 
Figura A.6 (a) Um elemento individual (b). A função aproximação usada para caracterizar 
a distribuição de temperatura ao longo do elemento. 
 
 
A.5. Aproximação Direta 
 
Como descrito no subitem 3.2.1, há uma variedade de aproximações para 
desenvolver as equações elementos. Nesse subitem, emprega-se duas delas. Primeiro, 
uma aproximação direta será usada em um caso simples, onde f(x)=0. Então, em função 
de sua aplicabilidade geral na Engenharia, devotar-se-á a maior parte da abordagem para 
o método dos resíduos ponderados. 
 
 
No caso onde f(x)=0, o método direto é empregado para gerar as equações 
elementos. O relacionamento entre o fluxo de calor e o gradiente de temperatura é assim 
representado pela lei de Fourier: 
 
dx
dT
kq ′−= 
 38
 
Nesta equação, q é fluxo [cal/(cm2. s)] e k´ é o coeficiente de condutividade térmica 
[cal/(s.cm.oC)]. Se uma função linear de aproximação é usada na caracterização da 
temperatura do elemento, o fluxo de calor para o elemento, através do nó 1, é 
representado por: 
 
12
21
1 xx
TT
kq
−
−
′= 
 
Onde q1 é o fluxo de calor no nó 1. Similarmente, com o nó 2: 
 
12
12
2 xx
TT
kq
−
−
′= 
 
Estas duas equações expressam o relacionamento da distribuição da temperatura 
interna do elemento (refletido pelas temperaturas dos nós) e o fluxo de calor nas suas 
extremidades. Assim, elas constituem as equações dos elementos desejadas, e serão 
simplificadas pelo reconhecimento da lei de Fourier, como útil para moldar o fluxo terminal 
em termos do gradiente de temperatura na fronteira. Ou seja: 
 
dx
xdT
kq
dx
xdT
kq
)(
 
)( 2
2
1
1 ′−=′−= 
 
Pode-se substitui-la dentro das equações dos elementos, resultando: 
 
��
�
�
�
��
�
�
�−
=
�
�
�
�
�
�
�
	
�
�
−
−
−
dx
xdT
dx
xdT
T
T
xx )(
)(
11
111
2
1
2
1
12
 (A.3) 
 
A Eq. (A.3) é moldada no formato da Eq. (3.54). Assim, obteve-se sucesso na 
geração da equação matricial onde o comportamento de um elemento típico no sistema é 
descrito. 
 
Outra forma de se chegar às mesmas equações elementos anteriores, é a partir da 
seguinte idéia : 
 
Se T’’(x)=0 então T´(x)= constante. 
 
Assim, pode-se ter: 
 
12
12
xx
TT
dx
dT
−
−
= (A.4) 
 
Aplicando (A.4) aos nós 1 e 2, tem-se: 
 39
�
�
�
�
�
�
�
=+−
−
�
−
−
=
−=−
−
�
−
−
=
2 nó 
)(
)(
)(
1)(
1 nó 
)(
)(
)(
1)(
2
21
1212
122
1
21
1212
121
dx
xdT
TT
xxxx
TT
dx
xdT
dx
xdT
TT
xxxx
TT
dx
xdT
 (A.5) 
 
Chegando-se, assim, a (A.3), e tem-se: 
 
��
�
�
�
��
�
�
�−
=
�
�
�
�
�
�
�
	
�
�
−
−
dx
xdT
dx
xdT
T
T
)(
)(
4,04,0
4,04,0
2
1
2
1 (A.6) 
 
Parte-se, então, para a etapa de montagem: 
 40
( )
( )
( )
( )
( )
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
�
−
=
�
�
�
�
��
�
�
�
�
�
�
�
��
�
�
�
�
�
�
�
�
�
	
�
�
�
�
�
�
�
−
−−
−−