O Método dos Elementos Finitos e os métodos “meshfree” na solução de Equações

Diferenciais Parciais

José Márcio MachadoUNESP - Campus de São José do Rio Preto

1) Ilustrando as idéias básicas do método com um exemplo simples:

Consideremos a equação de Laplace em domínios bidimensionais, que aparece comumente em vários tipos de problemas. Pede-se calcular o potencial u que é solução da e.d.p.

,02 u

em um domínio com uma fronteira C, dividida em subcontornos C1 e C2 onde u satisfaz respectivamente a condições de Neumann e de Dirichlet, e tais que

(1.1)

,21 CC.21 CCC

O princípio de mínima energia para o potencial requer que a energia do campo armazenada no domínio assuma um valor mínimo para a solução de (1.1). Em outras palavras, a solução procurada minimiza a funcional

.21 2

dSuuW (1.2)

Portanto, uma forma de resolver (1.1) é procurar uma expressão aproximada para W(u), supondo que u possa ser aproximado por uma expansão com funções simples adequadamenteescolhidas, e coeficientes indeterminados. A minimização da energia determinará então esses coeficientes, e portanto a aproximação para u. Esta “abordagem indireta” através deuma formulação integral equivalente e expansões de aproximação é, em essência, a idéia bá-sica seguida por todos os Métodos de Elementos Finitos e suas variantes: substituir o problema de encontrar diretamente a solução de uma dada equação diferencial parcial, pelo problema equivalente de encontrar uma solução minimizante para uma certa integral.

Ao construir essa expansão para o potencial u, inicialmente decompomos o domínio do problema em uma união de subdomínios topológicamente regulares (os elementos finitos),disjuntos dois a dois, não superpostos e não intersectantes (ou seja, construímos uma coberturado domínio), como ilustrado pela figura 1 abaixo. Seja então, um elemento individualdessa decomposição (fig. 2):

Fig. 1

Fig. 2No caso mais simples, uma aproximação U paraa solução u de (1.1) é dada no elemento da Fig. 2 simplesmente por

,, cybxayxU (1.3)com os coeficientes a, b e c satisfazendo

.111

33

22

11

3

2

1

cba

yxyxyx

UUU

(1.4)OBS.: Note na figura que os elementos po-dem possuir dimensões e orientações distin-tas, se adaptando às particularidades da geo-metria do domínio original.

Combinando (1.3) e (1.4), é imediato concluir que para um único elemento,

.111

1,

3

2

11

33

22

11

UUU

yxyxyx

yxyxU (1.5)

Com um pouco de Álgebra, a expansão (1.5) pode ser escrita na forma

3

1

,,i

ii yxUU (1.6)

onde, por exemplo, ,21

233223321 yxxxyyyxyxA

sendo A a área do elemento triangular. As demais funções alfa são fácilmente obtidas porpermutação cíclica de índices. Essas funções, chamadas na terminologia do método de fun-ções de base ou funções de forma, são interpolatórias nos vértices do triângulo, satisfazendoa propriedade do delta de Kronecker:

.,0,

,,1,

jiyx

jiyx

jji

jji

(1.7)

Substituindo agora a expansão (1.6) na equação (1.2), a energia associada a um único ele-mento é aproximada por

3

1

3

1

2 ..21

21

jjjii

i

e UdSUdSUW (1.8)

Se definirmos os elementos de matriz ,.

dSS jie

ij

a equação (1.8) será escrita então como uma forma matricial quadrática:

(1.9)

U,UT eS21

eW (1.10)

onde U representa o vetor coluna com os valores do potencial nos vértices (pontos nodais)do elemento triangular. Para determinar a aproximação que inclua todos os elementos deuma dada discretização (malha de elementos finitos), é suficiente impor a continuidadedo potencial em pontos nodais que são comuns a dois ou mais elementos.

Vamos considerar na figura ao lado, as numera-ções disjunta (a) e conjunta (b), que também são chamadas na terminologia corrente do método, respectivamente, de numerações local (a) e global (b).

Fig. 3

Para os dois elementos separados, o vetor das incógnitas se escreverá como

,654321 disUUUUUUTdisU (1.11)

e a energia total assocciada ao par de elementos será

,UU disTdis2

1disSW (1.12)

onde

.

000000000

000000000

266

265

264

256

255

254

246

245

244

133

132

131

123

122

121

113

112

111

SSSSSSSSS

SSSSSSSSS

disS (1.13)

A matriz S em (1.13) é conhecida como matriz de Dirichlet ou de rigidez. Se considerarmosagora a montagem conjunta com numeração global (elementos unidos), é necessário que ospotenciais variem contínuamente através das interfaces entre elementos. Em outras palavras,ao juntarmos dois elementos, os valores de potencial em vértices (pontos nodais) correspon-dentes devem ser iguais. Na figura 3.a, isto quer dizer que os potenciais em 1 e 6 são iguais,assim como os potenciais em 2 e 4.

Relacionando agora as numerações local e global por meio de uma matriz de conexões,

.,

000110000010010000100001

4

3

2

1

6

5

4

3

2

1

condis CUU

con

dis

UUUU

UUUUUU

(1.14)

Desta forma, substituindo (1.14) em (1.12), teremos

C,SS disTC S,UU con

Tcon2

1W (1.15)

.

00

255

254

256

133

132

131

245

123

244

122

246

121

265

113

264

112

266

111

SSSSSS

SSSSSSSSSSSS

S (1.16)

A eq. (1.16) define a matriz global para o problema, fornecendo juntamente com U uma aproximação para a energia do domínio constituído pelos dois elementos triangulares conectados.

O passo que falta é a minimização da energia W:

,0

p

l

pppl

lpllTp

Tl

l UU

SSSS

UUU kkU

W (1.17)

onde os subacritos l e p representam, respectivamente, os potenciais livres para variar (as in-cógnitas) e os valores prescritos (condições de contorno do tipo Dirichlet). Uma vez que a di-ferenciação óbviamente só é possível com respeito aos potenciais livres, temos

,0

p

llpll U

USS do que resulta por fim, .

p

plp-1ll

UUSS-

U (1.18)

Note que nesta formulação, as condições de Dirichlet são satisfeitas exatamente (são os valores prescritos de potencial). Pode-se mostrar que as condições de Neumann homogêneas(derivada normal nula) são implícitamente satisfeitas para este caso, não de forma precisa, mas com um valor aproximado que depende essencialmente da natureza da aproximação po-linomial que fizemos para U. Neste exemplo ela é linear nas coordenadas, mas poderia ser quadrática, cúbica, etc. Uma demonstração deste fato será feita na seção 2. Os pontos essenciais a mencionar agora são os seguintes:

1) Toda a técnica apresentada acima pode ser estendida sem dificuldade para qualquer número de elementos, aumentando apenas a dimensão dos vetores e matrizes. A formação da matriz global também pode ser feita diretamente, sem matrizes de conexão;

2) A matriz global em (1.16) é esparsa e simétrica. O padrão de distribuição dos elementos nulos depende da numeração global adotada, que pode ser inteiramente arbitrária;

3) O grau do polinômio em (1.3) pode ser aumentado para se obter uma aproximação tão boa quanto se queira. Na aplicação do MEF, o que se busca é um compromisso: soluções precisas com aproximações lineares requerem um grande número de elementos, mas aproximações quadráticas ou cúbicas aumentam a dimensão das matrizes de contribuição local (1.9). Aplicações correntes raramente empregam aproximações cúbicas ou de grau maior.

No exemplo anterior (P. P. Silvester & R. L. Ferrari, Finite Elements for Electrical Engineers, Cambridge University Press) estão presentes portanto os seguintes elementos característicos do MEF:A) Uma funcional, oriunda de princípios variacionais ou da integral de um resíduo ponderado(métodos de Galerkin, Bubnov-Galerkin, etc.) cuja solução minimizante é a solução da e.d.p.correspondente,B) Uma discretização do domínio segundo uma certa malha de elementos finitos;C) Um conjunto de funções de base, interpolantes nos elementos da discretização, satisfazendoas propriedades (1.7).

O ítem A) será discutido mais detalhadamente agora. É importante adiantar que, mesmo para problemas em que nào é possível encontrar uma formulação variacional correspondente, o método de Galèrkin e suas variantes permite a construção da integralprocurada.

2) Princípios Variacionais e o Método de Galerkin :Usando o Cálculo das Variações, é possível mostrar a validade do seguinte resultado:

yxfyxyxkyxk ,,,,,,, 21 Sejam funções dadas em uma região G bidimen- sional, e ss , funções definidas do segmento de arco s ao longo da fronteira C de G.

A função u(x,y) que torna a expressão integral

G

dydxuyxfuyxyuyxk

xuyxkI ,,

21,,

21 2

2

2

2

1

dsususC

2

21

estacionária sob a condição associada 1Csu em onde C1 é um segmento de C,será necessáriamente uma solução do problema de valores de contorno

yxfuyxyuyxk

yxuyxk

x,,,, 21

em G,

sob as condições de contorno de Dirichlet

e geral de Cauchy

su em C1

susnyukn

xuk yx

21 em C2, sendo

nx e ny os cosenos diretores da normal exterior relativamente ao contorno C, e sendo C2 o restodo contorno, de tal forma que ,21 CCC .21 CCClaramente temos uma vantagem aqui: a condição do tipo Cauchy não está presente naformulação integral de modo explícito. Apenas a condição de Dirhclet precisa ser levada em consideração, se formulamos o problema através de um princípio de extremo. (H. R. Schwarz, Finite Element Methods, Academic Press)

,0,02

2

2

2

fpara

yu

xu (equação de Laplace)

,0,,2

2

2

2

parayxf

yu

xu

(equação de Poisson)

.,0,02

2

2

2

fparau

yu

xu

(equação de Helmholtz)

Existem contudo, muitos problemas para os quais não é possível encontrar funcionais derivadas pelo Cálculo das Variações. Nestes casos, uma técnica bem mais abrangente para se encontrar a integral procurada é o chamado método de Galerkin, uma variante do Método dos Resíduos Ponderados. O método de Galerkin irá produzir matrizes globais formalmente idênticas às encontradas, por exemplo, pelo teorema anterior em problemas que admitem uma formulação variacional.

O método de Galerkin: Como já vimos, para um único elemento finito escrevemos a aproxima-ção

p

i

ei

ei

e yxUU1

,, onde supomos agora que a aproximação pode ser de ordem mais elevada, com p maior que 3.

Ora, a forma global da aproximação U em todo o domínio discretizado é uma combinação das a-proximações U(e)(x,y) sobre todos os elementos, ou seja, é a união das aproximações (2.1) sobretodos os elementos. Numerando todos os pontos nodais do domínio discretizado de 1 a n, a apro-ximação global U será dada então por

(2.1)

(2.2) .,1

n

kkk yxUU

Em (2.2), yxk , é a combinação (união) das funções de base individuais yxei ,

que tem o valor 1 no ponto nodal Pk, que é associado com a variável nodal Uk. Claramente afunção de base global índice k tem valores não nulos apenas em elementos que contém o ponto nodal Pk, um subdomínio muito restrito, seu suporte local. A expressão (2.2) é umaaproximação de Ritz, com funções que possuem apenas suporte local. Contribuições integrais globais serão dadas naturalmente pela soma das contribuições individuais, por exemplo,

yxk , são as funções de base globais.

. V Ve

dVdV Para exemplificar o método, consideremos como exemplo o problema de condução de calor em regime estacionário, dado pelas equações

,,,0

y

x

y

x

yx

qqem

yq

xq

Qkqy

kqx

uuA

da chamada forma mista, e em termos de operadores vetoriais, como

.0

,0

qn emqq

em

uB

(2.3)

(2.4)

Usando (2.2) para aproximar cada um dos componentes do vetor u em (2.3), podemos escrever então que

,,......,1,0 njdd

r

iii

Tj

r

iii

Tj uαBwuαAw

que é a integral da ponderação do resíduo tanto da equação (2.3) quanto de suas condições de contorno (2.4), resíduos estes obtidos ao substituirmos o valor exato da solução pela aproximação (2.2). Quando escolhemos as funções de ponderação wj como sendo as funções de base globais da aproximação (2.2), temos o método de Galerkin, ou de Bubnov-Galerkin.

(2.5)

Como podemos ver na equação (1.18), problemas de campo estacionário se reduzem em gerala uma equação matricial do tipo 0dSu sendo a matriz S bem esparsa.

Exemplo: discretização de um domínio, e matriz de Dirichlet global resultante:

Matriz de Dirichlet global para funções de base quadráticas nos elementos triangulares.

Numeração global arbitrária.Fig. 2.1

Fig. 2.2

OBS.: As regiões em branco na matriz são ocupadas por zeros.

Na técnica EFGM, aproximantes MLS são usadas para construir uma aproximação uh ( x ) para a solução exata ( x ) de uma dada equação diferencial parcial. Para tanto, um conjunto discreto de N pontos nodais é tomado sobre o domínio do problema, e uma função w de suporte compacto é associada a cada ponto nodal ou nó. Para cada nó existe um subdomínio, chamado domínio de influência, no qual o ponto nodal respectivo contribui para a solução. Este subdomínio consiste simplesmente da região em torno do ponto em que a função w (função pêso) é diferente de zero. Ao contrário do MEF, em que os elementos devem ser disjuntos dois a dois, domínios de influência na técnica EFGM podem se superpor. Para simplificar, consideremos um problema no plano e domínios de influência circulares ou retangulares, que correspondem respectivamente às duas funções pêso mais comumente encontradas, dadas por exponenciais gaussianas ou produtos tensoriais:

3) O Método EFGM (Element Free Galerkin Method) :

O domínio geométrico do problema é o contorno omega, e os círculos e retângulos superpostos são os domínios de influência. A solução procurada é avaliada nos pontos indicados, que continuaremos chamando de nodais mas que agora não fazem mais parte de uma malha de elementos conectados, mas são centrais a regiões que podem ser superpostas.

Uma aproximação local uhl (x) é construída como um polinomio de ordem m

com coeficientes não constantes, como segue:

.0

xxxxxx, apT

m

jjj

hl apu

A escolha mais simples é uma base linear pT, em uma dimensão:

. 1xT T

p p x

Para calcular os coeficientes a ( x ) e resolver a aproximação acima, uma norma discreta ponderada L2 (denotada por J) é construída de tal forma que

,1

2

1

2

n

llll

n

ll

hll

uw

uuwJ

xxxx

x,xxx l

apT

onde a soma é extendida a todos os nós para os quais a função pêso w é diferente de zero, e ul é um parâmetro nodal que controla a função u(x) no ponto xl . Note que para a técnica EFGM, os parâmetros nodais ul são diferentes dos valores uh(xl) da aproximação.

(3.1)

(3.2)

Calculando as derivadas parciais de J com respeito às componentes de a(x), resulta do critério de minimização que

,uBAa 1 xxx

onde as matrizes A e B e o vetor u, para uma base polinomial linear, são dados por

,1

...

...1

2

211

11

1

nn

nn

n

llll

xxx

xxw

xxx

xxw

w xxxx TppxA

,,11 nn xxxwxxxw ppB ... x

....u 21T

nuuu

(3.3)

Por substituição de (3.3) em (3.1), resulta

,1

uxxx

n

lllh uu

.1

l

m

jljjl p

BAp

xxxx

1T

BA 1

onde definimos

(3.4)

A abordagem EFGM consiste portanto em utilizar uma expansão como (3.4), e seus elementos constituintes, para a análise de uma dada EDP. Como um primeiro exemplo, consideremos a equação escalar de Helmholtz homogênea

,02 ukuem um domínio bidimensional que descreve a secção reta de uma guia de onda retangular. Na análise modal da equação acima, uma expansão como (3.4) é usada na funcional correspondente, que é modificada quando condições de Dirichlet são impostas.

Por exemplo, as soluções para modos TM na guia devem satisfazer condições de Dirichlet homogêneas no contorno do domínio. As soluções procuradas serão expressas por aproximantes EFGM, de modo que

u.xxx )()()(1

n

iii

h u

A funcional modificada incluirá multiplicadores de Lagrange para imposição das condições de contorno:

.

2

21

2

du

duu

duuF

iii

i jjiji

jii j

ji

As condições de contorno homogêneas requerem que seja zero na fronteira, de modo que o problema de auto-estados correspondente para a funcional F é dado por

(3.5)

,2

λλu

000Tu

0GGK

T (3.6)

Nas integrais acima, s é o comprimento de arco ao longo do contorno e as funções N(s) são interpolantes de Lagrange. Dependendo da geometria do problema, uma soma discreta também pode ser usada na imposição das CC:

onde

,

dK jiji

,

dT jiji

.

dsNG ikik

.

11

inniin

iii

u

uS

x

x1

...

...

Neste caso, a expressão para G é simplesmente

nnnn

n

n

xxx

xxxxxx

G

21

22212

12111

A equação (3.6) pode ser então resolvida para determinar o vetor u.

Para o caso de modos TE, condições de Dirichlet não são dadas explícitamente, o uso de multiplicadores de Lagrange é desnecessário e temos então que

.2 TuKu Conhecidos os valores de u, usamos novamente a expansão

uxxx )()()(1

n

iii

h u

para determinar os valores do campo escalar.

(3.7)

..

d

dk

hh

hh

Observe que (3.7) fornece os auto-estados físicos do problema, as frequências correspondentes serão

Note que o uso dos multiplicadores de Lagrange se deve ao fato de que as funções de forma vistas até aqui não são interpolantes, como no MEF, não satisfazendo a condição do delta de Kronecker.

A inconveniência disto é que os multiplicadores que aparecem em (5), embora essenciais para assegurar a satisfação das condições de contorno, não se relacionam diretamente com a solução física do problema. Isto representa um dispêndio extra de memória e processamento. Mesmo assim, resultados razoáveis podem ser conseguidos com esta abordagem. Como exemplo de desempenho, consideremos uma guia de onda retangular, com 10x20 mm, nós regularmente espaçados e integração por quadratura gaussiana 4x4 em cada célula. As soluções foram obtidas empregando-se funções de base bidimensionais lineares, tais como

.,1 yxyx xxpp TT

Como funções pêso, empregaram-se produtos tensoriais de splines cúbicas. Os gráficos seguintes ilustram as isolinhas de campo obtidas com o auxílio da técnica EFGM, e os êrros relativos no cálculos de alguns modos TE e TM (comparação com as soluções analíticas), feitos pela aplicação do MEF e da técnica EFGM, respectivamente. As figuras ilustrando as isolinhas demonstram a plena satisfação das condições de contorno.

1 2 3 4 5 6 7 8

2

4

6

8

10

12

14

Nod

e In

dex

J

mode TE 21

1 2 3 4 5 6 7 8

2

4

6

8

10

12

14

Node Index I

Nod

e In

dex

J

mode TE 32

Isolinhas para dois modos TE ..

0 500 1000 1500 2000degrees of freedom

0.00

0.04

0.08

0.12

relative error (%)

TE01 Mode

EFGM

first order FEM

second order FEM

( a )

0 500 1000 1500 2000degrees of freedom

0.00

0.40

0.80

1.20

relative error (%)

TE23 Mode

EFGM

first order FEM

second order FEM

( b )

0 500 1000 1500 2000degrees of freedom

0.00

0.50

1.00

1.50

2.00

relative error (%)

TE25 Mode

EFGM

first order FEM

second order FEM

( c )

0 500 1000 1500 2000degrees of freedom

0.00

1.00

2.00

3.00

4.00

5.00

6.00

relative error (%)

TE48 Mode

EFGM

first order FEM

second order FEM

( d )

Fig.1 Erros relativos nas frequencias, obtidos para uma guia retangular por EFGM e FEM em primeira e segunda ordens: (a) modo TE01, (b) modo TE23, (c) modo TE25, e (d) modo TE48 . 

0 500 1000 1500 2000 2500 3000degrees of freedom

0.0

0.5

1.0

1.5

2.0

2.5

Relative Error (%)

TM11 Mode

EFGM

first order FEM

( a )

0 500 1000 1500 2000 2500 3000degrees of freedom

0.0

1.0

2.0

3.0

4.0

5.0

Relative Error (%)

TM23 Mode

EFGM

first order FEM

( b )

0 500 1000 1500 2000 2500 3000degrees of freedom

0.0

1.0

2.0

3.0

4.0

5.0

6.0

Relative Error (%)

TM25 Mode

EFGM

first order FEM

( c )

0 500 1000 1500 2000 2500 3000degrees of freedom

0.0

1.0

2.0

3.0

4.0

5.0

6.0

7.0

8.0

Relative Error (%)

TM55 Mode

EFGM

first order FEM

( d )

Fig.2. Erros relativos nas frequencias, obtidos para uma guia retangular por EFGM e FEM em primeira ordem, como função dos graus de liberdade: (a) modo TM11, (b) modo TM23, (c) modo TM25, e (d) modo TM55.

Outro exemplo: a equação de Schrödingerunidimensional

Consideremos um domínio unidimensional dado por

,bxa e no qual um poço ou barreira de potencial quântico é descrito por uma função de potencial arbitrária V(x). Na chamada aproximação de massa efetiva, a função de onda para um elétron no domínio é dada pela equação de Schrödinger unidimensional

,0)())(()(2

22

xExVxmeff

onde é a função de onda do elétron, meff é a sua massa efetiva, E é o auto-valor da energia e

.2/ hComo antes, a funcional modificada por multiplicadoresde Lagrange se escreve como

(3.8)

b

a

ji

eff

oji

ji

dxdx

ddx

dmxuuF

2

2

b

ajiji

ji

dxxVuu

i j

b

ajiji dxuuE*

,

ibiib

iaiia buau

(3.9)Novamente, o problema de auto-valores se escreve como

,u

000Tu

0GGK

T

*E

já que a função de onda (estados confinados) deve se anularnos extremos do domínio.

b

aji

b

a

ji

eff

oji dxxVdx

dxd

dxd

mxK ,

2

2

b

ajiji dxT ,

bbbaaa

Gn

nT

......

......

21

21

Casos de teste:

1) Poço de potencial com GaAs/AlxGa1-xAs:

me= (0.067 + 0.083x)mo, x = 0.3

0 50 100 150 200 250 300 350W idth (Angstrons)

0

50

100

150

200

250

300

Energ

y (Me

V)

N = 1

N = 2

N = 3

N = 4

N = 5

N = 6

Resultados para U = 0.225 eV:

Auto-valores de energia para diferentes modos ligados, em função da largura do poço:

2) Poço quântico de GaAs/AlxGa1-xAs, com campo elétrico externo:

me= (0.0665 + 0.0835x)mo, x = 0.32,

d = 95 Angstrons,

db= 98 Angstrons,

U = 0.34 eV, com campos variando de 0.0 a 12x10**4 V/cm:

0 5 10 15Electric Field (10 V/cm)

-40

-20

0

20

40En

ergy

(meV

)

Electron

Electron

Heavy Hole

4

Estados fundamentais para elétrons e lacunas (N = 1) relatiamente aos potenciais no centro do poço quântico, como função do campo elétrico aplicado:

Análise de erros numéricos:

N

NN

e

ceN k

kk

Nek eNck

0 200 400 600nodes

0.0E+0

2.5E-3

5.0E-3

7.5E-3

1.0E-2

N = 1

N = 2

N = 3

N = 4

Erros dos auto-valores de energia para os modos N = 1, 2, 3, e 4 em função do número de pontos nodais.

Os resultados obtidos para estados ligados em estruturas com poços de potencial retangulares, e para estados quase-ligados em poços submetidos a campos elétricos externos apresentam boa concordância com dados obtidos na Literatura através de outros métodos. Isto justifica a continuação da pesquisa para simular dispositivos eletro-ópticos e componentes semicondutores pela técnica EFGM.. Possíveis melhorias incluem a utilização de aproximantes MLS interpolantes, que possibilitam a eliminação de multiplicadores de Lagrange nas funcionais das EDPs.

Conclusão:

Aproximantes MLS interpolantes:

Como já mencionamos, em geral aproximantes MLS do tipo (3.4) não interpolam dados, uma vez que a propriedade do delta Kronecker

jiiN jx

não é satisfeita para o conjunto de pontos nodais.

Contudo, já foi demonstrado que esta propriedade pode ser satisfeita tomando-se funções w singulares na definição das funções aproximantes (3.4). A construção dessas funções pêso singulares pode ser feita da seguinte forma:

I) Seleciona-se um conjunto original de funções pêso w tais que

ijiw jx (3.10)

i

ii w

ww

1

II) Constroi-se um novo conjunto de funções pêso singulares tomando-se o conjunto da etapa I), aplicando-se em seguida a definição

Claramente, isto resulta em

, ii xxx iw .jiiN jxpodendo-se mostrar que então,

Existe um processo recursivo para geração das funções de forma neste caso, conhecido na Literatura como abordagem de consistência (consistency approach).

(3.11)

Segundo os cálculos da abordagem de consistência, as funções de forma serão dadas agora por

Com as funções pêso modificadas, a expressão acima se escreve como

,,,,,,,,,

,,,,

,,,,

222,

2

2,

2

2,

2

rqDrqDrqpkjDikjk

kjDkjij

i

kjDikjk

kjDkjij

i

kjDikjk

kjDkjij

ii

wwwwww

www

wwd

wN

xxxxxxxxxxxx

xxxxxx

xxxxxx

ppii

i

i

.,,

,,,

,,,,,

2

22

,1,12,1

2,

2

kjkijkjiikijkjD

kjDkjnjk

inijni

kjDikjk

kjDkjij

ii

yxyxyxyxyxyx

wwwd

wwdw

N

xxx

xxx

xxxxxx

i

i

i

(3.12)

onde p, q e r são sempre diferentes de i. De acordo com as propriedades das funções pêso singulares,  0,,,, 22 rqDrqDrqp www xxxxxx pp ixx sempre que

e portanto .1

,,,,

,,,,

2,

2

2,

2

kjDikjk

kjDkjij

i

kjDikjk

kjDkjij

i

iiwww

www

Nxxxxxx

xxxxxxxx

ii

ii

Contudo, funções pêso singulares podem resultar em um mau condicionamento do algoritmo. Esse problema é comumente superado na implementação computacional tomando-se aproximantes MLS com funções pêso de Shepard normalizadas. Neste caso, define-se na vizinhança de um nó uma nova classe de funções pêso, Wi, dada por

 

.

11

1,

1

i i

i

j

j

j

j

jj

ww

ww

ww

W xx

Que tipo de funções podemos tomar e que satisfazem a condição (3.10)? Há um grande número de escolhas possíveis para funções pêso mencionadas na Literatura, mas as mais comuns em geral são splines e funções de Schwarz. Se os domínios de influência são retangulares por exemplo, podemos tomar produtos tensoriais de funções de Schwarz em que, para cada coordenada, a variável r é a norma da diferença entre o valor da coordenada para o ponto de observação e para o j-ésimo ponto nodal:

,1,0,10,12

2

rrerw r

r

(3.13)

e uma vez que ,, 21 jj yyrxxr

fica claro que o produto tensorial de funções como (3.13) (no plano, uma função para x e outra para y) satisfaz a condição do delta de Kronecker.

A Fig.1 mostra um problema de teste (C-type magnet) resolvido pelo MEF e pela técnica EFGM com funções de forma interpolantes. Os parâmetros são

3) Um teste da formulação interpolante:

Fe

Cu

Cu

Cuz yxJ

yxg,10

,,,

,0,

,0

3

0

Jz é a componente z do vetor densidade de corrente, e é a relutividade magnética do vácuo. Os pontos do domínio usados para a análise dos resultados obtidos com os dois métodos são mostrados na Fig. 2.

0

Fig. 1 Fig. 2

As equações de Maxwell para campos magnetostáticos são reduzidas neste caso a uma única equação de Poisson para o potential vetor magnético, que devido à simetria do problema só tem a componente z:

yxgyuyx

yxuyx

x,,,

. (3.14)

Uma vez que estamos usando aproximações EFGM interpolantes, multiplicadores de Lagrange são desnecessários para a imposição das condições de contorno. Obtemos agora o seguinte sistema de equações lineares:

,...,...,,, 21T

nuuu ufKu

.,

,,

dxdyyxgf

dxdyyxK

jj

jiij

(3.15)

O domínio foi discretizado em 3773 pontos nodais, e as integrações numéricas executadas com o uso de quadratura Gaussiana de quarta ordem com 1856 células. Resultados de ensaios preliminares em problemas de eletrostática mostraram que uma discretização menos refinada mas com uma ordem de integração mais elevada é preferível, ao invés de um refinamento maior na discretização mas com ordens de integração mais baixas. A figura 3 abaixo mostra as isolinhas do potencial obtidas pela técnica EFGM para o problema em questão.

Fig. 3

Fig. 4: Erros absolutos (tomando a solução FEM como referência) para os pontos de amostragem da Fig. 2. A aproximação FEM foi feita em 3779 pontos nodais e a EFGM, em 3773. Valores de erro absoluto em Wb/m2.

Notar a razoável concordância entre os dois métodos. Como esperado, o erro é maior em regiões onde o fluxo de campo é mais intenso (núcleo de ferro, perto dos enrolamentos).

Fig. 5: Erros absolutos (tomando a solução FEM como referência) para os pontos de amostragem da Fig. 2. A aproximação FEM foi feita em 7993 pontos nodais e a EFGM, em 3773. Valores de erro absoluto em Wb/m2.

Uma comparação dos erros relativos por sua vez, permite uma compreensão melhor do desempenho desses métodos:

Fig. 6: Erros relativos para os pontos da Fig. 2. A aproximação FEM foi feita em 3779 pontos nodais e a EFGM, em 3773.

Um número de conectividade em torno de 6-11 mostrou-se uma boa escolha neste caso, mas a técnica de truncamento do domínio, devido a descontinuidade do material, pode requerer outras escolhas em outros casos. Note que os erros relativos são menores (melhores) para a solução FEM com melhor discretização (Fig. 7, adiante).

Fig. 7:Erros relativos para os pontos da Fig. 2. A aproximação FEM foi feita em 7993 pontos nodais e a EFGM, em 3773.

A Fig. mostra a descontinuidade da norma do campo magnético, H, calculada sobre a linha de avaliação da Fig. 1, que concorda razoávelmente bem tanto qualitativa quanto quantitativamente com o resultado obtido via FEM.

Fig. 8: Descontinuidade na norma do campo segundo o método EFGM.

Outra aplicação: um problema de Magnetohidrodinâmica

,0

,2

V

EBJVpVV e

,,0o

eEE

,VBVEJ e

,0, BJB o

Fluxos magnetohidrodinâmicos em condição estacionária são normalmente descritos por um sistema de equações diferenciais parciais, composto pelas equações de Maxwell que descrevem o campo eletromagnético, e pelas equações de Navier-Stokes que descrevem o fluído condutor. Ou seja,

 

,,ayY

axX

,2

zpa

VV z

,2

o

z

zpa

BB

,aBM o

Após algumas manipulações algébricas, usando os índices 1 e 2 para denotar regiões de fluído e paredes do ducto respectivamente, e incluindo as variáveis adimensionais acima, o problema MHD será descrito por um conjunto de EDPs acopladas. Para a região de fluído:

,1121

2

21

2

XBM

YV

XV

,0121

2

21

2

XVM

YB

XB

Afora outras hipóteses adicionais, variáveis adimensionais são introduzidas para simplificar as equações diferenciais:

Para a região interna do ducto:

.022

2

22

2

YB

XB

As condições de contorno são

.,,0 21

12211 n

BnBBBV

Sob as condições da chamada aproximação de Shercliff, fazemos

,0,0 11

1

BnBV onde .21 ha

A figura a seguir ilustra as duas situações para condições de contorno:

h

a

X

Y

Y

X

Domínios computacionais: a) para o problema inteiro,b) usando a aproximação de Shercliff.

a)

b)

Nas simulações efetuadas, consideramos apenas as condições de contorno referentes a aproximação de Shercliff.

Aplicando o método de Galerkin aos sistemas de EDPs, temos:f,vK

,0

0,

00,

0

0,

0

0

y

xM

x

y

y

xI

I

I

I

I

I

I

I

IIII DCBA

,......11 NN BUBUTv

,

...

..................

...

...

,,,,,,1,,1,1,

,,,,,1,,1,,1,

,,1,,1,1,1,1,1,11,1

,1,,1,,11,1,1,1,1,1

yNyNxNxNNxNyyNxxNxN

NxNyNyNxNxNxNyyNxxN

yNyxNxNxyyxxx

NxyNyxNxxyyxx

MMMM

MMMM

M

Ω

,dDCBBAAMdK JTIJ

TIJ

TIJI

,0, II d

TII TTf

Como agora temos funções de forma interpolantes, condições decontorno de Dirichlet podem ser introduzidas diretamente no sistema f.vK Se usarmos contudo multiplicadores de Lagrange, esse sistema émodificado para

,

0fv

0GGK

T .

ds IKIK NG

As duas abordagens (com e sem multiplicadores de Lagrange)foram testadas neste problema, e os resultados comparados comaqueles obtidos pelo Método dos Elementos Finitos.

Testes comparativos:Em todas as simulações, tomaram-se 10406 nós para MEF em primeira ordem e EFGM, e 11801 nós para MEF em segunda ordem. Os valores de números de Hartman foram moderadamente altos (M=500) e as paredes do ducto são supostas como isolantes, uma vez que só existem soluções analíticas para este caso.

Arranjos dosnós:

0.00 0.40 0.80 1.20Normalized X Coordinate

0.0E+0

5.0E-4

1.0E-3

1.5E-3

2.0E-3

2.5E-3

Norm

alize

d Ve

loci

ty

Distribuição das velocidades ao longo do contorno (0X1,Y=0).

Campo de velocidades paratodo o domínio.

0.00 0.40 0.80 1.20Normalized X Coordinate

-2.0E-3

-1.5E-3

-1.0E-3

-5.0E-4

0.0E+0

5.0E-4

Nor

mal

ized

Indu

ced

Mag

netic

Fie

ld

 

Campo magnético induzido ao longo do percurso 0X1,Y=0).

Campo magnético induzidoem todo o domínio.

0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

0 . 9 9 0 0 . 9 9 2 0 . 9 9 4 0 . 9 9 6 0 . 9 9 8 1 . 0 0 0 N o r m a l i z e d X C o o r d i n a t e

0 . 0 E + 0

4 . 0 E - 4

8 . 0 E - 4

1 . 2 E - 3

1 . 6 E - 3

2 . 0 E - 3

N o r m a l i z e d V e l o c i t y

Resultados para EFGM com funções peso singulares (pontos) e solução analítica (curva) para um percurso dentro da camada de contorno (0.99X1, Y=0).

Isolinhas do campo induzido em todoo domínio.

Method CPU time (s) RC Number of non zero elements

EFGM with singular weight functions

154.31 7.43x10-4 989,543

EFGM with Lagrange multipliers 269.24 8.02x10-4 1,027,000

First order FEM 28.64 6.65x10-4 358,042

Second order FEM 94.86 1.70x10-4 705,104

Tabela comparativa de desempenho:

Erros relativos para EFGM com multiplicadores de Lagrange, em todo o domínio.

Erros relativos para funções peso singulares, em todo o domínio.

. Erros relativos em todo o domínio para o MEF em

primeira ordem.

0 . 1 0 0 . 3 0 0 . 5 0 0 . 7 0 0 . 9 0 0 . 0 0 0 . 2 0 0 . 4 0 0 . 6 0 0 . 8 0 1 . 0 0 Y C o o r d i n a t e ( X = 0 . 9 9 9 7 7 7 )

0 . 5 0

1 . 5 0

2 . 5 0

0 . 0 0

1 . 0 0

2 . 0 0

3 . 0 0

R e l a t i v e E r r o r ( % ) Erros relativos entre solução

analítica e EFGM usando funções peso singulares (linha sólida), EFGM usando multiplicadores de Lagrange (linha pontilhada) e MEF em primeira ordem (linha tracejada).

Algumas referências:

[1]     S. L. L.Verardi, J. M. Machado, Y. Shiyou, "The application of interpolating MLS approximations to the analysis of MHD flows", Finite Elements in Analysis and Design, vol. 39, pp. 1173-1187, 2003.[2]     S. L.Ho, Y. Shiyou, J. M. Machado, H. C. Wong, “Application of a Meshless Method in Electromagnetics”, IEEE Transactions on Magnetics, vol. 37, nr. 5, pp. 3198-3202, 2001.[3]     T. Belytschko, Y. Y. Lu, L. Gu, “Element Free Galerkin Methods”, International Journal for Numerical Methods in Engineering, vol. 37, pp. 229-256, 1994.[4]     P. Breitkopf, A. Rassineux, G. Touzot, P. Villon, “Explicit form and efficient computation of MLS shape functions and their derivatives”, International Journal for Numerical Methods in Engineering, vol. 48, pp. 451-466, 2000.[5]     L. W. Cordes, B. Moran, “Treatment of material discontinuity in the element-free Galerkin method”, Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 75-89, 1996.[6]     S. Beissel, and T. Belytschko, “An introduction to programming the meshless element-free Galerkin method”, Arch. Comp. Mechanics, vol. 5 (3), pp.207-241, 1998.[7]     Y. Krongauz, T. Belytschko, “EFG approximation with discontinuous derivatives”, Int. J. Num. Meth. Engn., v. 41, pp. 1215-1233, 1998.[8]     G. N. Marques, “O método element-free Galerkin: uma abordagem interpolante aplicada ao estudo de campos eletromagnéticos”. M.Sc. dissertation, Departament of Computer Science and Statistics – DCCE, IBILCE / Sao Paulo State University – UNESP, 2002.