144
IEEE LATIN AMERICA TRANSACTIONS, VOL. 6, NO. 2, JUNE 2008
Utilização do Método FDTD em Coordenadas
Não Ortogonais para Simular Estruturas de
Aterramento Elétrico
Rodrigo M. S. de Oliveira e Carlos Leonidas da S. S. Sobrinho
Resumo — Neste trabalho, o método FDTD em coordenadas
não ortogonais (LN-FDTD) foi implementado para a análise de
estruturas de aterramento com geometrias coincidentes ou não
com o sistema de coordenadas cartesiano. O método soluciona as
equações de Maxwell no domínio do tempo, permitindo a obtenção de dados a respeito da resposta transitória e de regime das
estruturas. O software desenvolvido foi testado e validado para o
caso de um sistema de aterramento em forma de prato e para um
eletrodo horizontal não alinhado com qualquer coordenada do
sistema cartesiano.
Palavras-chave — Método FDTD, Sistemas de Aterramento,
Equações de Maxwell, Sistema de Coordenadas Não Ortogonais.
S
I. INTRODUÇÃO
istemas de aterramento elétrico são indispensáveis para a
proteção e manutenção do funcionamento de sistemas
elétricos, bem como para a segurança de pessoas e de animais.
De uma forma geral, esses sistemas consistem em dispositivos
metálicos enterrados no solo que são conectados a algum
dispositivo que se quer proteger eletromagneticamente. O
principal objetivo é escoar correntes indesejadas para o solo
da maneira mais eficiente possível.
Existem diversas geometrias de sistemas de aterramento
possíveis, tais como eletrodos horizontais e verticais, eletrodos circulares, esféricos, malhas de terra, entre várias outras.
Deve-se ressaltar que o desempenho do sistema de aterramento depende também das características elétricas do solo (além,
é claro, da geometria citada anteriormente).
Para a obtenção da resposta completa de um dado sistema
de aterramento (transitório e regime), modelos eletromagnéticos que solucionem as equações de Maxwell no domínio do
tempo são desejáveis, de forma a representar com fidelidade a
realidade. Dessa forma, é possível desenvolver estruturas de
aterramento que respondam de forma próxima ao desejável
para um dado sistema a custos reduzidos e com excelente
precisão.
As primeiras análises de malhas de terra foram feitas por
Bewley [1][2] em 1934, trabalho no qual foram apresentadas
investigações teóricas e experimentais de cabos contrapeso.
Em 1943, eletrodos verticais foram investigados teoricamente
Artigo recebido em 06 de maio de 2007. Laboratório de Análises Numéricas em Eletromagnetismo (LANE), Universidade Federal do Pará, caixa
postal 8619, CEP 66075-907, Brasil; e-mails: [email protected], [email protected]. Este trabalho foi financiado pelo CNPq e Eletronorte.
por Bellaschi e Armingtom [3]. Mais tarde, Sunde [4] aplica
as Equações de Maxwell para demonstrar expressões para
cálculo de respostas DC de várias estruturas de aterramento.
A partir da década de 1980, a disponibilidade de computadores estimulou o uso de métodos numéricos na análise de
estruturas de aterramento. Entre esses métodos, podem ser
destacados: métodos de circuitos equivalentes [5]-[9], o método dos momentos [10]-[13], o método dos Elementos Finitos
[14][15] e mais recentemente o método das Diferenças Finitas
no Domínio do Tempo (FDTD) [16]-[18].
Em 2001, Tanabe publicou o artigo [17] no qual é demonstrado que o algoritmo de Yee [16] (o método FDTD) pode ser
aplicado para analisar estruturas de aterramento coincidentes
com o sistema cartesiano de coordenadas. Comparações dos
resultados gerados com o método FDTD com medidas experimentais mostraram que o método FDTD pode ser aplicado
com grande confiabilidade para esse tipo de problema, tanto
para respostas de regime quanto para transitórios eletromagnéticos. Eletrodos auxiliares foram utilizados para medir a tensão e para aplicar a corrente, como ilustrado pela Fig. 1.
Circuito de Corrente
Circuito de Tensão
V
I + -
Fig. 1. O Modelo experimental de Tanabe.
Em [18], Tuma et al. propõem um modelo computacional
baseado no algoritmo de Yee que elimina o eletrodo de potencial (o potencial é calculado em termos do campo elétrico) e o
eletrodo de corrente (que é introduzido na região absorvente
U-PML [18][19], simulando um canal de descarga de corrente
natural), como ilustrado na Fig. 2. A importância desse modelo é ressaltada pela eliminação do acoplamento eletromagnético entre o sistema de aterramento analisado e o circuito de
medição utilizado em [17].
MELO E SILVA DE OLIVEIRA AND LEONIDAS DA SILVA SOUZA : EMPLOYMENT OF THE FDTD
145
∂H x 1 ⎛ ∂E y ∂E z ⎞
⎟,
−
= ⎜⎜
∂t
∂y ⎟⎠
μ ⎝ ∂z
∂H y 1 ⎛ ∂E z ∂E x ⎞
= ⎜
−
⎟,
∂z ⎠
μ ⎝ ∂x
∂t
∂H z 1 ⎛ ∂E x ∂E y ⎞
⎟,
−
= ⎜⎜
∂x ⎟⎠
∂t
μ ⎝ ∂y
(3)
(4)
(5)
e
⎞
∂E x 1 ⎛ ∂H z ∂H y
−
− σE x ⎟⎟,
= ⎜⎜
∂t
∂z
ε ⎝ ∂y
⎠
∂E y 1 ⎛ ∂H x ∂H z
⎞
−
− σE y ⎟,
= ⎜
∂t
∂x
ε ⎝ ∂z
⎠
(6)
(7)
⎞
∂E z 1 ⎛ ∂H y ∂H x
−
− σE z ⎟⎟,
= ⎜⎜
∂t
∂y
ε ⎝ ∂x
⎠
Fig. 2. Modelo Computacional (Tuma et al.). O eletrodo horizontal representa a estrutura de terra. (a) Representação tridimensional, (b) Um corte bidimensional detalhando o caminho de integração do campo elétrico.
Dessa forma, neste trabalho, é proposta a aplicação do método FDTD em Coordenadas não ortogonais (LN-FDTD) [20]
para reduzir a restrição relativa à modelagem da geometria
dos sistemas de aterramento não compatíveis com o sistema
cartesiano de coordenadas. Inicialmente, o algoritmo de Yee é
apresentado no sistema de coordenadas retangulares. Posteriormente, as noções matemáticas relativas ao sistema de coordenadas não ortogonal são apresentadas, seguidas da sua aplicação ao algoritmo de Yee. São apresentados então os resultados obtidos e as considerações finais.
II. OS MÉTODOS FDTD E LN-FDTD
A. O algoritmo de Yee
As equações de Maxwell que governam a propagação de
ondas eletromagnéticas em meios isotrópicos, não dispersivos,
e com perda elétrica, em sua forma diferencial no domínio do
tempo, são dadas por:
r
r r
∂H
(1)
∇ × E = −μ
∂t
e
r
r r
∂E r
(2)
∇× H = ε
+ J,
∂t
r
r
nas quais E é o vetor intensidade de campo elétrico (V/m), H
é o vetor intensidade de campo magnético (A/m), ε e μ são,
respectivamente, permissividade elétrica (farad/m) e
r
permeabilidade magnética (henry/m) e J , o vetor densidade
2
de corrente elétrica de condução (A/m ). As equações (1) e
(2), quando expandidas, geram as respectivas equações
escalares
(8)
sendo Ex, Ey, Ez e Hx, Hy, Hz as componentes dos campos
r
r
elétrico E e magnético H ,respectivamente. Essas componentes
são, de uma forma geral, funções do tempo t e das três
coordenadas cartesianas x, y e z.
A equação (1), lei de Faraday, mostra que, quando há
r
r
variação no tempo do vetor B = μH (vetor densidade de fluxo
magnético), surgem componentes de campo elétrico
circulando em torno da direção dessa variação. Dualmente, a
lei de Ampère mostra que, quando há variação no tempo do
r
r
vetor D = εE numa certa direção, essas causam circulação de
campo magnético em torno da mesma direção. Essas duas
observações sobre as equações (1) e (2) são os pilares de toda
a teoria de propagação de ondas eletromagnéticas no universo
macroscópico [21], nas quais Yee se baseou ao definir seu
esquema de distribuição espacial e temporal das componentes
de campo de seu algorítimo numérico. A disposição espacial
dessas componentes está ilustrada na Figura 3: a célula
ortogonal de Yee.
z
x
y
Ez
Hx
∆z
Hy
Ey
(i,j,k )
Hz
Ex
∆x
∆y
r
r
Fig. 3. Distribuição espacial das componentes dos campos E e H para a
célula de Yee (i,j,k).
Dessa forma, as versões discretas das equações escalares
r
r
para as componentes x dos vetores E e H são dadas por
146
IEEE LATIN AMERICA TRANSACTIONS, VOL. 6, NO. 2, JUNE 2008
H
n+
1
2
1
1
x (i , j + ,k + )
2
2
=H
n−
1
2
1
1
x (i , j + ,k + )
2
2
+
⎡ En 1
− En 1 ⎤
y (i , j + ,k ) ⎥
Δ t ⎢ y (i , j + 2 ,k +1)
2
+ ⎢
⎥
Δz
μ
⎢⎣
⎥⎦
(9)
n
⎡ En
⎤
1 −E
1
z (i , j ,k + ) ⎥
Δ t ⎢ z (i , j +1, k + 2 )
2
− ⎢
⎥,
Δy
μ
⎢⎣
⎥⎦
e
Δ ⎞
⎛
⎜1−σ t ⎟
2
ε ⎟+
⎜
E 1
=E 1
x ( i + , j ,k )
x ( i + , j ,k ) ⎜
Δt ⎟
2
2
1
σ
+
⎟
⎜
2ε ⎠
⎝
1
n+
⎤
⎡ n+ 12
2
⎢ H z ( i + 1 , j + 1 ,k ) − H z ( i + 1 , j − 1 , k ) ⎥
Δt
2
2
2
2
⎥
⎢
+
⎥
Δ ⎞⎢
Δy
⎛
ε ⎜1 + σ t ⎟ ⎢
⎥
2ε ⎠ ⎣
⎝
⎦
n +1
Fig. 4. O sistema não ortogonal de coordenadas e os vetores unitários.
n
r
Um conjunto de três vetores complementares a l pode ser
definido de tal forma que cada um é normal a dois vetores
unitários simultaneamente. Esse conjunto é definido
matematicamente por (13).
(10)
1
n+
⎤
⎡ n+ 12
2
⎢ H y ( i + 1 , j ,k + 1 ) − H y ( i + 1 , j ,k − 1 ) ⎥
Δt
2 ⎥
2
2
2
⎢
,
−
⎥
Δt ⎞ ⎢
Δz
⎛
ε ⎜1 + σ ⎟ ⎢
⎥
2ε ⎠ ⎣
⎝
⎦
de tal forma que a posição indexada por (i, j , k ) é relacionada
a ( x, y , z ) por ( x = iΔ x , y = jΔ y e z = kΔ z ) e o tempo t é
relacionado ao índice n por t = nΔ t .
Para o esquema numérico ser estável, a condição de Courant
é utilizada para determinar o passo temporal Δ t . Tal condição
r r
r r
r r
r1 a2 × a3 r 2 a3 × a1 r 3 a1 × a2
,a =
,a =
.
a =
g
g
g
(13)
Os vetores definidos por (13) formam a chamada base
recíproca e são denominados de vetores recíprocos. Na
Equação (13), g representa o determinante do tensor métrico
[g] dado por
r r
⎛ a1 ⋅ a1
⎜r r
[ g ] = ⎜ a2 ⋅ a1
⎜ ar ⋅ ar
⎝ 3 1
r r
a1 ⋅ a2
r r
a2 ⋅ a2
r r
a3 ⋅ a2
r r
a1 ⋅ a3 ⎞ ⎛ g11
r r ⎟ ⎜
a2 ⋅ a3 ⎟ = ⎜ g 21
r r
a3 ⋅ a3 ⎟⎠ ⎜⎝ g 31
g12
g 22
g 32
g13 ⎞
⎟
g 23 ⎟.
g 33 ⎟⎠
(14)
é expressa por
1
Δt ≤
C
1
1
1
+ 2+ 2
2
Δx
Δy
Δz
De (14), nota-se que
,
Desenvolvendo o determinante de [g], Eq. (14), observa-se
que
r r r
(16)
g = a1 ⋅ (a2 × a3 ).
na qual C é a velocidade da luz no vácuo.
B. O sistema não ortogonal de coordenadas
A Equação (16) mostra que
Considerando um sistema de coordenadas curvilíneas
formado por (u1 , u 2 , u 3 ) , como ilustrado pela Fig. 4, o vetor
a3
comprimento diferencial de r , representado por dr , é dado
por
r
3
3
r
∂r
r
dr = ∑ l du l = ∑al du l .
l =1 ∂u
l =1
g é o volume do hexaedro
r r
r
formado pelos vetores a1 , a2 e a3 (veja a Figura 5).
r
r
(15)
g ll =| al |2 .
(11)
a1
(12)
a2
r
Os vetores al são chamados vetores unitários e os mesmos
formam a chamada base unitária. Como pode ser observado,
r
l
os vetores al são tangentes aos eixos u e podem ser escritos
como funções de x, y e z.
Fig. 5. Vetores unitários em uma célula não ortogonal (malha estruturada).
Usando a definição matemática (13), pode ser observado
que:
r r
a l ⋅ am = δ l ,m ,
na qual
δ l,m
é a Função Delta de Konecker.
(17)
MELO E SILVA DE OLIVEIRA AND LEONIDAS DA SILVA SOUZA : EMPLOYMENT OF THE FDTD
Uma definição métrica similar à introduzida pela Eq. (14)
pode ser feita para os vetores recíprocos da seguinte maneira:
r r
g lm = a l ⋅ a m .
147
C. O Método FDTD em Coordenadas Não-Ortogonais...
Para se obterem as componentes contravariantes do campo
(18)
r
H , representadas por h l ( l = 1,2,3 ), a equação vetorial (1)
Assim, com esses conceitos em mente, é possível expandir
pode ser decomposta em vetores paralelos aos vetores a l
usando-se a operação produto escalar. Dessa forma, sem perda
de generalidade, para l = 1 , tem-se:
r
um vetor V usando vetores unitários e recíprocos, como
mostrado em (19) (veja Eq.(12)),
3
r 3 r
r
V = ∑v l al = ∑vl a l ,
l =1
r
r
⎛
∂H
⎜⎜ − μ
∂t
⎝
(19)
⎞ r1 r r r 1
⎟⎟ ⋅ a = ∇ × E ⋅ a .
⎠
(
l =1
v l e vl são chamados, respectivamente, de l 'ésima
componente contravariante e l 'ésima componente co-variante
r
de V .
Para calcular a l 'ésima componente contravariante do vetor
r
r r
V , o produto escalar V ⋅ a l deve ser calculado. Dessa forma,
Empregando as Equações (20) e (13), obtém-se:
na qual
usando as equações (17) e (19), verifica-se que:
⎛
∂h1 ⎞ r r
⎜⎜ − μ
⎟ = ∇× E
∂t ⎟⎠
⎝
(
r r
é obtida por:
r r ⎛ 3 r ⎞ r
V ⋅ al = ⎜ ∑vl a l ⎟ ⋅ al = vl .
⎝ l =1
⎠
(21)
Pode-se obter uma relação para calcular as componentes cor
variantes de V a partir de suas componentes contravariantes
(e vice-versa). Se for calculado o produto escalar
r r
3
r r
V ⋅ am = (∑ v l al ) ⋅ am ( veja (21) e (14) ), obtém-se
l =1
3
vm = ∑g lm v l .
(22)
(
r r
r
× a3 ⎞⎟
.
g ⎟⎠
r r r r
r r r r
)(
) (
)(
)
r
Da Equação (19), nota-se que ∇ = ∑3 (∂/∂u l )ar l . Assim,
l =1
empregando (17)
que a primeira
contravariante de
e (21) ainda para o lado direto, verifica-se
equação para a primeira componente
r
H é dada por:
∂h1
∂e ⎞
⎛ ∂e
= −⎜ 32 − 23 ⎟ /( μ g ).
∂t
∂
∂
u
u ⎠
⎝
(25)
Desenvolvendo o mesmo procedimento para as demais
r
componentes de H , verifica-se que:
Seguindo procedimento similar, chega-se a:
∂h 2
⎛ ∂e ∂e ⎞
= −⎜ 13 − 31 ⎟ /( μ g )
∂t
⎝ ∂u ∂u ⎠
(26)
∂h 3
∂e ⎞
⎛ ∂e
= −⎜ 21 − 12 ⎟ /( μ g ).
∂t
⎝ ∂u ∂u ⎠
(27)
e que
v = ∑g vl .
lm
l =1
Finalmente, deve-se notar que os vetores
(23)
r
r
al e a l não
apresentam necessariamente módulo unitário. Isso significa
que as componentes
⎝
2
r r r r
r r r r
⎛
∂h1 ⎞ ∇ ⋅ a2 E ⋅ a3 − ∇ ⋅ a3 E ⋅ a2
⎜⎜ − μ
⎟=
.
∂t ⎟⎠
g
⎝
l =1
3
r
)⋅ ⎛⎜⎜ a
Utilizando a identidade (A× B) ⋅ (C× D) = (A⋅C)(B⋅ D) −(A⋅ D)(B⋅C)
no lado direito da equação anterior, chega-se a:
r r ⎛ 3 r⎞ r
V ⋅ a l = ⎜ ∑v l al ⎟ ⋅ a l = v l .
(20)
⎝ l =1
⎠
r
De forma similar, a l 'ésima componente co-variante de V
m
)
r
Para o campo E , chega-se a:
vl e v l apresentam magnitudes
dependentes do sistema de coordenadas. No entanto,
r
observando as expansões definidas por (19) e que g ll =| al |2
r
e g ll =| a l |2 , os valores fisicos das componentes Vl e V l são
calculados por:
Vl = vl g ll e V l = v l g ll .
(24)
ε
∂e1
⎛ ∂h ∂h ⎞
+ σe1 = ⎜ 32 − 23 ⎟ / g ,
∂t
∂u ⎠
⎝ ∂u
(28)
ε
∂e 2
⎛ ∂h ∂h ⎞
+ σe 2 = ⎜ 13 − 31 ⎟ / g
∂t
⎝ ∂u ∂u ⎠
(29)
ε
∂h ⎞
∂e 3
⎛ ∂h
+ σe3 = ⎜ 21 − 12 ⎟ / g .
∂t
⎝ ∂u ∂u ⎠
(30)
e
148
IEEE LATIN AMERICA TRANSACTIONS, VOL. 6, NO. 2, JUNE 2008
As equações (25)-(30) podem ser aproximadas de forma
usual por diferenças centradas. Desse modo, são obtidas as
seguintes equações de atualização:
1
1 ( n+ )
2
( i , j ,k )
=h
−
Δt ⎡e
⎢
μ g ⎢⎣
−e
Δu
(n)
3( i , j , k )
2
−
e
(n)
2( i , j , k +1)
−e
Δu
( n)
2( i , j , k )
3
⎤
⎥,
⎥⎦
(31)
−
1
2 ( n− )
2
( i , j ,k )
=h
(n)
(n)
(n)
(n)
Δ t ⎡ e1(i , j ,k +1) − e1(i , j ,k ) e3(i +1, j ,k ) − e3( i , j ,k ) ⎤
−
−
⎢
⎥,
Δu 1
Δu 3
μ g ⎢⎣
⎥⎦
1
3 ( n+ )
(32)
1
3 ( n− )
h( i , j ,k 2) = h(i , j ,k 2)
(n)
(n)
( n)
( n)
Δ t ⎡ e2(i +1, j ,k ) − e2(i , j ,k ) e1( i , j +1,k ) − e1(i , j ,k ) ⎤
−
−
⎥
⎢
Δu 1
Δu 2
μ g ⎢⎣
⎥⎦
e1(i(,nj +,k1))
⎛
⎜1−σ
= e1(i(,nj ,)k ) ⎜
⎜
⎜1+ σ
⎝
Δt
2ε
Δt
2ε
e(2i ,( nj ,+k1))
1
( n+ )
2
3( i , j , k )
⎛
⎜1−σ
= e(2i ,( nj ,)k ) ⎜
⎜
⎜1+σ
⎝
Δt
2ε
Δt
2ε
e
1
(n+ )
⎡ ( n + 12 )
⎤
⎢ h1( i , j ,k ) − h1( i , j2−1,k ) ⎥
⎢
⎥.
1
Δu 2
⎞
⎥
+ Δ tσ ⎟ g ⎢
⎣
⎦
2
⎠
Δt
⎛
⎜ε
⎝
Célula Secundária
h3
1
( n+ )
2
3( i , j −1, k )
2
e3
(34)
e2
c
h2
h1
e1
Célula Primária
Fig. 6.
Componentes Co-variantes em Células de Yee Não Ortogonais:
r
Célula Primária para o Campo Elétrico E e Célula Sencundária para o
r
Campo Magnético H .
⎞
⎟
⎟+
⎟
⎟
⎠
1
( n+ )
⎡ ( n+ 12 )
⎤
2
h
h
−
⎢ 1( i , j ,k ) 1( i , j ,k −1) ⎥
⎢
⎥
1
Δu 3
⎛
⎞
⎥
⎜ ε + Δ tσ ⎟ g ⎢⎣
⎦
2
⎠
⎝
1
( n+ )
⎡ ( n+ 12 )
⎤
2
−
h
h
Δt
⎢ 3(i , j ,k )
3( i , j −1,k ) ⎥
−
⎢
⎥,
1
Δu 2
⎞
⎛
⎥
⎜ ε + Δ tσ ⎟ g ⎢⎣
⎦
2
⎠
⎝
Δt
(33)
(36)
Como pode ser observado, as Equações (9)-(10) são
similares às Equações (31)-(36). Deve ser notado que as
componentes contravariantes dependem de componentes covariantes, as quais são calculadas através de (22). É
apropriado dizer que Δul (l = 1, 2 e 3) são considerados
unitários (iguais a um), já que as informações métricas da
célula estão contidas no termo g [20].
⎞
⎟
⎟+
⎟
⎟
⎠
⎡
⎤
−h
⎢h
⎥
⎢
⎥
1
u
Δ
⎛
⎞
⎥
⎜ ε + Δ tσ ⎟ g ⎢⎣
⎦
2
⎝
⎠
1
( n+ )
⎡ ( n + 12 )
⎤
2
Δt
⎢ h2( i , j ,k ) − h2( i , j ,k −1) ⎥
−
⎢
⎥,
1
Δu 3
⎛
⎞
⎢
⎥
ε
σ
g
+
Δ
⎜
t ⎟
⎣
⎦
2
⎝
⎠
Δt
⎞
⎟
⎟+
⎟
⎟
⎠
Δt
(n)
3( i , j +1, k )
h
Δt
2ε
Δt
2ε
1
( n+ )
⎡ ( n + 12 )
⎤
⎢ h2( i , j ,k ) − h2( i −1,2 j ,k ) ⎥
⎢
⎥
1
Δu 1
⎛
⎞
⎥
⎜ ε + Δ t σ ⎟ g ⎢⎣
⎦
2
⎝
⎠
1
1 ( n− )
2
( i , j ,k )
h
1
2 ( n+ )
2
( i , j ,k )
e(3i (, nj ,+k1))
⎛
⎜1−σ
= e(3i (, nj ,)k ) ⎜
⎜
⎜1+ σ
⎝
(35)
r
Malhas para as componentes co-variantes do campo E
devem ser construídas de modo que as mesmas contornem de
forma apropriada a geometria do problema a ser analisado.
Tais malhas são chamadas de malhas primárias, compostas
por células primárias (Fig. 6). Como pode ser observado pela
Fig.6, os cantos das células secundárias (ponto C) são
posicionados no centro das células primárias. Dessa forma, os
pontos da malha secundária são calculados a partir dos dados
da malha primária.
D. A Fonte de Excitação e Cálculos de Correntes e Tensões
Para garantir o significado físico dos cálculos, alguns
detalhes devem ser observados para implementar a fonte de
excitação e para calcular correntes e tensões. Tais detalhes são
discutidos neste tópico.
MELO E SILVA DE OLIVEIRA AND LEONIDAS DA SILVA SOUZA : EMPLOYMENT OF THE FDTD
1) Implementação da Fonte de Excitação: A função usada
como fonte de excitação neste trabalho segue [17]. Ela é
descrita pelas Equações (37) e (38), dadas em volts por:
e, sabendo-se que I 3 pode ser calculado por I 3 = J 3Δs3 , onde
r
I3 =
t
− a1
−e
t
− a2
)sen 2 (Ωt )
(37)
e se t ≥ 1.5T f :
Vs (t ) = (Vpeak /Ao )(e
nas
t
−e
t
− a2
a1 = 1,9315/T f ,
T
− a1 0
−e
T
− a2 0
,
Vs (t )
g 33( i , j ,k ) g (33i , j ,k )
.
I1 =
g11 ⎛ ∂h3 ∂h2 ⎞
2
⎜ 2 − 3 ⎟ g 22 g 33 − g 23 .
u
u
∂
∂
g ⎝
⎠
(43)
r
Br
VBA = − ∫ E ⋅ dl .
A
temos:
(39)
Assim, é possivel controlar a magnitude física da tensão da
fonte de excitação simplesmente pela especificação de Vs (t).
2) Cálculo de Correntes em Condutores Finos: Neste
trabalho, condutores elétricos são considerados perfeitos.
r
Dessa forma, o cálculo da densidade corrente J que é
conduzida por um fio pode ser realizado pela aplicação da lei
de Ampère em sua forma original
r r r
J = ∇× H.
(42)
Usando (19) e (17) e considerando A e B dois pontos
r
vizinhos, localizados nos dois pontos que definem o vetor a j ,
Dessa forma, tem-se:
e3 ( i , j ,k ) =
g 22 ⎛ ∂h1 ∂h3 ⎞
2
⎜ 3 − 1 ⎟ g11 g 33 − g13
g ⎝ ∂u ∂u ⎠
3) Cálculo de Tensões: A diferença de potencial entre dois
pontos A e B, VBA , pode ser calculada pela Equação:
que gera E3 ( i , j ,k ) (V/m), e, observando a equação Eq.(24),
g (33i , j ,k ) , para dar o sentido físico à componente.
I2 =
e
e3 (i , j ,k ) , Vs (t ) deve ser
dividido pelo comprimento da aresta da célula em que tal
componente está localizada, o qual é dado por g 33( i , j ,k ) , o
por
(40)
B⎛ 3
r ⎞ r
VBA = − ∫ ⎜ ∑ei a i ⎟ ⋅ a j = −e j .
A
⎠
⎝ i =1
r
J3 =
g 33 , é possível escrever:
1 ⎛ ∂h2 ∂h1 ⎞
⎜ 1 − 2 ⎟ g 33
g ⎝ ∂u ∂u ⎠
(44)
A Equação (44) não é válida para o caso no qual o ponto A
está na superfície de um fio fino, como ilustrado pela Fig.7.
Para esse caso, é necessário levar em conta o raio r do
condutor. Dessa forma, a maneira mais simples de solucionar
esse problema é calcular o valor físico Ej e a distância entre A
r
e B, d AB . Denominando esse potencial de VBA
, obtém-se (veja
(24) e (14)).
r
VBA
= − E j d AB = −e j g jj
(
)
g jj − r .
Para o cálculo da terceira componente co-variante da
r
corrente de condução I , o seguinte procedimento pode ser
empregado. Observando as Equações (19) e (30), pode ser
notado que:
r
1 ⎛ ∂h2 ∂h1 ⎞ r
J3 =
⎜ 1 − 2 ⎟a3 .
g ⎝ ∂u ∂u ⎠
Considerando que | a3 |=
(41)
De forma semelhante, chega-se a:
), a2 = 2,5584/Tt , Tt = 500 ×10 −6 , Ao = e
T0 = ln(a1/a2 )/(a1 − a2 ) e Ω = π/(3T f ) .
Para excitar, por exemplo,
g 33 ⎛ ∂h2 ∂h1 ⎞
2
⎜ 1 − 2 ⎟ g11 g 22 − g12 .
g ⎝ ∂u ∂u ⎠
(38)
T f = 0,058 ×10 −6 ,
quais
g 33 ⎛ ∂h2 ∂h1 ⎞ r r
⎜ 1 − 2 ⎟ | a1 × a2 | .
g ⎝ ∂u ∂u ⎠
Finalmente, após o desenvolvimento algébrico da Equação
acima, verifica-se que:
I3 =
− a1
r
a área Δs3 é dada por Δs3 =| a1 × a2 | , tem-se:
se t < 1.5T f :
Vs (t ) = (Vpeak /Ao )(e
149
Fig. 7. A geometria considerada para calcular
r
.
VBA
(45)
150
IEEE LATIN AMERICA TRANSACTIONS, VOL. 6, NO. 2, JUNE 2008
III. RESULTADOS
Este tópico mostra os resultados obtidos através do método
FDTD em coordenadas não ortogonais e a comparação da
resposta de regime com fórmulas disponíveis na literatura.
A. Eletrodo Horizontal
Neste tópico, são apresentados alguns resultados obtidos
com o software desenvolvido em coordenadas gerais. Um
eletrodo horizontal, com inclinação de 45° em relação ao eixo
x, foi simulado e os resultados obtidos, como esperado, foram
idênticos aos gerados com o FDTD convencional com o
eletrodo orientado paralelamente a x Ou seja, o problema é
independente do sistema de coordenadas, desde que o eletrodo
esteja nas mesmas condições em relação ao solo.
Em [22], um eletrodo horizontal de 1,5 metro de
comprimento, enterrado a uma profundidade de 0,1m, foi
simulado com um esquema de medição similar ao apresentado
em [17] (Fig.1). Todavia, não foi feita uma análise a respeito
da influência do eletrodo de injeção de corrente quando a
haste está localizada a profundidades maiores. Dessa forma,
foram simuladas várias situações de modo a se verificar o
funcionamento do método e os resultados obtidos foram
comparados com a Equação (46) [23], dada por:
1
[ln(2 L/a ) − 1] +
2πLσ
4h 2 + L2 ⎤
1 ⎡ ⎛⎜ L + 4h 2 + L2 ⎞⎟
⎢ln
⎥,
+ 2h/L −
⎟
L
2πLσ ⎢ ⎜⎝
2h
⎥⎦
⎠
⎣
TABELA I
COMPARAÇÃO ENTRE VALORES DE REGIME (EM OHMS) OBTIDOS POR
SIMULAÇÕES FDTD TRADICIONAIS E PELA EQUAÇÃO ANALÍTICA PARA VÁRIAS
PROFUNDIDADES DE ELETRODO HORIZONTAL
( L = 1,5m , a = 10mm E σ = 500 −1 S/m )
h (m)
1,50
1,25
1,00
0,75
0,50
0,25
LN-FDTD (Ω)
161,39
173,52
188,13
206,67
233,30
300,16
Eq.(46) (Ω)
276,06
270,49
270,33
274,70
284,52
307,96
Dif. Relativa
41,5 %
35,8 %
30,4 %
24,8 %
18,0 %
2,53 %
Para efeito de comparação com (46), a tensão e a corrente
foram calculadas em uma das extremidades do eletrodo
horizontal, como ilustrado pela Fig.8.
O campo elétrico foi integrado paralelamente ao solo,
abaixo da superfície e de forma normal ao eletrodo de injeção
(para o cálculo do potencial).
Considerando o caso no qual L = h = 1,5m, a = 10 mm e σ =
500-1 siemens por metro, o resultado obtido é mostrado pela
Fig.9, o qual concorda com a Eq. (46) no estado estacionário
(veja também a Tabela I).
R=
(46)
h
na qual L é o comprimento do eletrodo, h é a profundidade na
qual o eletrodo está enterrado, a é o raio do eletrodo e σ é a
condutividade elétrica do solo (veja Fig.8).
272 Ω
V/I
(Ω)
Iterações
(3us)
I
V
L
Fig. 8. Eletrodo horizontal e o cálculo da corrente e da tensão.
A Tabela I compara os resultados de regime obtidos
utilizando-se o método ilustrado pela Fig.2 com os obtidos
utilizando-se (46), para diversos valores de h (0,25 1,5 m).
Observa-se que quanto maior a profundidade, maior a
discordância entre os resultados. Isso se deve a dois motivos:
1) a Equação (46) foi desenvolvida para o eletrodo horizontal
isolado, sem o cabo ou eletrodo de injeção e 2) como
conseqüência do motivo citado anteriormente, o eletrodo de
injeção contribui para reduzir o valor da resistência total vista
pelo sistema de aterramento, que na verdade se comporta
como um “L” em posição vertical. Dessa forma, é esperada a
redução do valor V/I com a profundidade h.
Fig. 9. Relação tensão/corrente em uma das extremidades do eletrodo
horizontal colocado a 1,5 metro de profundidade (linha: método FDTD com
eletrodo alinhado com o eixo x ; círculos: método LN-FDTD com eletrodo
inclinado a 45 graus do eixo x).
Observa-se que os métodos FDTD e LN-FDTD podem ser
usados para identificar a influência de elementos como o cabo
de conexão do sistema de aterramento.
Deve-se ressaltar que os problemas simulados em [24]
foram novamente testados com o software aqui desenvolvido
utilizando-se o método LN-FDTD e resultados idênticos
àqueles foram obtidos.
B. Sistema de aterramento do tipo prato
Para testar um caso bastante geral, foi escolhido para
simulação o sistema de aterramento do tipo prato, ilustrado
pela Fig.10. De acordo com o trabalho de Yung [25], o valor
de regime para a relação tensão/corrente é dado pela Equação
(47).
MELO E SILVA DE OLIVEIRA AND LEONIDAS DA SILVA SOUZA : EMPLOYMENT OF THE FDTD
(
)(
151
)
13
a −0,818e −0,0093a 1 + e −1,427 h
(47)
100σ
Para esse caso, foi necessário construir uma malha circular
específica, cujo corte horizontal está ilustrado pela Fig.11.
Através dela, observa-se que foram utilizadas 22 células para
modelar o diâmetro 2a da estrutura.
As Figuras 12 e 13 mostram os resultados obtidos para a
relação tensão/corrente para dois pratos a meio metro da
superfície, num solo caracterizado pelos parâmetros
σ = 2.28 × 10 −3 , ε r = 50 e μ = μ0 , com diâmetros
R=
2a = 5,5 m e 2a = 8 m, respectivamente.
Fig. 12. Relação Tensão/Corrente para um sistema de aterramento do tipo
prato com h = 0,5m e diâmetro de 5,5m ( σ = 2.28 ×10 −3 e εr = 50).
Fig. 10. Sistema de aterramento do tipo prato.
Para esses dois casos, a Equação (47) fornece os valores
32,69 e 23,05 ohms, enquanto o método LN-FDTD resultou
em 33,42 e 25 ohms, conforme mostram as Figuras 12 e 13, o
que representa uma ótima aproximação numérica.
Dessa forma, nota-se que o método LN-FDTD é uma
poderosa ferramenta para análise de situações realísticas de
sistemas de aterramento, especialmente para as que não
possuem soluções analíticas e para a análise de transitórios.
Fig. 11. Corte horizontal da malha gerada para simular o sistema
de aterramento do tipo prato.
Fig.13. Relação Tensão/Corrente para um sistema de aterramento do tipo
prato com h = 0,5m e diâmetro de 8,0m ( σ = 2.28 ×10 −3 e εr = 50).
IV. CONCLUSÕES
Neste trabalho, foi desenvolvida uma nova metodologia para
a análise de sistemas de aterramento solucionando-se as
Equações de Maxwell no domínio do tempo em Coordenadas
Gerais pelo método LN-FDTD. A vantagem desta formulação
é que a mesma abre a possibilidade de modelar precisamente
estruturas que não coincidam com o sistema de coordenadas
cartesianos. Além disso, a formulação provê uma solução de
onda completa para as equações de Maxwell no domínio do
Tempo, na qual os fenômenos de reflexão, refração e difração
são todos considerados implicitamente, o que implica
soluções bastante realistas. Foi mostrado como realizar o
cálculo de tensões e correntes (além da implementação da
fonte) nesse sistema de coordenadas de forma a garantir a
consistência das grandezas físicas.
Foi demostrado que a técnica LN-FDTD pode ser aplicada
para a solução desse tipo de problema (aterramento) através
da simulação de estruturas simples, cuja solução é conhecida.
Foram simulados: 1) eletrodo horizontal, para o qual
inclinações em relação ao eixo x foram implementadas, com
obtenção de resultados idênticos aos gerados pelo método
152
IEEE LATIN AMERICA TRANSACTIONS, VOL. 6, NO. 2, JUNE 2008
FDTD convencional e 2) um prato perfeitamente circular
enterrado, sendo obtida uma excelente concordância com a
Equação (47). Para a simulação desse problema (diâmetro de
5,5m), foram geradas cinco malhas com os seguintes números
de pontos para representar o diâmetro do prato: 10, 15, 20, 30
e 40. Foi observado que, a partir de vinte pontos, os resultados
obtidos não sofreram mudanças significativas, de forma que a
convergência do método foi assegurada.
Deve-se ressaltar que, como os problemas de aterramento
são problemas abertos, há a necessidade de se restringir a
região de análise, o que é feito aqui através da técnica de
camadas perfeitamente casadas uniaxiais (técnica U-PML). A
formulação U-PML em Coordenadas Gerais foi estendida para
truncar meios condutivos (como o solo).
AGRADECIMENTOS
Os autores agradecem o suporte fornecido pela Eletronorte,
pelo CNPq e pela UFPA.
analysis of grounding systems in homogeneous and stratified soils using
the FDTD method”, in International Symposium on Lightning Protection
(SIPDA), Sao Paulo (Brazil), 2005.
[19] J. F. Almeida, R. O. dos Santos, and C. L. da S. S. Sobrinho, “Computational technique to UPML absorbing boundary conditions by FDTD: A
complete approach”, IEEE Latin America Transactions, vol. 3, no. 5, pp.
377–382, Dec. 2005.
[20] A. Taflove and S. C. Hagness, Computational Electrodynamics, The
Finite-Difference Time-Domain Method, 3rd ed. Artech House Inc.,
2005.
[21] J. D. Jackson, Classical Electrodynamics. Wiley, New York, 1962.
[22] M. Tsumura, Y. Baba, N. Nagaoka, and A. Ametani, “FDTD Simulation
of a Horizontal Grounding Electrode and Modeling of its Equivalent Circuit”, IEEE Transactions on Electromagnetic Compatibility, vol. 48, pp.
817–825, Nov. 2006.
[23] M. A. Mattos, Técnicas de Aterramento. Okime Eletromagnetismo
Aplicado, 2004.
[24] E. T. Tuma, R. O. dos Santos, R. M. S. de Oliveira, and C. L. da S. S.
Sobrinho, “Transient analysis of parameters governing grounding systems by the fdtd method”, IEEE Latin America Transactions, vol. 4, no.
1, pp. 55–61, March 2006.
[25] E. K.-N. Yung and W.-Y. Tam, “Analysis of a circular earthing plate”,
in IEE proceedings. Part C. Generation, transmission and distribution
(IEE proc. C), 1989.
REFERÊNCIAS
[1] L. V. Bewley, “Theory and tests of the counterpoise”, Elec. Engr.,
pp.1163–1172, Aug. 1934.
[2] L. V. Bewley, “The counterpoise”, G. E. Rev., vol. 37, pp. 73–81, 1934.
[3] L. P. Bellaschi, “Impulse and 60-cycle characteristics of driven grounds,
part III effect of lead in ground installation”, in AIEE Transactions, vol.
62, 1943, pp. 334–345.
[4] E. D. Sunde, Earth conduction effects in transmission systems. by Bell
Telephone Laboratories, 1949.
[5] A. P. Meliopoulos and M. G. Moharam, “Transient analysis of grounding
systems”, IEEE Trans. Power Apparatus and Systems, vol. PAS-02,
No.2, pp. 389–399, 1983.
[6] A. D. Papalexopoulos and A. P. Meliopulos, “Frequency dependent
characteristics of grounding systems”, IEEE Trans. Power Delivery, vol.
2, pp. 1073–1081, 1987.
[7] A. Geri, “Behaviour of grounding systems excited by high impulse
currents: the model and its validation”, IEEE Trans. Power Delivery, vol.
14, No.3, pp. 1008–1017, 1999.
[8] A. F. Otero, J. Cidras, and J. L. del Alamo, “Frequency-dependent grounding system calculation by means of a conventional nodal analysis technique”, IEEE Trans. Power Delivery, vol. 14, No.3, pp. 873–878, 1999.
[9] J. Cidras, A. F. Otero, and C. Garrido, “Nodal frequency analysis of
grounding systems considering the soil ionization effect”, IEEE Trans.
Power Delivery, vol. 15, No.1, pp. 107–107, 2000.
[10] L. Grcev and F. Dawalibi, “An electromagnetic model for transients in
grounding system”, IEEE Trans. Power Delivery, vol. 5, pp. 1773–1781,
1990.
[11] L. Grcev, “Computation of transient voltages near complex grounding
systems caused by lightning currents”, in Proceedings of IEEE International Symposium on EMC, 92CH3169-0, 1992, pp. 393–399.
[12] L. D. Grcev and F. E. Menter, “Transient electromagnetic fields near
large earthing systems”, IEEE Trans. Magnetics, vol. 32, No.3, pp. 1525–
1528, 1996.
[13] L. D. Grcev, “Computer analysis of transient voltages in large grounding
systems”, IEEE Trans. Power Delivery, vol. 11, pp. 815–823, 1996.
[14] B. Nekhoul, C. Cuerin, P. Labie, G. Meunier, and R. Feuillet, “A finite
element method for calculating the electromagnetic fields generated by
substation grounding systems”, IEEE Trans. Magnetics, vol. 31, No.3,
pp. 2150–2153, 1995.
[15] B. Nekhoul, P. Labie, F. X. Zgainski, and G. Meunier, “Calculating the
impedance of a grounding system”, IEEE Trans. Magnetics, vol. 33,
No.3, pp. 1509–1512, 1996.
[16] K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media”, IEEE Trans. Antennas and
Propagation, vol. 14, pp. 302–307, 1966.
[17] K. Tanabe, “Novel method for analyzing the transient behavior of grounding systems based on the finite-difference time-domain method”,
CRIEPI Report – Tokio, 2001.
[18] E. T. Tuma, R. M. S. de Oliveira, and C. L. da S. S. Sobrinho, “New
model of current impulse injection and potential measurement in transient
Rodrigo M. S. de Oliveira graduou-se em 2003 pela
Universidade Federal do Pará (Brasil), alcançou o título de
Mestre em 2004 e busca o grau de Doutor pela mesma
instituição (Engenharia Elétrica). Possui experiência em
calibração de grandezas elétricas, obtida nas Centrais
Elétricas do Norte do Brasil. Desde 2002, faz parte do
Laboratório de Análises Numéricas em Eletromagnetismo
(LANE). Atualmente, realiza pesquisas nas áreas de
aterramento elétrico, técnicas de otimização e antenas com
o método FDTD.
Carlos Leonidas da S.S. Sobrinho graduou-se em 1981
em Engenharia Elétrica pela Universidade Federal do Pará
(Brasil), obteve seu título de Mestre em 1989 pela Universidade Católica do Rio de Janeiro (PUC-RJ) e o grau de
Doutor em Engenharia Elétrica pela Universidade Estadual de Campinas (Unicamp), em 1992. Em 1999, concluiu
seu pós-doutorado na Inglaterra (Queen Mary Westfield –
University of London). Sua pesquisa envolve espalhamento eletromagnético, guias de onda, propagação eletromagnética, descargas atmosféricas, antenas, estruturas periódicas, métodos numéricos e sistemas de aterramento.