Calculo Numerico

g
8ν , indicando que K e´ pequeno para poros pequenos e alta viscosidade.
O valor da condutividade aproximada K e´ mostrada na tabela 8.5.
228 Cap´ıtulo 8 - Soluc¸a˜o Nume´rica de EDP’s
Tabela 8.5: Ordem de magnitude da condutividade K(m/s) de alguns materiais
material K(m/s)
argila < 109
areia fina 10−5 - 10−4
areia grossa 10−4 - 10−3
cascalho > 10−2
8.8.3 Problemas em aerodinaˆmica
O desenvolvimento de projetos de aerofo´lios e de asas ou he´lices de aeronaves en-
volve uma ana´lise do desempenho aerodinaˆmico dos mesmos. Em geral, esta ana´lise
e´ inicialmente feita atrave´s de me´todos computacionais que simulem a situac¸a˜o real
em que o corpo sera´ submetido, uma vez que experimentos reais envolvem um custo
mais elevado.
Por simplicidade, considere o escoamento potencial em toda a vizinhac¸a do perfil
aero- dinaˆmico. A equac¸a˜o que rege o escoamento potencial e´ de Laplace
∇2Ψ = 0.
A varia´vel Ψ representa o potencial de velocidades, sendo o vetor velocidade obtido
atrave´s da relac¸a˜o
~V = −~∇Ψ
que satisfaz ~∇× ~V = 0, ou seja, o escoamento e´ irrotacional.
Para a resoluc¸a˜o nume´rica destas equac¸o˜es foi necessa´ria a utilizac¸a˜o de um sistema
de coordenadas generalizadas (ξ, η, γ) que coincidem com a geometria da malha.
Estudos sobre as vantagens do uso destas coordenadas mostram que a utilizac¸a˜o de
sistemas de coordenadas ortogonais, como cartesiano, cil´ındrico e esfe´rico apresentam
Fund. de Ca´lculo Nume´rico para Engenheiros 229
dificuldades em problemas com geometria arbitra´ria, ocorrendo a necessidade de
interpolac¸o˜es nas fronteiras para a aplicac¸a˜o das condic¸o˜es de contorno, ocasionando
erros onde os gradientes sa˜o altos, ou seja, justamente nas regio˜es onde temos maior
interesse. Outra vantagem seria a versatilidade que o uso de coordenadas coincidentes
com a fronteira oferece.
Partindo da equac¸a˜o original para o caso bidimensional
∂2Ψ
∂x2
+
∂2Ψ
∂y2
= 0 (8.18)
atrave´s da regra da cadeia, resulta
∂Ψ
∂x
=
∂Ψ
∂η
ηx +
∂Ψ
∂ξ
ξx
∂Ψ
∂y
=
∂Ψ
∂η
ηy +
∂Ψ
∂ξ
ξy
Apo´s nova aplicac¸a˜o da regra da cadeia obte´m-se as expresso˜es para as derivadas de
segunda ordem
Ψxx = ξxxΨξ + ηxxΨη + ξ
2
xΨξξ + η
2
xΨηη + 2 ξx ηxΨξη
Ψyy = ξyyΨξ + ηyy Ψη + ξ
2
y Ψξξ + η
2
xΨηη + 2 ξy ηy Ψξη
sendo a expressa˜o final para a equac¸a˜o a ser implementada dada por
(ξxx + ξyy)Ψξ + (ηxx + ηyy)Ψη + (ξ
2
x + ξ
2
y)Ψξξ +
+(η2x + η
2
y)Ψηη + 2(ξx ηx + ξy ηy)ψξη = 0
Esta equac¸a˜o e´ resolvida numericamente, sendo suas derivadas aproximadas por
diferenc¸as finitas centrais. As me´tricas da transformac¸a˜o (ξx, ξy, ηx, ηy) sa˜o calcula-
230 Cap´ıtulo 8 - Soluc¸a˜o Nume´rica de EDP’s
dos atrave´s de sua interpretac¸a˜o geome´trica
∆xi =
xi+1,j − xi−1,j
2
∆xj =
xi,j+1 − xi,j−1
2
∆yi =
yi+1,j − yi−1,j
2
∆yj =
yi,j+1 − yi,j−1
2
∆ξ =
√
∆x2 +∆y2 ∆η =
√
∆x2 +∆y2
ξx =
∆ξ
∆x
ηx =
∆η
∆x
ξy =
∆ξ
∆y
ηy =
∆η
∆y
onde ∆ξ e´ a metade do comprimento de linha que liga os pontos (i+1, j) e (i−1, j).
As derivadas de segunda ordem ξxx, ξyy, ηxx, ηyy foram aproximadas por diferenc¸as
finitas centrais baseado nestas me´tricas.
Para a soluc¸a˜o da equac¸a˜o (8.18) utiliza-se o me´todo expl´ıcito de Jacobi, na qual
o valor da func¸a˜o Ψ calculada em uma iterac¸a˜o e´ baseada nos valores de Ψi−1,j,
Ψi+1,j, Ψi,j+1, Ψi,j−1, Ψi−1,j−1, Ψi+1,j+1 obtidos na iterac¸a˜o anterior. As relac¸o˜es
de discretizac¸a˜o em diferenc¸as finitas sa˜o dadas como segue
Ψξ =
Ψi+1,j −Ψi−1,j
2∆ξ
Ψη =
Ψi,j+1 −Ψi,j−1
2∆η
Ψξξ =
Ψi+1,j − 2Ψi,j +Ψi−1,j
∆ξ2
Ψηη =
Ψi,j+1 − 2Ψi,j +Ψi,j−1
∆η2
Ψξη =
Ψi+1,j+1 −Ψi+1,j−1 +Ψi−1,j−1 −Ψi−1,j+1
4∆ξ∆η
Apesar de uma certa instabilidade natural do me´todo de jacobi, este procedimento
convergiu para toleraˆncia 10−6. Na Fig. 8.8(a) apresenta-se a malha sobre a qual sa˜o
Fund. de Ca´lculo Nume´rico para Engenheiros 231
obtidas as soluc¸o˜es. As curvas equipotenciais resultantes da equac¸a˜o ∇2Ψ = 0 esta˜o
dispostas na figura 8.8(b). Na Fig. 8.9 mostra-se o campo de velocidades obtido a
partir da func¸a˜o corrente Ψ.
Figura 8.8: Malha de 98 × 41 pontos e linhas de corrente sobre o aerofo´lio
NACA 0012
Figura 8.9: Campo de velocidades do escoamento sobre o aerofo´lio GOE490
A seguir apresenta-se uma introduc¸a˜o ao me´todo de elementos finitos aplicados a
sistemas unidimensionais e a problemas gerais de engenharia.
232 Cap´ıtulo 9 - Int. ao me´todo de elementos finitos
9 INTRODUC¸A˜O AO ME´TODO DE
ELEMENTOS FINITOS
Assim como o me´todo de diferenc¸as finitas, elementos finitos e´ utilizado para aprox-
imar a soluc¸a˜o de problemas de valor inicial e de contorno. Ele foi originariamente
desenvolvido para utilizac¸a˜o em engenharia estrutural, mas atualmente e´ utilizado
para obter soluc¸o˜es de equac¸o˜es diferenciais parciais que surgem em todos os campos
da matema´tica aplicada; sendo mais vantajoso em problemas estruturais.
Segue uma breve introduc¸a˜o sobre como obter os sistemas matriciais em elementos
finitos. Diante da possibilidade de resolver os mesmos problemas por diferenc¸as
ou elementos finitos, os autores na˜o veˆem vantagens em utilizar o me´todo mais
complexo; por isso as aplicac¸o˜es anteriores foram resolvidas em diferenc¸as finitas.
9.1 Sistemas lineares unidimensionais
Considere sistemas massa-mola com duas massas, transfereˆncia de calor entre 2
pontos e de vaza˜o entre dois pontos num duto. No primeiro o deslocamento e´ pro-
porcional a` forc¸a aplicada. A troca de calor e´ dada pela lei de Fourier e a diferenc¸a
de pressa˜o e´ proporcional a` vaza˜o, de forma que
F = Kδ → δ = 1
k
F = CF
q = −K ∂T
∂x
→ q
∫
dx
K
= −dT
p1 − p2 = KQ
implicando em
F1 = Kδ1 − Kδ2
F2 = −Kδ1 + Kδ2
˛˛˛
˛ q1 = −K(T2 − T1)q2 = −K(T1 − T2)
˛˛˛
˛ Q1 = f(P1 − p2)Q2 = f(p2 − p1)
Fund. de Ca´lculo Nume´rico para Engenheiros 233
Assim, na forma matricial[
K −K
−K K
]{
δ1
δ2
}
=
{
F1
F2
} ∣∣∣∣
[
K −K
−K K
]{
T1
T2
}
=
{
q1
q2
}
∣∣∣∣
[
K −K
−K K
]{
p1
p2
}
=
{
Q1
Q2
}
ou ainda,
[K]{δ} = {F}
∣∣∣ [K]{T} = {q} ∣∣∣ [K]{p} = {Q}.
Observe a grande semelhanc¸a na obtenc¸a˜o destes sistemas matriciais.
9.2 Func¸o˜es de interpolac¸a˜o comuns para elementos
lineares, triangulares e tetrae´dricos
Seja a func¸a˜o interpoladora em coordenadas natuarais dada pela seguinte equac¸a˜o:
φ =
∑
Liφi
As func¸o˜es de interpolac¸a˜o para elementos lineares, triangulares e tetrae´dricos para
o me´todo de elementos finitos podem ser obtidos como segue (veja a figura 9.1):
a) Linear
234 Cap´ıtulo 9 - Int. ao me´todo de elementos finitos
Figura 9.1: Reta, triaˆngulo e tetraedro regular
Para uma dimensa˜o, conforme Fig. 9.1a, a func¸a˜o de interpolac¸a˜o e´ obtida fazendo
x = L1x1 + L2x2
1 = L1 + L2
Assim,
L1(x) =
x− x1
x2 − x1 L2(x) =
x2 − x
x2 − x1
b) Triangular
Utilizando coordenadas naturais, segundo a Fig. 9.1b, tem-se
x = L1x1 + L2x2 + L3x3
y = L1y1 + L2y2 + L3y3
1 = L1 + L2 + L3
Desta forma, aplicando o me´todo de Cramer obte´m-se
Fund. de Ca´lculo Nume´rico para Engenheiros 235
L1 =
˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x x2 x3
y y2 y3
˛˛˛
˛˛˛
˛˛˛
˛˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x1 x2 x3
y1 y2 y3
˛˛˛
˛˛˛
˛˛˛
˛
L2 =
˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x1 x x3
y1 y y3
˛˛˛
˛˛˛
˛˛˛
˛˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x1 x2 x3
y1 y2 y3
˛˛˛
˛˛˛
˛˛˛
˛
L3 =
˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x1 x2 x
y1 y2 y
˛˛˛
˛˛˛
˛˛˛
˛˛˛˛
˛˛˛
˛˛˛
˛
1 1 1
x1 x2 x3
y1 y2 y3
˛˛˛
˛˛˛
˛˛˛
˛
ou genericamente
Li =
1
2∆
(ai + bix+ ciy)
sendo o valor dos coeficientes
a1 = x2y3 − x3y2
a2 = x3y1 −