Método dos Elementos Finitos
Teoria Eletromagnética IIC
Prof. Antonio Lopes de Souza, Ph.D.
(Modelo estrutural em elementos finitos do Cubo d’Água de Pequim)
Soluções Para Problemas do
Eletromagnetismo
Solução Analítica: solução na forma de uma equação algébrica explícita, na
qual todos os valores dos parâmetros do problema podem ser substituídos.
Ex.: campo elétrico da carga pontual

E
Q
4 0 R 2

ar (V/m)
Vantagens das soluções analíticas: são exatas e tornam mais fácil observar o
comportamento da solução em função da variação dos parâmetros do
problema.
Figura 2 - Linhas de força do
campo elétrico de um sistema de
duas cargas pontuais
Desvantagem: somente são possíveis para problemas com configurações
simples e elevado grau de simetria.
Soluções Para Problemas do
Eletromagnetismo
Solução não Analítica: é usada quando a complexidade do problema torna
difícil a obtenção de uma solução analítica por métodos matemáticos
tradicionais. Os procedimentos não-analíticos incluem
Métodos gráficos
Métodos experimentais
Métodos numéricos
Os métodos gráficos e experimentais são aplicados a um número reduzido de
problemas e têm utilidade limitada.
Os métodos numéricos mais usados são:
Método das Diferenças Finitas (FDM - Finite Diference Method)
Método dos Momentos (MOM - Method of Moments)
Método dos Elementos Finitos (FEM – Finite Diference Method)
Soluções Para Problemas do
Eletromagnetismo
Solução Não-analítica: exemplo
Figura 3 - Visualização da deformação de um carro em um choque assimétrico usando
o método dos elementos finitos (imagem: Wikimedia Commons)
Método dos Elementos Finitos
O Método dos Elementos Finitos (FEM – Finite Element Method) é uma técnica
numérica usada para encontrar a solução de Equações Diferenciais Parciais (PDE).
Ele surgiu na década de 40 do século passado e era usado para resolver problemas
estruturais nas engenharias civil e aeronáutica.
Os primeiros trabalhos sobre o método foram publicados por Alexander Hrennikoff (1)
(1941) e Richard Courant (2) (1942). O uso do método para a soluções de problemas
de eletromagnetismo data de 1968 (3).
Figura 5: Solução bidimensional de um
problema de magetostática mostrando a
direção (as linhas) e a intensidade (cores)
da densidade de fluxo magnético. (imagem:
Wikimedia Commons).
Esquerda: Malha do problema
Direita: Solução gráfica do problema
(1) A. Hrennikoff; “Solution of Problems of Elasticity by the Frame-Work Method” (1941); ASME J. Appl. Mech. 8; A619–A715.
(2) Courant, R.L;“Variational Methods for the Solution of Problems of Equilibrium and Vibration”;Bulletin of the American Mathematical Society 49:1-23
(1943) .
(3) P.P.silvester e R.L.Ferrari; Finite Elements for Electrical Enginiers, Cambridge, England, Cambridge University Press, (1983).
Método dos Elementos Finitos
Etapas na modelagem de um problema através do Método dos Elementos
Finitos:
 Descrição geométrica da região;
Geração de uma malha de elementos interconectados por nós;
Definição das equações diferenciais parciais e respectivas
de contorno do problema;
condições
Solução numérica do sistema algébrico resultante;
Pós-processamento de resultados e visualização
Um modelo em CAD pode ser fornecido a um programa de geração de malhas que
vai gerar a entrada de outro programa que contém um núcleo de solução numérica
e, finalmente, o resultado é visualizado em uma ferramenta de pós-processamento
gráfico.
Método dos Elementos Finitos
Seqüência na apresentação do método:
1. Discretização do domínio em um número finito sub-regiões ou
elementos.
2. Obtenção das equações que regem um elemento típico.
3. Conexão de todos os elementos no domínio.
4. Resolução do sistema de equações obtido.
Método dos Elementos Finitos
1- Discretização
A representação do domínio: o domínio é dividido em regiões (elementos) que
não se sobrepõem. Buscamos uma aproximação para o potencial Ve dentro de
um elemento. Inter-relacionamos as distribuições de potencial em vários
elementos de tal modo que ela seja contínua através dos contornos entre os
elementos relacionados.
A solução aproximada do potencial para o domínio em estudo é do tipo:
N
V ( x, y)  Ve ( x, y)
e
Onde N é o número de elementos triangulares nos quais o domínio foi dividido
Figura 6 – elemento triangular
com a indicação dos nós e
tensões de nós. A orientação dos
nós segue o sentido anti-horário.
Método dos Elementos Finitos
1- Discretização
Representação do potencial no interior do elemento: quando o elemento
tem a forma triangular a forma mais comumente usada para representar o
potencial Ve no seu interior é a aproximação polinomial
Ve (x, y)  a  bx  cy
Quando o elemento for quadrangular a aproximação para representar Ve é dada
por:
Ve (x, y)  a  bx  cy  dxy
De um modo geral usamos elementos triangulares para problemas
bidimensionais porque esses permitem representar fronteiras irregulares com
maior precisão.
Supomos que o potencial varia linearmente dentro do elemento. Isso implica
em considerar o campo elétrico é uniforme dentro de um mesmo elemento.
Como os potenciais respeitam a equação de Laplace a distribuição de energia
do campo elétrico é mínima na região.
Método dos Elementos Finitos
2- Equações que regem os elementos
O elemento usado tem a forma triangular abaixo mostrada:
O potencial no seu interior é representado por: Ve ( x, y)  a  bx  cy
Ve1  a  bx1  cy1

Aplicando a equação aos três vértices do elemento temos: Ve 2  a  bx2  cy2
Ve  a  bx  cy
3
3
 3
Método dos Elementos Finitos
2- Equações que regem os elementos
Representando o sistema de equações anterior na forma matricial temos:
Ve1  a  bx1  cy1
Ve1  1 x1 y1  a 
  

  
Ve 2  a  bx2  cy2  Ve2   1 x2 y2   b 
Ve  a  bx  cy
Ve3  1 x3 y3  c 
3
3
 3
Onde a matriz
1 x1 y1 


1
x
y
2
2


1 x y 
 3 3
É denominada “matriz dos coeficientes”
Método dos Elementos Finitos
2- Equações que regem os elementos
Os coeficientes a,b,c podem ser calculados a partir da equação anterior como:
1
Ve1  1 x1 y1  a 
Ve   1 x y   b 
 2  2 2  
Ve3  1 x3 y3  c 
a  1 x1 y1  Ve1 
b   1 x y   Ve 
   2 2  2
c  1 x3 y3  Ve3 
E usando a regra de CRAMER
Ve1 x1 y1 
1
a  Ve 2 x2 y2 
D
Ve3 x3 y3 
1 Ve1 y1 
1
b  1 Ve 2 y2 
D
1 Ve3 y3 
1 x1 Ve1 
1
c  1 x2 Ve 2 
D
1 x3 Ve3 
Onde D é o determinante da matriz dos coeficientes dado por
1 x1 y1 
D  1 x2 y2    x2 y3  x3 y1  x1 y2  x2 y1  x1 y3  x3 y2
1 x3 y3 

Método dos Elementos Finitos
2- Equações que regem os elementos
É possível mostrar que a área A do triangulo que representa o elemento é dada por:
D 1
A    ( x2  x1 )( y3  y1 )  ( x3  x1 )( y2  y1 )
2 2

(Os nós devem ser numerados no
sentido anti-horário para evitar que a
área do elemento, acima calculada em
função das coordenadas dos vértices,
produza um valor negativo)
Ou seja:
D  2A
Método dos Elementos Finitos
2- Equações que regem os elementos
Substituindo o valor D=2A em a,b,e c e resolvendo os determinantes temos:
Ve 1 x1 y1 
1 
  1  Ve ( x y  x y )  Ve ( x y  x y )  Ve ( x y  x y )
a
Ve
x
y
2
2
2
1
2 3
3 2
2
3 1
1 3
3
1 2
2 1
2A 
2A
Ve 3 x 3 y 3 
1 Ve 1 y1 
1 
  1  Ve ( y  y )  Ve ( y  y )  Ve ( y  y )
b
1
Ve
y
2
2
1
2
3
2
3
1
3
1
2
2A 
2A
1 Ve 3 y 3 
1 x1 Ve 1
1
c  1 x 2 Ve 2
D
1 x 3 Ve 3

  1  Ve ( x  x )  Ve ( x  x )  Ve ( x  x )
1
3
2
2
1
3
3
2
1
 2A




Método dos Elementos Finitos
2- Equações que regem os elementos
As equações anteriores podem ainda ser colocadas na forma matricial como
Ve 1 
1
( x 2 y3  x 3 y 2 ) (x3 y1 - x1y3 ) (x1y 2 - x 2 y1 )Ve 2 
a
2A
Ve 3 
Ve 1 
1
( y 2  y3 ) (y3 - y1 ) (y1 - y 2 )Ve 2 
b
2A
Ve 3 
Ve 1 
1
( x3  x 2) (x1- x3) (x2- x1)Ve 2 
c
2A
Ve 3 
E lembrando a forma matricial da equação que descreve as distribuições de
potencial
a 
V ( x , y)  a  bx  cy  1, x , y  b 
c 
Método dos Elementos Finitos
2- Equações que regem os elementos
Substituindo as expressões matriciais para a, b e c na equação matricial
que descreve as distribuições de potenciais no elemento temos:
( x 2 y 3  x 3 y 2 )
1
1 x y (y2 - y1 )
Ve ( x , y) 
2A
(x3 - x 2 )
(x3 y1 - x1 y 3 ) (x1 y 2 - x 2 y1 ) Ve 1 
 Ve 
(y3 - y1 )
(y1 - y 2 )
  2
 Ve 3 
(x1 - x 3 )
(x 2 - x1 )
Expandindo a equação acima temos:
1
Ve ( x, y) 
Ve 1( x 2 y 3  x 3 y 2 )  ( y 2  y 3 ) x  ( x 3  x 2 ) y 
2A
1

Ve 2  ( x 3 y1  x1 y3 )  ( y3  y1 ) x  ( x1  x 3 ) y 
2A
1

Ve 3  ( x1 y 2  x 2 y1 )  ( y1  y 2 ) x  ( x 2  x1 ) y 
2A
Ou seja
Ve ( x, y)  1Ve1   2 Ve 2  3Ve 3 
3
  Ve
i
i 1
i
Método dos Elementos Finitos
2- Equações que regem os elementos
e
1 
1
(x 2 y3  x 3 y 2 )  ( y 2  y3 )x  (x 3  x 2 ) y
2A

1
 (x 3 y1  x1y3 )  ( y3  y1 )x  (x1  x 3 ) y 
2A
1
 ( x1y 2  x 2 y1 )  ( y1  y 2 ) x  ( x 2  x1 ) y 
3 
2A
1, 2 e 3 são funções lineares de interpolação que permitem a obtenção
dos potenciais Ve(x,y) dentro do elemento finito em termos dos potenciais
nos nós do mesmo elemento. Elas são denominadas “Funções de Forma
dos Elementos”. Desse modo a função
3
2 
Ve ( x, y)  1Ve1   2 Ve 2  3Ve 3 
  Ve
i
i
i 1
nos dá o potencial em qualquer ponto (x,y) dentro do elemento finito desde
que os potenciais nos vértices sejam conhecidos.
Método dos Elementos Finitos
2- Equações que regem os elementos
As funções de
propriedades:
forma
do
elemento
satisfazem
as
seguintes
1, i  j
 i ( x , y)  
0, i  j
Cada uma das três funções de forma se anula em todos os vértices com
exceção de um no qual assume o valor unitário.
Cálculo da Energia Armazenada no Elemento
Quando a quantidade analisada é o potencial eletrostático a energia do
campo elétrico é o funcional a ser utilizado porque o potencial
eletrostático minimiza esta energia. A energia do campo elétrico por
unidade de comprimento normal às duas dimensões (ou seja, na
direção z) pode ser obtida por:
2
1
We 
 e dS
2

Método dos Elementos Finitos
2- Equações que regem os elementos
A função de energia pode ser expandida como abaixo:
We 

A
2
1
1
1
2
 E dS 
 Ve dS   (Ve  Ve )dS
2
2A
2 A


Lembrando que
Ve ( x, y)  1Ve1   2 Ve 2  3Ve 3 
3
  Ve
i
i
i 1
A função de energia pode ser colocada na forma:
3
1
We   ( Ve1i 
2 A i 1

3
 Ve  )dS
1
i 1
i
1 3
We  
2 i1
3
 Ve (    )dS)Ve
i
j1
i
j
A

(e)
Onde denominamos Cij  ( i   j )dS
A
O termo Cij(e) pode ser considerado como o acoplamento entre os nós i e j
j
Método dos Elementos Finitos
2- Equações que regem os elementos
A equação da energia armazenada no campo elétrico
1 3
We  
2 i 1
3
 Ve (    )dS)Ve
onde
i
j1

i
j
j
A
Cij( e )  ( i   j )dS
é o acoplamento entre os nós i e j
A
pode ser escrita na forma matricial como: We 
 
1
T
Ve  C ( e ) Ve 
2
(e)
(e)
(e)


C
C
C
11
12
13
Ve1 


onde Ve   Ve 2 , Ve T  Ve 1 Ve 2 Ve 3  e C (e)  C 21( e ) C 22 ( e ) C 23 ( e ) 


 (e)
(e)
(e) 
Ve 3 
C31 C32 C33 
A matriz C(e) é denominada “Matriz dos Coeficientes do Elemento”
Método dos Elementos Finitos
2- Equações que regem os elementos
Cálculo dos Termos da Matriz dos Coeficientes do Elemento
Cálculo de C12(e)
1
(x 2 y3  x 3 y 2 )  ( y 2  y3 )x  (x 3  x 2 ) y 
1 
2A
1
 (x 3 y1  x1y3 )  ( y3  y1 )x  (x1  x 3 ) y 
2 
2A
Precisamos calcular
C12(e)   (1   2 )dS
A


( y2  y3 )a x  ( x3  x2 )a y
1

2A
1
 ( y3  y1 )ax  ( x1  x3 )a y
 2 ( x, x) 
2A
E resolvendo a integral temos que
 1 ( x , y ) 
C12
(e)



1
 ( y2  y3 )( y3  y1 )  ( x3  x2 )(x1  x3 )
2
4A

A
dS
Método dos Elementos Finitos
2- Equações que regem os elementos
1
 ( y2  y3 )( y3  y1 )  ( x3  x2 )( x1  x3 ) 
4A
Os outros coeficientes podem ser calculados de modo semelhante
1
(e)
1
(e)
2
2
 ( y3  y1 ) 2  ( x1  x3 ) 2
C

 ( y2  y3 )  ( x3  x2 )  22
C11 
4A
4A
1
(e)
 ( y1  y2 ) 2  ( x2  x1 ) 2
C33 
4A
1
(e)
(e)
 ( y2  y3 )( y1  y2 )  ( x3  x2 )( x2  x1 ) 
C13  C31 
4A
1
(e)
(e)
 ( y3  y1 )( y1  y2 )  ( x1  x3 )( x2  x1 ) 
C23  C32 
4A
P1  ( y2  y3 ) P2  ( y3  y1 ) P3  ( y1  y2 )
Usando a notação:
Q1  ( x3  x2 ) Q2  ( x1  x3 ) Q3  ( x2  x1 )
Logo: C12
(e)
 C21
(e)

Podemos escrever Cij
(e)

1
 Pi Pj  QiQ j
4A
 e A  1  P2Q3  P3Q2 
2

Método dos Elementos Finitos
2- Equações que regem os elementos
Em resumo: para a obtenção da função potencial no elemento é necessário conhecer as
funções de forma 1, 2 e 3 cujo cálculo depende do valor das coordenadas dos
3
vértices do elemento triangular.
Ve ( x, y)  1Ve1   2 Ve 2  3Ve 3 
Ve1 
Ve   Ve 2 , Ve T  Ve1 Ve 2 Ve 3  e C(e)
Ve 3 
i
i
i 1
1
T
Ve  C ( e ) Ve 
2
C11( e ) C12 ( e ) C13 ( e ) 
 (e)

(e)
(e)
 C 21 C 22 C 23 
 (e)
(e)
(e) 
C
C
C
32
33
 31

O cálculo da energia do elemento pode ser feito We 
Onde
 
  Ve
Como todos os elementos são escritos em função das diferenças entre as
coordenadas podemos definir novas variáveis P e Q como abaixo
Cij
(e)
1
 Pi Pj  QiQ j

4A
 e A  1  P2Q3  P3Q2
2
As novas variáveis satisfazem as
seguintes propriedades:

P1  ( y2  y3 ) P2  ( y3  y1 ) P3  ( y1  y2 )
Q1  ( x3  x2 ) Q2  ( x1  x3 ) Q3  ( x2  x1 )
P1  P2  P3  0  Q1  Q2  Q3
3
C
i 1
(e)
ij
3
  Cij
j 1
(e)
Método dos Elementos Finitos
3- Conexão de todos os elementos
Após a análise de um elemento típico o próximo passo é a conexão de todos
os elementos. A energia associada à conexão de todos os elementos da
malha é o somatório das energia armazenadas nos elementos:
V1 
V 2
N
1
T
 
W
We  V CV onde
V  . 
2
e 1
 
. 
Vn 

Nas equações acima “N” é o número de elementos e “n” é o número de nós.
A matriz [C ] é denominada Matriz de Rigidez Global e representa a conexão
das matrizes dos coeficientes dos elementos individuais. Não confundir [C]
com [C(e)]. Os coeficientes Cij são obtidos observando-se o fato de que os
potenciais devem ser contínuos através dos contornos dos elementos
Método dos Elementos Finitos
3- Conexão de todos os elementos
Considere a malha de três elementos abaixo:
Numeração global: 1,2,3,4,5 (nós da malha, qualquer ordem)
Numeração local: 1,2,3 (vértices do elemento, sempre no sentido antihorário para evitar que a área fique negativa)
Para o elemento 3 a numeração global 3,4,5 corresponde à numeração
local 1,2,3 respectivamente.
Método dos Elementos Finitos
3- Conexão de todos os elementos
C11 C12 C13 C14 C15 
C C C C C 
 21 22 23 24 25 
C  C31 C32 C33 C34 C35 


C
C
C
C
C
 41 42 43 44 45 
C C C C C 
 51 52 53 54 55 
Como a malha tem cinco nós globais a matriz de rigidez global é do tipo 5x5.
Note que a matriz de rigidez global é a superposição de todas a matrizes dos
elementos. Como muitos potenciais são iguais, vários termos da matriz se
superpõem de modo que a mesma fica reduzida a uma matriz de estrutura
nxn, onde n é o número de nós globais.
Para obter o termo C11 notamos que o nó global 1 pertence aos elementos1 e
2. Notamos, também que os nós locais que fazem parte do nó global 1 são
(1)
( 2)
os nós locais 1 de ambos os elementos, logo: C11  C11  C11
Já o coeficiente C14 é calculado lembrando que a fronteira global 1-4 contem
as fronteiras 1-2 do elemento 1 e 1-3 do elemento 2, então:
C14  C12  C13
(1)
( 2)
Método dos Elementos Finitos
3- Conexão de todos os elementos
C11 C12 C13 C14 C15 
C C C C C 
 21 22 23 24 25 
C  C31 C32 C33 C34 C35 


C
C
C
C
C
 41 42 43 44 45 
C C C C C 
 51 52 53 54 55 
Quando não há acoplamento entre os nós referenciados no termo o mesmo é
nulo, como no caso do termo C15, já que não há conexão direta entre os nós 1 e
5. Desse modo a matriz de rigidez global fica como:
(1)
(2)
C11(1)  C11( 2 ) C13 (1) C12 (2)

C12  C13
0
 (1)

(1)
(1)
C33
0
C32
0
C31

 (2)
(2)
( 3)
(2)
(3)
(3) 
C  C 21
0
C 22  C11
C 23  C13
C12 
(2)
(1)
(2)
(3)
(1)
(2)
(3)
(3) 
 (1)
C

C
C
C

C
C

C

C
C
21
31
23
32
31
22
33
33
32


(3)
(3)
(3) 
0
0
C
C
C
21
23
22


A matriz de rigidez global é simétrica (Cij =Cji) e o determinante formado por seus termos é nulo.
Método dos Elementos Finitos
4 - Resolução das equações resultantes
Do cálculo variacional temos que a equação de Laplace é satisfeita no domínio
quando a energia total no mesmo é mínima. Portanto é necessário que as
derivadas de W em relação ao potencial de cada nó sejam zero. Isso permite
estabelecer um sistema de equações que nos levará ao cálculo dos potenciais
dos nós globais V1, V2, V3 ...Vn, onde n é o número de nós globais.
W W W
W


 ......
0
V1 V2 V3
Vn
Método dos Elementos Finitos
4 - Resolução das equações resultantes
Temos que a energia elétrica armazenada em todo o domínio é dada por:
N

1
T
We  V CV
2
e 1
 V1 
V 
 2
Onde a matriz dos potenciais dos nós globais é dada por V    V3 
 
 V4 
 V5 
 
W
Desmembrando a equação da energia temos, para o exemplo tomado com
três elementos e cinco nós globais:
C11 C 22 C33 C 44 C55  V1 
C C C C C   V 
 21 22 23 24 25   2 
1
W  V1 V2 V3 V4 V5  C31 C32 C33 C34 C35  V3 
2

 
C
C
C
C
C
 41 42 43 44 45  V4 
C C C C C   V 
 51 52 53 54 55   5 
Método dos Elementos Finitos
4 - Resolução das equações resultantes
O produto das matrizes [C] e [V] gera uma matriz coluna como abaixo
mostrada:
C11V1  C12 V2  C13 V3  C14 V4  C15 V5 
C V  C V  C V  C V  C V 
22 2
23 3
24 4
25 5 
 21 1
1
W  V1 V2 V3 V4 V5  C31V1  C32 V2  C33 V3  C34 V4  C35 V5 
2


C
V

C
V

C
V

C
V

C
V
42 2
43 3
44 4
45 5 
 41 1
C V  C V  C V  C V  C V 
52 2
53 3
54 4
55 5 
 51 1
W

 C11V12  C12 V1V2  C13V1V3  C14 V1V4  C15 V1V5  C21V1V2  ...  C31V1V3  ...

V1 V1
...  C41V1V4  ...  C51V1V5

W
  2C11V1  C12 V1  C13V2  C14 V4  C15 V5  C21V2  ...  C31V3  ...  C41V4  ...  C51V5
V1

Método dos Elementos Finitos
4 - Resolução das equações resultantes
W
  2C11V1  C12 V1  C13V2  C14 V4  C15 V5  C21V2  ...  C31V3  ...  C41V4  ...  C51V5
V1
W
mas como C12=C21 temos
  2C11V1  2C12 V2  2C13 V3  2C14 V4  2C15 V5 
V1
como
W
 0,
V1
0  V1C11  V2C12  V3C13  V4C14  V5C15
n
W
nos leva a 0   Cik onde K  1,2,3,4....n
Ou no geral:
Vk
i 1
Dese modo obtemos um conjunto de equações que pode ser resolvido para
obter V1, V2, V3, V4, V5

Método dos Elementos Finitos
4 - Resolução das equações resultantes
FIM