Baixe Método dos Elementos Finitos Aplicado ao Eletromagnetismo e outras Notas de estudo em PDF para Eletromagnetismo, somente na Docsity!
Método dos Elementos Finitos Aplicado ao Eletromagnetismo
1. Introdução
Nesta apostila é apresentado os método do método dos elementos finitos de forma suscinta, baseado num exemplo de aplicação ao eletromagnetismo. Na primeira parte, constituida dos ítens 1-5, após uma breve introdução das Equações de Maxwell, é derivada destas a Equação de Poisson para o caso particular da magnetostática. Na segunda parte, constituida do ítem 6, é apresentada a solução desta pelo Método dos Elementos Finitos. Uma vez que Equação de Poisson descreve uma série de fenômenos físicos, a abordagem apresentada pode ser facilmente estendida, por analogia, a outras áreas da física tais como a distribuição de temperaturas, escoamento de fluidos em dutos, etc...
O Método dos Elementos Finitos tem suas origens nos anos 40, tendo sido entretanto vastamente utilizado apenas nos ultimos 20-30 anos, graças aos avanços experimentados pelo computadores. Ele consiste numa adaptação/modificação de métodos conhecidos já no início deste século, como por exemplo o Método de Ritz. Ele é atualmente definido como um Método Matemático para a solução de equações diferenciais parciais tal como a Equação de Poisson e Laplace. Devido às suas características de flexibilidade e estabilidade numérica, ele pode ser facilmente implementado em um sistema computacional, fato que explica a sua grande popularidade nos dias atuais. Um grande impulso para o seu desenvolivmento e aperfeiçoamento foi dado pela indústria aero- espacial, onde o método vem tendo larga aplicação desde os anos 50, sendo utilizado, entre outros para o projeto e análise estruturas complexas de aviões. As principais áreas de aplicação incluem: projeto e análise de estruturas, análise de escoamento de fluidos, distribuição de temperaturas e eletromagnetismo. Em muitos casos práticos, o Método dos Elementos Finitos é a única ferramenta capaz de fornecer uma solução aceitável, ainda que aproximada.
2. Equações de Maxwell
Os fenômenos eletromagnéticos são regidos pelas equações de Maxwell, as quais são dadas abaixo na sua forma diferencial:
∇ × = −
r r E
B
t
∇ × = +
r
r r H D t
∂ S ∂
r B 0 (3)
∇ • =
r D ρ (4)
r E - Intensidade campo elétrico (grandeza vetorial) r H - Intensidade campo magnético (grandeza vetorial) r D - Densidade campo elétrico (indução elétrica, grandeza vetorial) r B - Densidade campo magnético (indução magnética, grandeza vetorial) r S - Densidade superficial de corrente elétrica (grandeza vetorial)
ρ - Densidade volumétrica de cargas (grandeza escalar)
Nas expressões anteriores, ∇ representa o operador Nabla, dado pela seguinte relação em coordenadas cartesianas:
x ∂
i y
j z
k
r r r (5)
r i ,
r j e
r k são os vetores unitários nas direções x, z e y respectivamente. O uso deste operador é muito semelhante ao uso de um vetor; assim, operações ∇ × e ∇ • (rotacional e divergente) podem ser obtidas formalmente a partir do produto vetorial e do produto escalar entre ∇ e o vetor sobre o qual o operador efetua.
Entre as intensidades e as densidades de campo elétrico e magnético existem ainda as seguintes relações, válidas para materiais isotrópicos lineares: r r B = μ ⋅H (6)
D = ε ⋅E
r (7) r r S = σ ⋅E (8)
μ - permeabilidade magnética (Henry/metro)
ε - permissividade elétrica (Faraday/metro)
σ - condutividade elétrica ( 1 / (^) ( Ohm metro⋅ ))
μ , ε , e σ são grandeza escalares (constantes) para materiais lineares isotrópicos.
Para o caso especial de fenômenos estáticos as equações (1)-(4) se reduzem à seguinte forma:
∇ × =
r E 0 (9)
∇ × =
r r H S (10)
r B 0 (11)
r E ρ (12)
Por outro lado, as equações (6)-(8) seguem válidas. Neste caso existe uma independência entre o campo elétrico e o campo magnético. Para o estudo de campos magnéticos estáticos, os quais serão abordados mais adiante pelo método dos elementos finitos, necessita-se considerar apenas as seguintes equações:
∇ × =
r r H S (13)
r B 0 (14) r r B = μ ⋅H (15)
As equações acima são as equações fundamentais da magnetostática. Para a derivação do Método dos Elementos Finitos será utilizado ainda o Teorema de Green no Plano, o qual estabelece que para duas funções u(x,y) e v(x,y) vale a relação:
do potencial. Existe, desta forma, um determinado grau de liberdade de escolha para o potencia vetor. A fim de simplificar a expressão acima, pode-se optar por um potencial vetor que atenda a seguinte condição:
∇ • =
r A 0 (24)
Com esta condição a equação (23) se torna:
∇ × ∇ × = − ⋅
r r A μ S (25)
A equação acima é conhecida como Equação de Poisson no espaço, ela descreve não apenas os fenômenos eletromagnéticos, mas também muitos outros outros, tais como como a transmissão de calor, distribuição de temperaturas, escoamento de fluidos, etc...
Para caso o especial em que
r S é igual a zero, a equação assume a uma forma conhecida como Equação de Laplace :
∇ × ∇ × =
r A 0 (26)
A equação de Laplace pode portanto ser considerado um caso especial da equação de Poisson, sendo que os métodos numéricos em geral procuram a solução da equação de Poisson.
Conforme mostrado, por meio da introdução do potencial vetor chega-se a uma única equação que representa as equações (13) a (15). O processo de solução visa assim
determinar o vetor
r A x y z( , , ) , por meio do qual as grandezas eletromagnéticas de interesse
podem ser obtidas.
Para o caso particular em que o campo não variar segundo uma das variáveis (em geral a variável z), obtém-se um caso bi-dimensional. O vetor densidade de corrente é perpendicular ao plano em que o campo é descrito, conforme estabelece a equação (17). Sendo o campo dependente apenas das direções x e y, a densidade de corrente terá apenas componentes segundo o eixo z. O potencial vetor terá igualmente apenas componentes segundo o eixo z: r r S = S k⋅ (27) r r A = A ⋅ k (28)
A e S são as amplitudes da densidade de campo e do potencial vetor respectivamente, ambos funções escalares de x e y.
No caso bi-dimensional, a equação de Poisson para a magnetostática toma a seguinte forma:
∇ × ∇ × A = − μ⋅S
A última equação em geral é escrita da seguinte forma:
∇ 2 A = − μ ⋅S (29)
O operador ∇^2 é chamado de Laplaciano. Nas regiões do domínio em estudo em que S é zero, vale a equação de Laplace :
∇^2 A = 0 (30)
As expressão (28)-(29) representam uma equação diferencial parcial de segunda ordem, as quais descrevem problemas conhecidos como "Problemas de Valores de Contorno" ou ainda "Problemas de Potencial". Escrita de forma de uma forma mais explicita a equação
(28) toma a forma:
∂ ∂
μ
2 2
2 2
A
x
A
y
+ = − ⋅ S (31)
A determinação do campo nas direções x e y se reduz à determinação do potencial A x y( , )
segundo estas direções.
A solução da equação (30) na forma analítica só é possível para casos com geometrias muito simples e sob certas aproximações, as quais nem sempre são justificáveis na prática, fazendo com que a solução analítica, embora possível, não possua um valor inquestionável para a grande maioria dos casos práticos. A vantagem da solução analítica é, todavia, o fato de que a influência dos parâmetros físicos e geométricos aparece explicíta na solução, facilitando a sua análise. Por outro lado, o Método dos Elementos Finitos ( abreviadamente MEF) é um método de solução numérica da equação (30) que pode ser aplicado para qualquer espécie de domínio, ele não fornece, entretanto, uma solução na forma analítica. Uma análise da influência dos parâmetros físicos e geométricos precisa ser obtida por variações discretas de um grande número de casos semelhantes, fato que se torna cada vez mais irrelevante a medida que computadores cada vez mais potentes e sistemas de cálculo por elementos finitos mais sofisticados vão surgindo. Outra vantagem importante do MEF é a possibilidade de tratamento de casos não-lineares, o que
r
S
Fronteira externa Γ
∇ × =
r
H 0
∇ × =
r r
H S
Ferro
Ar
A = 0
A
n
= ∇ A= 0
Ferro
Linha de Simetria
Figura 1 - Exemplo da magnetostática: distribuição do campo magnético é regido pela equação de Poisson.
Baseado no exposto até aqui, a descrição dos problemas de potencial pode ser resumida na forma mostrada no quadro abaixo. Eles se compõem sempre da equação de Poisson, a qual é válida para o interior do domínio Ω e das condições de contorno impostas no contorno externo Γ. A condição de Dirichlet é imposta na parte do domínio designado por Γ 1 , enquanto que a condição de Neumann é imposta no restante do domínio designado
por Γ 2. O quadro 1 sumariza a descrição de um problema de potencial bi-dimensionnal.
6. Método dos Elementos Finitos:
Uma vez definido o problema de valor de contorno que se deseja resolver, pode-se passar para a solução do mesmo por meio do Método dos Elementos Finitos. O problema definido no quadro acima é um problema do tipo contínuo, uma vez que todos os pontos do domínio são incluidos tanto na descrição quanto na solução do problema. O MEF transforma este domínio contínuo num domínio discreto, onde a solução é conhecida em pontos discretos do dominio de cálculo, por exemplo em pontos de união de uma malha triangular (nós).
O MEF pode ser derivado basicamente por dois caminhos, conforme mostra a figura 4. No caminho mostrado à esquerda na figura as equações discretas são obtidas por meio da minimização de funções de energia (ou outra função equivalente), utilizando-se de princípios variacionais. O método variacional utilizado neste caso é o método de Ritz^2 com a particularidade de que as funções a serem minimizadas valerem apenas em subdomínios discreto (elementos da malha), por exemplo triângulos.
O caminho mais rápido para a derivação, e que será seguido aqui, é a aplicação de métodos residuais para a obtenção das equações discretas. Neste caso, o método empregado na derivação é o método de Galerkin, também com a particularidade de que as equações de aproximação ser referem a subdomínios discretos (elementos da malha).
6.1 Aproximação de Funções pelo Método de Galerkin
Multiplicando-se a equação (28) por uma função g x y( , ) (ainda a definir mais
precisamente) chamada de função teste (também chamada função peso ou função de amostragem) obtém-se a seguinte expressão:
∇ 2 A ⋅ g = − μ ⋅ S g⋅ (34)
Integrando-se ambos os lados da equação anterior sobre o domínio Ω resulta:
∇ ⋅ ⋅ = − ⋅ ⋅ ⋅ ∫∫ ∫∫
(^2) A g d S g d
Ω Ω
Ω μ Ω (35)
Usando-se o teorema de Green para o plano, enunciado anteriormente, chega-se à forma seguinte:
∇ ⋅ ∇ ⋅ − ∇ ⋅ ⋅ = ⋅ ⋅ ⋅ ∫∫ ∫ ∫∫
A g d A g d S g d Ω Γ Ω
Ω Γ μ Ω (36)
∇ ⋅ ∇ ⋅ − ∇ ⋅ ⋅ − ∇ ⋅ ⋅ = ⋅ ⋅ ⋅ ∫∫ ∫ ∫ ∫∫
A g d A g d A g d S g d Ω Γ Γ Ω
Ω Γ Γ Ω
1 2
μ (37)
(^2) Este método foi proposto originalmente num artigo em 1909, tendo sido desde então bastante utilizado
sobretudo em problemas de otimização e solução de equações diferenciais parciais.
Considerand0-se as condições de contorno definidas anteriormente, pode-se simplificar a última expressão para a que segue:
∂ ∂
A n
= ∇A = 0 sobre Γ 2
∇ = − ⋅
2 A μ S
Problema Discreto
Métodos Variacionais
Método de Ritz
Problema Contínuo
Leis Físicas
(eq. Maxwell)
Métodos Residuais
Método de Galerkin
Solução Discreta
(aproximada)
Discretização do Domínio Discretização do Domínio
K (^) [ ]A b G
[ ] (^) ⋅ =[ ]
Figura 4 - Representação esquemática da obtenção de uma solução aproximada pelo do Método dos Elementos Finitos
elemento. A forma mais simples de elemento utilizado na prática é o triângulo, ele pode aproximar domínios de formas quaisquer com boa precisão. Em cada elemento são definidos pontos característicos, nos quais a solução será determinada; a escolha das grandezas a determinar também determina parcialmente o tipo de função de aproximação que pode ser utilizado. No caso mais simples de elementos triangulares são escolhidos os três vértices dos mesmo como pontos característicos, os quais são chamados de nós. A figura 5 mostra a divisão de um domínio bastante simples em 6 subdomínios menores. Em sistemas de cálculo por elementos finitos modernos a discretização da região em estudo é feita de forma automática, a partir dos dados geométricos do domínio. O número de elementos utilizados na prática depende da natureza do domínio em estudo e do comportamento particular da solução. Em geral a precisão da solução aumenta com o número de elementos utilizados, havendo no entanto um limite para o número de elementos, a partir do qual os erros de arrendodamento se acumulam de tal forma que um aumento do número de elementos não traz uma melhora na precisão. Por outro lado, a precisão da solução depende também muito fortemente do tipo da função de aproximação utilizada (função linear, função quadrática, função cúbica, exponencial, etc...), do tipo de elemento utilizado (triangular, retangular, etc...).
O domínio discretizado é comumente designado por malha , sendo a mesma caraterizada pelo número de nós N e elementos M (subdomínios). A figura 6 mostra novamente o domínio da figura 1 após ter sido discretizado. Os elementos são em parte triângulos e em parte retângulos, havendo assim necessidade de um equacionamento diferente para cada um deles.
6.3 Equação Discretas
Uma vez que o domínio em estudo foi discretizado, pode-se obter as equações
Figura 5 - Estrutura mostrada na figura 1 após ter sido discretizada (malha resultante).
aproximadas válidas para cada subdomínio designado por Ωe , chamadas de equações
discretas. Considerando-se a subdivisão do domínio original em M elementos, pode-se inicialmente escrever a equação (40) na seguinte forma:
∇ ⋅ ∇ ⋅ = ⋅ ⋅ ⋅ ∑ (^) ∫∫ ∑ ∫∫ = =
A g d S g d
e e
e e
M e e
M
Ω Ω
Ω Ω 1 1
μ (41)
A última equação estabelece simplesmente que a integral foi divida numa soma de integrais parciais, abrangendo todos os M elementos que compõem o domínio. A fim de avaliar as integrais na equação (41), é preciso determinar a forma da solução aproximada e da função de teste g x y( , ). A fim de tornar a análise mais simples serão admitidas sempre funções lineares e elementos triangulares. Desta forma, de acordo com o Método de Galerkin, pode-se escrever a função de aproximação da seguinte forma:
A A x y ( ) A x y( ) R (^) j ( x y) j
N = ≅ = ⋅ j =
, , (^) ∑ , 1
φ (42)
Esta equação vale em princípio para todo o domínio em estudo, havendo N parâmetros R (^) j
a serem determinados para as N funções, a fim de que a solução seja conhecida. As funções φ (^) j ( ,x z ) são também conhecidas como funções de base e no caso mais simples são
funções lineares, as quais possuem em princípio um valor diferente de zero em todos os subdomínios. A solução aproximada é dada por uma combinação linear destas funções. Um procedimento semelhante é utilizado ao se aproximar funções periódicas por meio de séries de Fourier. A solução na forma (42) pode ser considerada uma generalização do método de Fourier, uma vez que além de funções senoidais, pode-se utilizar polinômios, exponenciais, funções lineares, quadráticas, etc...
O Método dos Elementos Finitos caracteriza-se pela forma particular como as funções φ (^) j ( ,x z ) e os parâmetros R (^) j são definidos:
- Os parâmetros R (^) j são escolhidos como sendo os potenciais nos nós da malha chamados de A (^) j , ou seja R (^) j = Aj. Este procedimento é muito semelhante ao que se utiliza no Método de Interpolação de Lagrange;
- As funções φj ( ,x y ) são definidas de tal forma que em cada subdomínio apenas algumas delas são diferentes de zero, a solução dentro do elemento uma torna-se combinação linear destas poucas funções, sendo assim mais fácil de ser determinada. Considerando-se, por exemplo, elementos triangulares em duas dimensões, apenas 3 delas tem valor diferente de zero em cada elemento.
A definição das funções na forma acima traz inúmeras vantagens práticas. Pode-se utilizar a mesma forma para todas as funções, deixando-se os seus parâmetros como incógnitas. Assim, pode-se definir a solução para apenas um triângulo e em seguida estendê-la para todos os demais, simplificando enormemente o equacionamento do problema. Considerando-se o triângulo genérico mostrado na figura 6, obtém-se para a solução aproximada:
A (^) j d S d j
j i e i e = e e
∑ ⋅^ ∫∫ ∇φ^ ⋅ ∇φ ⋅^ =^ ⋅^ ⋅^ ∫∫ ⋅ 1
3
Ω Ω
Ω μ φ Ω i = 1 2 3, , (54)
Considerando-se todo o domínio pode-se ainda estabelecer uma equação global:
e
M j j
A (^) j i d (^) e S (^) i d e = = e e
∑ ∑ ⋅^ ∫∫ ∇φ^ ⋅ ∇φ ⋅^ =^ ⋅^ ⋅^ ∫∫ ⋅ 1 1
3
Ω Ω
Ω μ φ Ω (55)
Dada a simplicidade das funções φj ( ,x y ) , a avaliação da integral se torna bastante fácil,
resultando tanto no lado esquerdo como no lado direito valores constantes.
Considerando-se a equação (54) e variando-se os índices i e j de 1 a 3 obtém-se um sistema local de equações, onde os potenciais dos nós aparecem como incógnitas:
A (^) j K b j
ij
e i =
∑ ⋅^ = 1
3 i = 1 2 3, , (56)
K (^) ije^ j i d e
e
= ∇φ ⋅ ∇φ ⋅ ∫∫ Ω
Ω i = 1 2 3, , j = 1 2 3, , (57)
bi e^ S (^) i d e
e
= ⋅ ⋅ ⋅ ∫∫
μ φ Ω
Ω i = 1 2 3, , (58)
Na forma matricial pode-se ainda estabelecer:
K K K K K K K K K
A A A
b b b
e e e e e e e e e
e e e
11 12 13 12 22 23 13 23 33
1 2 3
1 2 3
⋅
=
Considerando-se todos os elementos, obtém-se um sistema local semelhante ao (59) para cada elemento, onde a matriz de coeficientes locais é simétrica. A fim de se obter a matriz global do sistema, deve-se somar todos os sistemas locais, de acordo com a equação (55), obtendo-se desta forma um sistema único de equações chamado de sistema global. Os índices 1 , 2 e 3 nas equação locais de cada elemento possuem uma correspondência com a numeração global dos nós, a qual vai de 1 a M. Na montagem da matriz global deve-se avaliar a contribuição de cada um dos elementos e adicioná-los à matriz global. Por exemplo, supondo-se que o elemento 1 de uma malha contendo 10 nós esteja sendo considerado e que exista a seguinte correspondência entre os nós locais e globais:
1 ⇒ 5 , 2 ⇒ 7 , 3 ⇒ 9
Ao ser processado o elemento em questão, a sua contribuição será adicionada à matriz global da forma esquematizada abaixo.
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
11
1 12
1 13
1
12
1 22
1 23
1
13
1 23
1 33
1
1 2 3 4 5 6 7 8 9
10
K K K
K K K
K K K
A A A A A A A A A A
⋅
=
0 0 0 0 0 0 0
5
1
7
1
9
1
b
b
b
De acordo com a equação (60) as contribuições locais são colocadas nas linhas e colunas
da matriz global dadas pelo número dos nós globais; K 121 por exemplo será colocado na
linha 5 e coluna 7, K 131 na linha 5 e coluna 9, etc...
Supondo ainda que para o elemento 2 exista a seguinte correspondência entre os nós globais e locais:
1 ⇒ 3 , 2 ⇒ 5 , 3 ⇒ 10
Ao inserir-se o elemento 2, procede-se analogamente ao feito para o elemento 1, resultando então para a matriz global:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
11
2 12
2 13
2
12
2 11
1 22
2 12
1 13
1 23
2
12
1 22
1 23
1
13
1 23
1 33
1
13
2 23
2 33
2
1 2 3 4
K K K
K K K K K K
K K K
K K K K K K
A A A A
⋅ 5 6 7 8 9 10
3
2
5
1 5
2
7
1
9
1
10
2
0 0
0
0
0
A A A A A
b
b b
b
b b
=
Repetindo-se o procedimento ilustrado para todos os elementos da malha obtém-se um sistema de equações lineares do tipo:
[ A^ ] ⋅^ [^ K]^ g=^ [b^ ] (62)
Observa-se que a matriz de coeficientes (^) [ K] (^) g é simétrica e esparsa, com poucos elementos
fora da diagonal principal. O número de elementos fora da diagonal em cada linha equivale ao número de nós conectado ao nó relativo àquela linha. Estas características em geral são aproveitada em termos de armazenamento e solução do sistema.
An introduction to the mathematical theory of finite elements : J. T. Oden, J. N. Reddy; John Wiley & Sons, 1976.
The mathematical theory of finite element methods : Susanne C. Brenner, L. Ridgway Scott ; Springer Verlag, 1994.
Computer-aided analysis and design of electromagnetic devices : S. R. H. Holle ; Elsevier, 1989.
Finite Elements for electrical engineering : P. Silvester, R. L. Ferrari;Cambridge University Press, 1996.
The finite element method in electromagnetics , J. Jin, John Wiley & Sons, 1993.
