Introducción Al Metodo de Los Elementos Finitos

Livro desenvolvido no Departamento de Mecanica da PUC-Rio

INTRODUCAO AO METODO DEELEMENTOS FINITOS:

Aplicacao em Mecanica dos Fluidos

Marcio da Silveira Carvalho & Juliana Vianna Valerio

0

Conteudo

1 Introducao 1

2 Formulacao Integral 72.1 Formulacao forte e fraca de um problema . . . . . . . . . . . . . . 72.2 Metodo dos Resduos Ponderados . . . . . . . . . . . . . . . . . . . 10

3 Discretizacao e Funcoes Base 193.1 Ponto de Vista Global . . . . . . . . . . . . . . . . . . . . . . . . . 193.2 Ponto de Vista Elementar . . . . . . . . . . . . . . . . . . . . . . . 213.3 Elemento de Segunda Ordem . . . . . . . . . . . . . . . . . . . . . 273.4 Problema Unidimensional . . . . . . . . . . . . . . . . . . . . . . . 30

4 Problema Bi-Dimensional Linear 37

5 Problema Nao Linear 515.1 Solucao de problema nao-linear com uma variavel . . . . . . . . . . 525.2 Sistema de Equacoes Nao-Lineares . . . . . . . . . . . . . . . . . . 555.3 Metodo de Newton aplicado a uma Equacao Diferencial nao linear 575.4 Nocoes de Tecnicas de Continuacao . . . . . . . . . . . . . . . . . . 59

6 Formulacao Integral da Equacao de Navier-Stokes 61

7 Condicoes de Contorno 85

8 Problema de Valor Inicial - Regime Transiente 95

1

2 INTRODUCAO AO METODO DE ELEMENTOS FINITOS

2 Marcio Carvalho e Juliana Valerio

Captulo 1

Introducao

Varios fenomenos da natureza podem ser descritos, utilizando-se leis da fsica,em termos de equacoes diferenciais. Em algumas situacoes extremamente simplespode-se encontrar solucoes analticas, nas outras busca-se solucoes aproximadaspara o problema.

Os metodos para a obtencao de solucoes aproximadas de equacoes diferenciaispodem ser divididos em duas classes gerais:

(i)Metodo de Diferencas Finitas: As derivadas da funcao sao substitudas peladiferenca dos valores que a funcao assume em certos pontos do domnio. A equacaodiferencial da origem a um sistema de equacoes algebricas onde as incognitas saoos valores da funcao nos pontos escolhidos para representar as derivadas. As ideiasbasicas de Metodos de Diferencas Finitas 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 numeronito de pontos do domnio.

Se o domnio for dividido em N 1 intervalos atraves de N pontos (nos), deve-se obter o valor aproximado de u(x) em cada no, ou seja, um problema com Nincognitas. As N equacoes necessarias para resolver este problema sao obtidasaproximando a derivada em funcao dos valores nodais desconhecidos nas posicoesdos N nos.

Olhando nas proximidades do no 4:

du

dx

x4

=U5 U32x

.

1


Figura 1.1: Aproximacao da derivada da funcao por diferencas nitas.

A equacao diferencial no ponto x4 e aproximada como:

x4

[U5 U32x

]+ U4 = 0 x42xU3 + U4 +

x42x

U5 = 0.

O resultado deste procedimento sera um sistema linear N N , onde as incognitassao os valores de U nos N nos.

(ii)Metodos Variacionais: A solucao aproximada e escrita como combinacaolinear de funcoes base apropriadamente escolhidas:

u(x) =Ni=1

cii.

Os coecientes ci sao obtidos de forma que a equacao diferencial seja satisfeita.Para obter a solucao aproximada atraves de um metodo variacional, deve-se rescr-ever o problema em sua forma fraca (integral ou variacional). A necessidade do usoda formulacao fraca e ilustrada no Exemplo 2. Os diferentes metodos variacionaisutilizam diferentes formas integrais, funcoes base e funcoes peso.

O Metodo de Elementos Finitos se encaixa na classe de metodos variacionais,onde as funcoes base sao diferentes de zero apenas em uma pequena parte dodomnio, como sera mostrado a seguir.

A simples substituicao da aproximacao u(x) =N

i=1 cii na equacao diferen-cial nao garante a obtencao de N (numero de coecientes ci) equacoes algebricas


INTRODUCAO AO METODO DE ELEMENTOS FINITOS 3

linearmente independentes, por isso em metodos variacionais utiliza-se a formaintegral do problema, como mostrado no exemplo a seguir.

Exemplo 2:

ddx

(xdudx

) + u = 0 para 0 < x < 1,

sujeito a`s condicoes de contorno:

u(0) = 1 ; (xdudx

)x=1 = 0.

Solucao aproximada na forma:

u(x) = 0(x) +2

i=1

cii(x),

onde as funcoes base escolhidas sao:0(x) = 1,

1(x) = x2 2x,2(x) = x3 3x.

Precisamos encontrar os coecientes c1 e c2, o problema possui 2 incognitas.Entao substituindo 0, 1 e 2, a solucao aproximada torna-se:

u(x) = 1 + c1[x2 2x] + c2[x3 3x], derivandodu(x)dx

= 2c1x 2c1 + 3c2x2 3c2.Observe que as condicoes de contorno sao satisfeitas para qualquer valor de c1

e c2, que devem ser determinados de forma que a aproximacao satisfaca a equacaodiferencial em algum sentido.

Obrigando a solucao aproximada a satisfazer a equacao diferencial no sentidoexato, isto e, em todos os pontos do domnio:

ddx{x[c1(2x 2) + c2(3x2 3)]}+ 1 + c1(x2 2x) + c2(x3 3x) = 0,

2c1(x 1) 3c2(x2 1) 2c1x 6c2x2 + c1(x2 2x) + c2(x3 3x) + 1 = 0,1 + 2c1 + 3c2 + (6c1 3c2)x + (9c2 + c1)x2 + c2x3 = 0.

Para essa expressao ser valida para qualquer x, o coeciente de cada potenciade x deve ser nulo:

1 + 2c1 + 3c2 = 0,6c1 3c2 = 0,9c2 + c1 = 0,c2 = 0.



Esse sistema nao tem solucao! Vamos obrigar a solucao aproximada a satisfazera equacao diferencial num sentido integral ao inves de no sentido exato.

Podemos obriga-la a satisfazer a seguinte integral ponderada: 10

w(x){ ddx

(xdu

dx) + u}dx = 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 numerode funcoes peso utilizadas (tambem linearmente independentes). No exemplo es-colhemos 2 funcoes peso linearmente independentes: w = 1 e w = x , teremos:

10

1( d

dx{x[c1(2x 2) + c2(3x2 3)]}+ 1 + c1(x2 2x) + c2(x3 3x)

)dx = 0, 1

0

x

( d

dx{x[c1(2x 2) + c2(3x2 3)]}+ 1 + c1(x2 2x) + c2(x3 3x)

)dx = 0.

Entao: {2/3c1 + 5/4c2 = 1,3/4c1 31/20c2 = 1/2,

e nalmente, c1 = 222/23 e c2 = 100/23.

As funcoes bases utilizadas e a solucao aproximada estao ilustradas na gura(1.2):



Figura 1.2: Representacao das funcoes bases e solucao aproximada.


Captulo 2

Formulacao Integral

2.1 Formulacao forte e fraca de um problema

A formulacao forte de um problema de valor de contorno e denida 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 denida como:Determine u U , tal que

(W ) 10

dw

dx

du

dxdx =

10

fwdx + w(0)h ; para qualquer w V .U e o conjunto de funcoes teste denido como:

U = {u|u(1) = g e 10

(du

dx)2dx


Uma funcao f : R pertence ao espaco Ck() se djfdxj e contnua, com0 j k. C0 e o espaco das funcoes contnuas.

O espaco Sobolev de funcoes Hk() e denido como:

Hk() = {w|w L2; dwdx

L2; ...; dkw

dxk L2}, onde

L2() = {f |

f2d


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 +

10

wfdx, w V

Integrando por partes obtem-se:

10

wd2u

dx2dx + w(1)

du(1)dx

w(0)du(0)dx

= 10

wfdx + w(0)h

10

w[d2u

dx2+ f ]dx w(1)du(1)

dx+ w(0)[

du(0)dx

+ h] = 0.

Como w V w(1) = 0

10

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 (2.1) e valida para qualquer w V , vamos escolher um w deforma a concluirmos que d

2udx2 + 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 denicao) e [d2u

dx2+ f ]2 0,

entao a integral sera zero [d2u

dx2+ f ]2 = 0 d

2u

dx2+ f = 0.

Assim temos (i). Vamos ver (ii):

Usando (i) a expressao (2.1) ca 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) vericados, temos que se u e solucao de (W ) entao e tambemsolucao de (S).

Como mencionado anteriormente, o metodo de Elementos Finitos e baseado naformulacao fraca do problema.

A ideia basica e restringir os espacos U e V para espacos de dimensao nita.Isto implica que qualquer u U pode ser escrito como:

u(x) =Ni=1

cii,

onde os is formam uma base de U .

A m de simplicar, e comum seguir a notacao:

a(u,w) = 10

dw

dx

du

dxe (f, w) =

10

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 Resduos Ponderados

Neste curso vamos nos concentrar no Metodo dos Resduos Ponderados (WeightedResiduals).

Vamos novamente analisar o problema:

d2udx2 = f, em (0, 1);u(1) = g = 0 considerar g = 0, apenas para simplicar.du(0)dx = h.

Queremos encontrar uma aproximacao uh para a solucao do problema u.uh vai pertencer a um espaco de dimensao nita Uh, logo podemos escreve-lo

como:

uh(x) =Ni=1

cii,

onde N e a dimensao do espaco e is as funcoes base de Uh. Queremos tambemque Uh U logo, uh(1) = 0.



Como uh e apenas uma aproximacao de u, d2u

dx2 = fVamos chamar de resduo (R) da equacao, a medida de quanto uh aproxima

adequadamente o problema.

R =d2uhdx2

+ f,

Observe que se uh = u entao, R = 0.

Metodos dos Resduos Ponderados sao obtidos atraves da seguinte integral:

10

wR = 0

Observe que integrando o resduo com o peso w, estamos fazendo com que suamedia ponderada seja 0. Por isso a condicao e dita FRACA, pois nao estamosobrigando ponto a ponto que u = uh e sim na media (na integral).

Substituindo R, temos:

10

w[d2uhdx2

+ f ]dx = 0

integrando por partes:

10

dw

dx

duhdx

dx + 10

fwdx + w(0)h = 0 Equivalente a` formulacao variacional.

A integracao por partes alivia as restricoes impostas a` funcao aproximada uh.A expressao passa a depender da primeira derivada de uh, duhdx e nao mais da se-gunda d

2uhdx2 . Logo, a restricao sobre uh devera ser que sua derivada pertenca a

L2, duhdx L2, isto e, o quadrado de duhdx seja integravel. Observe que ate funcoescontnuas com descontinuidades nas derivadas podem satisfazer a estas condicoes,como ilustrada na gura(2.1):

Como uh Uh uh(x) =N

i=1 cii.

E entao, duhdx =N

i=1 cididx

Sendo o espaco Uh de dimensao N precisaremos de N equacoes linearmenteindependentes a m de determinar os coecientes ci.

As funcoes peso w devem pertencer a um espaco nito Vh V . Esse espacopode ser representado pela base i. Cada equacao esta relacionada com as funcoespeso i.



Deste modo as N equacoes L.I. para determinar os ci sao:

Ri = 10

didx

N

j=1

cjdjdx

dx + 1

0

fidx + i(0)h = 0

N

j=1

{ 10

didx

djdx

dx

}

Aij

cj = 10

fidx + i(0)h bi

; i = 1, ..., N.

Temos matricialmente Ac = b

Algumas referencias chamam A de matriz de rigidez.

O Metodo dos Resduos Ponderados recebe diferentes denominacoes de acordocom as funcoes 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 coecientes ci da aproximacao uh(x) =N

i=1 cii sao obtidos atraves daminimizacao da integral do quadrado dos resduos:

ci

10

R2dx = 0 10

R

ciRdx = 0.

i =R

ci.

METODO DE COLOCACAO

O resduo R e identicamente nulo em N pontos do domnio.

R(xi) = 0.

Tambem e um caso particular de metodo dos Resduos Ponderados.

wi = (x xi) 10

(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) =Ni=1

cii 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 Exerccio #1

R(x) = d2uhdx2

uh + x2.Queremos que: 1

0

w(x)R(x)dx = 0, onde w(x) e o vetor das funcoes peso i.

10

i[d2uhdx2

uh + x2]dx = 0

10

id2uhdx2

dx 10

i uh dx + 10

i x2 dx = 0



Integrando por partes, evitando assim a segunda derivada,temos:

duhdx

i

10

10

(didx

duhdx

)dx 10

i uh dx + 10

i x2 dx = 0

Usando que:

uh(1) = 1.i(0) = 0.

i(1) + 10

(didx

duhdx

)dx 10

i uh dx + 10

i x2 dx = 0

Como u(x) =N

i=1 cjj , entaodudx =

Ni=1 cj

jdx ; substituindo temos:

Ni=1

cj

10

[didx

djdx

ij]dx

Aij

= i(1) 10

i x2 dx

bi

Pelo Metodo de Galerkin i = i; entao,

Aij = 10

[didx

djdx ij

]dx e bi = i(1)

10i x

2 dx

Assim, usando1(x) = x ; d1dx = 12(x) = x2 ; d1dx = x3(x) = x3 ; d1dx = 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 .

c =

22801777

2031777

1757108

O graco mostra a solucao exata e a aproximada. Observe que a solucaoaproximada representa muito bem a solucao exata, utilizando-se apenas de tresfuncoes base.E evidente que a escolha de funcoes base apropriadas e fundamentalpara o bom desempenho do metodo.



Figura 2.1: Funcoes contnuas com derivadas descontnuas.



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


Captulo 3

Discretizacao e Funcoes Base

3.1 Ponto de Vista Global

Podemos perceber que para obter os elementos da matriz A, precisamos integraras funcoes teste e peso (e/ou suas derivadas) em todo o domnio.

O metodo dos elementos nitos consiste na escolha apropriada de funcoes i ei, de tal forma que as funcoes e suas derivadas so sao diferentes de zero em umapequena parte do domnio. Deste modo, a integral no domnio ca reduzida a umaintegral na parte do domnio onde a funcao e diferente de zero, como mostrado nagura(3.1) :

10

( )idx = x2x1

( )idx

A ideia e dividir o domnio em elementos, denidos por nos. Por exemplo: odomnio de 0 a 1 foi dividido em 8 elementos (9 nos).

Uma escolha simples que facilita a computacao e denir i e i de tal formaque:

i(xj) = ij i(xj) ={

1 se i = j,0 se i = j,

conforme mostrado na Fig.3.3 (Elementos Lagrangeanos).Desta forma o numero de funcoes base (dimensao do espaco) esta associado ao

numero de nos da malha. Quanto mais renada a malha, maior sera a dimensaodo espaco de funcoes e desta forma, mais precisa sera a solucao aproximada.

Se as funcoes teste i sao denidas da maneira mostrada anteriormente, oscoecientes ci passam a ter um signicado especial :

19


Figura 3.1: Exemplo de funcao base diferente de zero apenas em uma pequenaporcao do domnio.

u(xi) =N

j=1

cj j(xi) ij

=N

j=1

cjij = ci

ci = u(xi)

O coeciente ci e igual ao valor da funcao u no ponto xi.

Uma outra consequencia importante da denicao das funcoes teste e que varioselementos da matriz A sao iguais a zero, conforme ilustrado na Figura 3.4.

Por exemplo:

A26 = 10

d2dx

d6dx

dx = 0

A matriz A, do exemplo, passa a ser uma matriz esparsa:



Figura 3.2: 8 elem. 9 nos

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 nitos de uma forma global, isto e, cadafuncao i e denida em todo o domnio com a particularidade de ser diferente dezero somente em uma pequena parte dele.

Como cada funcao i so e diferente de zero em uma pequena parte do domnio,e mais eciente se a tratarmos nestas opcoes apenas. Isto leva a uma descricaolocal do metodo (elementar).

No exemplo mostrado anteriormente (Figura 3.3), existem duas funcoes difer-entes de zero em cada elemento e uma mesma funcao e diferente de zero em doiselementos, conforme mostrado na tabela a seguir:



Figura 3.3: i(xj) = ij

elemento no inicial no nal is = 01 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 globalA seja escrita 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

uma matriz 2 2.Por exemplo, a matriz do elemento 3, A(3), e igual a:

A(3) =

x4x3

d3dx

d3dx dx

x4x3

d3dx

d4dx dx x4

x3

d4dx

d3dx dx

x4x3

d4dx

d4dx dx

A forma de cada funcao nos elementos e sempre a mesma, independentementedo elemento em questao. Isto sugere uma troca de coordenadas para que possamos



Figura 3.4: Calculo de A26.

fazer uma descricao elementar das funcoes.

Descricao Elementar: O elemento se extende de = 1 a = +1 e possui doisnos como mostrado na gura:

As funcoes base elementares sao escritas como:

(e)1 () =

12

(e)2 () =

1+2

Assim d(e)1

d = 12 ;d

(e)2

d =12

Todos os elementos do domnio sao mapeados para esse elemento padrao.

Por exemplo, o mapeamento do elemento 3, que vai de x3 a x4, para o elementopadrao, de = 1 a = 1, pode ser escrito como:



Figura 3.5: Funcoes base lineares no sistema de coordenadas local.

x = x3 = 1,

(x) = 2x(x3+x4)x4x3 ,

x() = (x4x3) +x3+x42 ,

x = x4 = 1.

As integrais sao reescritas como:

A(3)12 =

x4x3

d3dx

d4dx

dx = 11

d(e)1

d

d

dx

d(e)2

d

d

dx

dx

dd =

= 11

d(e)1

d

d(e)2

d

d

dx= 2x4x3 =

2h(e)

d,

onde h(e) e o tamanho do elemento.

2h(3)

11

d(e)1

d

d(e)2

dd.

=para todos os elementos

A matriz do terceiro elemento entao ca:



A(3) =2

h(3)

11

d1d

d1d d

11

d1d

d2d d 1

1d2d

d1d d

11

d2d

d2d d

,

onde os A(3)ij sao:

A(3)11 =

2h(3)

11

(12

)(12

) d =2

h(3)142 =

1h(3)

A(3)12 =

2h(3)

11

(12

)(12) d =

1h(3)

= A(3)2,1

A(3)22 =

2h(3)

11

(12)(

12) d =

1h(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,precisamos saber a correspondencia entre a numeracao local dos nos no elementocom a numeracao global dos nos no domnio.

Logo, vamos denir a matriz DomNodeID (

#local dos nos ilnode , iele

#do elemento

) que fornece

a numeracao global do no local ilnode do elemento iele.

DomNodeID

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 locale a global.

Por exemplo, a matriz local do elemento 3 (2 2) e montada na seguinteposicao 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 nitos pode ser escrito como:

do iele = 1, NELE NELE # de elementoscall 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 ilnodedo jlnode = 1, NLOCALNODES

jgnode = DomNodeID (jlnode, iele)A(ignode, jgnode) = A(ignode, jgnode) + Aelem(ilnode, jlnode)

end

endend

Nota : Repare que a entrada A55 da matriz global tem contribuicao de 2 ele-mentos (elementos 4, 5).

A55 = A(4)22 + A

(5)11



Isto ocorre, pois:

A55 = 10

d5dx

d5dx

dx = x5x4

d5dx

d5dx

dx + x6x5

d5dx

d5dx

dx = 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 funcaoi foi denida tal que :

i(xj) = i,j ,

ou seja, 1(x1) = 1, 2(x1) = 0,1(x2) = 0, 2(x2) = 1.

dois pontos denem uma reta i e linear.

Figura 3.6: Elemento linear

Se desejamos ter uma aproximacao melhor podemos denir funcoes base quadraticasdentro de cada elemento. Como para denir uma funcao quadratica sao necessarios3 pontos, os elementos quadraticos possuem 3 nos.

Repare que i(xj) = i,j .



Figura 3.7: Elemento Quadratico

Nas coordenadas elementares, as funcoes base sao escritas como:

1() =( 1)

2,

2() =( + 1)

2,

3() = ( + 1)( 1).Logo, derivando-as temos:d1d

= 12

;d2d

= +12

;d3d

= 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)

11

d1d

d1d d

11

d1d

d2d d

11

d1d

d3d d 1

1d2d

d1d d

11

d2d

d2d d

11

d2d

d3d d 1

1d3d

d1d d

11

d3d

d2d d

11

d3d

d3d d



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

. . .

3.4 Problema Unidimensional

Vamos mostrar todos os passos para a solucao de um problema de conveccao/difusaousando o metodo dos elementos nitos.

d2udx2 Pe dudx = 0, 0 < x < 1;

du(0)dx = u(0) 1;

u(1) = 0.

Formulacao Fraca (Metodo de Galerkin)

Ri = 10

(d2u

dx2 Pe du

dx

)i dx = 0

10

du

dx

didx

dx +du

dx(1)i(1) du

dx(0)i(0) Pe

10

du

dxi dx = 0

10

(du

dx

didx

+ Pedu

dxi

)dx =

du

dx(1)i(1) du

dx(0)

u(0)1

i(0)

Ri = 10

(du

dx

didx

+ Pedu

dxi

)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) =Ni=1

Uii sem se preocupar com as condicoes de contorno.



Figura 3.9: Discretizacao do domnio atraves de 8 elementos lineares.

Ri =N

j=1

{ 10

(didx

djdx

+ Pedjdx

i)dx}

Uj = N

j=1

Ujj(0) 1i(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

{ 10

(didx

djdx

+ Pedjdx

i)dx}

Uj = [U1 1]i(0).

c.c. essencial u(1) = 0 Nj=1 Ujj(1) , lembrando que x = 1 e o noN UNN (1) = 0 UN = 0.

i = 1

R1 N

j=1

{ 10

(d1dx

djdx

+ Pedjdx

1)dx}

Uj = [U1 1]

i = 2, ..., N 1

Ri N

j=1

{ 10

(didx

djdx

+ Pedjdx

i)dx}

Uj = 0

i = N RN UN = 0



. . .

.... . .

0 0 0 0 0 0 0 0 1

U1U2U3U4U5U6U7U8U9

=

100000000

Ponto de Vista Local:

Formulacao Fraca: 10

(du

dx

didx

+ Pedu

dxi) dx = [u(0) 1]i(0)

condicao de contorno.

.

O ideal e montar as matrizes e vetores elementares sem se preocupar comas condicoes de contorno e depois modica-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=

2h(e)

.

A(e)ij =

2h(e)

11

djd

did

d + Pe 11

djd

id; f(e)i = 0.

(e)1 () =

1 2

d1d

=12

,

(e)2 () =

1 + 2

d2d

=12.

A(e)11 =

2h(e)

11

(12) (1

2) d + Pe

11

(12) (

1 2

) d =

=2

h(e)

11

14

d 11

Pe

4(1 ) d =

=1

h(e) Pe

2e assim por diante...



Apos obter as matrizes elementares, temos que nos preocupar com as condicoesde contorno.

Os termos correspondentes a`s condicoes de contorno se anulam em todos ostermos menos 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: 10

( ) = [U1 1]

N

j=1

( ) Uj + U1 = +1

A(1)1,1 = A(1)1,1 + 1f(1)1 = f

(1)1 + 1

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 aolongo do elemento. De um modo geral, temos que calcular: 1

1F ()d.

Para um problema simples, podemos calcular a integral analiticamente. Porem,e mais interessante faze-lo numericamente, pois engloba a resolucao nao so de prob-lemas simples como tambem complicados e ainda temos um programa generico.

Um metodo comumente utilizado e o da QUADRATURA GAUSSIANA, onde

11

F ()d =NGPi=1

F ()Wi

onde NGP e o numero de pontos utilizados para a integracao. Quanto maior onumero de pontos, 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.5553 0.0 0.888

+0.77459 0.555

PROJETO #1

d2udx2 Pe dudx = 0du(0)dx = u(0) 1

u(1) = 0

A solucao analtica 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 a seguinte funcao para gerar a malha concentrada:

xi = 1[(NNODES i)(NNODES 1)

]aa parametro que regula a concentracao da malha.

(d) Compare as solucoes obtidas com a solucao analtica em cada ponto damalha.

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 %

SUBROTINE SOLUTION

do iele = 1, NELE NELE # de elementoscall 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 ilnodedo jlnode = 1, NLOCALNODES

jgnode = DomNodeID (jlnode, iele)A(ignode, jgnode) = A(ignode, jgnode) + Aelem(ilnode, jlnode)

end

endendCALL 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 doend do

end do


Captulo 4

Problema Bi-DimensionalLinear

Vamos considerar o problema 2-D de conducao de calor. (kT ) = 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)

37


n T = q Tn |2 = q c.c Newman (uxo prescrito)

Formulacao Fraca : Ache T U tal que:

w2Td = 0 , w V ;onde:U = {T

|T |2d


w T d +

2w q d2 = 0

Discretizacao Uh U Uh de dimensao N .

Th =N

j=1

Cjj ,

onde C js sao as N incognitas.

As N equacoes linearmente independentes serao obtidas escolhendo-se N funcoespeso tambem linearmente independentes :

Ri i T d =

2

i q d2

Sendo Th =N

j=1 Cjj , entao, T =N

j=1 Cjj :

i (N

j=1

Cjj) d =2

i q d2

Nj=1

Cj

{

i j d}

Aij

=2

i q d2 bi

Aijcj = bi

Aij =

i j d

bi =2

i q d2

Metodo dos Elementos Finitos

Escolher i tal que (xj) = ij i(xj) ={

1 se i = j0 se i = j

Ponto de Vista Elementar

Elemento Bi-Linear



Figura 4.2: Funcoes base bi-dimensional (xj) = ij

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)

[ix

jx

+iy

jy

]d

A matriz A(e) e uma matriz 4 4 no caso de elementos bi-lineares.

Precisamos saber as derivadas das funcoes base, ix ;iy , para i = 1..4, ja

que os is sao denidos em termos de e .

Da mesma forma que no caso 1-D, precisamos de um mapeamento do sistemade coordenada local (, ) para o sistema de coordenadas global (x, y).

Uma vez obtida as matrizes elementares, temos que montar a matriz global.O processo de montagem e analogo ao caso unidimensional. Precisamos da matriz



Figura 4.3: Funcoes base de elemento bi-linear.

DomNodeID(ilnode, iele) = numero global do no local ilnode do elemento iele.

Esta matriz e obtida a partir da numeracao dos elementos e nos usados paradividir o domnio de calculo. A numeracao tanto dos nos como dos elementos ealeatoria. Como essa numeracao esta diretamente relacionada com a estruturada matriz global, e interessante usar uma numeracao otimizada, que minimize abanda da matriz global.

A Fig.4.4 ilustra um exemplo de numeracao dos elementos e nos. A matrizDomNodeID para 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



Figura 4.4: Exemplo de numeracao dos nos e elementos.

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 foi dito. Um esquema NAO OTIMIZADO pode levar a uma matriz globalcom banda alta, como mostrado na Fig.4.5. Existem diversos algoritmos paraotimizar a numeracao dos nos. Estes metodos nao serao tratados neste curso.

Figura 4.5: Dois esquemas de numeracao de nos e suas respectivas matriz global.

Para calcular elementos da matriz elementar, precisamos calcular as derivadasdas funcoes base com relacao a x e y:



A(e)ij =

(e)

[ix

jx

+iy

jy

]d

Mas, a funcao (, ) esta denida em termos de e , logo temos que mudarde coordenadas e para isso vamos usar um mapeamento que escreva x e y comofuncao de e :

x = x(, ) e y = y(, ).

MAPEAMENTO ISOPARAMETRICO:

Figura 4.6: Mapeamento isoparametrico

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 denido 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 inter-polacao ultilizada 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.

Uma vez conhecendo o mapeamento entre as coordenadas, pode-se determinar: ix e

iy .

i =

ix

x +

iy

y

i =

ix

x +

iy

y

i

i

=

x

y

x

y

Jacobiano da Transformacao

ix

iy

Logo:

ix

iy

= J1

i

i

; J1 = 1|J|

y yx x

Inverso do Jacobiano

|J| = x y y x .



Assim:

ix

=1

x

y y x

{y

i

y

i

}iy

=1

x

y y x

{x

i

+x

i

}

Para um mapeamento isoparametrico temos que:

x =

Ni=1 Xi

i ;

x =

Ni=1 Xi

i ; . . .

Logo, a integral elementare

[ix

jx +

iy

jy

]dxdy pode ser reescrita em

funcao das variaveis locais:

A(e)ij =

(e)

i j dxdy =(e)

J1i J1j |J| d d

A(e)ij =

11 11{(y

i y i )(y j y j )+

+(x i + x i )(x j + x j )}

1|J|dd

Como a expressao para o calculo da matriz elementar contem |J| no denom-inador, por isso e importante que |J | 0. Se o elemento car muito distorcidoisto pode acontecer e a precisao do calculo de A(e)ij ca comprometida, como nocaso mostrado a seguir:

Figura 4.7: Elemento distorcido.

dA = |J|dd = |dr1||dr2|sen



|J| = |dr1||dr2|sendd

Quando 0 ou 180 graus , |J| 0. Deve-se monitorar o valorde |J|. Se car negativo ou muito pequeno, os calculos devem ser interrompidos euma nova malha deve ser criada.

Integracao Numerica

Quadratura Gaussiana:

11 11 F (, )dd =

NGPIigp=1

NGPJjgp=1 F (i, j) wI wJ

Figura 4.8: Localizacao dos Pontos de Gauss.



ELEMENTO NUM.PONTOS PESOS COORDENADASINTEGRACAO 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

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 Xii

y =

NLNi=1 Yi

i

x =

NLNi=1 Xi

i

y =

NLNi=1 Yi

i

Calcular detJ = xy y x

Calcular ix = ...iy = ...

do ilnode = 1, NLOCALNODESdo jlnode = 1, NLOCALNODES

A(ilnode, jlnode) = A(ilnode, jlnode)+W

(ix jx + iy jy

) detJ

end do



end doend do

end do

Para executar os passos de calculo de forma eciente, pode-se utilizar umaestrutura matricial para armazenar e operar as variaveis. Um exemplo de tal es-trutura e descrito a seguir.

Armazenar as funcoes i(, ), i , i da seguinte forma:

PHI(i) =

1234

GradPHI(i, j) =

1

2

3

4

1

2

3

4

Armazenar as coordenadas dos pontos nodais de cada elemento da seguinteforma:

XY (i, j) =[

X1 X2 X3 X4Y1 Y2 Y3 Y4

] 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 Y1X2 Y2X3 Y3X4 Y4

detJ = JAC(1, 1) JAC(2, 2) JAC(1, 2) JAC(2, 1)

J1 =

JACINV (1, 1) = JAC(2, 2)/detJJACINV (1, 2) = JAC(1, 2)/detJJACINV (2, 1) = JAC(2, 1)/detJJACINV (2, 2) = JAC(1, 1)/detJ

As derivadas das funcoes base em relacao a`s coordenadas global sao ar-mazenadas na matriz



GradPHI XY (i, j) =

1x 2x 3x 4x

1y

2y

3y

4y

GradPHI XY (i, j) = J1 GradPHI

1|J|

y y

x

x

1

2

3

4

1

2

3

4

=

1x 2x 3x 4x

1y

2y

3y

4y

O algoritmo pode ser escrito como:

CALL BASISFUNC (XI,NETA,PHI,GradPHI)CALL MAPPING (XY,PHI,GradPHI, JAC, JACINV, detJ)GradPHI XY (i, j) = J1 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

end do

EXERCICIO # 3

Escreva a rotina BASISFUNC 1 para elementos bi-lineares e BASISFUNC 2para elementos quadraticos.

EXERCICIO # 4

Escreva a rotina MAPPING indicada no algoritmo acima.


Captulo 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 Resduos Ponderados (Galerkin),obtem-se:

Ri = 10

(udu

dx+

d2u

dx2+ f

)i dx = 0

10

du

dx

didx

dx +du

dx(1) i(1)

=0

dudx

(0) i(0) =0

+ 10

udu

dxidx +

10

fidx = 0,

Ri = 10

du

dx

didx

dx + 10

udu

dxidx +

10

fidx = 0.

Substituindo a expansao da funcao u(x) =N

j=1 Cjj(x):

Ri = 10

Nj=1

Cjdjdx

didx

dx+ 10

N

j=1

Cjj

( N

k=1

Ckdkdx

)idx+

10

fidx = 0.

51


Trocando a ordem da integracao e somatorio em j, a expressao pode ser ree-scrita como:

Nj=1

Cj

{ 10

djdx

didx

dx + 10

(N

k=1

Ckdkdx

)jidx

Aij

}= 10

fidx. bi

A expressao acima representa o produto matricial Ac = b. Porem, como oproblema e nao linear, os elementos da matriz A dependem dos coecientes Ck,solucao do problema. A notacao mais precisa de representar o problema acimaseria

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 deequacoes lineares. Muitos programas de simulacao numerica utilizam um metodoiterativo para resolver o sistema linear Ac = b oriundo de um problema linear. Avantagem de metodos iterativos e a economia de memoria, pois a matriz A nao earmazenada inteira na memoria, e as grandes desvantagens sao a nao garantia deconvergencia e o tempo computacional. O que estamos analisando neste captulosao metodos iterativos para resolver um sistema nao linear. A princpio, nestecurso, sistemas lineares sao resolvidos por metodos diretos (alguma variacao deeliminacao 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 comoMetodo de Picard.

METODO DE PICARD:

A funcao f(x) e reescrita como f(x) = x g(x). A solucao do problema, x,satisfaz a equacao x = g(x).



x0 chute inicialRepetir ate |xn+1 xn| <

xn+1 = g(xn)

Figura 5.1: Interpretacao geometrica do Metodo de Picard.

A interpretacao geometrica do metodo de Picard e apresentada na Figura 5.1.O metodo de Picard pode divergir ate quando o valor inicial e bem proximo dasolucao do problema, como no caso 3 da Figura 5.1 O criterio de convergencia dometodo:

CONVERGE |g(x)| < 1,



DIVERGE |g(x)| > 1.

METODO DE NEWTON:

O procedimento iterativo e determinado atraves da expansao da funcao f(x)em serie de 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 inicialRepetir ate |xn+1 xn| <

xn+1 = xn f(xn)f (xn)

A interpretacao geometrica do metodo de Newton e ilustrada na Figura 5.2.

Figura 5.2: Interpretacao Geometrica do Metodo de Newton.

O metodo de Newton possui uma propriedade de convergencia muito forte. Seo valor inicial 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 desenvolvida a seguir.

Prova

Vamos assumir que x e uma raiz simples: f (x) = 0.

Vamos denir o erro como n = xn x. Logo x = xn n.

Expandindo a funcao f(x) por serie de Taylor ao redor da solucao x:

f(x) = 0 = f(xn) + (x xn)f (xn) + 12(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 =122n

f ()f (xn)

.

Logo, tomando o limite quando x x:

n+1 =122n

f (x)f (x)

,

limn

n+12n

= C n+1 = C2n.

Por exemplo, se na quarta iteracao, a aproximacao x4 estiver perto da solucaodo problema e 4 = 103, pela prova acima, 5 = 106 e 6 = 1012.

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 serreescrita como

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

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. Oprocedimento iterativo converge se:

D(x) 1; Dij = gixj

.

METODO DE NEWTON:

Como no caso de uma equacao apenas, cada expressao do sistema deve ser de-senvolvida por serie de Taylor ao redor da solucao do problema x = (x1, x

2, ..., x

n):

f1(x) = f1(x(k)) + f1x1 (x1 x(k)1 ) + f1x2 (x2 x

(k)2 ) + ... +

f1xn

(xn x(k)n ) = 0;f2(x) = f2(x(k)) + f2x1 (x

1 x(k)1 ) + f2x2 (x2 x

(k)2 ) + ... +

f2xn

(xn x(k)n ) = 0;...fn(x) = fn(x(k)) + fnx1 (x

1 x(k)1 ) + fnx2 (x2 x

(k)2 ) + ... +

fnxn

(xn 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 denida como:

Jij =fixj

.

O algoritmo do metodo de Newton para um sistema de equacoes nao-linearespode ser escrito 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)

5.3 Metodo de Newton aplicado a uma EquacaoDiferencial nao linear

Como mostrado no incio do captulo, no caso da solucao de uma equacao diferen-cial nao linear pelo metodo de elementos nitos, a discretizacao leva a um sistemade equacoes algebricas nao linear:

Ri = Ri(C1, C2, ..., CN ) = 0 ; para i = 1, ..., N.

Cada equacao do sistema corresponde a um resduo ponderado:

Ri = 10

du

dx

didx

dx + 10

udu

dxidx +

10

fidx = 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 equacaoa cada incognita do problema.

Ji,j =RiCj

=

Cj

{Ri =

10

du

dx

didx

dx + 10

udu

dxidx +

10

fidx}

Ji,j = 10

didx

Cj

(du

dx

)dx +

10

{

Cj(u)

du

dx+ u

Cj

(du

dx

)}idx.



A derivada da funcao u e de sua derivada dudx em relacao aos coecientes Cjsao calculados da seguinte forma:

Cj(u) =

Cj

(N

k=1

Ckk

)=

Cj(C11 + ... + Cjj + ... + CNN ) = j ,

Cj

(du

dx

)=

Cj

(N

k=1

Ckdkdx

)=

Cj

(C1

d1dx

+ ... + Cjdjdx

+ ... + CNdNdx

)=

djdx

.

Logo, os elementos da matriz Jacobiana para este problema sao dados por:

Ji,j = 10

djdx

didx

dx + 10

{j

du

dx+ u

djdx

}idx.

O algoritmo de solucao de um problema nao linear pelo metodo de elementosnitos acoplado com o metodo de Newton pode ser escrito como mostrado a seguir:

Ler Chute inicial C(0)

do iter = 1, NITERCALL FORMRV (C(iter1))

Rotina para calculo do vetor resduo globalCALL FORMJM (C(iter1))

Rotina para calculo da matriz Jacobiana globalResolver o sistema JC = RC(iter) = C(iter1) +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 JC = Rdeve ser resolvido.

Se o problema for linear, o metodo de Newton converge em uma iteracao.

A visao elementar do problema ca inalterada. No lugar de calcular a matrize o vetor elementar A(e) e b(e), deve-se calcular a matriz jacobiana elementar J(e)

e o vetor resduo elementar R(e).



5.4 Nocoes de Tecnicas de Continuacao

O raio de convergencia do metodo de Newton e tipicamente menor do que o raiode convergencia do metodo de Picard. Isto signica que o metodo de Newtonsomente ira convergir se o chute inicial for proximo a` solucao do problema. Porem,se o chute inicial for adequado, o metodo de Newton converge quadraticamente,levando poucas iteracoes ate o calculo da solucao do problema. O metodo dePicard, apesar de possuir um raio de convergencia maior, tipicamente gasta umnumero muito maior de iteracoes para calcular a solucao do problema e em muitasvezes nao converge, independente do valor do chute inicial.

Logo, o procedimento mais seguro e eciente numericamente e utilizar o metodode Newton acoplado com tecnicas para procura de chutes iniciais adequados. Esteprocesso de calculo de chutes iniciais a` medida que um parametro do problema evariado e denominado de Continuacao.

Um sistema de equacoes algebricas originado a partir da discretizacao de umaequacao diferencial que descreve um determinado fenomeno fsico pode ser repre-sentado de uma forma geral como:

f(x, ) = 0,

fi sao as equacoes a serem resolvidas, xi representam a solucao do problema, eum parametro do problema (numero de Reynolds, por exemplo).

Geralmente a solucao do problema para uma determinada faixa do parametro e bem simples. Se for o numero de Reynolds, o problema f(x, 0) = 0 e linear,e o metodo de Newton converge em uma iteracao, independente do valor do chuteinicial.

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 chute inicial:

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 dasolucao em relacao ao parametro , como descrito a seguir. Este procedimento edenominado Continuacao de Primeira Ordem.

O sistema a ser resolvido e representado por

f(x, ) = 0.

Diferenciando esta expressao em relacao a , obtem-se:

fxJ

dxd

+f

= 0 Jdxd

= f



dxd

= J1( f

) x = J1 f

.

Logo, o valor do chute inicial para = i+1 e obtido pela seguinte expressao:

x(o)(i+1) = x(i)J1 f

.

Observe que se o metodo de Newton tiver sido utilizado para obter a solucaodo problema nao linear para = i, o custo computacional do calculo do chuteinicial e bem pequeno, uma vez que a matriz Jacobiana ja foi invertida para ocalculo de x(i).

A Figura 5.3 ilustra o chute inicial obtido com os dois procedimentos de con-tinuacao descritos nesta secao.

Figura 5.3: Chute inicial obtido com continuacao de ordem zero e primeira ordem.


Captulo 6

Formulacao Integral daEquacao de Navier-Stokes

A equacao de Navier-Stokes para uidos incompressveis :{ 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 pressaop.

Figura 6.1: Domnio para problema de escoamento.

61


As condicoes de contorno para a equacao de Navier-Stokes podem ser em ve-locidade ou forca agindo no contorno do domnio. A condicao de contorno utilizadadepende da fsica do problema. As fronteiras de um escoamento sao, em geral:

Entrada de uido. Sada de uido. Paredes solidas. Interface entre dois uidos.

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 a`s condicoes de im-permeabilidade e nao deslizamento.

u n = V n impermeabilidade,u t = V t nao deslizamento.

Para x 2, vamos prescrever a forca supercial

n T tracao

= f em x 2. (6.3)

Esta e uma condicao de contorno natural. Prescreve a forca na fronteira. Aforca prescrita f depende da fsica do problema.

Para resolver este sistema de equacoes diferenciais parciais pelo metodo dosresduos ponderados, temos que multiplicar o resduo da aproximacao de cadaequacao por um funcao peso e forcar a integral ao longo de todo o domnio a sernula.

Metodo dos Resduos Ponderados

Como a equacao de consevacao da quantidade de movimento (6.1) e umaequacao vetorial, o resduo ponderado correspondente a esta equacao sera calcu-lado atraves do produto interno do resduo da aproximacao com um funcao pesovetorial W. A equacao da continuidade e uma equacao escalar e desta forma, afuncao peso utilizada w e uma funcao escalar, como as utilizadas ate o captuloanterior.



Rm =

{ u u T} W d = 0;

Rc =

{ u}w d = 0.

Resduo associado a` equacao de quantidade de movimento:

Rm =

{ u u T} W d =

=

(u u) Wd

( T) W d.

O termo T apresenta segunda derivadas da velocidade (variavel primitivado problema). Vamos integrar por partes para transferir a derivada para a funcaopeso, conforme ja 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 resduo da equacao de conservacao da quantidadede movimento escrito em termos de suas componentes. No caso particular de sis-tema cartesiano, cada termo do resduo ponderado e escrito como:

1) (u u) W = W1(uux + v

uy

)+ W2

(u vx + v

vy

).



2) T : W = W1x[p + 2u

x

]

Txx

+ W1y

[(

u

y+

v

x)]

Txy

+

+ W2x

[(

u

y+

v

x)]

Tyx=Txy

+ W2y

[p + 2v

y

].

Tyy

3) (n T) W = fx W1 + fy W2.

Cada componente da funcao peso vetorial pode ser escrita como combinacaolinear de funcoes base escalar i. Se cada componente W1 e W2 pertence a umespaco vetorial de dimensao n, a funcao peso vetorial W pertence a um espaco dedimensao 2n que e gerado pela seguinte base de funcoes:

W =[

W1W2

]

Espacos de dim=2n

[

10

],

[20

], . . . ,

[n0

]

n funcoes

,

[01

],

[02

], . . . ,

[0n

].

n funcoes

Denindo o espaco das funcoes peso desta forma, as componentes da equacaode conservacao podem ser desacopladas, pois para as primeiras n funcoes baseaparecem termos correspondentes a` primeira componente da velocidade, ja queW2 = 0 e para as ultimas n funcoes base aparecem temos correspondentes a` se-gunda componente da velocidade, ja que para estas funcoes base W1 = 0.

Usando o fato do uido Newtoniano, ou seja, T = pI + [u + uT ], as2n equacoes algebricas associadas a` equacao de conservacao da quantidade demovimento sao escritas como:

Rimx =

i

[uu

x+ v

u

y

]+

ix

[p + 2u

x

]+

+iy

[(

u

y+

v

x)]d

i fx d ; i = 1, ..., n.

Rimy =

i

[uv

x+ v

v

y

]+

ix

[(

u

y+

v

x)]+

+iy

[p + 2v

y

]d

i fy d ; i = 1, ..., n.



Resduo associado a` equacao da continuidade:

Rc =

{ u}w d = 0.

Onde w esta no espaco gerado pelas seguintes funcoes peso: {1, 2, ..., m}.As m equacoes 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 funcoesbase. Os diferentes termos dos resduos associados a equacao de conservacao dequantidade de movimento contem pressao ou derivadas da velocidade. Como estestermos devem ter a mesma precisao de discretizacao, as funcoes bases utilizadaspara pressao nao precisam ser da mesma ordem que as usadas para o campo develocidade. Quando os dois campos de variaveis pertencem a espacos diferentes,diz-se que a Formulacao e Mista (Mixed Finite Element).

As variaveis do problema sao escritas em termos de funcoes base como:

u =[

uv

]=[ n

i=1 Uiini=1 Vii

] 2 n incognitas;

p =mi=1

Pi i m incognitas.

Repare que a funcao peso ultilizada na equacao de continuidade e a mesmafuncao ultilizada para interpolar a pressao (i) e assim o numero de incognitas(2n + m) e o mesmo que o numero de equacoes.

Como a formulacao nao apresenta derivadas da pressao ou da funcao pesoi, o espaco das funcoes i pode ser menor que o das funcoes i ultilizadas parainterpolar a velocidade (m < N) e as funcoes is podem ser descontnuas. Asfuncoes is devem ser contnuas, pois tem que satisfazer

|i|2 < . Nao ex-iste essa restricao para is.

Nao e qualquer combinacao arbitraria de is e is que leva a formulacoes combom desempenho 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 aseguir. Para facilitar a interpretacao geometrica, vamos analisar um problemade elasticidade linear de materiais incompressveis. As equacoes que descrevem oproblema sao identicas ao do problema de Stokes (Re = 0), sendo u o deslocamento.

T = 0 ;



u = 0.Vamos considerar a situacao de dois elementos triangulares, como mostrado na

Figura 6 e funcoes base i linear e i constante.

Figura 6.2: Exemplo de formulacao instavel - malha truncada.

(e)

( u) icte

d = 0 (e)

u d = 0

A condicao de incompressibilidade com i constante leva a uma condicao dearea constante 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 oproduto base 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.

Para manter a area II xa o no A so pode deslocar-se na horizontal.

uA = 0 Malha TrancadaIsso se deve a uma ma escolha da funcao peso. A tabela a seguir lista diferentes

combinacoes 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 discontnuo (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 defuncoes para a velocidade e pressao e valida CONDICAO DE BABUSKA-BREZZI.

A condicao de Babuska-Brezzi verica a consistencia das aproximacoes dasderivadas. A interpolacao de u e p nao podem ser escolhidas arbitrariamente.

i; isup i i d{ |i|2 d} 12 c {

2i d}

12

A vericacao e prova se uma determinada combinacao de i e i satisfazema condicao de Babuska-Brezzi sao complicados matematicamente e estao fora doescopo do curso, devendo ser tratado em um curso mais avancado de elementosnitos.

A formulacao mista apresentada pode ser interpretada como solucao de umaequacao diferencial (quantidade de movimento) com uma restricao imposta (con-tinuidade) pelos metodos dos Multiplicadores de Lagrange.

Na apresentacao do metodo de elementos nitos para equacao de Navier-Stokes,vamos nos concentrar no seguinte elemento : Bi-quadratico para a velocidade elinear discontnuo para pressao. Neste elemento, cada componente da velocidadepossui 9 graus de liberdade e o campo de pressao e representado por 3 graus deliberdade.

9i=1 Uii9i=1 Vii

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 sao21 21.

As funcoes base elementares i utilizadas para expandir o campo de velocidadeserao lagrangeanas, isto e:



Figura 6.3: Elemento bi-quadratico de 9 nos.

i(Xj) = ij ,

cada funcao base e igual a 1 no associado a funcao e nula nos demais. Cada co-eciente Uj e Vj representa o valor da velocidade u e v no no j. Em termos dascoordenadas locais, as funcoes base para a velocidade sao:

1(, ) =( 1)( 1)

4

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)

29(, ) = (1 2)(1 2)

A gura 6.4 ilustra a forma destas funcoes base.

Figura 6.4: Funcoes base bi-quadraticas.

As funcoes base utilizadas para expandir o campo de pressao nao sao La-grangeanas. Para este elemento, elas sao escolhidas de forma que o primeiro graude liberdade de pressao P1 represente o valor da pressao no centro do elemento;



o segundo grau de liberdade P2, a derivada da pressao na direcao ; e o terceirograu de liberdade P3, a derivada da pressao na direcao . Para isto, a variacao dapressao 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; dpd = P2 e dpd = P3.

Substituindo as expansoes para o campo de velocidade e pressao nas equacoesdos resduos ponderados, obtem-se um sistema de 2n+m equacoes algebricas naolineares. 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 matriz Jacobiana que representa a sensibilidade decada equacao em relacao a cada incognita:

Jij = Ricj Derivada do Resduo i em relacao ao coeciente cj .

Ate este captulo, a numeracao dos graus de liberdade dos elementos se confun-dia com a propria numeracao dos nos. Como agora estamos resolvendo um sistemade equacoes diferenciais parciais para os campos de velocidade (campo vetorial)e pressao (escalar), o numero de graus de liberdade de cada elemento e maior doque o numero de nos do elemento. O elemento que estamos trabalhando possui 9nos e 21 graus de liberdade. Sao eles: U1, ..., U9;V1, ..., V9;P1, P2, P3.

A numeracao dos graus de liberdade elementar e totalmente arbitraria. Anumeracao local adotada aqui e mostrada na tabela a seguir:



# grau de liberdade local grau de liberdade1 U12 U23

4...

5678 U89 U910 V111 V212

13...

14151617 V818 V919 P120 P221 P3

Figura 6.5: Numeracao local dos graus de liberdade.



Para a montagem da matriz jacobiana elementar e interessante construir umvetor que aponte o numero do grau de liberdade elementar associado a cada no:

ivxlocaldof(1) = 1(2) = 2...(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, ux ,

uy , v,

vx ,

vy.

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) = ...

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 umasmalha simples 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 numeracaodos elementos.

FIGS*AssMtrx

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 SOLUCAODA EQUACAO DE NAVIER-STOKES

METODO DE NEWTON

{JC = RC(k+1) = C(k) +C



PROGRAMA PRINCIPAL

Inicializar - Ler chute inicial ou inicializar o valor dos coeficientes c

Calcular o vetor de Resduo globalCALL FormRV R

Se R < stop. O chute inicial era soluc~ao do problema. Repetir enquanto R >

- Calcular a matriz Jacobiana globalCALL FormJM J

- Resolver JC = RCALL SOLVER (J,R,C)

- Calcular C = C+C- Calcular o vetor Resduo global

CALL FormRV R Gravar a soluc~ao C

FormRV Calculo do vetor de Resduo global

do iele = 1, NELECALL getelemRV ( elemRV )do ildof = 1, NDOFELE

igdof = ASSMtrx (ildof, iele)RV (igdof) = RV (igdof) + elemRV (ildof)

end doend do

FormJM Calculo da matriz jacobiana global

do iele = 1, NELECALL getelemJM ( elemJM)do ildof = 1, NDOFELE

igdof = ASSMtrx (ildof, iele)do jldof = 1, NDOFELE

jgdof = ASSMtrx (jldof, iele)JM(igdof, jgdof) = JM(igdof, jgdof)+elemJM(ildof, jgdof)

end doend do

end do



O que precisaremos para o calculo do Resduo elementar:Formulas:

Rimx =

Re i

[uu

x+ v

u

y

]+

ix

[p + 2u

x

]+

+iy

[(

u

y+

v

x)]d

Rimy =

Re i

[uv

x+ v

v

y

]+

ix

[(

u

y+

v

x)]+

+iy

[p + 2v

y

]d

Ric =

(u

x+

u

x) i d

Lembrar que vamos usar um mapeamento isoparametrico para transformaras coordenadas x, y em , .

F (x, y) dxdy =

F (, ) |J| dd.

Lembrar que vamos usar integracao de Gauss para calcular a integral:

F (, ) dd =igp

jgp

F (igp, jgp) WigpWjgp.

Estrutura de Dados:

Denir os seguintes vetores / matrizes:

PHI(i)i = 1, ..., 9 PHI

12...

9

GradPHI(i, j)i = 1, 2

j = 1, ..., 9GradPHI

1d . . .

9d

1d . . .

9d

PSI(i)i = 1, 2, 3 PSI

12

3

XY (i, j)i = 1, 2

j = 1, ..., 9XY

[X1 . . . X9Y1 . . . Y9

]



velocity(i, j)i = 1, 2

j = 1, ..., 9velocity

[U1 . . . U9V1 . . . V9

] pressure(i)i = 1, 2, 3 pressure

P1P2P3

GradPHI XY [

1dx . . .

9dx

1dy . . .

9dy

]= J1

[1d . . .

9d

1d . . .

9d

]

JAC [

dxd

dyd

dxd

dyd

]=

[1d . . .

9d

1d . . .

9d

]

X1 Y1... ...

X9 Y9

UV =[

UV

]=[

U1 . . . U9V1 . . . V9

]

1...

9

P = [ P ] = [ P1 P2 P3 ] 12

3

GradUV XY [ du

dxdvdx

dudy

dvdy

]=

[1dx . . .

9dx

1dy . . .

9dy

]

U1 V1... ...

U9 V9

getelemRV Calculo do Resduo 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 Tdo ilnode = 1, 9

ivx = ivxlocaldof(ilnode)ivy = ivylocaldof(ilnode)elemRV (ivx) = elemRV (ivx)+WdetJ(RePHI(ilnode)

(UV (1)GradUV XY (1, 1)+UV (2)GradUV XY (1, 2))+GradPHI XY (1, ilnode)(P (1)+2GradUV XY (1, 1))+GradPHI XY (2, ilnode)(GradUV XY (1, 2)+GradUV XY (2, 1)))

end do

end do

end do

Desta forma, as primeiras 9 linhas da matriz Jacobiana elementar correspon-dem aos 9 resduos da componente x da equacao de Navier-Stokes; as 9 linhasseguintes, aos 9 resduos da componente y; e as 3 ultimas linhas, aos 3 resduos daequacao de continuidade. Da mesma forma, as 9 primeiras colunas estao associ-adas aos coecientes Uj ; as 9 seguintes, aos coecientes Vj ; e as 3 ultimas colunas,aos 3 coecientes Pj . Seguindo a numeracao dos graus de liberdade apresentadana tabela, a matriz Jacobiana elementar possui a seguinte estrutura:

RimxUj

RimxVj

RimxPj

RimyUj

RimyVj

RimyPj

RicUj

RicVj

RicPj

.

A expressao de cada cada bloco desta estrutura a apresentada a seguir:

Jacobiano de Rimx :



RimxUj

=e

{ i

[j

u

x+ u

jx

+ vjy

]+

+ix

2 jx

+ iy

jy

}de;

RimxVj

=e

{ i j

u

y+

iy

jx

}de;

RimxPj

=e

ix

i de.

Jacobiano de Rimy :

RimyUj

=e

{ i j

v

x+

ix

jy

}de;

RimyVj

=e

{ i

[ujx

+ jv

y+ v

jy

]+

+ ix

jx

+ 2 iy

jy

}de;

RimyPj

=e

iy

i de.

Jacobiano da Continuidade :

RicUj

=e

jx

i de;

RicVj

=e

jy

i de;

RicPj

= 0.

O mapeamento da matriz elementar para a matriz global e feito de maneirasemelhante aos problemas anteriores. A diferenca e que agora cada no possui mais



de 1 grau de liberdade. As linhas e colunas da matriz elementar estao associadasao grau de liberdade e nao aos nos. Nos exemplos dos captulos anteriores, asduas numeracoes eram coincidentes. O mapeamento deixa de ser por nos (usandoa matriz DomNodeID) e passa a ser por graus de 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 globalca com a seguinte forma:

RimxUj

RimxVj

RimxPj

RimyUj

RimyVj

RimyPj

RicUj

RicVj

RicPj

= 0

=

K C

CT 0

.

O bloco de zeros na diagonal da matriz pode representar um problema a maisna solucao do sistema linear J C = R em cada passo do metodo de Newton.

Existem outras formulacoes alem da formulacao mista discutida anteriormenteutilizadas para resolver o problema de escoamento incompressvel.

Metodo da Penalidade

A ideia basica e escrever uma equacao que relacione a pressao com a velocidade.Se esta equacao de estado for quase incompressvel, entao a solucao obtida pelometodo da penalidade vai ser proxima da solucao do problema incompressvel.

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

Pode-se substituir a pressao na equacao de Navier-Stokes p = 1 u :

u u = [ 1 u + (u+uT )] .80 Marcio Carvalho e Juliana Valerio


A vantagem do Metodo da Penalidade e que a pressao foi eliminada e nao emais variavel do problema. O numero de variaveis do problema foi reduzido.

A proximidade da solucao obtida com a solucao do problema incompressveldepende do valor de . Deve-se utilizar um pequeno (o menor possvel). Porema medida que diminui (e consequentemente 1 aumenta), a contribuicao relativados termos viscosos da tensao diminui e a solucao obtida nao satisfaz corretamentea 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 ob-jetivo e evitar o problema descrito anteriormente . Como a formulacao e mista,resolve-se o campo de velocidade e pressao, porem a continuidade nao e satisfeitatotalmente: {

u u = [p + (u+uT )]; u = p.

Rm =

{ u u = [p + (u+uT )]} w d = 0;

Rc =

{ u p}m d = 0;

u =9

i=1

Uii ; v =9

i=1

Vii ; p =mi=1

Pi i.

Os resduos da equacao de quantidade de movimento sao identicos aos da for-mulcao mista. A unica diferenca e na equacao de continuidade.

Ric =

[u

x+

v

y p

]i d;

RicUj

=

jx

i d;

RicVj

=

jy

i d;

RicPj

=

i j d.

A matriz Jacobiana global ca 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 diagonale substituido por M.

Como M possui inversa, pode-se desacoplar a velocidade e pressao:

K C

CT M

U

P

=

Rm

Rc

.

K UC P = Rm

CT U+ M P = Rc P = 1 M1[RcCT U].

Substituindo-seP em funcao deU, pode-se chegar a um sistema de equacoessomente para o passo de velocidade:

K UC{ 1 M1(Rc CT U)} = Rm[K+ 1 C M

1 CT] U = Rm 1 C M1 Rc

O Problema desacoplado pode ser reescrito como:

K+ 1 C M1 CT 0

CT M

U

P

=

Rm 1 C M1 Rc

Rc

.



Calcula-se primeiramente a velocidade e depois a pressao:

p =1 u.


Captulo 7

Condicoes de Contorno

As condicoes de contorno de um problema em mecanica dos uidos podem ser de2 tipos:

Condicao de contorno em velocidade (Dirichlet) : Especica o valor da ve-locidade ou de um de seus componentes na fronteira do domnio.

Condicao de contorno em tracao (Newman) : Especica a forca ou um deseus componentes na fronteira do domnio.

As diversas situacoes fsicas 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 domnio considera-se um perl de velocidadeconhecido.

Por exemplo, perl parabolico:{

u = u(y),v = 0.

iii) Condicao de simetria

{Ft = 0 forca tangencial = 0,v = 0 velocidade vertical = 0.

85


y

x

Figura 7.1: perl parabolico.

y

x

Linha de simetria

Figura 7.2: linha de simetria.{Ft = 0 t (T n) = 0,v n = 0.

iv) Condicao de tracao nula.

T n = 0 {

Ft = t (T n) = 0,Fn = n (T n) = 0.

v) Condicao em interfaces.

n T =tensao supercial

dtds

+ npressao atm

pa

n (T n) =curvatura

K + pa,t (T n) = 0.



x

y

n

t

Figura 7.3: linha de simetria.

Figura 7.4: tracao nula.

As condicoes de contorno de Dirichlet sao geralmente impostas de umamaneira essencial. O resduo correspondente a` variavel e simplesmente substituidopor:

Ri = Ui U,onde U e o valor que deseja-se impor a Ui.

Por exemplo, vamos analisar o problema ja visto anteriormente:A condicao de contorno e:

velocidade horizontal u = 1 u 1 = 0,velocidade vertical v = 0.

Como o grau de liberdade local correspondente a` velocidade horizontal no nomarcado e #7 e o correspondente a` velocidade vertical e #16, vamos substituirestes respectivos resduos por:

R7 = U7 1,



liquido

ar

n

t

Figura 7.5: interfaces.

V = 1

Figura 7.6: no 7.

R16 = V7 0.Sendo Jij = RiUj , note que:

J7j = 0, para j = 7 e J77 = R7U7

= 1,

J16j = 0, para j = 16 e J16 16 = R16V7

= 1.

A procedimento usual e calcular os resduos e as jacobianas elementares in-dependentemente das condicoes de contorno. Se o elemento pertencer a` fronteirado domnio, impoe-se a condicao de contorno. Se a condicao for de Dirichlet, oresduo associado ao grau de liberdade e substituido por: Ri = Ui U e a linhada matriz jacobiana associada ao grau de liberdade e zerada a menos do elementoda diagonal, que e substituido por 1.

No caso do exemplo anterior, onde a condicao de contorno e de Dirichlet, parao elemento 2, a matriz jacobiana ca:



. . . . . . . . . 0 0 1 0 . . . 00 0 0 1 . . . 0...

. . ....

0 . . . 1

7 34

8 9 6

15

2

Figura 7.7: elemento 2.

As condicoes de Newman aparecem naturalmente na formulacao do problema.Lembrando a formulacao de resduos ponderados da equacao de Navier-Stokes,

tem-se:

Rm =

{Re u u T} W d =

=

Re (u u) W +T : W

n (T W) d.

Rimx =

Re i

[uu

x+ v

u

y

]+

ix

[p + 2u

x

]+

+iy

[(u

y+

v

x)]d

i (n T)x d.



Rimy =

Re i

[uv

x+ v

v

y

]+

ix

[(u

y+

v

x)]+

+iy

[p + 2v

y

]d

i (n T)y d.

Observe que a integral

i (n T)x/y d so sera calculada para os elementosda fronteira do domnio onde a condicao de contorno de Newman e imposta, poisnos elementos internos nao temos que nos preocupar com a integral ao longo dafronteira do domnio e nos elementos da fronteira onde a condicao de Dirichlete imposta o resduo ponderado e substituido da forma mostrada anteriormente.

Assim, quando a condicao e de Newman e o elemento esta na fronteira, o valordo integrando (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.

OBSERVACAO : A forma da equacao de Navier-Stokes a qual se aplica aformulacao fraca 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 simplicado usando que u = 0, pelaequacao da continuidade;

= [pI+u+uT]= p + (u)

2u

+ (uT) (u)

= p +2u+( u =0

)

= p +2u

Esta forma simplicada e muito utilizada na literatura de diferencas nitas,porem o termo que aparece na integral de linha nao possui signicado fsico, por



isso nos problemas onde as condicoes de contorno de Newman sao forcas atuandona fronteira do domnio nao e interessante usar a forma simplicada.

CONDICAO DE CONTORNO DE SAIDA DE ESCOAMENTO

Pode ocorrer condicoes de contorno na sada do escoamento que nao represen-tem uma forca, apenas sao a derivada da velocidade; como no caso do escoamentodesenvolvido, que a condicao de contorno pode ser escrita como n u = 0.

Em coordenadas cartesianas:

nxux + ny

uy = 0,

nxvx + ny

vy = 0.

***Fig de uma fronteira sintetica

A posicao da fronteira de sada do escoamento e arbitraria. Ela deve ser colo-cada de forma que o escoamento na regiao de interesse seja independente de suaposicao (uma maneira de saber se a fronteira esta o mais proximo possvel, porcausa de tempo computacional, mas ja com o escoamento desenvolvido, e avaliarum parametro especco).

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,

R =

( ) d

i [p n+ n u 0

+n uT] d.

Substitui-se a condicao de contorno, que e apenas um pedaco do argumento daintegral de 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 ux

+ nyv

x].



Rmy =

( ) d

i [p ny + nx uy

+ nyv

y].

Observe que os valores de p , ux ,uy

vx ,

vy na fronteira nao sao conheci-

dos e as expressoes em termos das funcoes base devem ser usadas no calculo daintegral de linha, isto e:

p =

Pjj ,ux =

jx Uj , ...

CALCULO DA INTEGRAL DE LINHA

i( ) d :

Considera-se o seguinte elemento, localizado na sada de um escoamento, ondea condicao de escoamento desenvolvido esta sendo imposta como condicao de con-torno.

**FIG

1) Calcular R sem se preocupar com a condicao de contorno.

2) Vericar se alguma condicao tem que ser imposta no elemento.

Sendo a condicao de contorno Natural e necessario acrescentar os termosi( ) d

nos graus de liberdade asso

Documents

Introducción Al Metodo de Los Elementos Finitos