Embed Size (px)
DESCRIPTION
Elementos finitos
Citation preview
Capıtulo 1
Introducao
Varios fenomenos da natureza podem ser descritos, utilizando-se leis da fısica, em termosde equacoes diferenciais. Em algumas situacoes extremamente simples pode-se encontrarsolucoes analıticas, nas outras busca-se solucoes aproximadas para o problema.
Os metodos para a obtencao de solucoes aproximadas de equacoes diferenciais podemser divididos em duas classes gerais:
(i)Metodo de Diferencas Finitas: As derivadas da funcao sao substituıdas pela diferencados valores que a funcao assume em certos pontos do domınio. A equacao diferencial daorigem a um sistema de equacoes algebricas onde as incognitas sao os valores da funcao nospontos escolhidos para representar as derivadas. As ideias basicas de Metodos de DiferencasFinitas sao mostradas no exemplo a seguir:
Exemplo 1:
Resolver: xdudx
+ u = 0, para 0 < x < 1,c.c: u(0) = 1.
A solucao sera aproximada atraves do calculo do valor da funcao em um numero finitode pontos do domınio.
Se o domınio for dividido em N − 1 intervalos atraves de N pontos (nos), deve-se obtero valor aproximado de u(x) em cada no, ou seja, um problema com N incognitas. As Nequacoes necessarias para resolver este problema sao obtidas aproximando a derivada emfuncao dos valores nodais desconhecidos nas posicoes dos N nos.
Olhando nas proximidades do no 4:
du
dx
∣∣∣∣x4
=U5 − U3
2∆x.
A equacao diferencial no ponto x4 e aproximada como:
x4
[U5 − U3
2∆x
]+ U4 = 0 ⇒ − x4
2∆xU3 + U4 +
x4
2∆xU5 = 0.
1
Figura 1.1: Aproximacao da derivada da funcao por diferencas finitas.
O resultado deste procedimento sera um sistema linear N × N , onde as incognitas sao osvalores de U nos N nos.
(ii)Metodos Variacionais: A solucao aproximada e escrita como combinacao linear defuncoes base apropriadamente escolhidas:
u(x) =N∑
i=1
ciφi.
Os coeficientes ci sao obtidos de forma que a equacao diferencial seja satisfeita. Paraobter a solucao aproximada atraves de um metodo variacional, deve-se rescrever o problemaem sua forma fraca (integral ou variacional). A necessidade do uso da formulacao fracae ilustrada no Exemplo 2. Os diferentes metodos variacionais utilizam diferentes formasintegrais, funcoes base e funcoes peso.
O Metodo de Elementos Finitos se encaixa na classe de metodos variacionais, onde asfuncoes base sao diferentes de zero apenas em uma pequena parte do domınio, como seramostrado a seguir.
A simples substituicao da aproximacao u(x) =∑N
i=1 ciφi na equacao diferencial naogarante a obtencao de N (numero de coeficientes ci) equacoes algebricas linearmente inde-pendentes, por isso em metodos variacionais utiliza-se a forma integral do problema, comomostrado no exemplo a seguir.
Exemplo 2:
− d
dx(x
du
dx) + u = 0 para 0 < x < 1,
sujeito as condicoes de contorno:
u(0) = 1 ; (xdu
dx)x=1 = 0.
2
Solucao aproximada na forma:
u(x) = φ0(x) +2∑
i=1
ciφi(x),
onde as funcoes base escolhidas sao:
φ0(x) = 1,
φ1(x) = x2 − 2x,
φ2(x) = x3 − 3x.
Precisamos encontrar os coeficientes c1 e c2, o problema possui 2 incognitas. Entaosubstituindo φ0, φ1 e φ2, a solucao aproximada torna-se:
u(x) = 1 + c1[x2 − 2x] + c2[x
3 − 3x], derivando
du(x)
dx= 2c1x − 2c1 + 3c2x
2 − 3c2.
Observe que as condicoes de contorno sao satisfeitas para qualquer valor de c1 e c2, quedevem ser determinados de forma que a aproximacao satisfaca a equacao diferencial emalgum sentido.
Obrigando a solucao aproximada a satisfazer a equacao diferencial no sentido exato, istoe, em todos os pontos do domınio:
− d
dxx[c1(2x − 2) + c2(3x
2 − 3)] + 1 + c1(x2 − 2x) + c2(x
3 − 3x) = 0,
− 2c1(x − 1) − 3c2(x2 − 1) − 2c1x − 6c2x
2 + c1(x2 − 2x) + c2(x
3 − 3x) + 1 = 0,
1 + 2c1 + 3c2 + (−6c1 − 3c2)x + (−9c2 + c1)x2 + c2x
3 = 0.
Para essa expressao ser valida para qualquer x, o coeficiente de cada potencia de x deveser nulo:
1 + 2c1 + 3c2 = 0,
−6c1 − 3c2 = 0,
−9c2 + c1 = 0,
c2 = 0.
Esse sistema nao tem solucao! Vamos obrigar a solucao aproximada a satisfazer aequacao diferencial num sentido integral ao inves de no sentido exato.
Podemos obriga-la a satisfazer a seguinte integral ponderada:∫ 1
0
w(x)− d
dx(x
du
dx) + udx = 0,
onde w(x) e a funcao peso.Esta integral pode representar uma medida do erro na aproximacao. O numero de
equacoes linearmente independentes que envolvam c1 e c2 sera igual ao numero de funcoes
3
peso utilizadas (tambem linearmente independentes). No exemplo escolhemos 2 funcoespeso linearmente independentes: w = 1 e w = x , teremos:
∫ 1
0
1
(− d
dxx[c1(2x − 2) + c2(3x
2 − 3)] + 1 + c1(x2 − 2x) + c2(x
3 − 3x)
)dx = 0,∫ 1
0
x
(− d
dxx[c1(2x − 2) + c2(3x
2 − 3)] + 1 + c1(x2 − 2x) + c2(x
3 − 3x)
)dx = 0.
Entao: 2/3c1 + 5/4c2 = 1,
3/4c1 − 31/20c2 = 1/2,
e finalmente, c1 = 222/23 e c2 = −100/23.
As funcoes bases utilizadas e a solucao aproximada estao ilustradas na figura (??):
Figura 1.2: Representacao das funcoes bases e solucao aproximada.
4
Capıtulo 2
Formulacao Integral
2.1 Formulacao forte e fraca de um problema
A formulacao forte de um problema de valor de contorno e definida como:Determine u : [0, 1] → R, tal que
(S)
d2udx2 + f = 0, em (0, 1);
u(1) = g → c.c. de Dirichlet;du(0)
dx= −h → c.c. de Newman.
A equacao diferencial e imposta em todos os pontos x ∈ (0, 1).A formulacao fraca equivalente a (S) e definida como:Determine u ∈ U , tal que
(W )
∫ 1
0
dw
dx
du
dxdx =
∫ 1
0
fwdx + w(0)h ; para qualquer w ∈ V .
U e o conjunto de funcoes teste definido como:
U = u|u(1) = g e
∫ 1
0
(du
dx)2dx < ∞
V e o conjunto de funcoes peso e e definido como:
V = w|w(1) = 0︸ ︷︷ ︸∗∗
e
∫ 1
0
(dw
dx)2dx < ∞
** A funcao peso vale zero quando aplicada no contorno onde a condicao e essencial.Isso porque neste ponto a funcao ja e conhecida.
Antes de continuarmos com a definicao de formulacao fraca e mostrarmos a equivalenciaentre as duas formulacoes precisaremos das definicoes:
• Uma funcao f : Ω → R pertence ao espaco Ck(Ω) se djfdxj e contınua, com 0 j k.
C0 e o espaco das funcoes contınuas.• O espaco Sobolev de funcoes Hk(Ω) e definido como:
5
Hk(Ω) = w|w ∈ L2;dw
dx∈ L2; ...;
dkw
dxk∈ L2, onde
L2(Ω) = f |∫
Ω
f 2dΩ < ∞.
Logo podemos escrever U e V como:
U = u|u ∈ H1e u(1) = g
V = w|w ∈ H1e w(1) = 0Notas:• Se f ∈ H1 ⇒ f e contınua e limitada.• A condicao de contorno de Dirichlet e imposta na definicao do espaco de funcoes teste
→ c.c. essencial.• A condicao de contorno de Newman e imposta naturalmente na formulacao → c.c.
natural.
Temos que mostrar que as duas formulacoes sao equivalentes:(i) Se u e solucao de (S) ⇒ u e solucao de (W ).(ii) Se u e solucao de (W ) ⇒ u e solucao de (S).
A demonstracao da parte (i), na realidade, corresponde ao metodo de achar a formulacaofraca equivalente a uma dada formulacao forte. Se u e solucao de (S), entao u satisfaz aequacao diferencial em todos os pontos x ∈ (0, 1) ⇒
⇒∫ 1
0
w[d2u
dx2+ f ]dx = 0 (w ∈ H1(0, 1))
Fazendo uma integracao por partes:
⇒ −∫ 1
0
dw
dx
du
dxdx + w(1)
du(1)
dx− w(0)
du(0)
dx+
∫ 1
0
wfdx = 0.
Como w ∈ V ⇒ w(1) = 0.
Como u e solucao de (S) ⇒ du(0)dx
= −h.
⇒∫ 1
0
dw
dx
du
dxdx = w(0)h +
∫ 1
0
wfdx, ∀w ∈ V
Logo, se u e solucao de (S), u satisfaz (W ).
Demonstracao da parte (ii):Se u e solucao de (W ), entao u ∈ U ⇒ u(1) = g e∫ 1
0
dw
dx
du
dxdx = w(0)h +
∫ 1
0
wfdx, ∀w ∈ V
6
Integrando por partes obtem-se:
−∫ 1
0
wd2u
dx2dx + w(1)
du(1)
dx− w(0)
du(0)
dx=
∫ 1
0
wfdx + w(0)h ⇒
⇒∫ 1
0
w[d2u
dx2+ f ]dx − w(1)
du(1)
dx+ w(0)[
du(0)
dx+ h] = 0.
Como w ∈ V ⇒ w(1) = 0
⇒∫ 1
0
w[d2u
dx2+ f ]dx + w(0)[
du(0)
dx+ h] = 0 ; ∀w ∈ V. (2.1)
Para provar que u e solucao de (S), temos que mostrar que:
(i)d2u
dx2+ f = 0 , para qualquer x ∈ (0, 1) e
(ii)du(0)
dx= −h.
Como a equacao (??) e valida para qualquer w ∈ V , vamos escolher um w de forma aconcluirmos que d2u
dx2 + f = 0:
w(x) = φ(x)[d2u
dx2+ f ]; com φ(0) = φ(1) = 0 e φ(x) > 0, x ∈ (0, 1)
Assim, w(0) = 0, pois,
w(0) = φ(0)[d2u
dx2+ f ], como φ(0) = 0 → w(0) = 0 e a equacao (2.1), torna-se:∫ 1
0
φ(x)[d2u
dx2+ f ]2dx = 0
sendo φ(x) > 0 (por definicao) e [d2u
dx2+ f ]2 ≥ 0,
entao a integral sera zero ⇔ [d2u
dx2+ f ]2 = 0 ⇒ d2u
dx2+ f = 0.
Assim temos (i). Vamos ver (ii):
Usando (i) a expressao (??) fica simplesmente:
w(0)[du(0)
dx+ h] = 0 ; ∀w ∈ V.
Como o espaco V nao restringe o valor de w(0), entao
w(0)[du(0)
dx+ h] = 0 ⇔ du(0)
dx+ h = 0 ⇒ du(0)
dx= −h
Com (i) e (ii) verificados, temos que se u e solucao de (W ) entao e tambem solucao de(S).
7
Como mencionado anteriormente, o metodo de Elementos Finitos e baseado na for-mulacao fraca do problema.
A ideia basica e restringir os espacos U e V para espacos de dimensao finita. Isto implicaque qualquer u ∈ U pode ser escrito como:
u(x) =N∑
i=1
ciφi,
onde os φ′is formam uma base de U .
A fim de simplificar, e comum seguir a notacao:
a(u,w) =
∫ 1
0
dw
dx
du
dxe (f, w) =
∫ 1
0
fwdx,
logo, a formulacao torna-se:Achar u ∈ U , tal que:
a(u,w) = (f, w) + w(0)h ; ∀w ∈ V.
2.2 Metodo dos Resıduos Ponderados
Neste curso vamos nos concentrar no Metodo dos Resıduos Ponderados (Weighted Resid-uals).
Vamos novamente analisar o problema:
d2udx2 = −f, em (0, 1);
u(1) = g = 0 considerar g = 0, apenas para simplificar.du(0)
dx= −h.
Queremos encontrar uma aproximacao uh para a solucao do problema u.uh vai pertencer a um espaco de dimensao finita Uh, logo podemos escreve-lo como:
uh(x) =N∑
i=1
ciφi,
onde N e a dimensao do espaco e φ′is as funcoes base de Uh. Queremos tambem que
Uh ⊂ U logo, uh(1) = 0.Como uh e apenas uma aproximacao de u, d2u
dx2 = −fVamos chamar de resıduo (R) da equacao, a medida de quanto uh aproxima adequada-
mente o problema.
R =d2uh
dx2+ f,
Observe que se uh = u entao, R = 0.
8
Metodos dos Resıduos Ponderados sao obtidos atraves da seguinte integral:∫ 1
0
wR = 0
Observe que integrando o resıduo com o peso w, estamos fazendo com que sua mediaponderada seja 0. Por isso a condicao e dita FRACA, pois nao estamos obrigando ponto aponto que u = uh e sim na media (na integral).
Substituindo R, temos:
⇒∫ 1
0
w[d2uh
dx2+ f ]dx = 0
integrando por partes:
⇒ −∫ 1
0
dw
dx
duh
dxdx +
∫ 1
0
fwdx + w(0)h = 0 → Equivalente a formulacao variacional.
A integracao por partes alivia as restricoes impostas a funcao aproximada uh. A ex-pressao passa a depender da primeira derivada de uh,
duh
dxe nao mais da segunda d2uh
dx2 .
Logo, a restricao sobre uh devera ser que sua derivada pertenca a L2,duh
dx∈ L2, isto e, o
quadrado de duh
dxseja integravel. Observe que ate funcoes contınuas com descontinuidades
nas derivadas podem satisfazer a estas condicoes, como ilustrada na figura(??):
Figura 2.1: Funcoes contınuas com derivadas descontınuas.
Como uh ∈ Uh ⇒ uh(x) =∑N
i=1 ciφi.
9
E entao, duh
dx=∑N
i=1 cidφi
dx
Sendo o espaco Uh de dimensao N precisaremos de N equacoes linearmente indepen-dentes a fim de determinar os coeficientes ci.
As funcoes peso w devem pertencer a um espaco finito Vh ⊂ V . Esse espaco pode serrepresentado pela base ψi. Cada equacao esta relacionada com as funcoes peso ψi.
Deste modo as N equacoes L.I. para determinar os ci sao:
Ri = −∫ 1
0
dψi
dx
(N∑
j=1
cjdφj
dx
)dx +
∫ 1
0
fψidx + ψi(0)h = 0 ⇒
⇒N∑
j=1
∫ 1
0
dψi
dx
dφj
dxdx
︸ ︷︷ ︸
Aij
cj =
∫ 1
0
fψidx + ψi(0)h︸ ︷︷ ︸bi
; i = 1, ..., N.
Temos matricialmente⇒ Ac = b
Algumas referencias chamam A de matriz de rigidez.
O Metodo dos Resıduos Ponderados recebe diferentes denominacoes de acordo com asfuncoes base φi e peso ψi.
METODO DE GALERKIN
O espaco das funcoes peso e identico ao espaco das funcoes teste → ψi = φi
⇒ Matriz A e simetrica.
METODO DE PRETROV-GALERKIN
ψi = φi ⇒ Matriz A nao e simetrica.
METODO DOS MINIMOS QUADRADOS
Os coeficientes ci da aproximacao uh(x) =∑N
i=1 ciφi sao obtidos atraves da minimizacaoda integral do quadrado dos resıduos:
∂
∂ci
∫ 1
0
R2dx = 0 ⇒∫ 1
0
∂R
∂ci
Rdx = 0.
φi =∂R
∂ci
.
METODO DE COLOCACAO
10
O resıduo R e identicamente nulo em N pontos do domınio.
R(xi) = 0.
Tambem e um caso particular de metodo dos Resıduos Ponderados.
wi = δ(x − xi) ⇒∫ 1
0
δ(x − xi)Rdx = 0 ⇒ R(xi) = 0.
EXERCICIO #1
Considere o problema:
−d2udx2 − u + x2 = 0, onde 0 ≤ x ≤ 1;
u(0) = 0,
u′(1) = 1.
Vamos assumir que uma solucao aproximada uh possa ser escrita como:
u(x) =N∑
i=1
ciφi onde, φi(0) = 0 ;
φ1(x) = x ,
φ2(x) = x2 ,
φ3(x) = x3 .
Detemine a solucao para o problema utilizando:
Metodo de Galerkin, ou seja , φi = ψi.
Compare a solucao aproximada com a exata do problema:
u(x) =2cos(−1 + x) − sen(x)
cos(1)+ x2 − 2.
Solucao do Exercıcio #1
R(x) = −d2uh
dx2− uh + x2.
Queremos que:∫ 1
0
w(x)R(x)dx = 0, onde w(x) e o vetor das funcoes peso ψi.
∫ 1
0
ψi[−d2uh
dx2− uh + x2]dx = 0
−∫ 1
0
ψid2uh
dx2dx −
∫ 1
0
ψi uh dx +
∫ 1
0
ψi x2 dx = 0
11
Integrando por partes, evitando assim a segunda derivada,temos:
duh
dxψi
∣∣∣∣10
−∫ 1
0
(dψi
dx
duh
dx)dx −
∫ 1
0
ψi uh dx +
∫ 1
0
ψi x2 dx = 0
Usando que:
•u′h(1) = 1.
•ψi(0) = 0.
−ψi(1) +
∫ 1
0
(dψi
dx
duh
dx)dx −
∫ 1
0
ψi uh dx +
∫ 1
0
ψi x2 dx = 0
Como u(x) =∑N
i=1 cjφj , entao dudx
=∑N
i=1 cjφj
dx; substituindo temos:
N∑i=1
cj
∫ 1
0
[dψi
dx
dφj
dx− ψiφj
]dx︸ ︷︷ ︸
Aij
= ψi(1) −∫ 1
0
ψi x2 dx︸ ︷︷ ︸bi
Pelo Metodo de Galerkin φi = ψi; entao,
Aij =∫ 1
0
[dφi
dx
dφj
dx− φiφj
]dx e bi = ψi(1) − ∫ 1
0ψi x2 dx
Assim, usandoφ1(x) = x ; ⇒ dφ1
dx= 1
φ2(x) = x2 ;⇒ dφ1
dx= x
φ3(x) = x3 ;⇒ dφ1
dx= 3x2
Calculamos a matriz A e o vetor b:
A =
23
34
45
34
1715
43
45
43
5835
b =
34
45
56
Fazendo Ac = b
Encontramos c(1) = 22801777
; c(2) = −2031777
; c(3) = −1757108
.
12
c =
22801777
−2031777
−1757108
0
0.2
0.4
0.6
0.8
1
1.2
0 0.2 0.4 0.6 0.8 1
line 1line 2
Figura 2.2: → uh e → u
O grafico mostra a solucao exata e a aproximada. Observe que a solucao aproximada rep-resenta muito bem a solucao exata, utilizando-se apenas de tres funcoes base.E evidente quea escolha de funcoes base apropriadas e fundamental para o bom desempenho do metodo.
13
Capıtulo 3
Discretizacao e Funcoes Base
3.1 Ponto de Vista Global
Podemos perceber que para obter os elementos da matriz A, precisamos integrar asfuncoes teste e peso (e/ou suas derivadas) em todo o domınio.
O metodo dos elementos finitos consiste na escolha apropriada de funcoes φi e ψi, detal forma que as funcoes e suas derivadas so sao diferentes de zero em uma pequena partedo domınio. Deste modo, a integral no domınio fica reduzida a uma integral na parte dodomınio onde a funcao e diferente de zero, como mostrado na figura(1.1) :
Figura 3.1: Exemplo de funcao base diferente de zero apenas em uma pequena porcao dodomınio.
∫ 1
0
( )φidx =
∫ x2
x1
( )φidx
A ideia e dividir o domınio em elementos, definidos por nos. Por exemplo: o domınio de0 a 1 foi dividido em 8 elementos (9 nos).
14
Figura 3.2: 8 elem. → 9 nos
Uma escolha simples que facilita a computacao e definir φi e ψi de tal forma que:
φi(xj) = δij ⇒ φi(xj) =
1 se i = j,0 se i = j,
conforme mostrado na Fig.1.3 (Elementos Lagrangeanos).Desta forma o numero de funcoes base (dimensao do espaco) esta associado ao numero
de nos da malha. Quanto mais refinada a malha, maior sera a dimensao do espaco defuncoes e desta forma, mais precisa sera a solucao aproximada.
Figura 3.3: φi(xj) = δij
Se as funcoes teste φi sao definidas da maneira mostrada anteriormente, os coeficientesci passam a ter um significado especial :
15
u(xi) =N∑
j=1
cj φj(xi)︸ ︷︷ ︸δij
=N∑
j=1
cjδij = ci
⇒ ci = u(xi)
O coeficiente ci e igual ao valor da funcao u no ponto xi.
Uma outra consequencia importante da definicao das funcoes teste e que varios elemen-tos da matriz A sao iguais a zero, conforme ilustrado na Figura 1.4.
Por exemplo:
A26 =
∫ 1
0
dφ2
dx
dφ6
dxdx = 0
Figura 3.4: Calculo de A26.
A matriz A, do exemplo, passa a ser uma matriz esparsa:
16
A =
× × 0 0 0 0 0 0 0× × × 0 0 0 0 0 00 × × × 0 0 0 0 00 0 × × × 0 0 0 00 0 0 × × × 0 0 00 0 0 0 × × × 0 00 0 0 0 0 × × × 00 0 0 0 0 0 × × ×0 0 0 0 0 0 0 × ×
3.2 Ponto de Vista Elementar
Ate agora vimos o metodo dos elementos finitos de uma forma global, isto e, cada funcaoφi e definida em todo o domınio com a particularidade de ser diferente de zero somente emuma pequena parte dele.
Como cada funcao φi so e diferente de zero em uma pequena parte do domınio, e maiseficiente se a tratarmos nestas opcoes apenas. Isto leva a uma descricao local do metodo(elementar).
No exemplo mostrado anteriormente (Figura 1.3), existem duas funcoes diferentes de zeroem cada elemento e uma mesma funcao e diferente de zero em dois elementos, conformemostrado na tabela a seguir:
elemento no inicial no final φ′is = 0
1 1 2 1-22 2 3 2-33 3 4 3-44 4 5 4-55 5 6 5-66 6 7 6-77 7 8 7-88 8 9 8-9
Podemos pensar em uma matriz elementar A(e), de modo que a matriz global A sejaescrita como:
A =⋃8
e=1 A(e)
A maneira com que essa “montagem” ou uniao e feita sera discutida a seguir.Como so existem dois φ′
is = 0 em cada elemento, a matriz elementar A(e) sera umamatriz 2 × 2.
Por exemplo, a matriz do elemento 3, A(3), e igual a:
A(3) =
∫ x4
x3
dφ3
dxdφ3
dxdx
∫ x4
x3
dφ3
dxdφ4
dxdx
∫ x4
x3
dφ4
dxdφ3
dxdx
∫ x4
x3
dφ4
dxdφ4
dxdx
17
A forma de cada funcao nos elementos e sempre a mesma, independentemente do ele-mento em questao. Isto sugere uma troca de coordenadas para que possamos fazer umadescricao elementar das funcoes.
Descricao Elementar: O elemento se extende de ξ = −1 a ξ = +1 e possui dois nos comomostrado na figura:
Figura 3.5: Funcoes base lineares no sistema de coordenadas local.
As funcoes base elementares sao escritas como:
φ(e)1 (ξ) = 1−ξ
2φ
(e)2 (ξ) = 1+ξ
2
Assimdφ
(e)1
dξ= −1
2;
dφ(e)2
dξ= 1
2
Todos os elementos do domınio sao mapeados para esse elemento padrao.
Por exemplo, o mapeamento do elemento 3, que vai de x3 a x4, para o elemento padrao,de ξ = −1 a ξ = 1, pode ser escrito como:
x = x3 ↔ ξ = −1,
→
ξ(x) = 2x−(x3+x4)x4−x3
,
x(ξ) = (x4−x3) ξ+x3+x4
2,
x = x4 ↔ ξ = 1.
18
As integrais sao reescritas como:
A(3)12 =
∫ x4
x3
dφ3
dx
dφ4
dxdx =
∫ 1
−1
dφ(e)1
dξ
dξ
dx
dφ(e)2
dξ
dξ
dx
dx
dξdξ =
=
∫ 1
−1
dφ(e)1
dξ
dφ(e)2
dξ
dξ
dx︸︷︷︸= 2
x4−x3= 2
h(e)
dξ,
onde h(e) e o tamanho do elemento.
⇒ 2
h(3)
∫ 1
−1
dφ(e)1
dξ
dφ(e)2
dξdξ.︸ ︷︷ ︸
=para todos os elementos
A matriz do terceiro elemento entao fica:
A(3) =2
h(3)
∫ 1
−1dφ1
dξdφ1
dξdξ
∫ 1
−1dφ1
dξdφ2
dξdξ
∫ 1
−1dφ2
dξdφ1
dξdξ
∫ 1
−1dφ2
dξdφ2
dξdξ
,
onde os A(3)ij sao:
A(3)11 =
2
h(3)
∫ 1
−1
(−1
2)(−1
2) dξ =
2
h(3)
1
42 =
1
h(3)
A(3)12 =
2
h(3)
∫ 1
−1
(−1
2)(
1
2) dξ =
−1
h(3)= A
(3)2,1
A(3)22 =
2
h(3)
∫ 1
−1
(1
2)(
1
2) dξ =
1
h(3)
A(3) =1
h(3)
[1 −1
−1 1
].
Montagem da Matriz Global
Precisamos obter a matriz global a partir das matrizes elementares. Para isso, precisamossaber a correspondencia entre a numeracao local dos nos no elemento com a numeracaoglobal dos nos no domınio.
Logo, vamos definir a matriz DomNodeID (
#local dos nos︷ ︸︸ ︷ilnode , iele︸︷︷︸
#do elemento
) que fornece a nu-
meracao global do no local ilnode do elemento iele.
DomNodeID
19
ilnodeiele 1 2 3 4 5 6 7 81 1 2 3 4 5 6 7 82 2 3 4 5 6 7 8 9
A matriz global e montada segundo a correspondencia entre a numeracao local e a global.Por exemplo, a matriz local do elemento 3 (2 × 2) e “montada” na seguinte posicao da
matriz global:
A =
× ×× ×
DomNodeID(1, 3) = 3,
→
A(3)11 → A33
A(3)12 → A34
A(3)21 → A43
A(3)22 → A44
DomNodeID(2, 3) = 4.
O processo de montagem tambem serve para o vetor f
f(3)1 → f3
f(3)2 → f4
Um algoritmo geral para a solucao por elementos finitos pode ser escrito como:
do iele = 1, NELE NELE → # de elementos
call getelem A → rotina para calcular a matriz A(iele) (local).
call getelem f → rotina para calcular o vetor local f (iele).
do ilnode = 1, NLOCALNODES
ignode = DomNodeID (ilnode, iele) → # global do no ilnode
do jlnode = 1, NLOCALNODES
jgnode = DomNodeID (jlnode, iele)A(ignode, jgnode) = A(ignode, jgnode) + Aelem(ilnode, jlnode)
end
end
end
20
Nota : Repare que a entrada A55 da matriz global tem contribuicao de 2 elementos(elementos 4, 5).
A55 = A(4)22 + A
(5)11
Isto ocorre, pois:
A55 =
∫ 1
0
dφ5
dx
dφ5
dxdx =
∫ x5
x4
dφ5
dx
dφ5
dxdx +
∫ x6
x5
dφ5
dx
dφ5
dxdx = A
(4)22 + A
(5)11
3.3 Elemento de Segunda Ordem
No exemplo anterior as funcoes base eram lineares. Cada elemento possuia 2 nos econsequentemente duas funcoes (uma associada a cada no). Isso porque a funcao φi foidefinida tal que :
φi(xj) = δi,j,
ou seja, φ1(x1) = 1, φ2(x1) = 0,
φ1(x2) = 0, φ2(x2) = 1.
⇒ dois pontos definem uma reta → φi e linear.
Figura 3.6: Elemento linear
Se desejamos ter uma aproximacao melhor podemos definir funcoes base quadraticasdentro de cada elemento. Como para definir uma funcao quadratica sao necessarios 3 pon-tos, os elementos quadraticos possuem 3 nos.
21
Figura 3.7: Elemento Quadratico
Repare que φi(xj) = δi,j.
Nas coordenadas elementares, as funcoes base sao escritas como:
φ1(ξ) =ξ(ξ − 1)
2,
φ2(ξ) =ξ(ξ + 1)
2,
φ3(ξ) = −(ξ + 1)(ξ − 1).
Logo, derivando-as temos:
dφ1
dξ= ξ − 1
2;
dφ2
dξ= ξ +
1
2;
dφ3
dξ= −2ξ.
Correspondencia entre descricao global e local:DomNodeID
ilnodeiele 1 2 3 4 5 61 1 3 5 7 9 112 3 5 7 9 11 133 2 4 6 8 10 12
As matrizes A(e) sao 3 × 3.
A3 =2
h(3)
∫ 1
−1dφ1
dξdφ1
dξdξ
∫ 1
−1dφ1
dξdφ2
dξdξ
∫ 1
−1dφ1
dξdφ3
dξdξ
∫ 1
−1dφ2
dξdφ1
dξdξ
∫ 1
−1dφ2
dξdφ2
dξdξ
∫ 1
−1dφ2
dξdφ3
dξdξ
∫ 1
−1dφ3
dξdφ1
dξdξ
∫ 1
−1dφ3
dξdφ2
dξdξ
∫ 1
−1dφ3
dξdφ3
dξdξ
22
Figura 3.8: Correspondencia da numeracao global e local.
A Matriz A(3) e montada na seguinte posicao:
A =
× × ×× × ×× × ×
A(3)11 → A55
A(3)12 → A57
A(3)13 → A56
. . .
23
3.4 Problema Unidimensional
Vamos mostrar todos os passos para a solucao de um problema de conveccao/difusaousando o metodo dos elementos finitos.
d2udx2 − Pe du
dx= 0, 0 < x < 1;
du(0)dx
= u(0) − 1;
u(1) = 0.
• Formulacao Fraca (Metodo de Galerkin)
Ri =
∫ 1
0
(d2u
dx2− Pe
du
dx
)φi dx = 0 ⇒
⇒ −∫ 1
0
du
dx
dφi
dxdx +
du
dx(1) φi(1) − du
dx(0)φi(0) − Pe
∫ 1
0
du
dxφi dx = 0
⇒∫ 1
0
(du
dx
dφi
dx+ Pe
du
dxφi
)dx =
du
dx(1) φi(1) − du
dx(0)︸ ︷︷ ︸
u(0)−1
φi(0)
Ri =
∫ 1
0
(du
dx
dφi
dx+ Pe
du
dxφi
)dx =
du
dx(1) φi(1)︸ ︷︷ ︸
=0
−(u(0) − 1)φi(0)
Tendo como c.c u(1) = 0, impoe - se φi(1) = 0.
• MEF
Usar elementos lineares: 8 elementos / 9 nos (N = 9)
u(x) =N∑
i=1
Uiφi → sem se preocupar com as condicoes de contorno.
Ri =N∑
j=1
∫ 1
0
(dφi
dx
dφj
dx+ Pe
dφj
dxφi)dx
Uj = −
[N∑
j=1
Ujφj(0) − 1
]φi(0); (i = 1, ..., N)
Observe que x = 0 e o no 1, entao φ1(0) = 1 e φj(0) = 0,∀j = 1. Assim, temos que:
Ri =N∑
j=1
∫ 1
0
(dφi
dx
dφj
dx+ Pe
dφj
dxφi)dx
Uj = −[U1 − 1]φi(0).
c.c. essencial ⇒ u(1) = 0 ⇒ ∑Nj=1 Ujφj(1) , lembrando que x = 1 e o no N ⇒
UNφN(1) = 0 ⇒ UN = 0.
24
Figura 3.9: Discretizacao do domınio atraves de 8 elementos lineares.
i = 1 ⇒
R1 ⇒N∑
j=1
∫ 1
0
(dφ1
dx
dφj
dx+ Pe
dφj
dxφ1)dx
Uj = −[U1 − 1]
i = 2, ..., N − 1 ⇒
Ri ⇒N∑
j=1
∫ 1
0
(dφi
dx
dφj
dx+ Pe
dφj
dxφi)dx
Uj = 0
i = N ⇒RN ⇒ UN = 0
× . . . ×
.... . .
× ×0 0 0 0 0 0 0 0 1
U1
U2
U3
U4
U5
U6
U7
U8
U9
=
100000000
• Ponto de Vista Local:
25
Formulacao Fraca:∫ 1
0
(du
dx
dφi
dx+ Pe
du
dxφi) dx = −[u(0) − 1] φi(0)︸ ︷︷ ︸
condicao de contorno.
.
→ O ideal e montar as matrizes e vetores elementares sem se preocupar com as condicoesde contorno e depois modifica-las em funcao destas.
A(e) =
[A
(e)11 A
(e)12
A(e)21 A
(e)22
]; f (e) =
[f
(e)1
f(e)2
];
dξ
dx=
2
h(e).
A(e)ij =
2
h(e)
∫ 1
−1
dφj
dξ
dφi
dξdξ + Pe
∫ 1
−1
dφj
dξφidξ; f
(e)i = 0.
φ(e)1 (ξ) =
1 − ξ
2→ dφ1
dξ=
−1
2,
φ(e)2 (ξ) =
1 + ξ
2→ dφ2
dξ=
1
2.
A(e)11 =
2
h(e)
∫ 1
−1
(−1
2) (−1
2) dξ + Pe
∫ 1
−1
(−1
2) (
1 − ξ
2) dξ =
=2
h(e)
∫ 1
−1
1
4dξ −
∫ 1
−1
Pe
4(1 − ξ) dξ =
=1
h(e)− Pe
2e assim por diante...
Apos obter as matrizes elementares, temos que nos preocupar com as condicoes decontorno.
Os termos correspondentes as condicoes de contorno se anulam em todos os termosmenos no primeiro e no ultimo.
φi(1) = 0 , so para i = N
φi(0) = 0 , so para i = 1.
No primeiro elemento, temos que acrescentar:∫ 1
0
( ) = −[U1 − 1] ⇒
⇒N∑
j=1
( ) Uj + U1 = +1
→A(1)1,1 = A
(1)1,1 + 1
f(1)1 = f
(1)1 + 1
26
No ultimo elemento temos uma condicao de contorno essencial.U9 = 0 → Essa e a equacao que se quer impor.
A(8)2,1 ← 0 A
(8)2,2 ← 1
f(8)2 ← 0.
se a c.c. fosse u(1) = g ⇒ f(8)2 ← g.
Para calcularmos os elementos de cada matriz, temos que integrar funcoes ao longo doelemento. De um modo geral, temos que calcular:∫ 1
−1
F (ξ)dξ.
Para um problema simples, podemos calcular a integral analiticamente. Porem, e maisinteressante faze-lo numericamente, pois engloba a resolucao nao so de problemas simplescomo tambem complicados e ainda temos um programa generico.
Um metodo comumente utilizado e o da QUADRATURA GAUSSIANA, onde
∫ 1
−1
F (ξ)dξ =NGP∑i=1
F (ξ)Wi
onde NGP e o numero de pontos utilizados para a integracao. Quanto maior o numero depontos, melhor a aproximacao.
A tabela a seguir fornece o valor de ξi e Wi para diferentes NGP .
NGP ξi Wi
1 0.0 2.0
2 −0.57735 1.0
+0.57735 1.0
−0.77459 0.555
3 0.0 0.888
+0.77459 0.555
PROJETO #1
d2udx2 − Pe du
dx= 0
du(0)dx
= u(0) − 1
u(1) = 0
27
A solucao analıtica do problema e dada por:
u(x) =ePex − ePe
1 − Pe − ePe
Escreva um programa para obter a solucao do problema com:
(a) Malha de 10 elementos lineares uniformemente espacados com Pe = 5, P e =10, P e = 30.
(b) Malha de 5 elementos quadraticos uniformemente espacados com Pe = 5, P e =10, P e = 30.
(c) Malha de 10 elementos lineares concentrados em x = 1, com Pe = 30. Utilize aseguinte funcao para gerar a malha concentrada:
xi = 1 −[
(NNODES − i)
(NNODES − 1)
]a
a → parametro que regula a concentracao da malha.(d) Compare as solucoes obtidas com a solucao analıtica em cada ponto da malha.
ESTRUTURA SUGERIDA PARA O PROGRAMA
CALL INPUT → Rotina para ler os dados do problema.
• Pe
• Tipo de Elemento
linear
Quadratico
• Numero de Elementos
• Malha
Uniforme
Concentrada
CALL GLOBALPOINTERS → rotina para:
•Calcular o numero de nos (NNODES).
•Gerar a matriz com o mapeamento da numerac~ao local → global
(DomNodeID) para cada tipo de elemento.
CALL MESH → Gerar a malha (coordenadas xi).
CALL SOLUTION → Calcular a soluc~ao por EF.
CALL EXACT → Calcular a soluc~ao exata em cada ponto da malha.
CALL OUTPUT → Criar uma tabela com :
xi | sol. exata | sol. aprox. | erro %
28
SUBROTINE SOLUTION
do iele = 1, NELE NELE → # de elementos
call getelemA → rotina para calcular a matriz A(iele) (local).
call getelemb → rotina para calcular o vetor local f (iele).
call BoundCond → rotina para impor condic~oes de contorno.
do ilnode = 1, NLOCALNODES
ignode = DomNodeID (ilnode, iele) → # global do no ilnode
do jlnode = 1, NLOCALNODES
jgnode = DomNodeID (jlnode, iele)A(ignode, jgnode) = A(ignode, jgnode) + Aelem(ilnode, jlnode)
end
end
end
CALL SOLVER
SUBROTINE getelmA
A(e) ← 0do igp = 1,NGP
XI = XGP (igp)CALL BASISFUNC (XI(igp), PHI,DPHIdXI︸ ︷︷ ︸
vetor
)
do ilnode = 1, NLOCALNODESdo jlnode = 1, NLOCALNODES
A(ilnode, jlnode) = A(ilnode, jlnode) + (w(igp) ∗ 2h(e)∗
dPHIdXI(ilnode) ∗ dPHIdXI(jlnode) + Pe ∗ PHI(ilnode)∗dPHIdXI(jlnode))
end do
end do
end do
29
Capıtulo 4
Problema Bi-Dimensional Linear
Vamos considerar o problema 2-D de conducao de calor.∇ · (k∇T ) = 0; x ∈ Ω,
T = 0; x ∈ Γ1,
n · ∇T = q; x ∈ Γ2.
Figura 4.1: Problema bi-dimensional.
k → condutividade termica. Vamos considerar k constante ⇒ ∇2T = 0.
T |Γ1 = 0 → c.c. Dirichlet (temperatura prescrita)n · ∇T = q ⇒ ∂T
∂n|Γ2 = q → c.c Newman (fluxo prescrito)
Formulacao Fraca : Ache T ∈ U tal que:
∫Ω
w∇2TdΩ = 0 , ∀w ∈ V ;
onde:
30
U = T ∣∣ ∫Ω|∇T |2dΩ < ∞ ; T (Γ1) = 0,
V = w ∣∣ ∫Ω|∇w|2dΩ < ∞ ; w(Γ1) = 0.
Desenvolvendo a expressao do resıduo ponderado:
Observe que:
d
dx
[w
dT
dx
]=
dw
dx
dT
dx+ w
d2T
dx2
⇒ wd2T
dx2= −dw
dx
dT
dx+
d
dx
[w
dT
dx
],
usando a mesma idea para duas variaveis:∫Ω
w
[∂2T
∂x2+
∂2T
∂y2
]dxdy = −
∫Ω
[∂w
∂x
∂T
∂x+
∂w
∂y
∂T
∂y
]dxdy +
∫Ω
∂
∂x
[w
∂T
∂x
]+
∂
∂y
[w
∂T
∂y
]dxdy
⇒∫
Ω
w∇2TdΩ = −∫
Ω
∇w · ∇TdΩ +
∫Ω
∇ · (w∇T )dΩ.
Usando o Teorema da Divergencia, que diz:∫Ω
∇ · fdΩ =
∫Γ
n · fdΓ
A formulacao fraca fica:∫Ω
w∇2TdΩ = −∫
Ω
∇w · ∇TdΩ +
∫Γ
n · (w∇T )dΓ
Dividindo o contorno em dois pedacos e como w|Γ1 = 0 (a funcao peso deve ser iden-ticamente nula ao longo do contorno onde e imposta uma condicao de contorno de Dirichlet)
⇒∫
Γ
wn · ∇TdΓ =
∫Γ1
wn · ∇T︸ ︷︷ ︸0
dΓ1 +
∫Γ2
w n · ∇T︸ ︷︷ ︸q
dΓ2
Formulacao Fraca:
− ∫Ω∇w · ∇T dΩ +
∫Γ2
w q dΓ2 = 0
Discretizacao → Uh ⊂ U Uh de dimensao N .
Th =N∑
j=1
Cjφj,
onde C ′js sao as N incognitas.
31
As N equacoes linearmente independentes serao obtidas escolhendo-se N funcoes pesotambem linearmente independentes :
Ri →∫
Ω∇ψi · ∇T dΩ =
∫Γ2
ψi q dΓ2
Sendo Th =∑N
j=1 Cjφj, entao, ∇T =∑N
j=1 Cj∇φj :
∫Ω
∇ψi · ∇(N∑
j=1
Cjφj) dΩ =
∫Γ2
ψi q dΓ2 ⇒
N∑j=1
Cj
∫Ω
∇ψi · ∇φj dΩ
︸ ︷︷ ︸
Aij
=
∫Γ2
ψi q dΓ2︸ ︷︷ ︸bi
Aijcj = bi
Aij =
∫Ω
∇ψi · ∇φj dΩ
bi =
∫Γ2
ψi q dΓ2
Metodo dos Elementos Finitos
Escolher φi tal que φ(xj) = δij ⇒ φi(xj) =
1 se i = j0 se i = j
Figura 4.2: Funcoes base bi-dimensional φ(xj) = δij
Ponto de Vista Elementar
Elemento Bi-Linear
32
Figura 4.3: Funcoes base de elemento bi-linear.
• φ1(ξ, η) =(η − 1)(ξ − 1)
4
• φ2(ξ, η) = −(η − 1)(ξ + 1)
4
• φ3(ξ, η) =(η + 1)(ξ + 1)
4
• φ4(ξ, η) = −(η + 1)(ξ − 1)
4
A(e)ij =
∫Ω(e)
∇φi · ∇φj dΩ
A(e)ij =
∫Ω(e)
[∂φi
∂x
∂φj
∂x+
∂φi
∂y
∂φj
∂y
]dΩ
A matriz A(e) e uma matriz 4 × 4 no caso de elementos bi-lineares.
Precisamos saber as derivadas das funcoes base, ∂φi
∂x; ∂φi
∂y, para i = 1..4, ja que os φ′
issao definidos em termos de ξ e η.
Da mesma forma que no caso 1-D, precisamos de um mapeamento do sistema de coor-denada local (ξ, η) para o sistema de coordenadas global (x, y).
Uma vez obtida as matrizes elementares, temos que “montar” a matriz global. O pro-cesso de montagem e analogo ao caso unidimensional. Precisamos da matriz
DomNodeID(ilnode, iele) = numero global do no local ilnode do elemento iele.
33
Esta matriz e obtida a partir da numeracao dos elementos e nos usados para dividir odomınio de calculo. A numeracao tanto dos nos como dos elementos e aleatoria. Como essanumeracao esta diretamente relacionada com a estrutura da matriz global, e interessanteusar uma numeracao otimizada, que minimize a banda da matriz global.
Figura 4.4: Exemplo de numeracao dos nos e elementos.
A Fig.4.4 ilustra um exemplo de numeracao dos elementos e nos. A matriz DomNodeIDpara este problema e dada por:
DomNodeID
ilnodeiele 1 2 3 4 5 6 7 8 91 1 2 3 5 6 7 9 10 112 5 6 7 9 10 11 13 14 153 6 7 8 10 11 12 14 15 164 2 3 4 6 7 8 10 11 12
34
A Matriz A(5) e montada na seguinte posicao:
A =
× × × ×× × × ×
× × × ×× × × ×
A(5)1 1 → A6 6
A(5)1 2 → A6 10
A(5)1 3 → A6 11
A(5)1 4 → A6 7
. . .
A estrutura da matriz global e funcao direta da numeracao global dos nos, como ja foidito. Um esquema NAO OTIMIZADO pode levar a uma matriz global com banda alta,como mostrado na Fig.4.5. Existem diversos algoritmos para otimizar a numeracao dos nos.Estes metodos nao serao tratados neste curso.
Para calcular elementos da matriz elementar, precisamos calcular as derivadas dasfuncoes base com relacao a x e y:
A(e)ij =
∫Ω(e)
[∂φi
∂x
∂φj
∂x+
∂φi
∂y
∂φj
∂y
]dΩ
Mas, a funcao φ(ξ, η) esta definida em termos de ξ e η, logo temos que mudar decoordenadas e para isso vamos usar um mapeamento que escreva x e y como funcao de ξ eη:
• x = x(ξ, η) e• y = y(ξ, η).
MAPEAMENTO ISOPARAMETRICO:
35
Figura 4.5: Dois esquemas de numeracao de nos e suas respectivas matriz global.
Requisitos para mapeamento:
→ Linhas de contorno do elemento sao mapeadas em linhas de contorno.
→ Nos globais sao mapeados em nos locais.
Um exemplo de mapeamento muito utilizado e definido por:x(ξ, η) =
∑4i=1 Xi φi(ξ, η)
y(ξ, η) =∑4
i=1 Yi φi(ξ, η)
Repare que :
• x1 = X1 φ1(−1,−1)︸ ︷︷ ︸=1
+X2 φ2(−1,−1)︸ ︷︷ ︸=0
+X3 φ3(−1,−1)︸ ︷︷ ︸=0
+X4 φ4(−1,−1)︸ ︷︷ ︸=0
= X1
e assim em diante...
• A funcao interpolacao ultilizada no mapeamento e a mesma funcao interpolacao ul-tilizada para aproximar a funcao T ⇒ elemento isoparametrico.
EXERCICIO #2:Mostre que o mapeamento isoparametrico satisfaz a condicao que a imagem da linha
η = 1 e mapeada no contorno do elemento entre os nos 3 e 4.
36
Figura 4.6: Mapeamento isoparametrico
Uma vez conhecendo o mapeamento entre as coordenadas, pode-se determinar : ∂φi
∂xe
∂φi
∂y.
∂φi
∂ξ= ∂φi
∂x∂x∂ξ
+ ∂φi
∂y∂y∂ξ
∂φi
∂η= ∂φi
∂x∂x∂η
+ ∂φi
∂y∂y∂η
⇒ ∂φi
∂ξ
∂φi
∂η
=
∂x
∂ξ∂y∂ξ
∂x∂η
∂y∂η
︸ ︷︷ ︸Jacobiano da Transformacao
∂φi
∂x
∂φi
∂y
Logo:
∂φi
∂x
∂φi
∂y
= J−1
∂φi
∂ξ
∂φi
∂η
; J−1 =
1
|J|
∂y
∂η−∂y
∂ξ
−∂x∂η
∂x∂ξ
⇒ Inverso do Jacobiano
|J| = ∂x∂ξ
∂y∂η
− ∂y∂ξ
∂x∂η
.
37
Assim:
∂φi
∂x=
1∂x∂ξ
∂y∂η
− ∂y∂ξ
∂x∂η
∂y
∂η
∂φi
∂ξ− ∂y
∂ξ
∂φi
∂η
∂φi
∂y=
1∂x∂ξ
∂y∂η
− ∂y∂ξ
∂x∂η
−∂x
∂η
∂φi
∂ξ+
∂x
∂ξ
∂φi
∂η
Para um mapeamento isoparametrico temos que:
∂x∂ξ
=∑N
i=1 Xi∂φi
∂ξ; ∂x
∂η=∑N
i=1 Xi∂φi
∂η; . . .
Logo, a integral elementar∫
Ωe
[∂φi
∂x
∂φj
∂x+ ∂φi
∂y
∂φj
∂y
]dxdy pode ser reescrita em funcao das
variaveis locais:
A(e)ij =
∫Ω(e)
∇φi · ∇φj dxdy =∫
Ω(e)J−1∇ξηφi · J−1∇ξηφj |J| dξ dη
A(e)ij =
∫ 1
−1
∫ 1
−1
(∂y
∂η∂φi
∂ξ− ∂y
∂ξ∂φi
∂η)(∂y
∂η
∂φj
∂ξ− ∂y
∂ξ
∂φj
∂η)+
+(−∂x∂η
∂φi
∂ξ+ ∂x
∂ξ∂φi
∂η)(−∂x
∂η
∂φj
∂ξ+ ∂x
∂ξ
∂φj
∂η)
1|J|dξdη
Como a expressao para o calculo da matriz elementar contem |J| no denominador, porisso e importante que |J | 0. Se o elemento ficar muito distorcido isto pode acontecer e a
precisao do calculo de A(e)ij fica comprometida, como no caso mostrado a seguir:
Figura 4.7: Elemento distorcido.
dA = |J|dξdη = |dr1||dr2|senα
|J| =|dr1||dr2|senα
dξdη
Quando α → 0 ou α → 180 graus , |J| → 0. Deve-se monitorar o valor de |J|. Seficar negativo ou muito pequeno, os calculos devem ser interrompidos e uma nova malha
38
deve ser criada.
Integracao Numerica
Quadratura Gaussiana:
∫ 1
−1
∫ 1
−1F (ξ, η)dξdη =
∑NGPIigp=1
∑NGPJjgp=1 F (ξi, ηj) wI wJ
Figura 4.8: Localizacao dos Pontos de Gauss.
ELEMENTO NUM.PONTOS PESOS COORDENADAS
INTEGRACAO w ξ η
linear 2 × 2 = 4 1.0 × 1.0 -0.57735 -0.577351.0 × 1.0 +0.57735 -0.577351.0 × 1.0 -0.57735 +0.577351.0 × 1.0 +0.57735 +0.57735
quadratico 3 × 3 = 9 0.555 × 0.555 -0.77459 -0.774590.888 × 0.555 0.0 -0.774590.555 × 0.555 +0.77459 -0.774590.555 × 0.888 -0.77459 0.00.888 × 0.888 0.0 0.00.555 × 0.888 +0.77459 0.00.555 × 0.555 -0.77459 +0.774590.888 × 0.555 0.0 +0.774590.555 × 0.555 +0.77459 +0.77459
39
A Localizacao dos pontos esta representada na Fig.4.8:
Algoritmo para o calculo de A(e)
do igp = 1, NGPXI = XGP (igp)WI = W (igp)do jgp = 1, NGP
NETA = XGP (jgp)WJ = W (jgp)W = WI ∗ WJCalcular func~oes base no ponto de Gauss
φi(ξ, η), ∂φi
∂ξ(ξ, η), ∂φi
∂η(ξ, η) i = 1, ..., NLOCALNODES
Calcular ∂x∂ξ
=∑NLN
i=1 Xi∂φi
∂ξ∂y∂ξ
=∑NLN
i=1 Yi∂φi
∂ξ∂x∂η
=∑NLN
i=1 Xi∂φi
∂η∂y∂η
=∑NLN
i=1 Yi∂φi
∂η
Calcular detJ = ∂x∂ξ
∂y∂η
− ∂y∂ξ
∂x∂η
Calcular ∂φi
∂x= ... ∂φi
∂y= ...
do ilnode = 1, NLOCALNODESdo jlnode = 1, NLOCALNODES
A(ilnode, jlnode) = A(ilnode, jlnode)+
W ∗(
∂φi
∂x∗ ∂φj
∂x+ ∂φi
∂y∗ ∂φj
∂y
)∗ detJ
end do
end do
end do
end do
Para executar os passos de calculo de forma eficiente, pode-se utilizar uma estrutura ma-tricial para armazenar e operar as variaveis. Um exemplo de tal estrutura e descrito a seguir.
• Armazenar as funcoes φi(ξ, η), ∂φi
∂ξ, ∂φi
∂ηda seguinte forma:
PHI(i) =
φ1
φ2
φ3
φ4
GradPHI(i, j) =
∂φ1
∂ξ∂φ2
∂ξ∂φ3
∂ξ∂φ4
∂ξ
∂φ1
∂η∂φ2
∂η∂φ3
∂η∂φ4
∂η
• Armazenar as coordenadas dos pontos nodais de cada elemento da seguinte forma:
XY (i, j) =
[X1 X2 X3 X4
Y1 Y2 Y3 Y4
]
40
• A matriz Jacobiana pode ser calculada por:
JAC(i, j) = GradPHI ∗ XY T ⇒
∂x
∂ξ∂y∂ξ
∂x∂η
∂y∂η
=
∂φ1
∂ξ∂φ2
∂ξ∂φ3
∂ξ∂φ4
∂ξ
∂φ1
∂η∂φ2
∂η∂φ3
∂η∂φ4
∂η
∗
X1 Y1
X2 Y2
X3 Y3
X4 Y4
detJ = JAC(1, 1) ∗ JAC(2, 2) − JAC(1, 2) ∗ JAC(2, 1)
J−1 =
JACINV (1, 1) = JAC(2, 2)/detJ
JACINV (1, 2) = −JAC(1, 2)/detJ
JACINV (2, 1) = −JAC(2, 1)/detJ
JACINV (2, 2) = JAC(1, 1)/detJ
• As derivadas das funcoes base em relacao as coordenadas global sao armazenadas namatriz
GradPHI XY (i, j) =
∂φ1
∂x∂φ2
∂x∂φ3
∂x∂φ4
∂x
∂φ1
∂y∂φ2
∂y∂φ3
∂y∂φ4
∂y
GradPHI XY (i, j) = J−1 ∗ GradPHI
1
|J|
∂y
∂η−∂y
∂ξ
∂x∂η
∂x∂ξ
∂φ1
∂ξ∂φ2
∂ξ∂φ3
∂ξ∂φ4
∂ξ
∂φ1
∂η∂φ2
∂η∂φ3
∂η∂φ4
∂η
=
∂φ1
∂x∂φ2
∂x∂φ3
∂x∂φ4
∂x
∂φ1
∂y∂φ2
∂y∂φ3
∂y∂φ4
∂y
O algoritmo pode ser escrito como:
CALL BASISFUNC (XI,NETA, PHI,GradPHI)CALL MAPPING (XY,PHI,GradPHI, JAC, JACINV, detJ)GradPHI XY (i, j) = J−1 ∗ GradPHIdo ilnode = 1, 4
do jlnode = 1, 4A(ilnode, jlnode) = A(ilnode, jlnode) + W ∗ detJ∗
(GradPHI XY (1, ilnode) ∗ GradPHI XY (1, jlnode)+GradPHI XY (2, ilnode) ∗ GradPHI XY (2, jlnode))
end do
41
end do
EXERCICIO # 3
Escreva a rotina BASISFUNC 1 para elementos bi-lineares e BASISFUNC 2 para ele-mentos quadraticos.
EXERCICIO # 4
Escreva a rotina MAPPING indicada no algoritmo acima.
42
Capıtulo 5
Problema Nao Linear
Vamos considerar o seguinte problema uni-dimensional:
ududx
+ d2udx2 = −f(x); x ∈ (0, 1),
u = 0; x = 0,
u = 0; x = 1.
Desenvolvendo a formulacao do Metodo dos Resıduos Ponderados (Galerkin), obtem-se:
Ri =
∫ 1
0
(udu
dx+
d2u
dx2+ f
)φi dx = 0 ⇒
⇒ −∫ 1
0
du
dx
dφi
dxdx +
du
dx(1) φi(1)︸ ︷︷ ︸
=0
−du
dx(0) φi(0)︸ ︷︷ ︸
=0
+
∫ 1
0
udu
dxφidx +
∫ 1
0
fφidx = 0,
Ri = −∫ 1
0
du
dx
dφi
dxdx +
∫ 1
0
udu
dxφidx +
∫ 1
0
fφidx = 0.
Substituindo a expansao da funcao u(x) =∑N
j=1 Cjφj(x):
Ri = −∫ 1
0
N∑j=1
Cjdφj
dx
dφi
dxdx +
∫ 1
0
(N∑
j=1
Cjφj
)(N∑
k=1
Ckdφk
dx
)φidx +
∫ 1
0
fφidx = 0.
Trocando a ordem da integracao e somatorio em j, a expressao pode ser reescrita como:
N∑j=1
Cj
∫ 1
0
−dφj
dx
dφi
dxdx +
∫ 1
0
(N∑
k=1
Ckdφk
dx
)φjφidx
︸ ︷︷ ︸Aij
=
∫ 1
0
−fφidx.︸ ︷︷ ︸bi
43
A expressao acima representa o produto matricial Ac = b. Porem, como o problema enao linear, os elementos da matriz A dependem dos coeficientes Ck, solucao do problema.A notacao mais precisa de representar o problema acima seria
A(c) c = b.
O sistema de equacoes nao lineares deve ser resolvido por um metodo iterativo.E importante perceber a diferenca entre metodos iterativos para resolver um sistema de
equacoes nao lineares e metodos iterativos para resolver um sistema de equacoes lineares.Muitos programas de simulacao numerica utilizam um metodo iterativo para resolver osistema linear Ac = b oriundo de um problema linear. A vantagem de metodos iterativos ea economia de memoria, pois a matriz A nao e armazenada inteira na memoria, e as grandesdesvantagens sao a nao garantia de convergencia e o tempo computacional. O que estamosanalisando neste capıtulo sao metodos iterativos para resolver um sistema nao linear. Aprincıpio, neste curso, sistemas lineares sao resolvidos por metodos diretos (alguma variacaode eliminacao gaussiana).
5.1 Solucao de problema nao-linear com uma variavel
Vamos considerar o seguinte problema de uma variavel:
x3 − 2x + 4︸ ︷︷ ︸f(x)
= 0
O metodo mais simples de determinar a raiz desta equacao e conhecido como Metodo dePicard.
METODO DE PICARD:
A funcao f(x) e reescrita como f(x) = x − g(x). A solucao do problema, x∗, satisfaz aequacao x∗ = g(x∗).
x0 ← chute inicial
Repetir ate |xn+1 − xn| < εxn+1 = g(xn)
A interpretacao geometrica do metodo de Picard e apresentada na Figura 3.1. O metodode Picard pode divergir ate quando o valor inicial e bem proximo da solucao do problema,como no caso 3 da Figura 3.1 O criterio de convergencia do metodo:
CONVERGE ⇒ |g′(x∗)| < 1,
44
Figura 5.1: Interpretacao geometrica do Metodo de Picard.
DIVERGE ⇒ |g′(x∗)| > 1.
METODO DE NEWTON:
O procedimento iterativo e determinado atraves da expansao da funcao f(x) em seriede Taylor ao redor da solucao x∗.
f(x∗) = f(x + ∆x) ≈ f(x) + ∆xf ′(x) = 0
⇒ ∆x = x∗ − x = − f(x)
f ′(x)⇒ x∗ = x − f(x)
f ′(x).
x0 ← chute inicial
Repetir ate |xn+1 − xn| < ε
xn+1 = xn − f(xn)f ′(xn)
45
A interpretacao geometrica do metodo de Newton e ilustrada na Figura 3.2.
Figura 5.2: Interpretacao Geometrica do Metodo de Newton.
O metodo de Newton possui uma propriedade de convergencia muito forte. Se o valorinicial do processo iterativo for proximo da solucao, o metodo converge quadraticamente,isto e, o valor do erro decai seguindo uma funcao quadratica, conforme a prova desenvolvidaa seguir.
Prova
Vamos assumir que x∗ e uma raiz simples: f ′(x∗) = 0.
Vamos definir o erro como εn = xn − x∗. Logo x∗ = xn − εn.
Expandindo a funcao f(x) por serie de Taylor ao redor da solucao x∗:
46
f(x∗) = 0 = f(xn) + (x∗ − xn)f ′(xn) +1
2(x∗ − xn)2f ′′(ξ); ξ ∈ (xn, x
∗) ⇒f(xn)
f ′(xn)+ (x∗ − xn) = −
12(x∗ − xn)2f ′′(ξ)
f ′(xn)⇒
x∗ −(
xn − f(xn)
f ′(xn)
)︸ ︷︷ ︸
xn+1
= −12(x∗ − xn)2f ′′(ξ)
f ′(xn)⇒
x∗ − xn+1 = −εn+1 = −12ε2nf
′′(ξ)f ′(xn)
⇒
εn+1 =1
2ε2n
f ′′(ξ)f ′(xn)
.
Logo, tomando o limite quando x → x∗:
εn+1 =1
2ε2n
f ′′(x∗)f ′(x∗)
,
limn→∞
εn+1
ε2n
= C ⇒ εn+1 = Cε2n.
Por exemplo, se na quarta iteracao, a aproximacao x4 estiver perto da solucao do pro-blema e ε4 = 10−3, pela prova acima, ε5 = 10−6 e ε6 = 10−12.
5.2 Sistema de Equacoes Nao-Lineares
Vamos agora considerar um sistema de equacoes nao-lineares, respresentado por
fi(x1, x2, x3, ..., xn) = 0; i = 1, ..., n.
METODO DE PICARD:
Como no caso de uma equacao apenas, cada expressao do sistema deve ser reescritacomo
xi = gi(x1, x2, x3, ..., xn); i = 1, ..., n.
O processo iterativo e similar ao caso de apenas uma variavel.
x(0) = (x(0)1 , x
(0)2 , ..., x
(0)n ) ← chute inicial
47
Repetir ate |x(k+1) − x(k)| < ε
x(k+1)i = gi(x
(k)1 , x
(k)2 , ..., x
(k)n ); i = 1, ..., n
As propriedades de convergencia sao similares as do caso com uma variavel. O procedi-mento iterativo converge se:
‖D(x)‖ ≤ 1; Dij =∂gi
∂xj
.
METODO DE NEWTON:
Como no caso de uma equacao apenas, cada expressao do sistema deve ser desenvolvidapor serie de Taylor ao redor da solucao do problema x∗ = (x∗
1, x∗2, ..., x
∗n):
f1(x∗) = f1(x
(k)) + ∂f1
∂x1(x∗
1 − x(k)1 ) + ∂f1
∂x2(x∗
2 − x(k)2 ) + ... + ∂f1
∂xn(x∗
n − x(k)n ) = 0;
f2(x∗) = f2(x
(k)) + ∂f2
∂x1(x∗
1 − x(k)1 ) + ∂f2
∂x2(x∗
2 − x(k)2 ) + ... + ∂f2
∂xn(x∗
n − x(k)n ) = 0;
...
fn(x∗) = fn(x(k)) + ∂fn
∂x1(x∗
1 − x(k)1 ) + ∂fn
∂x2(x∗
2 − x(k)2 ) + ... + ∂fn
∂xn(x∗
n − x(k)n ) = 0.
Reescrevendo a expressao em notacao matricial:
f(x) = 0 = f(x(k)) + J(x∗ − x(k)) ⇒J (x∗ − x(k))︸ ︷︷ ︸
∆x(k+1)
= −f(x(k)).
A matriz J e chamada de Matriz Jacobiana e definida como:
Jij =∂fi
∂xj
.
O algoritmo do metodo de Newton para um sistema de equacoes nao-lineares pode serescrito como:
x(0) = (x(0)1 , x
(0)2 , ..., x
(0)n ) ← chute inicial
Repetir ate |x(k+1) − x(k)| < εJ(∆x(k+1)) = −f(x(k))x(k+1) = x(k) + ∆x(k+1)
48
5.3 Metodo de Newton aplicado a uma Equacao Dife-
rencial nao linear
Como mostrado no inıcio do capıtulo, no caso da solucao de uma equacao diferencialnao linear pelo metodo de elementos finitos, a discretizacao leva a um sistema de equacoesalgebricas nao linear:
Ri = Ri(C1, C2, ..., CN) = 0 ; para i = 1, ..., N.
Cada equacao do sistema corresponde a um resıduo ponderado:
Ri = −∫ 1
0
du
dx
dφi
dxdx +
∫ 1
0
udu
dxφidx +
∫ 1
0
fφidx = 0.
A solucao pelo Metodo de Newton leva ao seguinte processo iterativo
J(∆C(k+1)) = −R(C(k)),
C(k+1) = C(k) + ∆C(k+1).
Os elementos da matriz Jacobiana, J, representam a sensibilidade da equacao a cadaincognita do problema.
Ji,j =∂Ri
∂Cj
=∂
∂Cj
Ri = −
∫ 1
0
du
dx
dφi
dxdx +
∫ 1
0
udu
dxφidx +
∫ 1
0
fφidx⇒
Ji,j = −∫ 1
0
dφi
dx
∂
∂Cj
(du
dx
)dx +
∫ 1
0
∂
∂Cj
(u)du
dx+ u
∂
∂Cj
(du
dx
)φidx.
A derivada da funcao u e de sua derivada dudx
em relacao aos coeficientes Cj sao calculadosda seguinte forma:
∂
∂Cj
(u) =∂
∂Cj
(N∑
k=1
Ckφk
)=
∂
∂Cj
(C1φ1 + ... + Cjφj + ... + CNφN) = φj,
∂
∂Cj
(du
dx
)=
∂
∂Cj
(N∑
k=1
Ckdφk
dx
)=
∂
∂Cj
(C1
dφ1
dx+ ... + Cj
dφj
dx+ ... + CN
dφN
dx
)=
dφj
dx.
Logo, os elementos da matriz Jacobiana para este problema sao dados por:
Ji,j = −∫ 1
0
dφj
dx
dφi
dxdx +
∫ 1
0
φj
du
dx+ u
dφj
dx
φidx.
O algoritmo de solucao de um problema nao linear pelo metodo de elementos finitosacoplado com o metodo de Newton pode ser escrito como mostrado a seguir:
Ler Chute inicial C(0)
49
do iter = 1, NITERCALL FORMRV (C(iter−1))
Rotina para calculo do vetor resıduo global
CALL FORMJM (C(iter−1))Rotina para calculo da matriz Jacobiana global
Resolver o sistema J∆C = −RC(iter) = C(iter−1) + ∆CErro = |R| + |∆C|if Erro < ε then STOP (Problema Convergiu)
end do
OBSERVACOES IMPORTANTES:
• A cada passo do metodo de Newton, o sistema linear de equacoes J∆C = −R deveser resolvido.
• Se o problema for linear, o metodo de Newton converge em uma iteracao.
• A visao elementar do problema fica inalterada. No lugar de calcular a matriz e o vetorelementar A(e) e b(e), deve-se calcular a matriz jacobiana elementar J(e) e o vetor resıduoelementar R(e).
5.4 Nocoes de Tecnicas de Continuacao
O raio de convergencia do metodo de Newton e tipicamente menor do que o raio deconvergencia do metodo de Picard. Isto significa que o metodo de Newton somente iraconvergir se o chute inicial for proximo a solucao do problema. Porem, se o chute inicialfor adequado, o metodo de Newton converge quadraticamente, levando poucas iteracoesate o calculo da solucao do problema. O metodo de Picard, apesar de possuir um raio deconvergencia maior, tipicamente gasta um numero muito maior de iteracoes para calculara solucao do problema e em muitas vezes nao converge, independente do valor do chuteinicial.
Logo, o procedimento mais seguro e eficiente numericamente e utilizar o metodo deNewton acoplado com tecnicas para procura de chutes iniciais adequados. Este processo decalculo de chutes iniciais a medida que um parametro do problema e variado e denominadode Continuacao.
Um sistema de equacoes algebricas originado a partir da discretizacao de uma equacaodiferencial que descreve um determinado fenomeno fısico pode ser representado de umaforma geral como:
f(x, α) = 0,
50
fi sao as equacoes a serem resolvidas, xi representam a solucao do problema, α e umparametro do problema (numero de Reynolds, por exemplo).
Geralmente a solucao do problema para uma determinada faixa do parametro α e bemsimples. Se α for o numero de Reynolds, o problema f(x, 0) = 0 e linear, e o metodo deNewton converge em uma iteracao, independente do valor do chute inicial.
O processo de continuacao consiste em, sabendo-se a solucao do problema para α = αi,determinar o chute inicial para o problema com α = αi+1.
O procedimento mais rudimentar e utilizar a solucao anterior (com α = αi) como chuteinicial:
x(o)(αi+1) = x(αi)
Este procedimento e conhecido como Continuacao de Ordem 0.
Um chute inicial bem melhor pode ser obtido usando o valor da derivada da solucaoem relacao ao parametro α, como descrito a seguir. Este procedimento e denominadoContinuacao de Primeira Ordem.
O sistema a ser resolvido e representado por
f(x, α) = 0.
Diferenciando esta expressao em relacao a α, obtem-se:
∂f
∂x︸︷︷︸J
·dxdα
+∂f
∂α= 0 ⇒ J
dx
dα= − ∂f
∂α⇒
dx
dα= J−1
(− ∂f
∂α
)⇒ ∆x = −∆αJ−1 ∂f
∂α.
Logo, o valor do chute inicial para α = αi+1 e obtido pela seguinte expressao:
x(o)(αi+1) = x(αi) − ∆αJ−1 ∂f
∂α.
Observe que se o metodo de Newton tiver sido utilizado para obter a solucao do problemanao linear para α = αi, o custo computacional do calculo do chute inicial e bem pequeno,uma vez que a matriz Jacobiana ja foi invertida para o calculo de x(αi).
A Figura 3.3 ilustra o chute inicial obtido com os dois procedimentos de continuacaodescritos nesta secao.
51
Figura 5.3: Chute inicial obtido com continuacao de ordem zero e primeira ordem.
52
Capıtulo 6
Formulacao Integral da Equacao deNavier-Stokes
A equacao de Navier-Stokes para fluidos incompressıveis :ρ u · ∇u = ∇ · T, (em regime permanente);
∇ · u = 0;(6.1)
onde para um Fluido Newtoniano ; T = −pI + µ[∇u + ∇uT ].
As variaveis do problema sao o campo de velocidade u e o campo de pressao p.
Figura 6.1: Domınio para problema de escoamento.
As condicoes de contorno para a equacao de Navier-Stokes podem ser em velocidade ouforca agindo no contorno do domınio. A condicao de contorno utilizada depende da fısicado problema. As fronteiras de um escoamento sao, em geral:
• Entrada de fluido.• Saıda de fluido.• Paredes solidas.
53
• Interface entre dois fluidos.
Vamos assumir que para x ∈ Γ1, prescrevemos a velocidade; i.e.
u = V em x ∈ Γ1. (6.2)
Esta e uma condicao de contorno essencial. Corresponde as condicoes de impermeabili-dade e nao deslizamento.
u · n = V · n → impermeabilidade,u · t = V · t → nao deslizamento.
Para x ∈ Γ2, vamos prescrever a forca superficial
n · T︸ ︷︷ ︸tracao
= f em x ∈ Γ2. (6.3)
Esta e uma condicao de contorno natural. Prescreve a forca na fronteira. A forca pre-scrita f depende da fısica do problema.
Para resolver este sistema de equacoes diferenciais parciais pelo metodo dos resıduos pon-derados, temos que multiplicar o resıduo da aproximacao de cada equacao por um funcaopeso e forcar a integral ao longo de todo o domınio a ser nula.
Metodo dos Resıduos Ponderados
Como a equacao de consevacao da quantidade de movimento (4.1) e uma equacao veto-rial, o resıduo ponderado correspondente a esta equacao sera calculado atraves do produtointerno do resıduo da aproximacao com um funcao peso vetorial W. A equacao da con-tinuidade e uma equacao escalar e desta forma, a funcao peso utilizada w e uma funcaoescalar, como as utilizadas ate o capıtulo anterior.
Rm =
∫Ω
ρ u · ∇u −∇ · T · W dΩ = 0;
Rc =
∫Ω
∇ · uw dΩ = 0.
• Resıduo associado a equacao de quantidade de movimento:
54
Rm =
∫Ω
ρ u · ∇u −∇ · T · W dΩ =
=
∫Ω
ρ (u · ∇u) · WdΩ −∫
Ω
(∇ · T) · W dΩ.
O termo ∇ · T apresenta segunda derivadas da velocidade (variavel primitiva do prob-lema). Vamos integrar por partes para transferir a derivada para a funcao peso, conformeja feito anteriormente, usando a seguinte igualdade tensorial:
T : ∇W = ∇ · (T · W) − (∇ · T) · W
⇒∫
Ω
(∇ · T) · W dΩ = −∫
Ω
T : ∇W +
∫Ω
∇ · (T · W) dΩ.
Usando o teorema de Gauss podemos escrever:∫
Ω∇ · (T · W) dΩ =
∫Γn · (T · W) dΓ
logo,
⇒∫
Ω
(∇ · T) · W dΩ = −∫
Ω
T : ∇W +
∫Γ
n · (T · W) dΓ.
Finalmente:
Rm =
∫Ω
ρ (u · ∇u) · W +
∫Ω
T : ∇W −∫
Γ
(n · T) · W dΓ.
A funcao peso vetorial W pode ser escrita em termos de suas componentes: W =[W1,W2] e cada termo do resıduo da equacao de conservacao da quantidade de movimentoescrito em termos de suas componentes. No caso particular de sistema cartesiano, cadatermo do resıduo ponderado e escrito como:
1) (u · ∇u) · W = W1
(u∂u
∂x+ v ∂u
∂y
)+ W2
(u ∂v
∂x+ v ∂v
∂y
).
2) T : ∇W = ∂W1
∂x
[−p + 2µ
∂u
∂x
]︸ ︷︷ ︸
Txx
+ ∂W1
∂y
[µ(
∂u
∂y+
∂v
∂x)
]︸ ︷︷ ︸
Txy
+
+ ∂W2
∂x
[µ(
∂u
∂y+
∂v
∂x)
]︸ ︷︷ ︸
Tyx=Txy
+ ∂W2
∂y
[−p + 2µ
∂v
∂y
].︸ ︷︷ ︸
Tyy
3) (n · T) · W = fx W1 + fy W2.
Cada componente da funcao peso vetorial pode ser escrita como combinacao linear defuncoes base escalar ψi. Se cada componente W1 e W2 pertence a um espaco vetorial de
55
dimensao n, a funcao peso vetorial W pertence a um espaco de dimensao 2n que e geradopela seguinte base de funcoes:
W =
[W1
W2
]︸ ︷︷ ︸Espacos de dim=2n
⇒[
ψ1
0
],
[ψ2
0
], . . . ,
[ψn
0
]︸ ︷︷ ︸
n funcoes
,
[0ψ1
],
[0ψ2
], . . . ,
[0ψn
].︸ ︷︷ ︸
n funcoes
Definindo o espaco das funcoes peso desta forma, as componentes da equacao de con-servacao podem ser desacopladas, pois para as primeiras n funcoes base aparecem termoscorrespondentes a primeira componente da velocidade, ja que W2 = 0 e para as ultimas nfuncoes base aparecem temos correspondentes a segunda componente da velocidade, ja quepara estas funcoes base W1 = 0.
Usando o fato do fluido Newtoniano, ou seja, T = −pI + µ[∇u + ∇uT ], as 2n equacoesalgebricas associadas a equacao de conservacao da quantidade de movimento sao escritascomo:
Rimx =
∫Ω
ρ ψi
[u∂u
∂x+ v
∂u
∂y
]+
∂ψi
∂x
[−p + 2µ
∂u
∂x
]+
+∂ψi
∂y
[µ(
∂u
∂y+
∂v
∂x)
]dΩ −
∫Γ
ψi fx dΓ ; i = 1, ..., n.
Rimy =
∫Ω
ρ ψi
[u
∂v
∂x+ v
∂v
∂y
]+
∂ψi
∂x
[µ(
∂u
∂y+
∂v
∂x)
]+
+∂ψi
∂y
[−p + 2µ
∂v
∂y
]dΩ −
∫Γ
ψi fy dΓ ; i = 1, ..., n.
• Resıduo associado a equacao da continuidade:
Rc =
∫Ω
∇ · uw dΩ = 0.
Onde w esta no espaco gerado pelas seguintes funcoes peso: χ1, χ2, ..., χm. As mequacoes algebricas associadas a equacao da continuidade sao escritas como:
Ric =
∫Ω
(∂u
∂x+
∂v
∂y) χi dΩ ; i = 1, ...,m.
Os campos de velocidade e pressao devem ser escritos em termos de funcoes base. Osdiferentes termos dos resıduos associados a equacao de conservacao de quantidade de movi-mento contem pressao ou derivadas da velocidade. Como estes termos devem ter a mesma
56
precisao de discretizacao, as funcoes bases utilizadas para pressao nao precisam ser damesma ordem que as usadas para o campo de velocidade. Quando os dois campos devariaveis pertencem a espacos diferentes, diz-se que a Formulacao e Mista (Mixed FiniteElement).
As variaveis do problema sao escritas em termos de funcoes base como:
u =
[uv
]=
[ ∑ni=1 Uiφi∑ni=1 Viφi
]⇒ 2 n incognitas;
p =m∑
i=1
Pi χi ⇒ m incognitas.
Repare que a funcao peso ultilizada na equacao de continuidade e a mesma funcao ul-tilizada para interpolar a pressao (χi) e assim o numero de incognitas (2n + m) e o mesmoque o numero de equacoes.
• Como a formulacao nao apresenta derivadas da pressao ou da funcao peso χi, o espacodas funcoes χi pode ser menor que o das funcoes φi ultilizadas para interpolar a velocidade(m < N) e as funcoes χ′
is podem ser descontınuas. As funcoes φ′is devem ser contınuas,
pois tem que satisfazer∫ |∇φi|2 < ∞. Nao existe essa restricao para χ′
is.
• Nao e qualquer combinacao arbitraria de φ′is e χ′
is que leva a formulacoes com bomdesempenho numerico. Uma ma escolha de φ′
is e χ′is pode levar a metodos instaveis.
Um exemplo simples de um dos problemas que podem ocorrer e discutido a seguir. Parafacilitar a interpretacao geometrica, vamos analisar um problema de elasticidade linearde materiais incompressıveis. As equacoes que descrevem o problema sao identicas ao doproblema de Stokes (Re = 0), sendo u o deslocamento.
∇ · T = 0 ;
∇ · u = 0.
Vamos considerar a situacao de dois elementos triangulares, como mostrado na Figura4 e funcoes base φi linear e χi constante.
⇒∫
Ω(e)
(∇ · u) χi︸︷︷︸cte
dΩ = 0 ⇒∫
Ω(e)
∇ · u dΩ = 0
A condicao de incompressibilidade com χi constante leva a uma condicao de area con-stante em todos os elementos.∫
I
∇ · u dΩ = 0. ⇔ Area do elemento I e constante.
Observe que em se tratando de um triangulo, se a area e constante , entao o produtobase por altura nao muda, logo o no A so pode deslocar-se na vertical.
Por outro lado,∫II
∇ · u dΩ = 0 ⇔ Area do elemento II e constante.
57
Figura 6.2: Exemplo de formulacao instavel - malha truncada.
Para manter a area II fixa o no A so pode deslocar-se na horizontal.
uA = 0 ⇒ Malha Trancada
Isso se deve a uma ma escolha da funcao peso. A tabela a seguir lista diferentes com-binacoes de funcoes base φi e ψi que levam a metodos estaveis.
Elementos Retangulares φi χi
bilinear (n = 4) constante (m = 1)biquadratico (n = 9) bilinear (m = 4)biquadratico (n = 9) linear discontınuo (m = 4)
bicubico (n = 16) biquadratico (m = 9)
Elementos Triangulares quadratico (n = 4) linear (m = 3)cubico (n = 4) quadratico (m = 6)
Existe uma condicao que determina se uma certa combinacao de espacos de funcoes paraa velocidade e pressao e valida ⇒ CONDICAO DE BABUSKA-BREZZI.
A condicao de Babuska-Brezzi verifica a consistencia das aproximacoes das derivadas.A interpolacao de u e p nao podem ser escolhidas arbitrariamente.
∀χi; φisup
∫Ω χi ∇·φi dΩ
∫Ω |∇φi|2 dΩ 12≥ c ∫
Ωχ2
i dΩ 12
A verificacao e prova se uma determinada combinacao de φi e χi satisfazem a condicaode Babuska-Brezzi sao complicados matematicamente e estao fora do escopo do curso, de-vendo ser tratado em um curso mais avancado de elementos finitos.
A formulacao mista apresentada pode ser interpretada como solucao de uma equacaodiferencial (quantidade de movimento) com uma restricao imposta (continuidade) pelos
58
metodos dos Multiplicadores de Lagrange.
Na apresentacao do metodo de elementos finitos para equacao de Navier-Stokes, vamosnos concentrar no seguinte elemento : Bi-quadratico para a velocidade e linear discontınuopara pressao. Neste elemento, cada componente da velocidade possui 9 graus de liberdadee o campo de pressao e representado por 3 graus de liberdade.
Figura 6.3: Elemento bi-quadratico de 9 nos.
∑9i=1 Uiφi∑9i=1 Viφi
9 graus de liberdade para cada componente da velocidade,
p =∑m
i=1 Pi χi ⇒ 3 graus de liberdade para pressao.
Cada elemento possui 21 graus de liberdade. As matrizes elementares sao 21 × 21.As funcoes base elementares φi utilizadas para expandir o campo de velocidade serao
lagrangeanas, isto e:
φi(Xj) = δij,
cada funcao base e igual a 1 no associado a funcao e nula nos demais. Cada coeficiente Uj
e Vj representa o valor da velocidade u e v no no j. Em termos das coordenadas locais, asfuncoes base para a velocidade sao:
φ1(ξ, η) =ξ(ξ − 1)η(η − 1)
4
59
φ2(ξ, η) =ξ(ξ + 1)η(η − 1)
4
φ3(ξ, η) =ξ(ξ + 1)η(η + 1)
4
φ4(ξ, η) =ξ(ξ − 1)η(η + 1)
4
φ5(ξ, η) =(1 − ξ2)η(η − 1)
2
φ6(ξ, η) =ξ(ξ + 1)(1 − η2)
2
φ7(ξ, η) =(1 − ξ2)η(η + 1)
2
φ8(ξ, η) =ξ(ξ − 1)(1 − η2)
2
φ9(ξ, η) = (1 − ξ2)(1 − η2)
A figura 4.4 ilustra a forma destas funcoes base.As funcoes base utilizadas para expandir o campo de pressao nao sao Lagrangeanas. Para
este elemento, elas sao escolhidas de forma que o primeiro grau de liberdade de pressao P1
represente o valor da pressao no centro do elemento; o segundo grau de liberdade P2, aderivada da pressao na direcao η; e o terceiro grau de liberdade P3, a derivada da pressaona direcao ξ. Para isto, a variacao da pressao em cada elemento deve ser escrita como:
p = P1 + P2 η + P3 ξ, consequentemente, as funcao base ψi sao:
ψ1(ξ, η) = 1
ψ2(ξ, η) = η
ψ3(ξ, η) = ξ
Observe que:• p(ξ = 0, η = 0) = P1;• dp
dη= P2 e
• dpdξ
= P3.
Substituindo as expansoes para o campo de velocidade e pressao nas equacoes dosresıduos ponderados, obtem-se um sistema de 2n + m equacoes algebricas nao lineares.As nao linearidades vem dos termos convectivos da equacao de Navier-Stokes. Neste curso,o sistema nao linear sera resolvido pelo metodo de Newton. Para isto, deve-se calcular a ma-triz Jacobiana que representa a sensibilidade de cada equacao em relacao a cada incognita:
Jij = ∂Ri
∂cj⇒ Derivada do Resıduo i em relacao ao coeficiente cj.
Ate este capıtulo, a numeracao dos graus de liberdade dos elementos se confundia coma propria numeracao dos nos. Como agora estamos resolvendo um sistema de equacoesdiferenciais parciais para os campos de velocidade (campo vetorial) e pressao (escalar),
60
Figura 6.4: Funcoes base bi-quadraticas.
o numero de graus de liberdade de cada elemento e maior do que o numero de nos doelemento. O elemento que estamos trabalhando possui 9 nos e 21 graus de liberdade. Saoeles: U1, ..., U9; V1, ..., V9; P1, P2, P3.
A numeracao dos graus de liberdade elementar e totalmente arbitraria. A numeracaolocal adotada aqui e mostrada na tabela a seguir:
61
# grau de liberdade local grau de liberdade1 U1
2 U2
3
4...
5678 U8
9 U9
10 V1
11 V2
12
13...
14151617 V8
18 V9
19 P1
20 P2
21 P3
Figura 6.5: Numeracao local dos graus de liberdade.
Para a montagem da matriz jacobiana elementar e interessante construir um vetor queaponte o numero do grau de liberdade elementar associado a cada no:
ivxlocaldof(1) = 1(2) = 2
62
...(9) = 9
ivylocaldof(1) = 10(2) = 11...(9) = 18
ipresslocaldof(1) = 19ipresslocaldof(2) = 20ipresslocaldof(3) = 21
Algoritmo para calculo de A(e)
do igp = 1, NGPXI = XGP (igp)WI = WGP (igp)do jgp = 1, NGP
NETA = XGP (jgp)WJ = WGP (jgp)W = WI ∗ WJCALL BASISFUNC (XI,NETA, PHI,GradPHI, PSI)CALL MAPPING (XY,PHI,GradPHI, JAC, JACINV, detJ)GradPHI XY = JACINV ∗ GradPHICalcular u, ∂u
∂x, ∂u
∂y, v, ∂v
∂x, ∂v
∂y.
do ilnode = 1, 9ivx = ivxlocaldof(ilnode)ivy = ivylocaldof(ilnode)do jlnode = 1, 9
jvx = jvxlocaldof(jlnode)jvy = jvylocaldof(jlnode)JacElem(ivx,jvx) = ...
JacElem(ivx,jvy) = ...
JacElem(ivy,jvx) = ...
JacElem(ivy,jvy) = ...
end do
do jpress = 1, 3jpress = ipresslocaldof(jlpress)JacElem(ivx,jpress) = ...
JacElem(ivy,jpress) = ...
63
end do
do jpress = 1, 3jpress = jpresslocaldof(ilpress)do jlnode = 1, 9
jvx = ivxlocaldof(jlnode)
jvy = ivylocaldof(jlnode)
JacElem(ipress, jvx) = ...
JacElem(ipress,jvy) = ...
end do
end do
Final do algoritmo ⇒ JacElem (21 × 21)
Precisamos agora montar a matriz global. Vamos fazer esse processo para umas malhasimples com 4 elementos.
EXEMPLO
*FIG
Malha com 4 elementos e 25 nos.
• Nos → A numeracao global segue de acordo com a numeracao DomNodeID (ilnode, iele)dos elementos.
DomNodeID
ilnode iele 1 2 3 41 1 4 2 32 2 3 16 173 3 10 17 224 4 11 3 105 5 7 18 206 6 12 19 237 7 13 20 248 8 14 6 129 9 15 21 25
Numeracao local dos nos.
• Graus de liberdade → A numeracao global segue de acordo com a numeracao doselementos.
FIGS*AssMtrx
64
ildofiele 1 2 3 41 1 4 2 32 2 3 37 383 3 22 38 524 4 23 3 225 5 7 39 416 6 24 40 537 7 25 41 548 8 26 6 249 9 27 42 55
10 10 13 11 1211 11 12 43 4412 12 28 44 5613 13 29 12 2814 14 16 45 4715 15 30 46 5716 16 31 47 5817 17 32 15 3018 18 33 48 5919 19 34 49 6020 20 35 50 6121 21 36 51 62
ALGORITMO DETALHADO PARA SOLUCAO
DA EQUACAO DE NAVIER-STOKES
METODO DE NEWTON
J∆C = −R
C(k+1) = C(k) + ∆C
PROGRAMA PRINCIPAL
• Inicializar - Ler chute inicial ou inicializar o valor dos coeficientes
→ c• Calcular o vetor de Resıduo global
CALL FormRV → R• Se ‖R‖ < ε → stop. O ‘‘chute inicial’’ era soluc~ao do problema.
• Repetir enquanto ‖R‖ > ε- Calcular a matriz Jacobiana global
65
CALL FormJM → J- Resolver J∆C = −R
CALL SOLVER (J,R,∆C)
- Calcular C = C + ∆C- Calcular o vetor Resıduo global
CALL FormRV → R• Gravar a soluc~ao C
FormRV → Calculo do vetor de Resıduo global
do iele = 1, NELECALL getelemRV (→ elemRV )
do ildof = 1, NDOFELEigdof = ASSMtrx (ildof, iele)RV (igdof) = RV (igdof) + elemRV (ildof)
end do
end do
FormJM → Calculo da matriz jacobiana global
do iele = 1, NELECALL getelemJM (→ elemJM)
do ildof = 1, NDOFELEigdof = ASSMtrx (ildof, iele)do jldof = 1, NDOFELE
jgdof = ASSMtrx (jldof, iele)JM(igdof, jgdof) = JM(igdof, jgdof) + elemJM(ildof, jgdof)
end do
end do
end do
O que precisaremos para o calculo do Resıduo elementar:Formulas:
Rimx =
∫Ω
Re ψi
[u∂u
∂x+ v
∂u
∂y
]+
∂ψi
∂x
[−p + 2µ
∂u
∂x
]+
+∂ψi
∂y
[µ(
∂u
∂y+
∂v
∂x)
]dΩ
66
Rimy =
∫Ω
Re ψi
[u
∂v
∂x+ v
∂v
∂y
]+
∂ψi
∂x
[µ(
∂u
∂y+
∂v
∂x)
]+
+∂ψi
∂y
[−p + 2µ
∂v
∂y
]dΩ
Ric =
∫Ω
(∂u
∂x+
∂u
∂x) χi dΩ
• Lembrar que vamos usar um mapeamento isoparametrico para transformar as coorde-nadas x, y em ξ, η. ∫
Ω
F (x, y) dxdy =
∫Ω
F (ξ, η) |J| dξdη.
• Lembrar que vamos usar integracao de Gauss para calcular a integral:∫Ω
F (ξ, η) dξdη =∑igp
∑jgp
F (ξigp, ηjgp) WigpWjgp.
Estrutura de Dados:
Definir os seguintes vetores / matrizes:
• PHI(i)i = 1, ..., 9
PHI ⇒
φ1
φ2...
φ9
• GradPHI(i, j)i = 1, 2
j = 1, ..., 9GradPHI ⇒
φ1
dξ. . . φ9
dξ
φ1
dη. . . φ9
dη
• PSI(i)i = 1, 2, 3
PSI ⇒ χ1
χ2
χ3
• XY (i, j)
i = 1, 2j = 1, ..., 9
XY ⇒[
X1 . . . X9
Y1 . . . Y9
]
• velocity(i, j)i = 1, 2
j = 1, ..., 9velocity ⇒
[U1 . . . U9
V1 . . . V9
] • pressure(i)i = 1, 2, 3
pressure ⇒ P1
P2
P3
• GradPHI XY ⇒[ φ1
dx. . . φ9
dxφ1
dy. . . φ9
dy
]= J−1 ∗
[φ1
dξ. . . φ9
dξφ1
dη. . . φ9
dη
]
• JAC ⇒[
dxdξ
dydξ
dxdη
dydη
]=
[φ1
dξ. . . φ9
dξφ1
dη. . . φ9
dη
]∗
X1 Y1
......
X9 Y9
67
• UV =
[UV
]=
[U1 . . . U9
V1 . . . V9
]∗
φ1
...φ9
• P =[
P]
=[
P1 P2 P3
] ∗ χ1
χ2
χ3
• GradUV XY ⇒[
dudx
dvdx
dudy
dvdy
]=
[ φ1
dx. . . φ9
dxφ1
dy. . . φ9
dy
]∗
U1 V1
......
U9 V9
getelemRV → Calculo do Resıduo elementar
do igp = 1, NGPXI = XGP (igp)WI = WGP (igp)do jgp = 1, NGP
NETA = XGP (jgp)WJ = WGP (jgp)W = WI ∗ WJCALL BASISFUNC (XI,NETA, PHI,GradPHI, PSI)CALL MAPPING (XY,PHI,GradPHI, JAC, JACINV, detJ)CALL VELOCITY (velocity, PHI, UV )
CALL PRESSURE (pressure, PSI, P)
GradPHI XY = JACINV ∗ GradPHIGradUV XY = GradPHI XY ∗ UV T
do ilnode = 1, 9ivx = ivxlocaldof(ilnode)ivy = ivylocaldof(ilnode)elemRV (ivx) = elemRV (ivx) + W ∗ detJ ∗ (Re ∗ PHI(ilnode)∗
(UV (1) ∗ GradUV XY (1, 1) + UV (2) ∗ GradUV XY (1, 2))+GradPHI XY (1, ilnode) ∗ (−P (1) + 2 ∗ GradUV XY (1, 1))+GradPHI XY (2, ilnode)∗(GradUV XY (1, 2)+GradUV XY (2, 1)))
end do
end do
end do
68
Desta forma, as primeiras 9 linhas da matriz Jacobiana elementar correspondem aos 9resıduos da componente x da equacao de Navier-Stokes; as 9 linhas seguintes, aos 9 resıduosda componente y; e as 3 ultimas linhas, aos 3 resıduos da equacao de continuidade. Damesma forma, as 9 primeiras colunas estao associadas aos coeficientes Uj; as 9 seguintes,aos coeficientes Vj; e as 3 ultimas colunas, aos 3 coeficientes Pj. Seguindo a numeracao dosgraus de liberdade apresentada na tabela, a matriz Jacobiana elementar possui a seguinteestrutura:
∂Rimx
∂Uj
∂Rimx
∂Vj
∂Rimx
∂Pj
∂Rimy
∂Uj
∂Rimy
∂Vj
∂Rimy
∂Pj
∂Ric
∂Uj
∂Ric
∂Vj
∂Ric
∂Pj
.
A expressao de cada cada bloco desta estrutura a apresentada a seguir:
• Jacobiano de Rimx :
∂Rimx
∂Uj
=
∫Ωe
ρ ψi
[φj
∂u
∂x+ u
∂φj
∂x+ v
∂φj
∂y
]+
+∂ψi
∂x2 µ
∂φj
∂x+ µ
∂ψi
∂y
∂φj
∂y
dΩe;
∂Rimx
∂Vj
=
∫Ωe
ρ ψi φj
∂u
∂y+ µ
∂ψi
∂y
∂φj
∂x
dΩe;
∂Rimx
∂Pj
=
∫Ωe
− ∂ψi
∂xχi dΩe.
• Jacobiano de Rimy :
69
∂Rimy
∂Uj
=
∫Ωe
ρ ψi φj
∂v
∂x+ µ
∂ψi
∂x
∂φj
∂y
dΩe;
∂Rimy
∂Vj
=
∫Ωe
ρ ψi
[u∂φj
∂x+ φj
∂v
∂y+ v
∂φj
∂y
]+
+ µ∂ψi
∂x
∂φj
∂x+ 2 µ
∂ψi
∂y
∂φj
∂y
dΩe;
∂Rimy
∂Pj
=
∫Ωe
− ∂ψi
∂yχi dΩe.
• Jacobiano da Continuidade :
∂Ric
∂Uj
=
∫Ωe
∂φj
∂xχi dΩe;
∂Ric
∂Vj
=
∫Ωe
∂φj
∂yχi dΩe;
∂Ric
∂Pj
= 0.
O mapeamento da matriz elementar para a matriz global e feito de maneira semelhanteaos problemas anteriores. A diferenca e que agora cada no possui mais de 1 grau de liber-dade. As linhas e colunas da matriz elementar estao associadas ao grau de liberdade e naoaos nos. Nos exemplos dos capıtulos anteriores, as duas numeracoes eram coincidentes. Omapeamento deixa de ser por nos (usando a matriz DomNodeID) e passa a ser por grausde liberdade.
Defne-se a matriz:
AssMtx(ildof, iele) = # global do grau de liberdade local ildof do elemento iele
Apos a montagem de todas as matrizes elementares, a matriz jacobiana global fica coma seguinte forma:
∂Rimx
∂Uj
∂Rimx
∂Vj
∂Rimx
∂Pj
∂Rimy
∂Uj
∂Rimy
∂Vj
∂Rimy
∂Pj
∂Ric
∂Uj
∂Ric
∂Vj
∂Ric
∂Pj= 0
=
K −C
CT 0
.
70
O bloco de “zeros” na diagonal da matriz pode representar um problema a mais nasolucao do sistema linear J ∆C = −R em cada passo do metodo de Newton.
Existem outras formulacoes alem da formulacao mista discutida anteriormente utilizadaspara resolver o problema de escoamento incompressıvel.
⇒Metodo da Penalidade
A ideia basica e escrever uma equacao que relacione a pressao com a velocidade. Seesta equacao de estado for quase incompressıvel, entao a solucao obtida pelo metodo dapenalidade vai ser proxima da solucao do problema incompressıvel.
A formulacao do problema pelo Metodo da Penalidade e:ρ u · ∇u = ∇ · T ; T = −pI + µ[∇u + ∇uT ];
∇ · u = ε p ; onde ε e pequeno.
Quando ε → 0, recupera-se a formulacao incompressıvel.
Pode-se substituir a pressao na equacao de Navier-Stokes ⇒ p = −1ε∇ · u :
ρ u · ∇u = ∇ · [−1ε∇ · u + µ (∇u + ∇uT )
].
A vantagem do Metodo da Penalidade e que a pressao foi eliminada e nao e mais variaveldo problema. O numero de variaveis do problema foi reduzido.
A proximidade da solucao obtida com a solucao do problema incompressıvel dependedo valor de ε. Deve-se utilizar um ε pequeno (o menor possıvel). Porem a medida que εdiminui (e consequentemente 1
εaumenta), a contribuicao relativa dos termos viscosos da
tensao diminui e a solucao obtida nao satisfaz corretamente a equacao de Navier-Stokes.Esta e uma grande desvantagem do metodo.
⇒Metodo da Penalidade com Interpolacao Mista.
E uma combinacao da Formulacao de Penalidade e Formulacao Mista. O objetivo eevitar o problema descrito anteriormente . Como a formulacao e mista, resolve-se o campode velocidade e pressao, porem a continuidade nao e satisfeita totalmente:
ρ u · ∇u = ∇ · [−p + µ(∇u + ∇uT )];
∇ · u = ε p.
Rm =
∫Ω
ρ u · ∇u = ∇ · [−p + µ(∇u + ∇uT )]
· w dΩ = 0;
Rc =
∫Ω
∇ · u − ε pm dΩ = 0;
71
u =9∑
i=1
Uiφi ; v =9∑
i=1
Viφi ; p =m∑
i=1
Pi χi.
Os resıduos da equacao de quantidade de movimento sao identicos aos da formulcaomista. A unica diferenca e na equacao de continuidade.
Ric =
∫Ω
[∂u
∂x+
∂v
∂y− ε p
]χi dΩ;
∂Ric
∂Uj
=
∫Ω
∂φj
∂xχi dΩ;
∂Ric
∂Vj
=
∫Ω
∂φj
∂yχi dΩ;
∂Ric
∂Pj
=
∫Ω
ε χi χj dΩ.
A matriz Jacobiana global fica com a seguinte estrutura:
K −C
CT ε M
; onde Mij = −
∫Ω
χi χj dΩ.
A diferenca para a formulacao mista usual e que o bloco de “zeros” na diagonal e sub-stituido por ε M.
Como M possui inversa, pode-se desacoplar a velocidade e pressao:
K −C
CT ε M
∆U
∆P
=
−Rm
−Rc
.
K ∆U − C ∆P = −Rm
CT ∆U + ε M ∆P = −Rc ⇒ ∆P = 1εM−1[−Rc − CT ∆U].
Substituindo-se ∆P em funcao de ∆U, pode-se chegar a um sistema de equacoes so-mente para o passo de velocidade:
72
K ∆U − C
1εM−1(−Rc − CT ∆U)
= −Rm
[K + 1εC M−1 CT] ∆U = −Rm − 1
εC M−1 Rc
O Problema desacoplado pode ser reescrito como:
K + 1εC M−1 CT 0
CT ε M
∆U
∆P
=
−Rm − 1εC M−1 Rc
−Rc
.
Calcula-se primeiramente a velocidade e depois a pressao:
p =1
ε∇ · u.
73
Capıtulo 7
Condicoes de Contorno
As condicoes de contorno de um problema em mecanica dos fluidos podem ser de 2 tipos:
• Condicao de contorno em velocidade (Dirichlet) : Especifica o valor da velocidadeou de um de seus componentes na fronteira do domınio.
• Condicao de contorno em tracao (Newman) : Especifica a forca ou um de seus com-ponentes na fronteira do domınio.
As diversas situacoes fısicas de interesse se enquadram nestas categorias:
i) Paredes solidas → condicao de nao deslizamento e impermeabilidade.
impermeabilidade ⇒ n · v = n · V ⇒
u = Vx,v = Vy;
nao deslizamento ⇒ t · v = t · V;
onde t e o vetor tangencial, n o normal a parede e V a velocidade da parede.
ii) Condicao de entrada do domınio → considera-se um perfil de velocidade conhecido.
Por exemplo, perfil parabolico:
u = u(y),v = 0.
y
x
Figura 7.1: perfil parabolico.
74
iii) Condicao de simetria
Ft = 0 → forca tangencial = 0,v = 0 → velocidade vertical = 0.
y
x
Linha de simetria
Figura 7.2: linha de simetria.Ft = 0 ⇒ t · (T · n) = 0,v · n = 0.
x
y
n
t
Figura 7.3: linha de simetria.
iv) Condicao de tracao nula.
Figura 7.4: tracao nula.
T · n = 0 ⇒
Ft = t · (T · n) = 0,Fn = n · (T · n) = 0.
75
v) Condicao em interfaces.
liquido
ar
n
t
Figura 7.5: interfaces.
n · T =
tensao superficial︷︸︸︷γ
dt
ds+ n
pressao atm︷︸︸︷pa
⇒ n · (T · n) =
curvatura︷︸︸︷K γ + pa,
t · (T · n) = 0.
• As condicoes de contorno de Dirichlet sao geralmente impostas de uma maneiraessencial. O resıduo correspondente a variavel e simplesmente substituido por:
Ri = Ui − U∗,
onde U∗ e o valor que deseja-se impor a Ui.
Por exemplo, vamos analisar o problema ja visto anteriormente:
V = 1
Figura 7.6: no 7.
A condicao de contorno e:
76
velocidade horizontal ⇒ u = 1 ⇒ u − 1 = 0,
velocidade vertical ⇒ v = 0.
Como o grau de liberdade local correspondente a velocidade horizontal no no marcadoe #7 e o correspondente a velocidade vertical e #16, vamos substituir estes respectivosresıduos por:
R7 = U7 − 1,
R16 = V7 − 0.
Sendo Jij = ∂Ri
∂Uj, note que:
J7j = 0, para j = 7 e J77 =∂R7
∂U7
= 1,
J16j = 0, para j = 16 e J16 16 =∂R16
∂V7
= 1.
A procedimento usual e calcular os resıduos e as jacobianas elementares independente-mente das condicoes de contorno. Se o elemento pertencer a fronteira do domınio, impoe-sea condicao de contorno. Se a condicao for de Dirichlet, o resıduo associado ao grau deliberdade e substituido por: Ri = Ui − U∗ e a linha da matriz jacobiana associada ao graude liberdade e zerada a menos do elemento da diagonal, que e substituido por 1.
No caso do exemplo anterior, onde a condicao de contorno e de Dirichlet, para o ele-mento 2, a matriz jacobiana fica:
× × . . . ×× × . . . ×× × . . . ×0 0 1 0 . . . 00 0 0 1 . . . 0...
. . ....
0 . . . 1
• As condicoes de Newman aparecem naturalmente na formulacao do problema.Lembrando a formulacao de resıduos ponderados da equacao de Navier-Stokes, tem-se:
77
7 34
8 9 6
15
2
Figura 7.7: elemento 2.
Rm =
∫Ω
Re u · ∇u −∇ · T · W dΩ =
=
∫Ω
Re (u · ∇u) · W + T : ∇W −∫
Γ
n · (T · W) dΓ. ⇒
Rimx =
∫Ω
Re φi
[u∂u
∂x+ v
∂u
∂y
]+
∂φi
∂x
[−p + 2
∂u
∂x
]+
+∂φi
∂y
[(∂u
∂y+
∂v
∂x)
]dΩ −
∫Γ
φi (n · T)x dΓ.
Rimy =
∫Ω
Re φi
[u
∂v
∂x+ v
∂v
∂y
]+
∂φi
∂x
[(∂u
∂y+
∂v
∂x)
]+
+∂φi
∂y
[−p + 2
∂v
∂y
]dΩ −
∫Γ
φi (n · T)y dΓ.
Observe que a integral∫
Γφi (n · T)x/y dΓ so sera calculada para os elementos da
fronteira do domınio onde a condicao de contorno de Newman e imposta, pois nos elementosinternos nao temos que nos preocupar com a integral ao longo da fronteira do domınio Γ enos elementos da fronteira onde a condicao de Dirichlet e imposta o resıduo ponderado esubstituido da forma mostrada anteriormente.
Assim, quando a condicao e de Newman e o elemento esta na fronteira, o valor dointegrando (n ·T)x ou (n ·T)y vem da condicao de contorno a ser imposta. Por exemplo,no caso da tracao nula:
(n · T)x = 0, (n · T)y = 0;
no caso de interface:
(n · T)x = γdtxds
+ nx pa, (n · T)y = γdtyds
+ ny pa.
78
• OBSERVACAO : A forma da equacao de Navier-Stokes a qual se aplica a formulacaofraca e importante na imposicao da condicao de contorno de Newman.
A forma utilizada no curso foi a do divergente da tensao:
∇ · T = ∇ · [−pI + ∇u + ∇uT].
Observe que esse termo pode ser simplificado usando que ∇ · u = 0, pela equacao dacontinuidade;
= ∇ · [−pI + ∇u + ∇uT]
= −∇p + ∇ · (∇u)︸ ︷︷ ︸∇2u
+∇ · (∇uT)︸ ︷︷ ︸∇(∇·u)
= −∇p + ∇2u + ∇(∇ · u︸ ︷︷ ︸=0
)
= −∇p + ∇2u
Esta forma simplificada e muito utilizada na literatura de diferencas finitas, porem otermo que aparece na integral de linha nao possui significado fısico, por isso nos problemasonde as condicoes de contorno de Newman sao forcas atuando na fronteira do domınio naoe interessante usar a forma simplificada.
CONDICAO DE CONTORNO DE SAIDA DE ESCOAMENTO
Pode ocorrer condicoes de contorno na saıda do escoamento que nao representem umaforca, apenas sao a derivada da velocidade; como no caso do escoamento desenvolvido, quea condicao de contorno pode ser escrita como n · ∇u = 0.
Em coordenadas cartesianas:
nx∂u∂x
+ ny∂u∂y
= 0,
nx∂v∂x
+ ny∂v∂y
= 0.
***Fig de uma fronteira sintetica
A posicao da fronteira de saıda do escoamento e arbitraria. Ela deve ser colocada deforma que o escoamento na regiao de interesse seja independente de sua posicao (umamaneira de saber se a fronteira esta o mais proximo possıvel, por causa de tempo computa-cional, mas ja com o escoamento desenvolvido, e avaliar um parametro especıfico).
Esta condicao sera imposta de uma forma natural:
R =
∫Ω
( ) dΩ −∫
Γ
φi (n · T)︸ ︷︷ ︸(?)
dΓ,
n · T = n · [−p I + ∇u + ∇uT] = −p n + n · ∇u + n · ∇uT,
79
R =
∫Ω
( ) dΩ −∫
Γ
φi [−p n + n · ∇u︸ ︷︷ ︸0
+n · ∇uT] dΓ.
Substitui-se a condicao de contorno, que e apenas um pedaco do argumento da integralde linha, com isso, permanecem incognitas.
R =
∫Ω
( ) dΩ −∫
Γ
φi [−p n + n · ∇uT] dΓ,
onde, p e u sao as incognitas.
Assim, Rmx e Rmy sao:
Rmx =
∫Ω
( ) dΩ −∫
Γ
φi [−p nx + nx∂u
∂x+ ny
∂v
∂x].
Rmy =
∫Ω
( ) dΩ −∫
Γ
φi [−p ny + nx∂u
∂y+ ny
∂v
∂y].
Observe que os valores de p , ∂u∂x
, ∂u∂y
∂v∂x
, ∂v∂y
na fronteira Γ nao sao conhecidos e as ex-pressoes em termos das funcoes base devem ser usadas no calculo da integral de linha, isto e:
p =∑
Pjχj , ∂u∂x
=∑ ∂φj
∂xUj , ...
CALCULO DA INTEGRAL DE LINHA
∫Γ
φi( ) dΓ :
Considera-se o seguinte elemento, localizado na saıda de um escoamento, onde a condicaode escoamento desenvolvido esta sendo imposta como condicao de contorno.
**FIG
1) Calcular R sem se preocupar com a condicao de contorno.
2) Verificar se alguma condicao tem que ser imposta no elemento.
Sendo a condicao de contorno Natural e necessario acrescentar os termos∫
Γφi( ) dΓ nos
graus de liberdade associados aos nos da fronteira. Assim, acrescenta-se a integral de linhanos resıduos 2, 3, 6, 11, 12 e 15.
R2 = R2 +∫
Γφ2 (n · T)x dΓ, R11 = R11 +
∫Γφ2 (n · T)y dΓ,
R3 = R3 +∫
Γφ3 (n · T)x dΓ, R12 = R12 +
∫Γφ3 (n · T)y dΓ,
R6 = R6 +∫
Γφ6 (n · T)x dΓ, R15 = R15 +
∫Γφ6 (n · T)y dΓ.
80
Utilizando para a integral de linha 3 pontos de Gauss, de forma geral, tem-se:∫Γ
F(x, y) dΓ =
∫ 1
−1
F(ξ = 1, η)dΓ
dηdη,
onde,dΓ
dη=
√(dx
dη
)2
+
(dy
dη
)2
.
Como o valor de ξ esta fixo, a funcao so varia com η, assim, a intergracao em coorde-nadas locais e ao londo da direcao η.
Logo,
∫Γ
φ2 (n · T)x dΓ =
=3∑
IGP=1
φi(ηGP ) −pnx + nx∂u
∂x+ ny
∂v
∂xηGP
dΓ
dη∗ WGP .
CONDICAO DE CONTORNO DE PRESSAO
Observe que nas equacoes so aparecem as derivadas de pressao e portanto o nıvel depressao deve ser fixado nas condicoes de contorno ou artificialmente.
Reu · ∇u = −∇P + . . .
Se P satisfaz a equacao de Navier-Stokes, P1 = P +P0, com P0 constante, tambem satisfaz,pois, ∇P1 = ∇P + ∇P0︸︷︷︸
0
⇒ ∇P1 = ∇P.
Exemplos:
1) Para problemas com superfıcie livre (interface), o nıvel de pressao e especificado au-tomaticamente na condicao de contorno da interface:
**Pode ter uma fig para ilustrar...
n · T = γdt
ds+ nPatm.
2) A pressao pode ser tambem especificada nas fronteiras de entrada e/ou saıda doescoamanto. Neste caso, o termo de pressao que aparece na integral de linha nao deve sercalculado e sim especificado. Por exemplo, se a pressao na fronteira for especificada P ∗, atracao e igual a:
n · T = −P ∗n + n · ∇u + n · ∇uT ,
81
∫Γ
φi (n · T)x dΓ =
∫Γ
φi
(−P ∗n + n · ∇u + n · ∇uT)
xdΓ.
**Pode ter uma fig para ilustrar...
3) No caso de um escoamento sem entrada e nem saıda (como o projeto 2), deve-seespecificar o nıvel de pressao de um ponto do escoamento. Os resıduos associados aos grausde liberdade de pressao correspondem aos resıduos de continuidade. Logo, para especificaro valor da pressao de um dado elemento, R19 passa a ser R19 = P1 − P ∗ e J19j = 0, paraj = 19 e J1919 = 1.
**Pode ter uma fig para ilustrar...
82
Capıtulo 8
Problema de Valor Inicial - RegimeTransiente
Vamos comecar analisando o problema de conducao em regime transiente.
∂u∂t
= ∂2u∂x2 ,
u(x, 0) = u0(x), → condicao inicial.
∂u∂x
(0, t) = ∂u∂x
(1, t) = 0. → condicao de contorno.
Vamos usar a expansao por elementos finitos somente para descrever a variacao de u noespaco, enquanto que na discretizacao do tempo usaremos diferencas finitas.
u(x, t) =N∑
i=1
Ui(t)φi(x) ⇒ ∂u
∂t(x, t) =
N∑i=1
Ui(t)φi(x) ;∂u
∂x(x, t) =
N∑i=1
Ui(t)∂φi(x)
∂x.
Pode-se tambem utilizar o metodo de elementos finitos para discretizar no tempo, poremexistem poucas evidencias de alguma vantagem.
Ri =
∫ 1
0
[ut − uxx] φi dx = 0 ⇒
Integrando por partes:
⇒∫ 1
0
[utφi + uxφix] dx − [uxφi]10︸ ︷︷ ︸
=0 pela c.c.
= 0.
Ri =
∫ 1
0
[utφi + uxφix] dx = 0.
Substituindo as expansoes de u temos:
83
Ri =
∫ 1
0
N∑
j=1
Uj φj φi +N∑
j=1
Uj φjx φix
dx,
=N∑
j=1
(∫ 1
0
φj φi dx
)︸ ︷︷ ︸
Mij
Uj +
(∫ 1
0
φjx φix dx
)︸ ︷︷ ︸
−Kij
Uj,
=N∑
j=1
MijUj − KijUj,
M U = K U ⇒ U = M−1 K U.
M e comumente denotada matriz Massa.
A representacao geral para um problema de valor inicial e :
U = f(t,U) ; U(0) = U0.
comentario : Talvez seja interessante colocar uma secao para os metodos...
Os metodos para a solucao desse tipo de sistema podem ser divididos em duas categoriasgerais :
• METODOS EXPLICITOS → f e calculada no tempo anterior, logo e conhecida(metodo bastante simples, porem instavel).
• METODOS IMPLICITOS → f e calculada no instante presente, logo e desconhecida.
Uma outra classificacao usual e devido as informacoes utilizadas na integracao do tempo:
• Passo-Unico → Toda a informacao utilizada esta no intervalo tn, tn+1.
• Multi-Passo → Usa-se mais de uma solucao em tempos anteriores para aproximar asolucao em tn+1.
84
METODOS DE PASSO UNICO
Metodos Explıcitos
• Diferenca a frente
du
dt≈ un+1 − un
∆t= f(un, tn) ⇒ un+1 = un + ∆t f(un, tn).
Observe que todas as variaveis sao conhecidas explicitamente sem ser preciso calcular nen-hum sistema.
• Runge - Kuta
A ideia e calcular f(u, t) em varios pontos entre tn e tn+1 e combinar os resultados paraoblter uma boa aproximacao de u em tn+1.
Metodos Implıcitos
• Diferenca a Re (Backward Difference).
du
dt≈ un+1 − un
∆t= f(un+1, tn+1) ⇒ un+1 = un + ∆t f(un+1, tn+1)
Depois e necessario resolver o sistema de equacoes, que pode ou nao ser linear.
Observe que na diferenca a frente o problema e montado no tempo n e o instante seguinte(n+1) e usado no calculo da derivada, ja na diferenca a re o problema e montado no instanten + 1 sendo o tempo anterior (n) usado no calculo da derivada.
• Crank-Nickolson (Trapezio)
un+1 − un
∆t=
1
2[f(un+1, tn+1) + f(un, tn)].
METODOS MULTI-PASSO
Formula Geral: un+1 =∑k1
i=0 αi un−i + ∆tn∑k2
i=0 βi f(un+1−i, tn+1−i).
Metodos Explıcitos
• Adams-Bashforth (ordem q)
un+1 = un + ∆tn+1
q−1∑i=1
βi f(un+1−i, tn+1−i).
85
• Adams-Moulton
un+1 = un + ∆tn+1
q−1∑i=0
βi f(un+1−i, tn+1−i).
ESTABILIDADE DOS METODOS
Problema modelo:y = −λy,y(0) = y0.
Solucao analıtica: y = y0e−λt.
Usando quey(t) = y0e
−λt,y(t + 1) = y0e
−λ(t+1),
e fazendo y(t+1)y(t)
temos a solucao exata :
y(t + 1) = y(t)e−λ∆t. (8.1)
Resolvendo por diferenca a frente para comparar:yn+1 − yn = −∆t λyn
⇒ yn+1 = (1 − λ∆t) yn. (8.2)
Erro : εn = y(tn) − yn
εn+1 = y(tn+1) − yn+1
Substituindo y(tn+1) e yn+1 em εn+1 e ainda somando e subtraindo (1 − λ∆t) y(tn), temos:
εn+1 = e−λ∆t y(tn) − (1 − λ∆t) yn + (1 − λ∆t) y(tn) − (1 − λ∆t) y(tn)
⇒ εn+1 = (1 − λ∆t) [y(tn) − yn]︸ ︷︷ ︸εn
+[e−λ∆t − (1 − λ∆t)] y(tn)
⇒ εn+1 = (1 − λ∆t)εn + [e−λ∆t − (1 − λ∆t)] y(tn)
n → ∞ ⇒ |εn+1||εn| = |1 − λ∆t| < 1.
Se |1 − λ∆t| > 1 → o erro cresce com o tempo, e o metodo e dito instavel.
Para o metodo ser estavel ⇒ ∆t < 1λ, essa condicao impoe um limite no passo do tempo.
Pode-se, assim, exigir passos de tempo muito pequenos.
86
FORMULACAO POR ELEMENTOS FINITOS E DIFERENCA A RE
• Problema Linear:
dudt
= d2udx2 ,
u(x, t = 0) = u0(x),
dudx
(x = 0, t) = dudx
(x = 1, t) = 0.
Ri =
∫ 1
0
[du
dt− d2u
dx2] φi dx =
∫ 1
0
du
dtφi +
du
dx
dphiidx
dx.
u(x, t) =N∑
j=1
Uj(t)φj(x).
⇒N∑
j=1
∫ 1
0
φj φi dx︸ ︷︷ ︸Mij
[U(n+1)j − U
(n)j ]︸ ︷︷ ︸
∆t
+N∑
j=1
∫ 1
0
dφj
dx
dφi
dxdx︸ ︷︷ ︸
−Kij
U(n+1)j = 0.
1
∆tM[U(n+1) − U(n)] = K U(n+1) ⇒
⇒ [ 1∆t
M − K]U(n+1) = 1∆t
M U(n) (∗)A cada intervalo de tempo resolver o sistema (∗).
• Problema Nao-Linear:
dudt
+ u dudx
= d2udx2 ,
u(x, t = 0) = u0(x),
dudx
(x = 0, t) = dudx
(x = 1, t) = 0.
Ri =
∫ 1
0
[du
dt+ u
du
dx− d2u
dx2] φi dx
87
=
∫ 1
0
[du
dtφi + u
du
dxφi +
du
dx
dφi
dx] dx.
u(x, t) =N∑
j=1
Uj(t)φj(x).
Ri =
∫ 1
0
N∑j=1
[U(n+1)j − U
(n)j ]︸ ︷︷ ︸
∆t
φj φi +N∑
j=1
U(n+1)j φj
N∑k=1
U(n+1)k
dφk
dxφi +
N∑j=1
U(n+1)j
dφj
dx
dφi
dx
Sendo um sistema na-linear de equacoes em U(n+1) utilizar o metodo de Newton, ouseja, para cada instante de tempo resolver:
∂Ri
∂U(j+1)j
∆U(j+1)j = −Rj,
U(n+1) = U(n+1) + ∆U(n+1).
∂Ri
∂U(j+1)j
=1
∆t
∫ 1
0
φj φi dx︸ ︷︷ ︸M
+
∫ 1
0
u(n+1) ∂φj
∂xφi +
∂u(n+1)
∂xφj φi +
∂φj
∂x
∂φi
∂xdx︸ ︷︷ ︸
J
Para cada instante de tempo t(n+1):
• Chutar U(n+1).• Se R (U(n+1)) = 0 → convergiu → Passa para o proximo instante de tempo.
→ Se nao: ( 1
∆tM + J) ∆U(j+1) = −R (U(n+1),U(n)),
U(n+1) = U(n+1) + ∆U(n+1).
88
