MÉTODOS DE ELEMENTOS FINITOS PARA AS EQUAÇÕES DE …saulopo/maxwell/MiniCursoBelem_1.pdf · Saulo P. Oliveira M´etodos de Elementos Finitos para as Equac¸ oes de Maxwell 3/61˜

Universidade Federal do Para

Coordenacao do Programa de Pos-Graduacao em Geofsica

METODOS DE ELEMENTOSFINITOS PARA AS EQUACOES

DE MAXWELL

Saulo Pomponet Oliveira

Departamento de Matematica, Universidade Federal do Parana, Curitiba-PR

Cronograma

22/08 : Equacoes de MaxwellMetodos de elementos finitos para a equacao do potencialEquacoes de Maxwell no regime harmonicoMetodos de elementos finitos nodais

23/08 : Metodos de elementos finitos de arestaPanorama da pesquisa na area

Saulo P. Oliveira Metodos de Elementos Finitos para as Equacoes de Maxwell 2/61

Sumario - Quarta-feira 22/08Equacoes de Maxwell

Interpretacao fsica (meios homogeneos)Equacao do Potencial Eletrico

Formulacao variacionalMetodo de Galerkin

Metodos de Elementos Finitos (2D)Geracao da malhaElemento de referenciaAlgoritmo de montagemExemplos

Equacao de Maxwell no regime harmonicoFormulacao VariacionalFormulacao Variacional no plano

Metodos de Elementos Finitos NodaisExemplo


Equacoes de Maxwell

Lei de Ampere-Maxwell

Lei de Faraday

Lei de Gauss

Lei de Gauss do magnetismo

Grandezas envolvidas: E: campo eletrico H: campo magnetico D: inducao eletrica B: inducao magnetica J : densidade de corrente eletrica : densidade de carga eletrica


Equacoes de Maxwell:

D

tH + J = 0

B

t+E = 0

D =

B = 0

Equacoes de MaxwellMotivacao

As equacoes de Maxwell servem como modelo para umavasta gama de fenomenos eletromagneticos

A solucao destas equacoes e relevante em Comunicacao sem fio

Exames medicos nao-invasivos

Levantamento nao-destrutivo do subsolo


Equacoes de Maxwell:

D

tH + J = 0

B

t+E = 0

D =

B = 0

Equacoes de MaxwellEquacoes constitutivas

Vamos assumir relacoes constitutivas lineares:

D = E (x) : tensor de permissividade eletricaB = H (x) : tensor de permeabilidade magnetica

Substituindo as eqs. constitutivas:





Caso Isotropico: (x) =

(x) 0 00 (x) 00 0 (x)

Substituindo as

eqs. constitutivas:





Caso Isotropico: (x) =

(x) 0 00 (x) 00 0 (x)

Substituindo

as eqs. constitutivas:






D

tH + J = 0

B

t+E = 0

D =

B = 0






H Et

= J

H

t+E = 0

(E) =

(H) = 0


Equacoes de MaxwellCalculo Vetorial 3D

Sejam u =

u1u2u3

, v = v1v2v3

e w = w1w2w3

Produto vetorial: u v =

u2v3 u3v2u3v1 u1v3u1v2 u2v1

=i j ku1 u2 u3v1 v2 v3

Rotacional: v =

yv3 zv2zv1 xv3xv2 yv1

=i j kx y zv1 v2 v3

Produto misto: (w,u,v) = w (u v) =

w1 w2 w3u1 u2 u3v1 v2 v3



Operadores compostos:

( v)

= 0

Analogamente, () =

i j kx y zx y z

= 0 = 0Em domnios simplesmente conexos (Bossavit, 1998): w = v w = 0 w = w = 0

Teorema de Helmholtz: Dado um campo vetorial w, v, tais que w = v +

Expressamos w em termos de dois potenciais:v: potencial vetorial (termo solenoidal)q: potencial escalar (termo irrotacional)




( v) = (,,v)

= 0

Analogamente, () =

i j kx y zx y z







( v) =

x y zx y zv1 v2 v3

= 0

Analogamente, () =

i j kx y zx y z







( v) = 0

Analogamente, () =

i j kx y zx y z







( v) = 0

Analogamente, () =

i j kx y zx y z

() =i j kx y zx y z







( v) = 0

Analogamente, () =

i j kx y zx y z

= 0 = 0

Em domnios simplesmente conexos (Bossavit, 1998): w = v w = 0 w = w = 0






( v) = 0

Analogamente, () =

i j kx y zx y z

= 0 = 0OBS: estas operacoes sao formalizadas e generalizadas nocontexto de formas diferenciais (Bossavit, 1988)







( v) = 0

Analogamente, () =

i j kx y zx y z

= 0 = 0Consequencias: Se w = v, entao w = 0 Se w = , entao w = 0







( v) = 0

Analogamente, () =

i j kx y zx y z





Equacoes de MaxwellMeios homogeneos

Para meios homogeneos,

H Et

= J

H

t+E = 0

(E) =

(H) = 0

H Et

= J

H

t+E = 0

E = /

H = 0 H = 0 =H = A, logo

A = J + Et

E = At

E = /




H Et

= J

H

t+E = 0

E =

H = 0

H Et

= J

H

t+E = 0

E = /

H = 0 H = 0 =H = A, logo

A = J + Et

E = At

E = /




H Et

= J

H

t+E = 0

E = /

H = 0

H = 0 =H = A, logo

A = J + Et

E = At

E = /




H Et

= J

H

t+E = 0

E = /

H = 0 H = 0 =H = A, logo

(A) = J + Et

E = (A)t

E = /

A = J + Et

E = At

E = /




H Et

= J

H

t+E = 0

E = /

H = 0

H = 0 =H = A, logo

A = J + Et

E = At

E = /



Regime dependente do tempo (dinamico):

A = J + Et

E = At

E = / (H = A)

Por convencao, E = , de modo que

A = J = / (E = , H = A)

Interpretacao (S. Guerreiro) :Estatico: as fontes de E (H) sao cargas (correntes) eletricasOndas: a oscilacao de E gera fonte em H, e vice-versa




A = J + Et

E = At

E = / (H = A)Considerando o regime estatico,

A = JE = 0 E = /

Por convencao, E = , de modo queA = J

= / (E = , H = A)Interpretacao (S. Guerreiro) :Estatico: as fontes de E (H) sao cargas (correntes) eletricasOndas: a oscilacao de E gera fonte em H, e vice-versa




A = J + Et

E = At


A = JE =

E = /






A = J + Et

E = At


A = J () = /






A = J + Et

E = At

E = / (H = A)

Por convencao, E = , de modo que

A = J = / (E = , H = A)

Interpretacao (S. Guerreiro) :Estatico: as fontes de E (H) sao cargas (correntes) eletricasOndas: a oscilacao de E gera fonte em H, e vice-versa


Equacao do Potencial EletricoDey & Morrison (1979)

Da lei de Gauss,

(E) =

(E

t

)=

t

Da lei de Ampere-Maxwell, H Et = J , logo

J = t

Substitucao: lei de Gauss conservacao de carga:

H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


J = t


H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


(H J) = t

J = t


H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


(H) J = t

J = t


H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


J = t


H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


J = t


H Et

= J

H

t+E = 0

(E) =

(H) = 0

H Et

= J

H

t+E = 0

J = t

(H) = 0



Da lei de Gauss,

(E

t

)=

t


J = t


H Et

= J

H

t+E = 0

J = t

(H) = 0Saulo P. Oliveira Metodos de Elementos Finitos para as Equacoes de Maxwell 11/61


Corrente eletrica (Ward e Hohmann, 88): I(t) =J(x, t) n d

I(t) =

t(x, t) d

Assumindo t = Ioxo(x) = f(x) e baixas frequencias J = f(x)

Vamos considerar a lei de Ohm J = E, assim,

((x)(x)) = f(x)

Vamos considerar o seguinte problema: ((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2, 1 2 = .



Corrente eletrica (Ward e Hohmann, 88): I(t) =

J(x, t) d

I(t) =

t(x, t) d

Assumindo t = Ioxo(x) = f(x) e baixas frequencias J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d

Assumindo t = Ioxo(x) = f(x) e baixas frequencias

J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d

Fonte pontual: (x, t) = o(t)xo(x) =o(t)t = I(t)


J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d


H Et

= J

H

t+E = 0

J = t

(H) = 0

J = f(x)Vamos considerar a lei de Ohm J = E, assim,

((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d

Assumindo t = Ioxo(x) = f(x) e baixas frequenciasH = JE = 0 J = f(x)

(H) = 0


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d


E = 0 J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d


E = J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d

Assumindo t = Ioxo(x) = f(x) e baixas frequenciasE =


E = (E) = f(x)

((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .




t(x, t) d


E = J = f(x)


((x)(x)) = f(x)


(x) = 0, x 1(x) n = 0, x 2, 1 2 = .


Formulacao VariacionalExemplo 1D

Considere o seguinte problema de valores de contorno:{u(x) = f(x), 0 < x < 1u(0) = u(1) = 0

Para f C([0, 1]), u V = {u C2([0, 1]) | u(0) = u(1) = 0}.Dado v V , 1

0u(x)v(x) dx =

10f(x)v(x) dx

Integrando por partes: 10u(x)v(x) dx =

10f(x)v(x) dx




Para f C([0, 1]), u V = {u C2([0, 1]) | u(0) = u(1) = 0}.

Dado v V , 10u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx




Para f C([0, 1]), u V = {u C2([0, 1]) | u(0) = u(1) = 0}.Dado v V ,

10u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx





u(x) = f(x)

10u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx





u(x)v(x) = f(x)v(x)

10u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx





0u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx





0u(x)v(x) dx =

10f(x)v(x) dx

Integrando por partes:

u(1)v(1) + u(0)v(0) + 1

0u(x)v(x) dx =

10f(x)v(x) dx

10u(x)v(x) dx =

10f(x)v(x) dx





0u(x)v(x) dx =

10f(x)v(x) dx


10f(x)v(x) dx



A solucao do problema de valores de contorno{u(x) = f(x), 0 < x < 1u(0) = u(1) = 0

satisfaz 10u(x)v(x) dx =

10f(x)v(x) dx

Novo espaco das solucoes:

V =

{v L2([0, 1])

10v(x)2 + v(x)2 dx

Equacao do Potencial EletricoEspaco vetorial da solucao

Dado u : IR, vamos definir

u0 =(

|u(x)|2 d

)1/2|u|1 =

(u(x) u(x) d

)1/2

|u|1 =(

|u(x)|2 d

)1/2u1 =

(u20 + |u|21

)1/2Usaremos os seguintes espacos:

L2() = {v : IR | u0



u0 =(

|u(x)|2 d

)1/2|u|1 =

(|u(x)|2 d

)1/2u1 =

(u20 + |u|21

)1/2

Usaremos os seguintes espacos:

L2() = {v : IR | u0



u0 =(

|u(x)|2 d

)1/2|u|1 =

(|u(x)|2 d

)1/2u1 =

(u20 + |u|21

)1/2Usaremos os seguintes espacos:

L2() = {v : IR | u0

Equacao do Potencial EletricoFormulacao Variacional

((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2

Seja V = {v H1() | v|1 = 0}.

Multiplicando a equacao por v V e integrando em :

1a. Identidade de Green: escolhendo F = v

((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2

Seja V = {v H1() | v|1 = 0}.Multiplicando a equacao por v V e integrando em :


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x)) = f(x)


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) = f(x)v(x)


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d

Integrar por partes ?


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d

Teorema da divergencia de Gauss: F (x) d =

F (x) n d


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d

1a. Identidade de Green: escolhendo F = v ((x)(x)v(x)) d =

((x)(x)v(x)) n d

((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d

1a. Identidade de Green: escolhendo F = v((x)(x))v(x)+(x)(x)v(x)d =

v(x)(x)(x)nd

((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d


((x)(x))v(x)d =

v(x)(x)(x)nd

(x)(x)v(x)d

((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd

((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2


((x)(x))v(x) d =

f(x)v(x) d


((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd



((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd

Como v V , temos v |1= 0:

v(x)(x)(x) n d =

2

v(x)(x)(x) n d

Da condicao de contorno em 2, n |2= 0, logov(x)(x)(x) n d = 0

Assim,

((x)(x))v(x) d =

(x)(x) v(x) d



((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd


v(x)(x)(x)nd =

1

v(x)(x)(x)nd

2

v(x)(x)(x)nd

v(x)(x)(x) n d =

2

v(x)(x)(x) n d


Assim,

((x)(x))v(x) d =

(x)(x) v(x) d



((x)(x))v(x)d =

(x)(x)v(x)d

v(x)(x)(x)nd


v(x)(x)(x) n d =

2

v(x)(x)(x) n d


Assim,

((x)(x))v(x) d =

(x)(x) v(x) d



Vimos que ((x)(x)) = f(x), x

(x) = 0, x 1(x) n = 0, x 2

(x)(x) v(x) d =

f(x)v(x) d

Definindo

a(u, v) =

(x)u(x)v(x)d e F (v) =

f(x)v(x)d,

podemos escrever a nova formulacao como:

Formulacao Variacional: encontrar V tal quea(, v) = F (v) v V

V = {v H1() | v|1 = 0}



Vimos que

((x)(x))v(x) d =

f(x)v(x) d

(x)(x) v(x) d =

f(x)v(x) d

Definindo

a(u, v) =


f(x)v(x)d,



V = {v H1() | v|1 = 0}



Vimos que(x)(x) v(x) d =

f(x)v(x) d

Definindo

a(u, v) =


f(x)v(x)d,



V = {v H1() | v|1 = 0}


Equacao do Potencial EletricoMetodo de Galerkin

a(, v) = F (v) v V

Dados v1, v2, . . . vn V , L.I., defina Vh V como sendo

Vh = span{v1, . . . vn} = {u V | u(x) = u1v1(x) + . . .+ unvn(x)}

Metodo de Galerkin: encontrar h Vh tal quea(h, vh) = F (vh) vh Vh

h Vh = h(x) =nj=1

jvj(x) (n incognitas)

Escolhendo vh = v1, . . . , vn, obtemos as seguintes n equacoes:

Au = b



a(, v) = F (v) v V

Dados v1, v2, . . . vn V , L.I., defina Vh V como sendo



h Vh = h(x) =nj=1



a(h, vi) = F (vi), 1 i n

Au = b



a(, v) = F (v) v VDados v1, v2, . . . vn V , L.I., defina Vh V como sendoVh = span{v1, . . . vn} = {u V | u(x) = u1v1(x) + . . .+ unvn(x)}


h Vh = h(x) =nj=1



a

nj=1

jvj , vi

= F (vi), 1 i n

Au = b



a(, v) = F (v) v VDados v1, v2, . . . vn V , L.I., defina Vh V como sendo



h Vh = h(x) =nj=1


Escolhendo vh = v1, . . . , vn, obtemos as seguintes n equacoes:nj=1

j a(vj , vi) = F (vi) 1 i n

Au = b



a(, v) = F (v) v VDados v1, v2, . . . vn V , L.I., defina Vh V como sendo



h Vh = h(x) =nj=1


Escolhendo vh = v1, . . . , vn, obtemos as seguintes n equacoes:nj=1

a(vj , vi)j = F (vi) 1 i n

Au = b



Relacao entre e h ?

a(, v) = F (v) v Va(h, vh) = F (vh) vh Vh

Assim,

a( h, h) =Sabendo que a(u, u) C1u21 e a(u, v) C2u1v1,

h1 C2C1 vh1

Como vh Vh foi arbitrario,

h1 minvhVh

C2C1 vh1



Relacao entre e h ?

a(, vh) = F (vh) vh Vha(h, vh) = F (vh) vh Vh

a(, vh) a(h, vh) = 0 vh Vh

a( h, vh) = 0 vh Vh

Assim,


h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?


a( h, vh) = 0 vh Vh

Assim,


h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?


a( h, vh) = 0 vh Vh

Assim,

a( h, h) = a( h, vh + vh h)

Sabendo que a(u, u) C1u21 e a(u, v) C2u1v1,

h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?


a( h, vh) = 0 vh Vh

Assim,

a( h, h) = a( h, vh) + a( h, vh h)


h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?


a( h, vh) = 0 vh VhAssim,

a( h, h) = a( h, vh) + a( h, vh h)


h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?



a( h, h) = a( h, vh)


h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?



a( h, h) = a( h, vh)Sabendo que a(u, u) C1u21 e a(u, v) C2u1v1,

C1 h21 a( h, h)

h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?




C1 h21 a( h, vh)

h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?




C1h21 a(h, vh) C2h1vh1

h1 C2C1vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?




C1 h21 C2 h1 vh1

h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?




C1 h1 C2 vh1

h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



Relacao entre e h ?




h1 C2C1 vh1


h1 minvhVh

C2C1 vh1



h1 minvhVh

C2C1 vh1

Vh

V!

Vh

V!

!h



h1 minvhVh

C2C1 vh1

Vh

V!

!h


Metodo de Elementos Finitos

Os metodos de elementos finitos escolhem as funcoes debase N1, . . . Nn com base numa discretizacao do domnio

Seja = [0, 1] [0, 1].

Associamos uma funcao de base a cada vertice (exceto em 1)


Metodo de Elementos FinitosFuncoes de base

Dados os vertices x1, . . . ,xn, escolhemos N1, . . . Nn tais que: Ni(x) e contnua;

Ni(x) |Ke= Ni(x, y) |Ke= ae + bex+ cey (Ni |Ke tem grau 1)

Ni(xj) ={

1, i = j0, i 6= j

0

0.5

1 0

0.5

1

0

0.5

1

yx



Dada f(x), sejam f1 = f(x1), . . . , fn = f(xn) e

f(x) = f1N1(x) + f2N2(x) + . . .+ fn1Nn1(x) + fnNn(x)

Temos que f Vh. Alem disso,

f(x2) = f(x2)

Em geral, f(xi) = f(xi), ou seja, f interpola f em x1, . . .xn

0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x1) = f1N1(x1)+f2N2(x1)+. . .+fn1Nn1(x1)+fnNn(x1)

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x1) = f1(1) + f2(0) + . . .+ fn1(0) + fn(0)

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x1) = f1

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x1) = f(x1)

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x2) = f1N1(x2)+f2N2(x2)+. . .+fn1Nn1(x2)+fnNn(x2)

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x2) = f1(0) + f2(1) + . . .+ fn1(0) + fn(0)

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x2) = f2

f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1





Temos que f Vh. Alem disso,f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1

0

0.5

1 0

0.5

1

1

0.5

0

0.5

1






f(x2) = f(x2)


0

0.5

1 0

0.5

1

1

0.5

0

0.5

1



Lembrando que

h1 minvhVh

C2C1 vh1,

temos em particular que

h1 C2C1 1

Assim, 1 0 = h1 0


Metodo de Elementos FinitosImplementacao

Passos do metodo de elementos finitos:

1 Geracao da malha

2 Calculos no elemento de referencia

3 Algoritmo de montagem

4 Solucao do sistema linear

5 Visualizacao, pos-processamento


Geracao da malhaPacote utilizado

Utilizaremos o pacote triangle :

riangleA Two-Dimensional Quality MeshGenerator and DelaunayTriangulator.

Jonathan Richard ShewchukComputer Science DivisionUniversity of California atBerkeleyBerkeley, California 94720-1776

Winner of the 2003 James Hardy Wilkinson Prize in Numerical Software.Created at Carnegie Mellon University as part of the Quake project (tools for large-scale earthquakesimulation).Supported by an NSERC 1967 Science and Engineering Scholarship and NSF Grant CMS-9318163.

Triangle generates exact Delaunay triangulations, constrained Delaunay triangulations, conformingDelaunay triangulations, Voronoi diagrams, and high-quality triangular meshes. The latter can begenerated with no small or large angles, and are thus suitable for finite element analysis.

Triangle (version 1.6, with Show Me version 1.6) is available as a .zip file (159K) or as a .shar file(829K) (extract with sh) from Netlib in the voronoi directory. Please note that although Triangle isfreely available, it is copyrighted by the author and may not be sold or included in commercial productswithout a license.

New features in Version 1.6 (released July 28, 2005). Improved handling of domains with smallangles (thanks to an algorithm of Gary Miller, Steven Pav, and Noel Walkington). In particular, Trianglenow offers a no-large-angle guarantee even for domains that have lots of tiny input angles (which makea no-small-angle guarantee impossible). Meshes sometimes have fewer triangles than in previousversions, thanks to two changes. First, Triangle now uses Paul Chew's Delaunay refinement algorithm,which is more conservative about splitting segments than previous versions of Triangle when the anglebound is under 30 degrees. (Ruppert's algorithm is still available through the -D switch, offeringall-Delaunay meshes.) Second, a change in the priority queue of bad triangles (suggested by Alperngr) yields fewer triangles when the angle bound is large. Many bugs are fixed, including three bugsthat were causing segmentation faults. (If you use Triangle version 1.5, I urge you to replace itimmediately. Earlier versions are stable, though.)

Of special interest. Two papers about Triangle are available. These and related papers are availablefrom the Research Credit page. Triangle's robust geometric predicates are available separately from theRobust Predicates page. These geometric predicates are in the public domain (though Triangle is not).

Instructions for using Triangle

A demonstration of what Triangle can doHow fast is Triangle?

Triangle: A Two-Dimensional Quality Mesh Generator and D... http://www.cs.cmu.edu/~quake/triangle.html

1 of 2 06/02/12 13:11


Geracao da malhaDados relevantes / meios homogeneos

O gerador de malha deve fornecer: Numero de vertices Coordenadas dos vertices

Numero de elementos (no caso, triangulos) Matriz de conectividade dos elementos

Identificacao dos vertices na fronteira


Geracao da malhaDados relevantes / arquivo gerador da malha do Exemplo 1

# quad.poly: geracao de uma malha em um dominio quadrado## Numero de vertices, dimensao (2), numero de atributos e # e numero de ident. de fronteira (1: Dirichlet; 2: Neumann)#4 2 0 1## Vertices#1 0 0 12 0 1 13 1 1 14 1 0 1## Linhas (sobretudo do contorno)#4 11 1 2 12 2 3 13 3 4 14 4 1 1## Numero de furos no dominio (seguido pelas coords)#0## Numero de sub-regioes (seguido pelas coords)#0


Geracao da malhautilizacao do pacote triangle

Sintaxe (ha mais opcoes):

./triangle -qANG -aA INPUT.poly

ANG: angulo mnimo dos triangulos (ex: 30 (30o)) A: area maxima dos triangulos (ex: 0.01 (0.01u.a.)) INPUT: nome do arquivo de entrada (extensao .poly)

Aplicativo adicional para visualizar a malha gerada:

./showme INPUT.poly


Geracao da malhaArquivos gerados

quad.1.node: quad.1.ele:9 2 0 1 1 0 0 1 2 0 1 1 3 1 1 1 4 1 0 1 5 0.5 0.5 0 6 0 0.5 1 7 0.5 0 1 8 1 0.5 1 9 0.5 1 1# Generated by ./triangle -q -a0.2 quad.poly

8 3 0 1 8 3 5 2 6 1 5 3 5 9 2 4 5 2 6 5 5 1 7 6 7 4 5 7 5 4 8 8 5 3 9# Generated by ./triangle -q -a0.2 quad.poly

Do arquivo quad.1.node: Numero de vertices: 9

Coordenadas dos vertices:[0, 0], [0, 1], . . . [0.5, 1]

Identificacao dos verticesna fronteira:[1, 1, 1, 1, 0, 1, 1, 1, 1]


Geracao da malhaExtracao dos dados relevantes

Do arquivo quad.1.ele: Numero de elementos: 8

Matriz de conectividade:

8 6 5 53 1 9 35 5 2 9


Elemento de referenciaHa uma transformacao do triangulo K para qualquer triangulo:

O triangulo K e definido por x1 = (1, 0), x2 = (0, 1), x3 = (0, 0)

Vamos escrever a transformacao na forma{x(x, y) = x1N1(x, y) + x2N2(x, y) + x3N3(x, y)

y(x, y) = y1N1(x, y) + y2N2(x, y) + y3N3(x, y)

Vetorialmente, x(x) = x1N1(x) + x2N2(x) + x3N3(x)


Elemento de referenciaFuncoes de base locais

x(x) = x1N1(x) + x2N2(x) + x3N3(x)

Para que x(x1) = x1, vamos impor

N1(x1) = 1, N2(x1) = 0, N3(x1) = 0

Analogamente,

x(x2) = x2 = N1(x2) = 0, N2(x2) = 1, N3(x2) = 0x(x3) = x3 = N1(x3) = 0, N2(x3) = 0, N3(x3) = 1



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja . b1 = 1

c1 = 0

a1 = 0

Portanto, N1(x, y) = x. Ao final, encontramos:

N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.

OBS: A transformacao x = x(x) e linear, i.e.,

x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .

b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .N1(x1) = 1

N1(x2) = 0

N1(x3) = 0

b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .N1(1, 0) = 1

N1(0, 1) = 0

N1(0, 0) = 0

b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .a1 + b1(1) + c1(0) = 1

a1 + b1(0) + c1(1) = 0

a1 + b1(0) + c1(0) = 0

b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .a1 + b1 = 1

a1 + c1 = 0

a1 = 0

b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j

Seja N1(x) = N1(x, y) = a1 + b1x+ c1y .b1 = 1

c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j


c1 = 0

a1 = 0

Portanto, N1(x, y) = x.

Ao final, encontramos:

N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j


c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j


c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = x1N1(x) + x2N2(x) + x3N3(x)

x(x) = JKx+ x0



Assim,

Ni(xj) =

{1, i = j0, i 6= j


c1 = 0

a1 = 0


N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


x(x) = JKx+ x0


Elemento de referenciaRelacao entre Ni e Ni

Sejam = K, = 2 e x(x) = JKx+ x0.

Neste caso, Vh = span{N1, N2, N3}.

Temos que

Ni(x(xj)) =

{1, i = j0, i 6= j , Ni(xj) =

{1, i = j0, i 6= j

Ni(x(xj)) = Ni(xj), j = 1, 2, 3

Como Ni e Ni sao lineares, Ni(x(x)) = Ni(x).

Pode-se mostrar (Jin, 2003) que Ni(x(x)) = J1K Ni(x)





Temos que

Ni(xj) =

{1, i = j0, i 6= j , Ni(xj) =

{1, i = j0, i 6= j , x(xj) = xj

Ni(x(xj)) =

{1, i = j0, i 6= j , Ni(xj) =

{1, i = j0, i 6= j

Ni(x(xj)) = Ni(xj), j = 1, 2, 3







Temos que

Ni(x(xj)) =

{1, i = j0, i 6= j , Ni(xj) =

{1, i = j0, i 6= j

Ni(x(xj)) = Ni(xj), j = 1, 2, 3




Elemento de referenciaCalculo das Integrais

Vamos aproximar f(x) por f(x) =3

j=1 fjNj(x),

Kf(x)Ni(x) d

3j=1

fj

KNj(x)Ni(x) |JK |d

Por outro lado, tomando (x) |K K ,K(x)u(x)v(x)d K

K

(J1K Nj(x))(J1K Ni(x))|JK |d




j=1 fjNj(x),Kf(x)Ni(x) d

3j=1

fj

KNj(x)Ni(x) d

Kf(x)Ni(x) d

3j=1

fj

KNj(x)Ni(x) |JK |d


K






3j=1

fj

KNj(x(x))Ni(x(x)) |JK |d

Kf(x)Ni(x) d

3j=1

fj

KNj(x)Ni(x) |JK |d


K






3j=1

fj

KNj(x)Ni(x) |JK |d


K






3j=1

fj

KNj(x)Ni(x) |JK |d

N1(x, y) = x, N2(x, y) = y, N3(x, y) = 1 x y.


K






3j=1

fj

KNj(x)Ni(x) |JK |d


K



Elemento de referenciaIntegracao numerica

Visando automacao, calculamos as integrais com quadraturas:K

Nj(x)Ni(x) |JK |d nintl=1

Nj(pl)Ni(pl) |JK |wlK

(J1K Nj(x)) (J1K Ni(x)) |JK |d

nintl=1

(J1K Nj(pl) (J1K Ni(pl)) |JK |wl

OBS: Para funcoes de base de ordem mais alta, |JK | = |JK(pl)|. E possvel usar pl = xl

Ex: elementos espectrais GLL (Komatitsch e Tromp, 1999)


Elemento de referenciaIntegracao numerica

Visando automacao, calculamos as integrais com quadratura

Documents

MÉTODOS DE ELEMENTOS FINITOS PARA AS EQUAÇÕES DE …saulopo/maxwell/MiniCursoBelem_1.pdf · Saulo P. Oliveira M´etodos de Elementos Finitos para as Equac¸ oes de Maxwell 3/61˜