MÉTODO FDTD INCONDICIONALMENTE ESTÁVEL BASEADO EM POLINÔMIOS DE LAGUERRE PARA SIMULAÇÕES ELETROMAGNÉTICAS EM REGIME TRANSITÓRIO
Matheus Soares da Silva
Projeto de Graduação apresentado ao Curso de Engenharia Elétrica da Escola Politécnica, Universidade Federal do Rio de Janeiro, como parte dos requisitos necessários à obtenção do título de Engenheiro.
Orientadores: Wescley Tiago Batista de Sousa Rubens de Andrade Júnior
Rio de Janeiro Setembro de 2017
MÉTODO FDTD INCONDICIONALMENTE ESTÁVEL BASEADO EM POLINÔMIOS DE LAGUERRE PARA SIMULAÇÕES ELETROMAGNÉTICAS EM REGIME TRANSITÓRIO
Matheus Soares da Silva
PROJETO CURSO DA DOS
DE
DE
GRADUAÇÃO ENGENHARIA
UNIVERSIDADE REQUISITOS
SUBMETIDO ELÉTRICA
FEDERAL
DO
NECESSÁRIOS
RIO
PARA
AO
CORPO
DA
ESCOLA
DE
JANEIRO
A
OBTENÇÃO
DOCENTE
DO
POLITÉCNICA COMO DO
PARTE
GRAU
DE
ENGENHEIRO ELETRICISTA.
Examinado por:
Prof. Rubens de Andrade Júnior, D.Sc.
Prof. Antonio Carlos Siqueira de Lima, D.Sc.
Eng. Bárbara Maria Oliveira Santos,
RIO DE JANEIRO, RJ BRASIL SETEMBRO DE 2017
Soares da Silva, Matheus Método FDTD Incondicionalmente Estável Baseado em Polinômios de Laguerre para Simulações Eletromagnéticas em Regime Transitório/Matheus Soares da Silva. Rio de Janeiro: UFRJ/ Escola Politécnica, 2017. XV, 81 p.: il.;
29, 7cm.
Orientadores: Wescley Tiago Batista de Sousa Rubens de Andrade Júnior Projeto de Graduação UFRJ/ Escola Politécnica/ Curso de Engenharia Elétrica, 2017. Referências Bibliográ cas: p. 56 58. 1.
WLP-FDTD.
Computacional.
et al.
2.
FDTD.
3.
Eletromagnetismo
I. Tiago Batista de Sousa, Wescley
II. Universidade Federal do Rio de Janeiro, Escola
Politécnica, Curso de Engenharia Elétrica. III. Título.
iii
Agradecimentos Gostaria de agradecer:
Ao meu orientador, Wescley, por ter me dado um tema desa ador e interessante para pesquisar enquanto Aluno de Iniciação Cientí ca.
Creio que esta foi
uma das melhores IC's que eu poderia ter feito e a experiência me ajudou a decidir pela carreira acadêmica. Agradeço também pelos conselhos variados que me deu e pelo bom humor que sempre teve. Obrigado, te desejo sucesso e felicidades pra toda a família.
Ao Professor Rubens, primeiro pela excelente disciplina de Eletromagnetismo II, uma das que mais gostei durante a graduação e o motivo pelo qual, em última análise, estou escrevendo este Projeto Final. Em seguida, pela paciência que teve comigo nestes últimos meses enquanto escrevia este trabalho, além das sugestões e dicas. Te levo como exemplo de acadêmico para o futuro.
A todos os professores do departamento que eventualmente me orientaram e deram dicas pro ssionais:
Antônio Lopes, Antônio Carlos Ferreira, Helói
Moreira, Maurício Aredes, Robson Dias, Sérgio Hazan e Walter Suemitsu.
Aos colegas de LASUP por terem sempre criado um ambiente agradável de trabalho, devo agradecer dentre eles especialmente à Bárbara pelas diversas dicas sobre como navegar esse período de m de curso e por ter me cedido momentaneamente o seu computador, muito obrigado! Também devo agradecer ao Felipe Costa por ter me recomendado muito material sobre diferenças nitas e ao Professor Elkin por estar sempre me dando incentivo e perguntando como estava andando o trabalho.
A Felipe Cabral por ter disponibilizado o padrão LaTeX que utilizei pra escrever este trabalho.
A minha mãe, Teresa Cristina, pelo suporte e pela con ança que depositou em mim.
Aos meus amigos de faculdade pelo companheirismo, especialmente ao Felipe
iv
Dicler pelas diversas conversas sobre os mais variados assuntos, dicas acadêmicas que me deu em diversos momentos e por ter me apresentado o Inkscape! Também devo agradecer especialmente a Mike Mattos por ter sido meu bom amigo desde os tempos de CEFET e a Erick Gama, meu amigo desde o primeiro período. Desejo sucesso e felicidades a todos vocês!
Ao meu amigo, Fred Guimarães, pelas dicas que me deu relativas a este trabalho.
v
Resumo do Projeto de Graduação apresentado à Escola Politécnica/ UFRJ como parte dos requisitos necessários para a obtenção do grau de Engenheiro Eletricista.
MÉTODO FDTD INCONDICIONALMENTE ESTÁVEL BASEADO EM POLINÔMIOS DE LAGUERRE PARA SIMULAÇÕES ELETROMAGNÉTICAS EM REGIME TRANSITÓRIO
Matheus Soares da Silva
Setembro/2017
Orientadores: Wescley Tiago Batista de Sousa Rubens de Andrade Júnior Curso: Engenharia Elétrica
O método das Diferenças Finitas no Domínio do Tempo (FDTD, do inglês:
Dierence Time Domain ),
Finite
desenvolvido por K. Yee em 1966, é uma forma larga-
mente utilizada e relativamente simples de, numericamente, resolver as Equações de Maxwell. O algoritmo encontra um vasto ramo de aplicações em engenharia elétrica, do projeto de antenas à análise de transitórios eletromagnéticos. No entanto, esta técnica pode se tornar instável na análise de alguns sistemas físicos devido a condição de Courant-Friedrichs-Lewy, que exige que tempo e espaço estejam discretizados de forma acoplada. Trabalhos tem sido desenvolvidos desde os anos 2000 no sentido de desenvolver versões incondicionalmente estáveis do algoritmo. Uma destas versões se baseia nas propriedades dos polinômios de Laguerre. Este trabalho contém o desenvolvimento do algoritmo FDTD baseado em polinômios de Laguerre e sua aplicação a dois casos de eletrodinâmica; um deles unidimensional e o outro bidimensional. Ambos os casos são validados por comparação com resultados presentes na literatura.
vi
Abstract of Undergraduate Project presented to POLI/UFRJ as a partial ful llment of the requirements for the degree of Engineer.
UNCONDITIONALLY STABLE FDTD SCHEME BASED ON LAGUERRE POLYNOMIALS FOR SIMULATIONS OF ELECTROMAGNETIC TRANSIENTS
Matheus Soares da Silva
September/2017
Advisors: Wescley Tiago Batista de Sousa Rubens de Andrade Júnior Course: Electrical Engineering
The Finite-Dierence Time-Domain method based on the Yee Algorithm is a traditional and relatively simple numerical technique to solve the Maxwell Equations.
This algorithm nds many applications throughout electrical engineering,
from antenna design to electromagnetic transients analysis. However, this technique becomes unstable in the analysis of some physical systems because of the CourantFriedrichs-Lewy Condition, that introduces a coupling requirement between the discretization of time and space. Work on the development of unconditionally stable versions of the algorithm is ongoing since the 2000's. One of these versions is based on the Laguerre polynomials properties. The present work presents the development of the Laguerre Polynomials based FDTD algorithm and application to two eletrodynamics problems. The results are validated by comparison with those available in the scienti c literature.
vii
Sumário Lista de Figuras
x
Lista de Tabelas
xii
Lista de Símbolos
xiii
Lista de Abreviaturas
xv
1 Introdução
1
1.1
Motivação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Objetivos
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3
Organização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2 O Método FDTD
3
2.1
Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2
Diferenças Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.3
Equações de Maxwell . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3.1
Modos Transversais Eletromagnéticos (Modos TEM)
. . . . .
8
2.3.2
Modo Transversal Elétrico (Modo TE)
. . . . . . . . . . . . .
10
2.3.3
Modo Transversal Magnético (Modo TM)
2.4
2.5
. . . . . . . . . . .
11
Algoritmo de Yee . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4.1
15
Condições de Fronteira . . . . . . . . . . . . . . . . . . . . . .
Estabilidade e a Condição Courant-Friedrichs-Lewy
. . . . . . . . . .
3 O Método FDTD baseado em Polinômios de Laguerre
17
21
3.1
Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.2
Polinômios de Laguerre . . . . . . . . . . . . . . . . . . . . . . . . . .
22
3.3
Ortogonalidade
23
3.4
O Método FDTD Baseado em Polinômios de Laguerre
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
3.4.1
Modo TE
3.4.2
Modo TEM
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
3.4.3
Condições de Fronteira . . . . . . . . . . . . . . . . . . . . . .
35
viii
4 Simulações e Resultados
38
4.1
Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
4.2
Estrutura do Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . .
38
4.3
Simulação Bidimensional: Sistema Físico Apresentado por Chung et al. 40 4.3.1
4.4
Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
Simulação Unidimensional: Propagação de Onda Transversal Eletromagnética 4.4.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
5 Conclusão e Trabalhos Futuros
54
5.1
Conclusão
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
5.2
Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
Referências Bibliográ cas
56
A Critério de Von Neumann para Estabilidade Numérica
59
B Correlação Linear
62
C Código Desenvolvido para o Caso Bidimensional utilizando WLPFDTD 63 D Código Desenvolvido para o Caso Unidimensional utilizando WLPFDTD 77
ix
Lista de Figuras 2.1
Aumento no número de artigos publicados em IEEE
Microwaves
utilizando FDTD a partir dos anos 90.
Transactions on . . . . . . . . . .
2.2
Diferenças centradas (DC), regressivas (DR) e progressivas (DP).
2.3
Malha de Yee para o modo
2.4
Malha de Yee para o modo
2.5
Visualização da Malha de Yee no espaço-tempo, o círculo vermelho
TEz . . . . . . . . . . TEM em relação a x
4
. .
6
. . . . . . . . . . . .
12
polarizado em
y.
. .
13
mostra o ponto onde o operador diferencial é discretizado para que se obtenha a Equação de Atualização para Campos Elétricos (Equação 2.43). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6
14
Visualização da Malha de Yee no espaço-tempo, o círculo vermelho mostra o ponto onde o operador diferencial é discretizado para que se obtenha a Equação de Atualização para Campos Magnéticos (Equação 2.44).
2.7
Simulação de Campo Elétrico
1.0. 2.8
Ey
utilizando método FDTD com
Simulação de Campo Elétrico
utilizando método FDTD com
Ey
utilizando método FDTD com
18
Sc =
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulação de Campo Elétrico
1.005.
Ey
14
Sc =
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.005. 2.8
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
Sc =
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3.1
Polinômios de Laguerre de ordem 0 a 4. . . . . . . . . . . . . . . . . .
23
3.2
Bases de Laguerre até ordem 4, escala temporal
3.3
Síntese de um sinal utilizando diferentes números de bases de Laguerre. 27
3.4
Síntese de um sinal utilizando 100 bases de Laguerre.
4.1
Fluxograma do Algoritmo Implementado para o Método WLP-FDTD. 40
4.2
Simulação bidimensional para validação do método. . . . . . . . . . .
4.3
Forma da densidade de corrente elétrica inserida no problema bidi-
4.4
s = 1.
. . . . . . . .
. . . . . . . . .
26
27
41
mensional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
Resultado da simulação bidimensional mostrada na Figura 4.2. . . . .
43
x
4.5
Resultados para a simulação bidimensional utilizando o método WLPFDTD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5
45
Resultados para a simulação bidimensional utilizando o método WLPFDTD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
4.6
Estrutura do problema unidimensional. . . . . . . . . . . . . . . . . .
48
4.7
Forma de Onda da Densidade Super cial de Corrente
4.8
Sobreposição dos resultados da simulação unidimensional utilizando
Ky . .
. . . . . .
o método WLP-FDTD e dos apresentados por Feynman et al. [1]. . . 4.8
50
Sobreposição dos resultados da simulação unidimensional utilizando o método WLP-FDTD e dos apresentados por Feynman et al. [1]. . .
4.8
49
51
Sobreposição dos resultados da simulação unidimensional utilizando o método WLP-FDTD e dos apresentados por Feynman et al.[1]. . . .
xi
52
Lista de Tabelas 2.1
Algumas das Grandezas Físicas Eletromagnéticas. . . . . . . . . . . .
4.1
Comparação entre os valores obtidos por Chung et al. e os obtidos por algoritmo próprio, para os pontos
4.2
a, b , c
e
d.
. . . . . . . . . . . .
7
44
Comparação entre os valores analíticos e os obtidos pelo método WLP-FDTD para campo elétrico no sistema físico unidimensional. . .
xii
52
Lista de Símbolos A
Matriz de coe cientes utilizada no método WLP-FDTD, p. 33
Bf
Máxima frequência, em Hertz, presente em uma forma de onda., p. 28
H E Cx,y , Cx,y
Constantes utilizadas no método WLP-FDTD, p. 32
Ln
Polinômio de Laguerre de ordem
n,
NL
Número necessário de bases de Laguerre para sintetizar uma
p. 22
certa função., p. 28
O(hn )
Notação
Big-Oh
denotando a ordem
n
de um certo método de
diferenças nitas., p. 5
RP earson Tf U ∆t
Coe ciente de Correlação de Pearson, p. 38 Duração, em segundos, de uma forma de onda., p. 28 Função qualquer descrita como uma série., p. 29 Discretização Temporal, p. 13
∆x
Discretização Espacial, coordenada
x,
p. 13
∆y
Discretização Espacial, coordenada
y,
p. 13
Φ0,1,2... β q−1
Conjunto ortogonal in nito de funções, p. 24 Vetor coluna contendo os somatórios dos termos de campos elétricos e magnéticos de todas as iterações anteriores do método WLP-FDTD, p. 33
δmn
Delta de Kronecker, p. 24
Permissividade Elétrica, Unidade SI:
µ
Permeabilidade Magnética, Unidade SI:
xiii
F/m,
p. 7
H/m,
p. 7
φp,q σ
Base de Laguerre de ordem
p
ou
q,
p. 25
Condutividade Elétrica, Unidade SI:
S/m,
p. 7
~ B
Densidade de Fluxo Magnético, Unidade SI:
~ D
Densidade de Fluxo Elétrico, Unidade SI:
~ E
Intensidade de Campo Elétrico, Unidade SI:
~ H
Intensidade de Campo Magnético, Unidade SI:
J~ ~ K c0,1,2...
Densidade de Corrente, Unidade SI:
W b/m2 ,
C/m2 ,
A/m2 ,
p. 7
p. 7
V /m,
p. 7
A/m,
p. 7
p. 7
Densidade Super cial de Corrente, p. 46 Coe cientes a serem determinados., p. 24
f
Função qualquer., p. 4
h
Raio de um certo intervalo., p. 4
i
Índice horizontal de um certo ponto na malha de Yee, p. 14
j
Índice vertical de um certo ponto na malha de Yee, p. 14
k
Um número natural qualquer., p. 22
p1
Ponto de medição de campo elétrico., p. 41
s
Fator de escala temporal, p. 26
t
Tempo, p. 12
u, v, w v
Vetores quaisquer., p. 23 Velocidade qualquer de propagação de uma onda eletromagnética, dada em
w(t) x0
m/s,
p. 16
Função peso., p. 25 Ponto central de um intervalo, ponto este no qual se deseja avaliar o operador diferencial., p. 4
L xN max
Maior zero dentre os polinômios de Laguerre utilizados para sintetizar uma certa função., p. 28
xiv
Lista de Abreviaturas ABC
Absorbing Boundary Conditions, p.
CAE
Computer Assisted Engineering, p.
CFL
Courant-Friedrichs-Lewy, p. 1
FDTD IEEE LASUP
Finite Dierence Time Domain, p.
PML
Perfectly Matched Layer, p.
WLP-FDTD
vi 3
Laboratório de Aplicações de Supercondutores, p. 2
Perfect Electric Conductor, p.
TM
1
Institute of Electrical and Electronic Engineers, p.
PEC
TE
15
15
16
Transversal Elétrica(o), p. 9 Transversal Magnética(o), p. 9
Weighted Laguerre Polynomials Finite Dierence Time Domain, p. 2
xv
Capítulo 1 Introdução 1.1
Motivação
Métodos computacionais encontram vasta aplicação na solução de problemas em engenharia.
A utilização destes nas últimas décadas permitiu o crescimento de
uma área conhecida como CAE (do inglês:
Computer Assisted Engineering ).
Esta
popularidade advém da possibilidade de se fazer análise física de sistemas complexos de maneira mais rápida e simples, sem necessidade de se criar um modelo analítico para cada novo problema, além de se conseguir tratar problemas para os quais não existe modelo analítico.
Com estes programas, em muitos casos, só se torna
necessário representar corretamente o problema em termos geométricos, de materiais e de condições de contorno. No entanto, a utilização de
softwares comerciais na simulação de materiais super-
condutores do tipo II se revela um problema em aberto. Estes materiais apresentam acoplamento entre temperatura, densidade de corrente e campo magnético, além de um comportamento não-linear na relação entre campo elétrico e densidade de corrente. Tais di culdades motivam a criação de novos modelos de simulação, utilizando códigos próprios [2] [3] [4] [5] [6]. Em eletrodinâmica, o método das Diferenças Finitas no Domínio do Tempo (ou
Finite-Dierence Time Domain ) é largamente empregado e possui diversas implementações na forma de softwares, comerciais e gratuitos. Este método
FDTD, do inglês
se baseia na discretização tanto do espaço quanto do tempo em termos de diferenças nitas. No entanto, para que o método mantenha sua estabilidade, a discretização do espaço e do tempo deve estar acoplada. Este acoplamento é dado pela condição de Courant-Friedrichs-Lewy (ou condição CFL) [7] [8] [9]. Com isso, o método se torna impróprio para as chamadas simulações multiescala, onde se deseja observar o comportamento de um sistema em escalas diferentes em uma mesma simulação, ou seja, em casos onde a discretização do espaço de
1
simulação é bastante na em relação ao espaço, ou tempo, total de simulação. Nestas situações o método se torna muito dispendioso computacionalmente [10].
Este se
revela o caso em aplicações de supercondutores em sistemas de potência, onde a malha de discretização do espaço deve ser na e, ao mesmo tempo, as frequências são relativamente baixas e um tempo relativamente longo de simulação se faz necessário [5]. O método das Diferenças Finitas no Domínio do Tempo usando Polinômios de Laguerre (também conhecido como Laguerre-FDTD ou WLP-FDTD, do inglês
Weigh-
ted Laguerre Polynomials - Finite-Dierence Time Domain ) utiliza as propriedades dos polinômios de Laguerre para modi car o método FDTD tradicional e torná-lo incondicionalmente estável, permitindo que tempo e espaço possam ser discretizados de maneira desacoplada e que as simulações sejam mais e cientes [11] [10]. Utilizando o método WLP-FDTD, este trabalho foi desenvolvido no Laboratório de Aplicações de Supercondutores da UFRJ (LASUP) com o objetivo de dar uma contribuição inicial a uma linha de pesquisa que busca avançar simulações multifísicas de supercondutores, adicionando a descrição das contribuições da parcela eletromagnética a métodos baseados em diferenças nitas que já realizam, com sucesso, simulações eletrotérmicas [5] [12].
1.2
Objetivos
O objetivo é apresentar as bases teóricas do método FDTD baseado em polinômios de Laguerre, desenvolver um algoritmo próprio e, por comparação com resultados da literatura, demonstrar que, utilizando o domínio das bases de Laguerre é possível construir um método preciso e incondicionalmente estável de solução numérica das Equações de Maxwell para regimes transitórios.
1.3
Organização
A organização deste trabalho é feita da seguinte forma:
O Capítulo 1 contém a
motivação, os objetivos e esta organização. O Capítulo 2 trata do Método FDTD baseado no Algoritmo de Yee, passando pelas Equações de Maxwell e por alguns conceitos de métodos numéricos. O Capítulo 3 trata dos Polinômios de Laguerre e do Método das Diferenças Finitas no Domínio do Tempo utilizando estes Polinômios. O Capítulo 4 traz a metodologia, a estrutura geral dos algoritmos desenvolvidos, a descrição das simulações feitas e os resultados e análise destes. Por m, o Capítulo 5 traz a conclusão e sugestões para trabalhos futuros.
2
Capítulo 2 O Método FDTD 2.1
Introdução
Desenvolvidas no século XIX, as Equações de Maxwell são de fundamental importância em todas as áreas da ciência que utilizam teoria eletromagnética. Até meados do século XX, a modelagem de sistemas eletromagnéticos era feita, principalmente, utilizando o domínio da frequência em estado estacionário com soluções normalmente encontradas analiticamente, utilizando séries in nitas ou formas fechadas. Tais abordagens possuem algumas limitações no tratamento de materiais com geometrias complexas, ou que possuam não linearidades. Em 1966, K. Yee desenvolve um método numérico no domínio do tempo que, embora simples conceitualmente, é capaz de resolver as Equações de Maxwell [7]. A popularidade de soluções diretamente no domínio do tempo, discretizando o espaço-tempo em forma de malha, aumenta durante as décadas de 80 e 90, impulsionada pelo desenvolvimento dos computadores modernos.
Também durante as décadas de 80 e 90 o método de
Yee é re nado por outros pesquisadores como Ta ove, Mur, Berenger e Sullivan, adicionando funcionalidades ao método e estudando seus critérios de estabilidade. O grá co da Figura 2.1 mostra o crescimento na quantidade de artigos publicados no periódico IEEE
Transactions on Microwaves
80 e 90 [13].
3
tratando do assunto entre os anos
Publicações Utilizando FDTD
Número de Publicações
500
400
300
200
100
0 1965
1970
1975
1980
1985
1990
1995
Ano Figura 2.1: Aumento no número de artigos publicados em IEEE
Microwaves 2.2
Transactions on
utilizando FDTD a partir dos anos 90.
Diferenças Finitas
Métodos numéricos de Diferenças Finitas para solução de Equações Diferenciais se baseiam fundamentalmente na discretização das variáveis de interesse.
Isso é
feito aplicando-se aproximações por diferenças nitas para operadores diferenciais. O método FDTD utiliza, especi camente, a chamada aproximação por diferenças centradas tanto para o espaço quanto para o tempo, fazendo deste um dos métodos conhecidos como métodos raio
h
ao redor do ponto
leapfrog
x0 ,
[14]. Tomamos uma função
f
e um intervalo de
usar expansão por Série de Taylor para representar a
função no extremo positivo do intervalo resulta na Equação (2.1) [9] [14] [15]:
h2 ∂ 2 f (x) h3 ∂ 3 f (x) ∂f (x) + + + ··· f (x0 + h) = f (x0 ) + h ∂x x=x0 2! ∂x2 x=x0 3! ∂x3 x=x0 Rearranjando termos e dividindo por
h
(2.1)
os dois lados da Equação (2.1) obtemos:
∂f (x) h ∂ 2 f (x) h2 ∂ 3 f (x) f (x0 + h) − f (x0 ) = + + + ··· h ∂x x=x0 2! ∂x2 x=x0 3! ∂x3 x=x0
(2.2)
Truncando o lado direito desta equação no segundo termo conseguimos uma aproximação pra derivada conhecida como aproximação por diferenças progressivas (Equação (2.3)):
4
∂f (x) f (x0 + h) − f (x0 ) + O(h) = ∂x x=x0 h O termo
O(h)
truncamento. a
h,
está representado na notação
Big Oh
(2.3)
e simboliza o erro devido ao
Como neste caso o primeiro termo do truncamento é proporcional
o erro de truncamento cresce linearmente com o passo de discretização
h.
Aproximações que apresentam esta relação linear entre erro e discretização são ditas aproximações de primeira ordem.
Por outro lado, realizar a expansão por Séries de Taylor para representar a função no extremo inferior do intervalo
[x0 − h, x0 + h]
resulta na Equação (2.4):
h2 ∂ 2 f (x) h3 ∂ 3 f (x) ∂f (x) + − + ··· f (x0 − h) = f (x0 ) − h ∂x x=x0 2! ∂x2 x=x0 3! ∂x3 x=x0
(2.4)
Seguindo o mesmo procedimento anterior chegamos na chamada aproximação por diferenças regressivas, apresentada na Equação (2.5). Este tipo de aproximação é, assim como a aproximação por diferenças progressivas, de primeira ordem.
∂f (x) f (x) − f (x0 − h) + O(h) = ∂x x=x0 h
(2.5)
É possível, a partir das Equações (2.1) e (2.4), chegar em mais um tipo de aproximação, conhecida como aproximação por diferenças centradas (ou centrais). Começamos subtraindo as duas equações:
∂f (x) 2 h3 ∂ 3 f (x) f (x0 + h) − f (x0 − h) = h + + ··· ∂x x=x0 3! ∂x3 x=x0 Dividindo por
h
(2.6)
nos dois lados da equação obtemos:
f (x0 + h) − f (x0 − h) ∂f (x) 2 h2 ∂ 3 f (x) = + + ··· h ∂x x=x0 3! ∂x3 x=x0 Importante notar que, neste caso, o termo a ser truncado depende de
(2.7)
h2
e, por-
tanto, o erro de truncamento cai com o quadrado da discretização. A aproximação por diferenças centradas (Equação (2.8)) é dita de segunda ordem e, como já dito anteriormente, é a base do método FDTD.
5
f (x0 + h) − f (x0 − h) ∂f (x) = + O(h2 ) h ∂x x=x0
(2.8)
A aproximação por diferenças centradas também é aplicável, naturalmente, às funções de mais de uma variável e suas derivadas parciais, como mostrado na Equação (2.9) [15].
∂f (x, y) f ((x0 + h), y) − f ((x0 − h), y) = + O(h2 ) h ∂x x=x0
(2.9)
A Figura 2.2 ilustra os três tipos de aproximação apresentados. É notável que, para uma aproximação razoável da derivada da função
f (x)
no ponto
x0 ,
as apro-
ximações por diferenças regressivas e progressivas necessitam de um intervalo
h
de
discretização consideravelmente menor do que a aproximação por diferenças centradas.
DC
DP DR
Figura 2.2: Diferenças centradas (DC), regressivas (DR) e progressivas (DP).
6
2.3
Equações de Maxwell
As Equações de Maxwell na forma diferencial mostradas nas Equações (2.10), (2.11), (2.12) e (2.13) governam fenômenos eletromagnéticos em um meio sem cargas livres descrito por coordenadas cartesianas [1] [7] [16] [17] :
~ ~ = J~ + ∂ D ∇×H ∂t ~ =− ∇×E
(2.10)
~ ∂B ∂t
(2.11)
~ =0 ∇·D
(2.12)
~ =0 ∇·B
(2.13)
Supondo a presença exclusiva de materiais lineares, isotrópicos e não dispersivos, podemos usar as relações constitutivas (2.14) e (2.15).
~ = E ~ D
(2.14)
~ = µH ~ B
(2.15)
Na presença de materiais condutores é preciso adicionar também a Equação (2.16).
~ J~ = J~independente + σ E
(2.16)
A Tabela 2.1 descreve as grandezas físicas utilizadas até agora. Tabela 2.1: Algumas das Grandezas Físicas Eletromagnéticas.
Representação
Grandeza Física
Unidade SI
~ E ~ D
Campo Elétrico
~ H ~ B
Densidade de Fluxo Magnético
J~
Densidade de Corrente Elétrica
Permissividade Elétrica
µ
Permeabilidade Magnética
V m C m2 A m Wb m2 A m2 F m H m
σ
Condutividade Elétrica
S
Densidade de Fluxo Elétrico Campo Magnético
7
Utilizando as equações constitutivas (2.14) e (2.15), podemos reescrever as Equações de Maxwell em função somente dos campos
~ E
e
~. H
~ ∂E 1 ~ − J) ~ = (∇ × H ∂t
(2.17)
~ 1 ∂H ~ = − (∇ × E) ∂t µ
(2.18)
~ =0 ∇ · E
(2.19)
~ =0 ∇ · µB
(2.20)
Escrevendo os rotacionais de maneira explícita em termos de derivadas espaciais nas Equações (2.17) e (2.18) obtemos as seguintes seis equações acopladas:
∂Hx 1 ∂Ey ∂Ez = − ∂t µ ∂z ∂y
(2.21)
1 ∂Ez ∂Ex ∂Hy = − ∂t µ ∂x ∂z
(2.22)
∂Hz 1 ∂Ex ∂Ey = − ∂t µ ∂y ∂x
(2.23)
1 ∂Hz ∂Hy ∂Ex = − − Jx ∂t ∂y ∂z 1 ∂Hx ∂Hz ∂Ey = − − Jy ∂t ∂z ∂x
(2.24)
(2.25)
∂Ez 1 ∂Hy ∂Hx = − − Jz ∂t ∂x ∂y
(2.26)
É possível modi car as Equações (2.21), (2.22), (2.23), (2.24), (2.25) e (2.26) de modo a reduzir a ordem do espaço de simulação e descrever o comportamento de sistemas eletromagnéticos bidimensionais e unidimensionais.
2.3.1
Modos Transversais Eletromagnéticos (Modos TEM)
Considerando que a geometria estudada, os campos e eventuais excitações de corrente se estendem in nitamente sem sofrer qualquer variação na direção de duas das três coordenadas espaciais, é possível anular duas das três derivadas espaciais presentes nas Equações (2.21), (2.22), (2.23), (2.24), (2.25) e (2.26). O desenvolvimento neste trabalho é feito anulando as derivadas em relação a
8
y
e
z.
Este procedimento
resulta nas Equações (2.31), (2.32), (2.33) e (2.34):
∂Ez 1 ∂Hy ~ = − Jz ∂t ∂x
(2.27)
∂Hy 1 ∂Ez = ∂t µ ∂t
(2.28)
∂Ey 1 = ∂t
∂Hz − − Jy ∂x
∂Hz 1 = ∂t µ
∂Ey − ∂x
(2.29)
(2.30)
Estas quatro equações podem ser organizadas em dois conjuntos compostos de duas equações e independentes entre si. Estes conjuntos são denominados modos. O conjunto formado pelas Equações (2.31) e (2.32) é chamado de modo transversal eletromagnético polarizado em
z
[7] ou de modo transversal magnético unidimensi-
onal (1D TM) [18]. Neste trabalho é feita a opção de utilizar a nomenclatura usada por Ta ove em [7]. Este modo será mencionado como modo
TEMx polar. z .
1 ∂Hy ∂Ez = − J~z ∂t ∂x
(2.31)
1 ∂Ez ∂Hy = ∂t µ ∂t
(2.32)
O conjunto formado pelas Equações (2.33) e (2.34) é chamado de modo transversal eletromagnético polarizado em
y
[7] ou de modo transversal elétrico unidi-
mensional (1D TE) [18]. Neste trabalho é feita a opção de utilizar a nomenclatura usada por Ta ove em [7]. Este modo será mencionado como modo
∂Ey 1 = ∂t
∂Hz − − Jy ∂x
∂Hz 1 = ∂t µ
Os dois modos apresentados,
∂Ey − ∂x
TEMx polar. z
e
TEMx polar. y .
(2.33)
TEMx polar. y ,
(2.34)
constituem duas for-
mas possíveis de se realizar simulações unidimensionais no eixo
x.
É evidente que
estes modos somente podem ser usados para se obter resultados precisos quando o problema é composto de geometrias e grandezas físicas que possam ser aproximadas
9
como invariantes em duas das três coordenadas cartesianas.
2.3.2
Modo Transversal Elétrico (Modo TE)
Considerando que a geometria estudada, os campos e eventuais excitações de corrente se estendem in nitamente em apenas uma das direções sem sofrer qualquer variação, as derivadas em relação a esta direção se tornam nulas e podemos construir equações bidimensionais. Como exemplo, é suposto que não há qualquer variação na direção
z
e que, portanto, as derivadas em relação a
z
se tornam nulas e as
Equações (2.21), (2.21), (2.22), (2.23), (2.24), (2.25) e (2.26) podem ser um pouco simpli cadas e resultam em outro conjunto de seis equações:
∂Hx 1 = ∂t µ
∂Ez − ∂y
(2.35)
∂Hy 1 ∂Ez = ∂t µ ∂x
(2.36)
1 ∂Ex ∂Ey ∂Hz = − ∂t µ ∂y ∂x
(2.37)
1 ∂Hz ∂Ex = − Jx ∂t ∂y ∂Ey 1 ∂Hz = − Jy − ∂t ∂x 1 ∂Hy ∂Hx ∂Ez = − − Jz ∂t ∂x ∂y
(2.38)
(2.39)
(2.40)
As Equações (2.35), (2.36), (2.37), (2.38), (2.39) e (2.40) podem dar origem a dois conjuntos de equações independentes entre si, denominados modos, de acordo com a orientação dos campos elétricos e magnéticos no espaço [18]. No modo transversal elétrico o campo magnético somente possui componente não-nula na direção na qual não há qualquer variação (no caso, a direção
z ),
enquanto os campos elétricos só
possuem componentes não-nulas nas direções perpendiculares a esta, no caso, no plano
xy .
Para de nir este modo, basta agrupar as Equações (2.37), (2.38) e (2.39).
z que não há qualquer variação, este modo é dito modo transversal elétrico em relação a z , ou modo TEz [7]. Como é na direção
1 ∂Hz ∂Ex = − Jx ∂t ∂y ∂Ey 1 ∂Hz = − − Jy ∂t ∂x
10
(2.41)
(2.42)
∂Hz 1 ∂Ex ∂Ey = − ∂t µ ∂y ∂x
2.3.3
(2.43)
Modo Transversal Magnético (Modo TM)
Já no modo transversal magnético, o campo elétrico somente possui componente não-nula na direção na qual não há variação (no caso, a direção
z ),
enquanto os
campos magnéticos só possuem componentes não-nulas nas direções perpendiculares (no caso, o plano
xy ).
Para de nir este modo basta agrupar as Equações (2.35),
(2.36) e (2.40). Este modo é dito Modo
modo transversal magnético em relação a z ,
ou
TMz : ∂Hx 1 = ∂t µ
∂Ez − ∂y
(2.44)
∂Hy 1 ∂Ez = ∂t µ ∂x
(2.45)
1 ∂Hy ∂Hx ∂Ez = − − Jz ∂t ∂x ∂y
(2.46)
Os dois modos apresentados,TMz e
TEz
constituem as duas formas possíveis
de se realizar simulações bidimensionais no plano
xy .
É evidente que estes modos
somente podem ser usados para se obter resultados precisos quando o problema é composto de geometrias e grandezas físicas que possam ser aproximadas como invariantes na direção
z.
11
2.4
Algoritmo de Yee
O método FDTD discretiza o tempo (t) e o espaço em diferenças centrais e desloca as grandezas eletromagnéticas no espaço, fazendo com que cada componente de campo elétrico
~ E
esteja cercada por quatro componentes de campo magnético
~, H
e vice-versa [9]. A Figura 2.3 mostra como exemplo a malha de Yee para estudos
TEz ) e a Figura 2.4 para estudos unidimensionais, modo TEM x, polarizado em y .
bidimensionais (modo em relação a
Figura 2.3: Malha de Yee para o modo
12
TEz .
Figura 2.4: Malha de Yee para o modo
Os índices
i
e
j
TEM
em relação a
x
polarizado em
y.
são os índices de discretização horizontal e vertical, respecti-
vamente. Estes índices são números inteiros de nidos tais que, sendo
∆x
e
∆y
os
passos de discretização espacial, um ponto no espaço localizado na malha de Yee e com coordenadas cartesianas da forma uma expressão da forma
(i∆x, j∆y)
(x, y)
possa ter sua posição descrita com
[7].
É importante observar que campos elétricos e magnéticos estão deslocados também no tempo. Como exemplo trataremos do caso unidimensional. Após aplicar a aproximação por diferenças centradas para as derivadas temporais e espaciais das Equações (2.33) e (2.34) do modo transversal eletromagnético polarizado na direção
y,
existirão duas equações conhecidas como Equações de Atualização. Uma destas
equações, (2.47), atualizará os valores para os campos elétricos ção (2.48)) atualizará os valores para os campos magnéticos
Ez e a outra (EquaHy . Estas equações
serão aplicadas de maneira intermitente de forma a avançar a simulação no tempo. As Figuras 2.5 e 2.6 ilustram o processo.
t Ey |t+1 i = Ey |i −
∆t Hz |ti −Hz |ti−1 ∆x
(2.47)
t Hz |t+1 i = Hz |i −
∆t Ey |ti+1 −Ey |ti µ∆x
(2.48)
13
Figura 2.5:
Visualização da Malha de Yee no espaço-tempo, o círculo vermelho
mostra o ponto onde o operador diferencial é discretizado para que se obtenha a Equação de Atualização para Campos Elétricos (Equação 2.43).
Figura 2.6:
Visualização da Malha de Yee no espaço-tempo, o círculo vermelho
mostra o ponto onde o operador diferencial é discretizado para que se obtenha a Equação de Atualização para Campos Magnéticos (Equação 2.44).
14
Adicionando
Jy (t)
como função forçante obtemos:
∆t t+1 ∆t Hz |ti − Hz |ti−1 + Jy |i ∆x ∆t Hz |t+1 = Hz |ti − Ey |ti+1 − Ey |ti i µ∆x
Ey |t+1 = Ey |ti − i
(2.49)
(2.50)
As Equações (2.49) e (2.50) são as Equações de Atualização do método FDTD para excitação com densidade de corrente. A adição de corrente não modi ca a lógica demonstrada nas Figuras 2.5 e 2.6, em cada iteração temporal do método, campos elétricos e magnéticos são recalculados utilizando estas equações e atualizados antes de avançar para a próxima iteração, este processo é conhecido como processo de marcha no tempo ou esquema
2.4.1
leapfrog.
Condições de Fronteira
O espaço discretizado pela simulação é, assim como a memória disponível para descrevê-lo, nito. Portanto, é fundamental determinar o comportamento das fronteiras da simulação.
Condições do tipo PEC Em alguns casos deseja-se determinar um espaço de simulação cercado por condutores perfeitos. Com isso, haverá re exão de ondas eletromagnéticas que atinjam seus limites. Este tipo de condição é útil no caso de, por exemplo, câmaras ressonantes. Fronteiras deste tipo são conhecidas como fronteiras PEC (do inglês:
tric Conductor ).
Perfect Elec-
De nir uma condição deste tipo é bem simples, só é necessário
determinar os campos elétricos tangentes às fronteiras como zero [9].
Condições Absorventes Em alguns casos se deseja simular uma malha de simulação in nita, isto é, deseja-se que ondas que atinjam a fronteira do espaço de simulação o abandonem de nitivamente, sem qualquer tipo de re exão. Criar um espaço de simulação muito maior do que o necessário em geral não é uma opção, pois exigiria uma quantidade não disponível de memória. A solução é, portanto, encerrar o espaço de simulação arti cialmente, usando condições que simulem uma absorção perfeita de energia. Em 1981, Gerrit Mur [19] de niu uma formulação para condições de fronteira que absorvem energia. Muitas outras formas se seguiram mas trataremos neste trabalho exclusivamente das chamadas condições de Mur de Primeira Ordem, que tem como base a Equação (2.51) [18]. condição ABC (do inglês:
Este tipo de condição de fronteira é conhecida como
Absorbing Boundary Conditions ). 15
sendo
v
∂ 1 ∂ ± ∂x v ∂t
Ey (x, y, t) = 0
(2.51)
a velocidade de propagação da onda eletromagnética no meio estudado.
A Equação (2.51) pode ter duas soluções, a saber:
1 ∂ ∂ + ∂x v ∂t
∂ 1 ∂ − ∂x v ∂t
Ey (x, t) = 0 → Ey (x, t) = f (x − vt)
(2.52)
Ey (x, t) = 0 → Ey (x, t) = f (x + vt)
(2.53)
Note que as soluções descrevem ondas que se propagam em somente uma direção, a Equação (2.51) é a chamada Equação da Onda Unidirecional e pode ser discretizada utilizando diferenças centrais. O desenvolvimento é feito, como exem-
x = 0, equivalente a i = 1. O ponto de 1 discretização espacial utilizado é normalmente i = , de modo a caracterizar o início 2
plo, para uma fronteira vertical presente em
da fronteira no ponto médio entre dois campos elétricos, o ponto de discretização temporal é
t =
1 , ou seja, também no ponto médio entre dois campos elétricos, 2
mas no tempo. Estas escolhas para realizar a aproximação por diferenças centrais resultam na Equação (2.54) [18]:
t Ey |t+1 i=1+1/2 − Ey |i=1+1/2
v∆t
t+1/2
t+1/2
Ey |i=2 − Ey |i=1 = ∆x
(2.54)
Como o campo elétrico não está efetivamente de nido para os pontos intermediários de tempo e espaço que foram escolhidos, uma interpolação linear é feita, resultando em:
t+1 t+1 t t t t Ey |t+1 Ey |t+1 i=2 +Ey |i=2 − Ey |i=1 +Ey |i=1 i=1 +Ey |i=2 − (Ey |i=1 +Ey |i=2 ) = 2v∆t 2∆x
(2.55)
Finalmente, rearranjando para o campo elétrico na fronteira, conseguimos chegar à Equação de Atualização (2.56):
Ey |t+1 i=1 =
Ey |ti=2
+
v∆t − ∆x v∆t + ∆x
t Ey |t+1 i=2 −Ey |i=1
(2.56)
Esta é uma formulação que funciona bem para situações onde se espera que qualquer incidência de onda eletromagnética ocorra de forma aproximadamente perpendicular às fronteiras da simulação. posteriores tais como o
Caso este não seja o caso, outros métodos
Perfectly-Matched Layer
16
(PML) são recomendados [8].
2.5
Estabilidade e a Condição Courant-FriedrichsLewy
Estabilidade é de nida como a característica de métodos numéricos que não possuem crescimento ilimitado de erro, sendo uma condição necessária mas não su ciente para que as simulações entreguem resultados precisos [20]. O método FDTD, como dito anteriormente, é um método explícito e está sujeito a instabilidade de acordo com os passos de discretização espaço-temporais. Métodos numéricos explícitos de resolução de alguns tipos de Equações Diferenciais Parciais (nas quais se encaixam as Equações de Maxwell) estão sujeitos à condição CourantFriedrichs-Lewy (ou condição CFL) para que se mantenham estáveis. Para o caso do método FDTD bidimensional descrito no plano
xy ,
a condição CFL se traduz na
Inequação (2.57) [7] [18].
1 ∆t ≤ q v ∆x1 2 +
(2.57)
1 ∆y 2
Esta condição para a relação entre discretização espacial e temporal pode ser deduzida utilizando o chamado Critério de Von Neumann, apresentado no Anexo A. Mais informação sobre esta linha de análise de estabilidade também pode ser encontrada nas referências [18] e [21]. Esta relação também pode ser justi cada por um argumento físico: a discretização da malha não pode ser feita de forma a violar a causalidade. Como em cada passo temporal informações sobre grandezas físicas passam de um ponto da malha para outro, a discretização tem que ser feita de tal forma que esta transmissão de informação ocorra em uma velocidade próxima da prevista analiticamente para ondas eletromagnéticas,
v =
√1 [14] [22]. Em outras palavras, a de nição do valor µ
assumido por um grandeza física qualquer em um ponto no espaço-tempo deve ser em função somente de pontos dos quais a própria solução analítica dependa.
Ou
seja, é natural esperar que a condição CFL especi que claramente que a máxima discretização espacial não pode ser superior a propagação da luz em um passo temporal, caso contrário o método FDTD estaria fazendo uso de informação ainda não disponível, de acordo com o modelo analítico. Em simulações unidimensionais a razão mostrada na Equação (2.58) é útil, esta é conhecida como Número de Courant.
Sc =
v ∆t ∆x
A Figura 2.7 mostra resultados para campo elétrico
(2.58)
Ey
de uma simulação uni-
dimensional da propagação de uma onda transversal eletromagnética utilizando o
17
método FDTD com [9]).
Sc = 1.0
(valor considerado ideal para este tipo de simulação
Note que, utilizando Número de Courant unitário, a onda eletromagnética
percorre a geometria sem sofrer distorções.
Distribuição de Ey ao longo da geometria em diferentes instantes 1 t = 0.93 ns t = 1.06 ns t = 1.19 ns t = 1.33 ns
0.8 0.6 0.4
Ey (V/m)
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
Figura 2.7: Simulação de Campo Elétrico
0.3
0.35
0.4
Ey utilizando método FDTD com Sc = 1.0.
A mesma simulação é agora apresentada para o caso onde a condição de Courant foi desrespeitada, ou seja:
Sc > 1.0.
Mais especi camente:
2.8 mostra a intensidade do campo elétrico
Ey
Sc = 1.005.
A Figura
para diferentes tempos de simula-
ção, de forma a ressaltar a instabilidade na simulação. Note que, com o passar do tempo de simulação, o erro cresce sem limite. Este comportamento, como já dito anteriormente, caracteriza a instabilidade numérica.
18
Distribuição de Ey ao longo da geometria em um certo instante.
Distribuição de Ey ao longo da geometria em um certo instante.
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2 Ey (V/m)
1
Ey (V/m)
1
0
0
−0.2
−0.2
−0.4
−0.4
−0.6
−0.6
−0.8
−0.8
−1
0
0.05
0.1
0.15
(a)
0.2 0.25 Eixo x (cm)
0.3
0.35
−1
0.4
0
0.05
0.1
0.15
(b)
t = 0.134 ns
Distribuição de E ao longo da geometria em um certo instante.
0.2 0.25 Eixo x (cm)
0.3
0.35
0.4
t = 0.335 ns
Distribuição de E ao longo da geometria em um certo instante.
y
y
1
2
0.8
1.5
0.6 1 0.4 0.5 Ey (V/m)
Ey (V/m)
0.2 0 −0.2
0
−0.5
−0.4 −1 −0.6 −1.5
−0.8 −1
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
(c)
t = 0.469 ns
0.3
0.35
−2
0.4
0
2
2
1.5
1.5
1
1
0.5
0.5
0
−0.5
−1
−1
−1.5
−1.5
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
(e)
t = 0.603 ns
0.3
0.35
0.15
0.2 0.25 Eixo x (cm)
0.3
(d)
t = 0.536 ns
0.35
0.4
0
−0.5
−2
0.1
Distribuição de Ey ao longo da geometria em um certo instante.
Ey (V/m)
Ey (V/m)
Distribuição de Ey ao longo da geometria em um certo instante.
0.05
−2
0.4
Figura 2.8: Simulação de Campo Elétrico
1.005.
19
Ey
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
(f)
t = 0.670 ns
0.3
0.35
0.4
utilizando método FDTD com
Sc =
Distribuição de Ey ao longo da geometria em um certo instante. 50
8
40
6
30
4
20
2
10 Ey (V/m)
Ey (V/m)
Distribuição de Ey ao longo da geometria em um certo instante. 10
0
0
−2
−10
−4
−20
−6
−30
−8
−40
−10
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
0.3
(g)
t = 0.804 ns
0.35
−50
0.4
Figura 2.8: Simulação de Campo Elétrico
1.005.
20
Ey
0
0.05
0.1
0.15
0.2 0.25 Eixo x (cm)
0.3
(h)
t = 0.938 ns
0.35
0.4
utilizando método FDTD com
Sc =
Capítulo 3 O Método FDTD baseado em Polinômios de Laguerre 3.1
Introdução
A partir dos anos 2000 surgiram novas formulações do método FDTD com o objetivo de eliminar,
ou reduzir,
a limitação imposta pela condição CFL. A
primeira destas versões foi o método FDTD Alternado-Direto Implícito, ou ADI-
Alternate-Direct Implicit FDTD ).
FDTD (do inglês:
Embora o método ADI seja
incondicionalmente estável, ou seja, não dependa da condição CFL para que o erro não cresça monotonicamente, ele ainda não garante precisão para qualquer razão entre a discretização espacial e temporal [11].
Em outras palavras, ainda haverá
erro numérico para razões de discretização que se distanciem da condição CFL, tudo que o método ADI-FDTD garante é que este erro não crescerá sem limites. Outras tentativas surgiram na última década (tais como o método
Cole's Non-Standard
FDTD ) mas um tratamento destas está além do objetivo deste trabalho.
Em um artigo publicado em 2003 por Chung, Sarkar, Jung e Salazar-Palma [11], uma versão incondicionalmente estável utilizando polinômios de Laguerre é sugerida. A ideia fundamental é abandonar o esquema de marcha no domínio do tempo mostrado no Capítulo 2 e adotar um outro domínio, que tem como bases os polinômios de Laguerre. Este método é abreviado na literatura como LaguerreFDTD ou WLP-FDTD (do inglês:
Weighted Laguerre Polinomials-FDTD)
21
3.2
Polinômios de Laguerre
A Equação Diferencial (3.1) é chamada a Equação de Laguerre.
x sendo
x≥0
e
k
∂ 2y ∂y + ky = 0 + (1 − x) 2 ∂x ∂x
(3.1)
um número natural, a solução em termos de séries in nitas desta
equação diferencial é descrita pela Equação (3.2).
k(k − 1) 2 k(k − 1)(k − 2) 3 x − x 22 22 · · · 32 k(k − 1)(k − 2)(k − 3) 4 (−1)n k! n x + x + ... 22 · · · 32 · · · 42 (n! )2 (k − n)!
yk (x) = a0 1 − kx +
Aplicar
k = 0, 1, 2, 3...
e
a0 = k!
(3.2)
na Equação (3.2) resulta na família de polinômios
mostrada na Equação (3.3)
L0 (x) = 1 L1 (x) = 1 − x L2 (x) = 2 − 4x + x2 L3 (x) = 6 − 18x + 9x2 − x3
(3.3)
... 2
Lk (x) = (k! )
k X n=0
(−1)n xn 2 (n! ) (k − n)!
Estes polinômios são denominados Polinômios de Laguerre e podem ser escritos de uma forma recursiva, ou seja, o polinômio
Ln ,
se
n > 2,
pode ser escrito em
função dos dois polinômios anteriores, como mostrado na Equação (3.4) [23].
L0 (x) = 1 L1 (x) = 1 − x Lk (x) =
(3.4)
(2k − 1 − x) Lk−1 (x) − (k − 1) Lk−2 (x) k
22
para
k≥2
A Figura 3.1 mostra os cinco primeiros polinômios de Laguerre.
20 Ordem 0 Ordem 1 Ordem 2 Ordem 3 Ordem 4
15
L(x)
10
5
0
−5
−10
0
5
10 x
15
20
Figura 3.1: Polinômios de Laguerre de ordem 0 a 4.
A partir daqui passaremos a tratar a variável da qual dependem os polinômios de Laguerre como a variável tempo,
3.3
t.
Ortogonalidade
Os conceitos vetoriais de produto interno escalar e ortogonalidade podem ser generalizados para funções. Consideremos dois vetores espaço vetorial
IR
2
u
e
v,
de nidos, por exemplo, no
. O produto interno entre esses dois vetores, representado
u·v,
é de nido pelas seguintes propriedades [24]: i. ii.
u·v
=
ku · v
v·u
=
k(u · v),
iii.
u·u=0
se
iv.
u+v·w
=
sendo
u=0 u·w
+
e
k
escalar.
u·u>0
v · w,
sendo
se
w
u 6= 0 vetor no
IR2 .
Uma integral de nida do produto entre duas funções
f1
e
f2
no intervalo
[a, b]
apresenta as propriedades (i-iv) apresentadas. Com isso, é possível de nir o produto interno entre funções como (3.5):
23
b
Z f1 · f2 =
f1 (t)f2 (t)dt
(3.5)
a As duas funções
f1
e
f2
serão ditas ortogonais se o resultado da integral em
(3.5) for zero.
Prosseguindo com a analogia entre vetores e funções: um conjunto formado por vetores mutualmente ortogonais forma uma base do espaço de vetores. Isso signi ca dizer que qualquer outro vetor neste mesmo espaço pode ser representando como uma combinação linear destes. Por exemplo, se mensional
IR
2
, qualquer vetor
u
e
v
formam base no espaço bidi-
w também contido neste espaço pode ser representado
pela combinação linear (3.6).
w = c1 u + c2 v
(3.6)
É possível fazer mais uma vez um paralelo com funções: um conjunto ortogonal in nito de funções
[Φ0 , Φ1 , Φ2 ...]
Φn
de nido no intervalo
t = [a, b], ou seja, um conjunto da forma
que respeita a condição (3.7).
Z
b
Φm (t) Φn (t) dt = δmn
(3.7)
a
sendo
δmn
o delta de Kronecker, de nido por:
δmn =
1,
se
m = n,
0,
se
m 6= n.
(3.8)
constitui uma base do chamado espaço de funções e, portanto, pode representar funções
f (t)
de nidas no intervalo
[a, b]
como mostrado na Equação (3.9):
f (t) = c0 Φ0 + c1 Φ1 + c2 Φ2 ... + cn Φn onde
cn
(3.9)
são coe cientes a serem determinados.
Alguns conjuntos de funções são ditos ortogonais em relação a uma função peso
w(t) > 0:
24
Z
b
w(t) Φm (t) Φn (t) dt = δmn
(3.10)
a
O conjunto de polinômios de Laguerre se encaixa no caso (3.10). Estas tem uma
e−t . Observando também t = [a, b] está de nido para este
relação de ortogonalidade para uma função peso igual a que o intervalo anteriormente representado por caso como
t = [0, ∞].
O conceito de ortogonalidade pode ser descrito, para o
caso especí co dos polinômios de Laguerre, pela Equação (3.11).
É interessante
observar que enquanto em análise vetorial é possível dar um sentido geométrico de perpendicularidade à palavra ortogonal, o mesmo não ocorre no presente caso, sendo esta uma simples de nição advinda da generalização do produto interno entre vetores.
Z
∞
e−t Lp Lq dt = δpq
(3.11)
0
Isto signi ca dizer que é possível de nir bases do espaço de funções utilizando polinômios de Laguerre, como mostrado na Equação (3.12). A Figura 3.2 mostra as cinco primeiras bases de Laguerre.
φp (t, s) = e−
25
s.t 2
Lp (s.t)
(3.12)
1 Ordem 0 Ordem 1 Ordem 2 Ordem 3 Ordem 4
φp (t)
0.5
0
−0.5 0
5
10 t
15
20
Figura 3.2: Bases de Laguerre até ordem 4, escala temporal
s = 1.
Importante notar que na Equação (3.12) zemos uso de um fator de escala temporal
s > 0.
Isso é justi cável ao observar a Figura 3.2, note que as bases de
Laguerre decaem a zero após a passagem de um tempo
t
su cientemente grande,
na ordem de segundos. Mas o comportamento das grandezas físicas em simulações não é xo e muito menos xo na escala de segundos.
Em algumas simulações as
variáveis de interessem podem tender a zero em escalas de tempo muito inferiores e, neste caso, os polinômios de Laguerre não seriam capazes de sintetizá-las. Isso signi ca que, para sintetizar sinais em escalas temporais diferentes, é necessário utilizar uma escala de modo a ajustar as bases ao problema especí co a ser estudado [10].
Em termos práticos, dado um número su cientemente grande de bases de Laguerre, é possível representar (ou sintetizar) qualquer função transitória como um somatório ponderado dessas bases, um conceito similar ao da Série de Fourier. A Figura 3.3 mostra a síntese do pulso
e−(
t−15 2 ) 5
fazendo uso de
5, 8
e
12
bases de
Laguerre. Fica evidente que com o número crescente de bases a descrição da função ca cada vez mais próxima da função original.
26
1
Sinal, f(t)
0.5
0 Sinal Original Sintese com 5 bases de Laguerre Sintese com 8 bases de Laguerre Sintese com 12 bases de Laguerre −0.5
0
5
10
15
20 t
25
30
35
40
Figura 3.3: Síntese de um sinal utilizando diferentes números de bases de Laguerre.
A Figura 3.4 mostra, por m, a síntese deste mesmo sinal utilizando
100
de Laguerre.
1
Sinal, f(t)
0.5
0
Sinal Original Sintese com 100 bases de Laguerre −0.5
0
5
10
15
20 t
25
30
35
40
Figura 3.4: Síntese de um sinal utilizando 100 bases de Laguerre.
27
bases
As escolhas para o número temporal
s
NL
de bases de Laguerre e para o valor da escala
a serem utilizados variam e se baseiam, respectivamente, nas Equações
(3.13) e (3.14) [25].
NL = 2Bf Tf + 1
(3.13)
sendo:
NL
o número necessário de bases de Laguerre.
Bf
a máxima frequência, em Hertz, presente na forma de onda a ser descrita.
Tf
a duração em segundos da forma de onda.
s=
L xN max Tf
(3.14)
sendo: L xN max
o maior zero dentre os polinômios de Laguerre utilizados.
Uma pergunta razoável é: Por que utilizar especi camente polinômios de Laguerre? Existem diversas outras classes de funções que apresentam ortogonalidade. Além dos já citados senos e cossenos poderíamos citar os polinômios de Chebyshev e Legendre. Alguns dos motivos estão apresentados a seguir [10]:
Os polinômios de Laguerre formam uma classe de funções ortogonais de nidas de 0 a
∞
(os polinômios de Chebyshev, por exemplo, estão de nidos de -1
a 1). Isto faz com que polinômios de Laguerre sejam uma boa escolha para simulação de sistemas físicos em função do tempo, pois é neste mesmo intervalo que o tempo está de nido.
Como a função peso dos polinômios de Laguerre é uma exponencial decrescente, todas as bases de Laguerre eventualmente decaem a zero. Para simulações transitórias, isto é ideal, pois garante estabilidade incondicional, ou seja, o método numérico jamais retornará um sinal que cresça sem limite e, portanto, jamais resultará em um erro que cresça sem limite.
A relação de recursividade entre os polinômios de Laguerre faz com que seja simples computacionalmente obtê-los em tempo de execução, sem necessidade de ocupar memória armazenando as bases a priori.
28
3.4
O Método FDTD Baseado em Polinômios de Laguerre
Utilizando as propriedades dos polinômios de Laguerre se torna possível representar grandezas físicas de comportamento transitório em termos de somas de bases de Laguerre ponderadas por coe cientes.
A Equação (3.15) mostra, como exemplo,
o caso do campo elétrico. A mesma equação se aplica para campos magnéticos e correntes elétricas.
Ex (r, t) =
∞ X
Exp (r) φp (t¯)
(3.15)
p=0 sendo:
Ex (r, t)
o campo elétrico na direção
variável tempo
t¯ = s.t r φp
x em um determinado ponto r,
função da
t.
simplesmente uma representação para o tempo escalado pelo fator
s.
a variável posição no espaço. a base de Laguerre de ordem
Exp (r)
p.
o coe ciente de ponderação da base Laguerre de ordem
sentação do campo elétrico
Ex
em uma posição especí ca
p
para a repre-
r.
É possível demonstrar que a derivada temporal de uma função
U
qualquer escrita
na forma da Equação (3.15), será escrita como na Equação (3.16) [11]:
p−1 ∞ X X ∂U (r, t) =s 0.5 Up (r) + Uk (r) φp (t¯) ∂t p=0 k=0,p>0
(3.16)
Partindo das Equações (3.15) e (3.16) é possível desenvolver o método para os diferentes modos já apresentados (TE e
3.4.1
TEM).
Modo TE
Primeiro, reapresentamos as equações que representam o Modo
∂Ex 1 ∂Hz = − Jx ∂t ∂y ∂Ey 1 ∂Hz = − − Jy ∂t ∂x 29
TEz : (3.17)
(3.18)
∂Hz 1 ∂Ex ∂Ey = − ∂t µ ∂y ∂x
(3.19)
O primeiro passo é utilizar as Equações (3.15) e (3.16) para reescrever as grandezas desconhecidas (campos elétricos e magnéticos) e suas derivadas temporais presentes em (3.17), (3.18) e (3.19) utilizando formulações em termos de polinômios de Laguerre.
s
∞ X
0.5Exp (x, y)
p−1 X
+
p=0
∞ X
∞ X p=0
(3.20)
0.5Eyp (x, y)
p−1 X
+
Eyk (x, y)
φp (t¯)
k=0
=−
(3.21)
∞ 1 X ∂H p
∞ X
p=0
Jy (x, y, t) z (x, y) φp (t¯) − ∂y
0.5Hzp (x, y)
+
p=0
1 =− µ
φp (t¯)
∂Hzp Jx (x, y, t) (x, y) φp (t¯) − ∂y
p=0
s
k=0
1 =
s
Exk (x, y)
p−1 X
Hzk (x, y)
φp (t¯)
k=0 ∞ X p=0
(3.22)
∂Eyp ∂Exp (x, y) − (x, y) φp (t¯) ∂y ∂x
É importante lembrar que o objetivo do método é abandonar o domínio do tempo e, no entanto, as bases
φp (t¯)
e as densidades de corrente ainda apresentam termos
com dependência em relação ao tempo. Para eliminar esta dependência utilizamos um procedimento conhecido como Procedimento de Teste Temporal de Galerkin [11] que consiste em, fazendo uso das propriedades ortogonais das bases de Laguerre, multiplicar os dois lados da Equação (3.20) por uma base especí ca em função de
t¯ no
intervalo
[0, Tf ),
sendo
Tf
φq
e integrar
um período de tempo su cientemente
grande, de modo a garantir que todas as formas de onda já tenham decaído a praticamente zero. Aplicar este procedimento elimina a dependência explicita em relação ao tempo das seguintes formas:
Como as bases são ortogonais entre si e o procedimento consiste em realizar o produto
φp φq
dentro de uma integral, todos os termos dos somatórios presentes
30
nas Equações (3.20), (3.21) e (3.22) com
desaparecem, e quando
q = p,
1.
o termo resulta simplesmente em
q 6= p
Os termos relativos a densidade de corrente passam a ser:
Jxq (x, y)
TF
Z
J(x) (x, y, t) φq (t¯) dt¯
= 0
Jyq (x, y)
Z
(3.23)
TF
J(y) (x, y, t) φq (t¯) dt¯
= 0
ou seja, determinam-se assim os termos do somatório de bases de Laguerre que sintetizam a densidade de corrente (que já é conhecida para todo o tempo de simulação).
Aplicando estes dois resultados do procedimento às Equações (3.20), (3.21) e (3.22) temos:
q−1 X J q (x, y) 1 ∂Hzq q k s 0.5Ex (x, y) + (x, y) − x Ex (x, y) = ∂y k=0
(3.24)
q−1 X Jyq (x, y) 1 ∂Hzq q k (x, y) − s 0.5Ey (x, y) + Ey (x, y) = − ∂y k=0
(3.25)
q−1 X ∂Eyq 1 ∂Exq q k s 0.5Hz (x, y) + Hz (x, y) = − (x, y) − (x, y) µ ∂y ∂x k=0
(3.26)
Feito o procedimento de passagem do domínio do tempo para o domínio das bases de Laguerre, faz-se necessário discretizar o espaço, o método WLP-FDTD não difere do método FDTD neste sentido e utiliza a mesma malha deslocada de Yee, ou seja, campos elétricos e magnéticos são avaliados em pontos diferentes, como já mostrado para o caso bidimensional na Figura 2.3.
Aplicando a discretização espacial em
diferenças centrais, seguida de alguma manipulação algébrica, nas Equações (3.24), (3.25) e (3.26) obtemos (3.27), (3.28) e (3.29).
Exq |i,j
Eyq |i,j
=
=
CyE
[Hzq |i,j −Hzq |i,j−1 ]
−CxE
q−1 X 2 q − Jx |i,j − 2 Exk |i,j s k=0
[Hzq |i,j −Hzq |i−1,j ]
31
q−1 X 2 q − Jy |i,j − 2 Eyk |i,j s k=0
(3.27)
(3.28)
Hzq |i,j
−CxH
=
Eyq |i+1,j −Eyq |i,j
+
CyH
q−1 X q q Ex |i,j+1 −Ey |i,j − 2 Hzk |i,j
(3.29)
k=0
sendo:
CxE =
2 . s ∆x
CyE =
2 . s ∆y
CxH =
2 . sµ ∆x
CyH =
2 . sµ ∆y
Observando as Equações (3.27), (3.28) e (3.29) é possível notar que existe uma relação implícita entre campos elétricos e magnéticos. Em outras palavras, é possível substituir a Equação (3.29) nas Equações (3.27) e (3.28), o que resultará no conjunto composto pelas Equações (3.30) e (3.31), apresentado a seguir:
1 H H + Cy + Cy + Exq |i,j−1 −CyH + Exq |i,j+1 −CyH + E Cy −CxH + Eyq |i,j−1 CxH + Eyq |i+1,j CxH + Eyq |i+1,j−1 −CxH
Exq |i,j Eyq |i,j
q−1
= −∆y Jxq |i,j −
(3.30)
q−1
X X 2 X k k E | + 2 H | − 2 Hzk |i,j i,j z i,j−1 CyE k=0 x k=0 k=0
1 H H q H q H + C + C + E | −C + E | −C + i−1,j i+1,j x x y x y x CxE −CyH + Exq |i−1,j CyH + Exq |i−1,j+1 −CyH + Exq |i,j+1 CyH
Eyq |i,j Exq |i,j
q−1
q−1
= −∆x Jyq |i,j −
q−1
(3.31)
q−1
X X 2 X k k E | − 2 H | + 2 Hzk |i,j i,j z i−1,j CxE k=0 y k=0 k=0
O método se baseia em descobrir os termos da série que de nem o campo elétrico na Equação (3.15) utilizando as Equações (3.30) e (3.31). Para isso, são necessárias tantas iterações quanto termos da série. Em outras palavras, serão necessárias tantas iterações quanto bases de Laguerre utilizadas para sintetizar as grandezas físicas. A relação descrita pelas Equações (3.30) e (3.31) valerá para todos os pontos da geometria, exceto os pontos nas fronteiras, como será mostrado adiante.
32
Isto dá
origem a um sistema linear como o mostrado na Equação (3.32). Este sistema será resolvido em cada uma das iterações.
Cada uma das linhas da matriz descreve a
equação para os campo elétricos de cada ponto da geometria. No caso do modo
TE,
existem dois campos elétricos a serem determinados para cada ponto da geometria, por isso, a matriz terá como número de linhas o dobro do número de pontos utilizados para discretizar a geometria.
[A]2m × 2m [E q ]2m × 1 = [J q ]2m × 1 + [β q−1 ]2m × 1
para
q = 1, 2...
(3.32)
sendo:
m,
o número de pontos utilizados para discretizar a geometria.
[A]2m × 2m ,
a matriz dos coe cientes que multiplicam os termos de campo elé-
trico.
[E q ]2m × 1
o vetor coluna contendo todos os termos de Laguerre para os campos
elétricos.
[J q ]2m × 1
o vetor coluna contendo os termos de Laguerre para todas as corren-
tes.
[β q−1 ]2m × 1
o vetor coluna contendo os somatórios dos termos de campos elé-
tricos e magnéticos de todas as iterações anteriores. Portanto, este vetor não existe para
q = 0.
É importante observar que os termos do campo elétrico da iteração atual
q
cam
dependendo somente de termos já descobertos dos campos magnéticos. Descobrir os próprios campos magnéticos da iteração atual não é um problema, isto pode ser feito em termos de campos elétricos através da Equação (3.29). Também não há problema no fato dos termos do campo elétrico da iteração atual
Eq
terem dependência em
relação aos termos de densidade de corrente de mesma iteração
J q,
pois, no método
WLP-FDTD, a densidade de corrente atua como conhecida a priori e, obviamente, todos os seus termos já são também conhecidos. Como não existem termos anteriores na primeira iteração,
0
0
[A] [E ] = [J ].
[β q−1 ]
é inexistente e a Equação (3.32) se torna simplesmente
Uma observação nal é que a matriz
A
só precisa ser calculada
uma vez no início do problema, a partir daí ela se mantém a mesma para todas as iterações.
Isto permite a utilização de métodos de aceleração de resolução do
sistema linear, tais quais fatoração LU [11].
33
3.4.2
Modo TEM
O modo transversal eletromagnético tem um desenvolvimento similar ao do modo transversal elétrico mostrado anteriormente, começamos reapresentando o modo
TEMx polar. y
e suas Equações, (3.33) e (3.34).
∂Ey 1 ∂Hz = − Jy ∂t ∂x
(3.33)
∂Hz 1 ∂Ey = ∂t µ ∂x
(3.34)
Aplicando as Equações (3.15) e (3.16) obtemos as Equações (3.35) e (3.36):
s
∞ X
0.5Eyp (x)
+
p=0
s
q−1 X
Eyk (x)
k=0 ∞ X
0.5Hzp (x)
+
p=0
q−1 X
1 φp (t¯) = −
Hzk (x)
k=0
∞ X ∂H p
Jy (x, t) (x) φp (t¯) − ∂x z
p=0
1 φp (t¯) = − µ
∞ X ∂Eyp p=0
∂x
(x) φp (t¯)
(3.35)
(3.36)
Aplicando o Procedimento de Teste Temporal de Galerkin:
q−1 X Jyq (x) 1 ∂Hzq q k s 0.5Ey (x) + Ey (x) = − (x) − ∂x k=0
(3.37)
q−1 X 1 ∂Eyq q k (x) s 0.5Hz (x) + Hz (x) = − µ ∂x k=0
(3.38)
Aplicando a discretização espacial por diferenças centradas e rearranjando termos camos com as Equações (3.39) e (3.40):
Eyq |i
=
−CxE
Hzq |i
=
(Hzq |i
−CxH
−
Hzq |i−1 )
Eyq |i+1
−
q−1 X 2 q Eyk |i − Jy |i − 2 s k=0
Eyq |i
− 2
q−1 X k=0
sendo:
CxE =
2 . s ∆x
CxH =
2 . sµ ∆x
34
Hzk |i
(3.39)
(3.40)
Assim como no caso do modo TE, existe uma relação implícita entre campos elétricos e magnéticos. Substituindo a Equação (3.40) em (3.39) se torna possível escrever:
1 H H + Cx + Cx + Eyq |i−1 −CxH + Eyq |i+1 −CxH E Cx q−1 q−1 q−1 X X 2 X k q k = −∆x Jy |i − E Ey |i − 2 Hz |i−1 + 2 Hzk |i,j Cx k=0 k=0 k=0
Eyq |i
(3.41)
Assim como no caso anterior, chega-se a um sistema linear. O sistema para o
TEMx polar. y
modo
é da forma mostrada na Equação (3.42).
[A]m × m [E q ]m × 1 = [J q ]m × 1 + [β q−1 ]m × 1
para
q = 1, 2...
(3.42)
sendo:
m,
o número de pontos utilizados para discretizar a geometria.
[A]m × m , a matriz dos coe cientes que multiplicam os termos de campo elétrico. [E q ]m × 1
o vetor coluna contendo todos os termos de Laguerre para os campos
elétricos.
[J q ]m × 1 o vetor coluna contendo os termos de Laguerre para todas as correntes. [β q−1 ]m × 1
o vetor coluna contendo os somatórios dos termos de campos elétri-
cos e magnéticos de todas as iterações anteriores. Este vetor não existe para
q = 0. Observe que nesse caso, ao contrário do modo TE, só existe um campo elétrico a ser determinado para cada ponto da geometria.
Com isso, enquanto o sistema
de um problema bidimensional tem tantas linhas quanto o dobro de pontos usados pra discretizar a geometria, o sistema linear de um problema unidimensional terá simplesmente tantas linhas quanto pontos usados pra discretizar a geometria. Todas as observações feitas anteriormente sobre a relação entre iterações e ordens dos polinômios, além da estrutura do algoritmo, se mantém.
3.4.3
Condições de Fronteira
As condições de fronteira apresentadas na Seção 2.4.1 também podem ser implementadas no método WLP-FDTD. Estas condições são adicionadas exatamente da mesma forma tanto para problemas bidimensionais (modos unidimensionais (modo
TEM). 35
TE
e
TM)
quanto para
Condição PEC (Perfect
Electric Conductor )
A condição do tipo PEC é adicionada de maneira razoavelmente simples modi cando o sistema linear (3.32) em dois passos :
As linhas da matriz
[A]2m × 2m
correspondentes aos pontos presentes no limite
da geometria devem ter, excetuando-se a diagonal principal, somente elementos nulos.
As linhas dos vetores
[β q−1 ]2m × 1
e
[J q ]2m × 1
devem ser zeradas.
Condição Absorvente de Mur (ABC) A condição de Mur de primeira ordem para limites laterais (por exemplo,
x=L
para uma geometria de nida no intervalo
[0, L])
x=0
ou
é dada pela Equação (2.51),
apresentada na Seção 2.4.1 [11]. Note que somente a componente do campo elétrico tangente à fronteira é tratada (no caso,
Ey ).
Para inserir esta condição no algoritmo do método WLP-FDTD é necessário escrevê-la em termos das bases de Laguerre e, para isso, as Equações (3.15) e (3.16) são utilizadas. Substituindo estas em (2.51) obtemos a Equação (3.43):
∂Eyq (x, y) s + ∂x v
! q−1 Eyq (x, y) X k + Ey (x, y) = 0 2 k=0
Aplicando a discretização espacial por diferenças centrais no ponto
(3.43)
x = 0 (Equa-
ção (3.44)) e, fazendo alguma manipulação algébrica, obtém-se a condição de Mur de primeira ordem aplicada ao WLP-FDTD para o ponto
x = 0,
mostrada na Equação
(3.45):
Eyq |2,j −Eyq |1,j ∂Ey |1+ 1 ,j = 2 ∂x 2
Eyq |1,j
(3.44)
q−1 s 1 s 1 s X k q + + Ey |2,j − = − Ey |2,j + Eyk |1,j 4v ∆x 4v ∆x 2v k=0
Aplicando o mesmo procedimento no outro extremo da geometria (x
(3.45)
= L),
ob-
temos a relação dada pela Equação (3.46), semelhante à Equação (3.45).
Eyq |L,j
q−1 s 1 s 1 s X k q + + Ey |L−1,j − = − Ey |L,j + Eyk |L−1,j 4v ∆x 4v ∆x 2v k=0
(3.46)
As Equações (3.45) e (3.46) podem ser inseridas nas linhas referentes às fron-
36
teiras do sistema linear da Equação (3.32). O sistema linear terá sua matriz seu vetor coluna
β
A
e
modi cados nas linhas referentes aos extremos da geometria.
O desenvolvimento para fronteiras verticais se dá da mesma forma, apenas se faz necessária a substituição da coordenada utilizada na derivada espacial da condição de Mur (x por
y ),
além da utilização da componente tangente de campo elétrico
correspondente.
37
Capítulo 4 Simulações e Resultados 4.1
Metodologia
Este trabalho tem, como dito anteriormente, o objetivo de validar as características de estabilidade incondicional e precisão do método WLP-FDTD. Para isso, foram desenvolvidos algoritmos próprios utilizando linguagem MATLAB para os casos
TEMx polar. y
e
TEz .
Utilizando estes algoritmos, foram realizadas simulações de
dois sistemas físicos, ambos com respostas eletromagnéticas presentes na literatura. A validação é feita, para o caso
TEz ,
por comparação entre o grá co obtido pelo
algoritmo próprio e o presente na referência [11]. Para o caso
TEMx polar. y ,
a validação é feita por comparação com o resultado
analítico disponível na referência [1]. Grá cos baseados na solução analítica foram gerados para comparação com os grá cos obtidos pelo algoritmo próprio. disso, o coe ciente de correção linear
RP earson
entre as duas respostas é obtido. Este
coe ciente está limitado entre -1 e 1, inclusive. positiva completa entre as grandezas, completa e
RP earson = 0
Além
RP earson
RP earson = 1 denota correlação = −1 denota correlação negativa
denota a completa ausência de correlação [26]. O Anexo B
trata brevemente deste conceito.
4.2
Estrutura do Algoritmo
Tanto o algoritmo unidimensional quanto o bidimensional se baseiam, fundamentalmente, em quatro
loops :
Loop de de nição das bases de Laguerre (NL iterações):
3.4 e 3.12 são utilizadas para gerar as bases de Laguerre. são obtidos os coe cientes de Laguerre escolhida, utilizando a Equação (3.23).
38
Jxq
e
Jyq
loop, as equações Neste mesmo loop já
Neste
para a densidade de corrente
Loop
A (2m
de de nição da matriz
iterações para casos bidimensionais,
m
iterações para casos unidimensionais. Sendo m o número de pontos escolhido para discretizar a geometria): Durante estas iterações, as equações (3.30) e (3.31) (para o caso
TEz )
ou a Equação (3.41) (para o caso
utilizadas para construir a matriz (3.32) (para o caso Existem
TEz )
ou na
TEMx polar. y )
são
A do sistema linear apresentado na Equação Equação (3.42) (para o caso TEMx polar. y ).
loops internos para tratar os extremos do espaço de simulação e aplicar
as condições de fronteira.
Loop
principal (NL iterações): Neste
o vetor
β,
loop,
três processos ocorrem; Primeiro,
dependente do somatório dos termos já descobertos, é atualizado
com os coe cientes da última iteração, se esta é a primeira iteração, este vetor simplesmente não existe.
Jxq e
Em seguida, utilizando o vetor
Jyq descobertos no primeiro
Loop,
e a matriz
A
Loop, o sistema linear ((3.32) ou (3.42)) é resolvido. coe cientes
E
q
os termos
construída no segundo
Em seguida, utilizando os
e a Equação (3.29) (bidimensional), ou a Equação (3.40) (uni-
dimensional), os coe cientes
Hq
são descobertos. Estes termos são utilizados
na próxima iteração, de modo a atualizar o vetor
β,
β.
Loop de retorno ao domínio do tempo (2m iterações para o caso bidimensional, m iterações para o caso unidimensional): Neste loop nal, os termos descobertos no Loop principal são aplicados na Equação (3.15) para todos os pontos da geometria, de modo a fazer com que as grandezas físicas retornem ao domínio do tempo.
A Figura 4.1 mostra um uxograma com a visão geral do algoritmo implementado neste trabalho para o algoritmo WLP-FDTD.
39
Loop de definição das bases de Laguerre
Loop de retorno ao domínio do tempo
Loop de definição da Matriz A
Loop principal
Gerador de Gráficos
Dados de Geometria e Discretização
Figura 4.1: Fluxograma do Algoritmo Implementado para o Método WLP-FDTD.
4.3
Simulação Bidimensional: Sistema Físico Apresentado por Chung et al.
Para validação do algoritmo próprio para o caso bidimensional, decidiu-se replicar os resultados apresentados na referência [11]. A con guração simulada é a mostrada na Figura 4.2.
40
Camada de Material Condutor
Camada de Material Condutor
Figura 4.2: Simulação bidimensional para validação do método.
Como mostrado na Figura 4.2, uma corrente elétrica uindo na direção
y
é es-
tabelecida entre duas películas de material condutor separadas por uma camada de ar.
y (Ey ) é medida exatamente no ponto p1. Esta con guração exige a
A componente de campo elétrico na direção
centro da geometria, denominada aqui como utilização do Modo
TEz
e dos dois tipos de condição de fronteira apresentados ante-
riormente: condições absorventes (nas extremidades laterais) e condições PEC (nas extremidades superior e inferior). A geometria tem e
0.1 m
na direção da coordenada
y.
1 m na direção da coordenada x
A discretização foi escolhida de modo a fazer
com que a condição de Courant (Equação (2.57)) fosse desrespeitada. Dessa forma, decidiu-se por
∆x = 0.01, ∆y = 0.001 m
e
∆t = 0.0293 ns.
A forma de onda da densidade de corrente utilizada é a mostrada na Figura 4.3 e descrita pela Equação (4.1). Os valores de do artigo e são
s = 6.07 × 1010
e
NL = 150.
x = 0.3 m.
41
s
e
NL
foram replicados diretamente
A corrente elétrica é estabelecida em
Forma de onda de J (t) y
1 0.8 0.6 0.4
A/m2
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
0
0.2
0.4
0.6 0.8 Tempo (segundos)
1
1.2 −8
x 10
Figura 4.3: Forma da densidade de corrente elétrica inserida no problema bidimensional.
J(t) = e Sendo
fc = 1 GHz , Td =
1 2fc
2 c − t−T T d
= 50 ns
e
sin(2πfc (t − Tc ))
(4.1)
Tc = 3Td = 150 ns.
O código do algoritmo próprio escrito em linguagem MATLAB para esta simulação é apresentado no Anexo C.
4.3.1
Resultados
O campo elétrico
Ey
no ponto
p1
da Figura 4.2 foi obtido através de simulação
utilizando o algoritmo desenvolvido para o método WLP-FDTD. A forma obtida é apresentada na Figura 4.4.
42
E (t) no ponto p1 y
2 b 1.5
1
d
V/m
0.5
0
−0.5 a −1
−1.5 c −2
0
0.2
0.4
0.6 0.8 Tempo (segundos)
1
1.2 −8
x 10
Figura 4.4: Resultado da simulação bidimensional mostrada na Figura 4.2.
Os pontos
a, b , c
e
d
foram tomados como referências visuais para comparação
entre os valores apresentados na Figura 4.4 e os apresentados na referência [11]. A Tabela 4.1 apresenta os valores de tempo e intensidade da componente de campo elétrico referência.
Ey
nos pontos
a, b , c
e
d
obtidos tanto neste trabalho quanto na
Os valores retirados do artigo são uma estimativa visual e, portanto,
têm uma incerteza associada de 0.1
V /m para o campo elétrico e 0.2 ns para o tempo.
43
Tabela 4.1: Comparação entre os valores obtidos por Chung et al. e os obtidos por algoritmo próprio, para os pontos
Ponto
a b c d
a, b , c
e
d.
Grandeza
Chung et al.
Algoritmo Próprio
Ey
−0.3 ± 0.1 V /m
−0.28 V /m
t
1.5 ± 0.2 ns
1.5 ns
Ey
1.5 ± 0.1 V /m
1.53 V /m
t
2.0 ± 0.2 ns
1.93 ns
Ey
−1.5 ± 0.1 V /m
−1.54 V /m
t
2.4 ± 0.2 ns
2.35 ns
Ey
0.3 ± 0.1 V /m
0.28 V /m
t
2.9 ± 0.2 ns
2.8 ns
A Figura 4.5, apresentada na próxima página, mostra uma sequência de pontos de tempo de simulação. Nesta gura é possível observar a propagação da componente
Ey
do campo elétrico causada pela passagem de corrente.
44
Ey(t)
Ey(t)
V/m
9
V/m
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(a)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(b)
t = 0.878 ns E (t)
9
80
90
t = 1.17 ns E (t)
V/m
y
70
V/m
y
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(c)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(d)
t = 1.46 ns Ey(t)
80
90
t = 1.75 ns Ey(t)
V/m
9
70
V/m
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(e)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(f)
t = 2.05 ns
70
80
90
t = 2.34 ns
Figura 4.5: Resultados para a simulação bidimensional utilizando o método WLPFDTD.
45
Ey(t)
Ey(t)
V/m
9
V/m
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(g)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(h)
t = 2.63 ns Ey(t)
80
90
t = 3.22 ns Ey(t)
V/m
9
70
V/m
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(i)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(j)
t = 3.51 ns Ey(t)
80
90
t = 3.80 ns Ey(t)
V/m
9
70
V/m
9
3
8
3
8 2
2
7
7 1
5
1
6 Eixo y (cm)
Eixo y (cm)
6
0
4
5
0
4 −1
−1
3
3 −2
−2
2
2 −3
1
−3
1
0
0 10
20
30
40 50 60 Eixo x (cm)
(k)
70
80
90
10
20
30
40 50 60 Eixo x (cm)
(l)
t = 4.10 ns
70
80
90
t = 4.39 ns
Figura 4.5: Resultados para a simulação bidimensional utilizando o método WLPFDTD.
46
Dos resultados apresentados é possível concluir que, primeiro, o método se manteve numericamente estável mesmo atuando fora da condição CFL. Além disso, observando a Tabela 4.1, é possível a rmar que os resultados apresentados estão em concordância com os encontrados na referência [11].
Por m, as condições de
fronteira do tipo Mur de Primeira Ordem implementadas no método WLP-FDTD se comportaram como o previsto, as ondas eletromagnéticas incidentes nas fronteiras são absorvidas e retiradas do espaço de simulação sem re exão perceptível.
4.4
Simulação
Unidimensional:
Propagação
de
Onda Transversal Eletromagnética Como dito anteriormente na Seção 2.3.1, uma das suposições necessárias para que se chegue a uma simulação unidimensional é a de que não exista variação da geometria em relação a dois dos eixos cartesianos.
Mais especi camente, no caso modelado
pelas Equações (2.33) e (2.34), não há variação em relação às coordenadas
y
e
z.
Dentro dessas condições, caso a excitação do sistema físico discretizado seja feita com densidade de corrente, esta situação equivale à de uma densidade super cial de corrente
~ K(t)
que se estende in nitamente no plano yz. No caso a ser tratado a
partir de agora utilizaremos o modo
TEM
polarizado na direção
dizer que esta densidade super cial estará na direção 4.6.
47
y,
y,
o que signi ca
como mostrado na Figura
Figura 4.6: Estrutura do problema unidimensional.
Esta con guração dá origem, como mostrado na Figura 4.6, a propagação de onda eletromagnética, mais especi camente ao modo de propagação transversal eletromagnético, ou simplesmente onda transversal eletromagnética.
A solução ana-
lítica para os valores de campos elétricos e magnéticos causados por tal densidade super cial é dada pela Equação (4.2) [1]:
Ey (t) = ∓c µHz (t) = − sendo
K
Ky (t ± xc ) 20 c
a densidade super cial de corrente, dada no SI por
Os termos
±
(4.2)
A . m
da Equação (4.2) são decididos de acordo com o sentido de propa-
gação da onda. Caso esta esteja se propagando no sentido positivo da direção x, a Equação (4.2) se torna:
48
Ky (t − xc ) Ey (t) = c µHz (t) = − 20 c
(4.3)
Caso a propagação seja no sentido negativo:
Ky (t + xc ) 20 c
(4.4)
O desenvolvimento do método WLP-FDTD para o caso
TEM polarizado na dire-
Ey (t) = −c µHz (t) = −
ção
y já foi demonstrado na Seção 3.4.2.
O algoritmo próprio desenvolvido utilizando
este desenvolvimento é aplicado à simulação de uma geometria unidimensional, consistindo de vácuo, com dimensão de
10cm e limitada nos dois extremos por fronteiras
do tipo ABC. Um plano in nito percorrido por uma densidade super cial de corrente
~ K
na direção
y
é colocado em
de corrente super cial
~ K
x = 5 cm.
A forma escolhida para a densidade
é a descrita pela Equação (4.5) e mostrada na Figura 4.7.
Forma de onda de Ky(t) 1.2
1
0.8
A/m
0.6
0.4
0.2
0
−0.2
0
0.5
1
1.5 2 Tempo (segundos)
2.5
3 −7
x 10
Figura 4.7: Forma de Onda da Densidade Super cial de Corrente
49
Ky .
J =e
t−5Tk 2 − T
(4.5)
k
Tk = 10 ns.
sendo
O número de bases de Laguerre
NL
necessário para sintetizar esta forma de onda
s foi obtido utilizando a 9 Equação (3.14). Os seguintes resultados foram obtidos: s = 1.13 × 10 e NL = 91. A discretização foi feita da seguinte maneira: ∆t = 0.037 ns e ∆x = 0.01 m. foi obtido utilizando a Equação (3.13), e o fator de escala
Discretizações estas que também desrespeitam a condição de Courant, o que pode ser visto facilmente pela Equação (2.58).
1.11 > 1.0.
Utilizando estes valores, obtemos
Sc =
O código do algoritmo próprio escrito em linguagem MATLAB para
esta simulação é apresentado no Anexo D.
4.4.1
Resultados
O algoritmo para simulação do caso unidimensional foi desenvolvido e seus resultados comparados com a solução analítica dada pela Equação (4.2). A Figura 4.8 mostra, em sequência temporal, os valores de campo elétrico
Ey
ao longo da geometria
previstos pela referência [1] sobrepostos aos valores retornados pelo método WLPFDTD. Distribuição do campo eletrico Ey ao longo da geometria
Distribuição do campo eletrico Ey ao longo da geometria
0
0
−0.5 −5 −1 −10
E (V/m)
−2
−15
y
Ey (V/m)
−1.5
−2.5
−20
−3 −25 −3.5 Feynman et. al WLP−FDTD −4 0
2
4
6
8
Feynman et. al WLP−FDTD
−30 10
0
2
4
Eixo x (cm)
(a)
6
8
10
Eixo x (cm)
(b)
t = 30.0 ns
t = 36.0 ns
Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando o método WLP-FDTD e dos apresentados por Feynman et al. [1].
50
Distribuição do campo eletrico Ey ao longo da geometria
Distribuição do campo eletrico Ey ao longo da geometria
0
0
−10
−20
−20
−40
−30
−60 −80
−50
E (V/m)
−60
−100
y
Ey (V/m)
−40
−120
−70 −140 −80 −160
−90
−180
−100
Feynman et. al WLP−FDTD
−110 0
2
4
6
8
Feynman et. al WLP−FDTD
−200 10
0
2
4
Eixo x (cm)
(c)
(d)
t = 42.0 ns
−20
−20
−40
−40
−60
−60
−80
−80
−100
−100
y
−120
−140
−160
−160
−180
−180
2
4
6
8
t = 48.0 ns
−200
Feynman et. al WLP−FDTD 0
Feynman et. al WLP−FDTD
−220 10
0
2
4
Eixo x (cm)
(e)
10
−120
−140
−220
8
Distribuição do campo eletrico Ey ao longo da geometria 0
E (V/m)
Ey (V/m)
Distribuição do campo eletrico Ey ao longo da geometria 0
−200
6 Eixo x (cm)
6
8
10
Eixo x (cm)
(f)
t = 54.0 ns
Distribuição do campo eletrico Ey ao longo da geometria
t = 60.0 ns
Distribuição do campo eletrico Ey ao longo da geometria
0
0
−20
−20
−40 −40 −60 −60
−100
E (V/m)
−120
−80
y
Ey (V/m)
−80
−100 −140 −120
−160 −180
−140
−200
Feynman et. al WLP−FDTD
Feynman et. al WLP−FDTD
−160
−220 0
2
4
6
8
10
0
2
4
Eixo x (cm)
(g)
6
8
10
Eixo x (cm)
(h)
t = 66.0 ns
t = 72.0 ns
Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando o método WLP-FDTD e dos apresentados por Feynman et al. [1].
51
Distribuição do campo eletrico Ey ao longo da geometria
Distribuição do campo eletrico Ey ao longo da geometria
0
0 −1
−10
−2 −3
−20
E (V/m)
−30
y
Ey (V/m)
−4 −5 −6 −7
−40
−8 −50
−9
−60
Feynman et. al WLP−FDTD
−10
Feynman et. al WLP−FDTD
−11 0
2
4
6
8
10
0
2
4
Eixo x (cm)
(i)
6
8
10
Eixo x (cm)
(j)
t = 78.0 ns
t = 84.0 ns
Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando o método WLP-FDTD e dos apresentados por Feynman et al.[1].
A Tabela 4.2 mostra, para cada um dos instantes mostrados na Figura 4.8, o maior valor absoluto de erro (denonimado pontos da geometria, o ponto (em módulo)
|Ey pico |
xerromax
|erro|max )
encontrado dentre todos os
no qual este erro ocorre, o valor de pico
do campo elétrico retornado pelo método WLP-FDTD dentre
todos os pontos da geometria e o coe ciente de correlação linear
RP earson
entre as
duas respostas.
Tabela 4.2: Comparação entre os valores analíticos e os obtidos pelo método WLPFDTD para campo elétrico no sistema físico unidimensional.
Instante
|erro|max
|Ey pico |
RP earson
t = 30.0 ns
0.025 V /m
10 cm
3.38 V /m
1.000
t = 36.0 ns
0.075 V /m
10 cm
26.36 V /m
1.000
t = 42.0 ns
0.052 V /m
7.24 cm
98.18 V /m
1.000
t = 48.0 ns
0.034 V /m
8.18 cm
180.72 V /m
1.000
t = 54.0 ns
0.103 V /m
0.01 cm
188.45 V /m
1.000
t = 60.0 ns
0.213 V /m
0.01 cm
188.45 V /m
1.000
t = 66.0 ns
0.518 V /m
0.01 cm
188.83 V /m
1.000
t = 72.0 ns
0.949 V /m
0.01 cm
142.76 V /m
1.000
t = 78.0 ns
1.240 V /m
0.01 cm
53.44 V /m
0.9996
t = 84.0 ns
1.358 V /m
0.01 cm
10.71 V /m
0.9782
Ponto
xerromax
Como não há crescimento monotônico do erro, pode-se a rmar que a simulação não resultou em instabilidade numérica mesmo com discretizações que desrespeitam
52
a condição CFL. É importante observar, no entanto, que há um crescimento do erro absoluto máximo com o decorrer da simulação exatamente em uma das fronteiras da simulação e, consequentemente, uma queda no coe ciente de correlação linear entre as duas respostas. A hipótese aqui adotada é a de que, embora a fronteira utilizando condição ABC do tipo Mur de Primeira Ordem tenha funcionado razoavelmente bem e tenha efetivamente absorvido a onda incidente e a retirado do espaço de simulação, este erro crescente a partir de
t = 54.0 ns
em uma das extremidades poderia ser
mitigado com o uso de outras formulações mais recentes para a condição de fronteira ABC, como por exemplo a formulação PML [8]. Outro ponto a ser observado é que a discretização feita resulta em uma quantidade par de pontos na geometria, ou seja, se torna impossível impor a densidade de corrente verdadeiramente no centro da simulação.
Esta é uma explicação possível para o surgimento assimétrico do
erro. Por m, a conclusão é a de que houve boa concordância entre a simulação e os resultados analíticos.
53
Capítulo 5 Conclusão e Trabalhos Futuros 5.1
Conclusão
Este trabalho apresentou a formulação teórica que permite a construção de um método FDTD baseado em Polinômios de Laguerre, tanto para casos bidimensionais quanto unidimensionais. Utilizando esta formulação, algoritmos foram escritos para dois casos, um deles unidimensional e o outro bidimensional. Os resultados obtidos foram comparados com o resultado derivado analiticamente para o caso unidimensional e com o publicado na literatura para o caso bidimensional. Os resultados se revelaram próximos e, portanto, é possível concluir que o algoritmo desenvolvido utilizando o método WLP-FDTD de fato apresenta estabilidade incondicional e entrega resultados precisos para simulações eletromagnéticas mesmo com passos de discretização que não respeitam a condição CFL. Isto faz com que o método seja uma opção para simulações FDTD que atualmente sejam feitas com discretização desnecessariamente pequena motivada somente por cuidados com estabilidade, o que resulta em simulações custosas computacionalmente.
5.2
Trabalhos Futuros
Sendo este um trabalho mais próximo da vertente teórica, existem diversos trabalhos futuros que podem ser realizados levando-se em conta uma abordagem mais prática: Pode-se tentar utilizar o método WLP-FDTD em um caso prático de engenharia elétrica, particularmente problemas de transitórios eletromagnéticos para os quais existam simetrias. É possível também desenvolver a formulação deste método para as Equações das Linhas de Transmissão. Também é possível continuar na vertente teórica e estudar trabalhos mais recentes que aperfeiçoam o método WLP-FDTD, expandindo-o para simulações de longa duração no domínio do tempo e o integrando com a Análise Nodal Modi -
54
cada [10]. Outra opção é explorar as outras versões incondicionalmente estáveis do método FDTD. Por m, também há a possibilidade de se tentar expandir o método WLP-FDTD para simulações eletrotérmicas.
55
Referências Bibliográ cas [1] FEYNMAN, R., LEIGHTON, R., SANDS, M.
sics, v. 2: [2] KASAL, R.
The Feynman Lectures on Phy-
Mainly Electromagnetics and Matter. Addison-Wesley, 1964.
Simulação de Supercondutores pelo Modelo do Estado Crítico.
Dis-
sertação de mestrado, COPPE/UFRJ, 2006.
Modelagem de Supercondutores Aplicada ao Projeto de Mancais Magnéticos. Tese de doutorado, COPPE/UFRJ, 2007.
[3] SOTELO, G.
Modelagem do Comportamento de Mancais Magnéticos Utilizando Fitas e Blocos Maciços Supercondutores. Tese de doutorado, COPPE/UFRJ,
[4] SASS, F.
2015.
Simulações de Regime Transitório de Limitadores de Corrente de Curto-Circuito Supercondutores. Tese de doutorado, COPPE/UFRJ,
[5] DE SOUSA, W.
2015.
Simulação de Sistemas Eletromagnéticos usando Método de Elementos Finitos e Método de Diferenças Finitas no Tempo.
[6] OLIVEIRA SANTOS, B.
Trabalho de conclusão de curso, Universidade Federal do Rio de Janeiro, 2017.
Computational Electrodynamics: The FiniteDierence Time-Domain Method. Artech House, 2000.
[7] TAFLOVE, A., HAGNESS, S.
Implementação do Método FDTD para as Equações de Maxwell em Ambientes de Programação Paralela. Dissertação de mestrado, Uni-
[8] VELOSO, L.
versidade Federal do Rio de Janeiro, 2015. [9] SCHNEIDER, J.
Understanding the Finite-Dierence Time-Domain Method.
2015.
EM Simulation Using the Laguerre-FDTD Scheme for Multiscale 3D Interconnections. Tese de doutorado, Georgia Institute of Technology,
[10] HA, M.
2011.
56
[11] CHUNG, Y., SARKAR, T., JUNG, B., et al.
An Unconditionally Stable
IEEE Transactions on Microwave Theory and Techniques, VOL 51. NO. 3, 2003. Scheme for the Finite-Dierence Time-Domain Method,
[12] DE SOUSA, W., POLASEK, A., DIAS, R., et al.
Thermalelectrical ana-
logy for simulations of superconducting fault current limiters,
Cryogenics,
2014. [13] SHLAGER, K., SCHNEIDER, J.
A Survey of the Finite-Dierence Time
Advances in Computational Electrodynamics: The Finite-Dierence Time-Domain, Artech House.
Domain Literature. In:
Ta ove, A. (Ed.),
Discretização de Equações Diferenciais Parciais - Técnicas de Diferenças Finitas. Sociedade Brasileira de
[14] CUMINATO, J., MENEGUETTE, M.
Matematica, 2013. [15] FARLOW, S.
Partial Dierential Equations for Scientists and Engineers. Dover
Publications, Inc., 1993. [16] GRIFFITHS, D.
Introduction to Electrodynamics.
[17] KRAUS, J., CARVER, K.
Electromagnetics.
[18] INAN, U., MARSHALL, R.
Prentice-Hall, 1999.
McGraw-Hill, 1973.
Numerical Electromagnetics: The FDTD Method.
Cambridge University Press, 2011. [19] MUR, G.
Absorbing Boundary Conditions for the Finite-Dierence Appro-
IEEE Transactions on Electromagnetic Compatibility VOL. EMC 23. NO 4., ximation of the Time-Domain Electromagnetic-Field Equations,
1981. [20] NAJM, F.
Circuit Simulation.
Wiley-IEEE, 2010.
Modeling of Electromagnetic Wave Propagation in Printed Circuit Board and Related Structures. Tese de doutorado, Multimedia University,
[21] KUNG, F.
Malasia, 2003.
Finite-Dierence Time-Domain Studies of the Optical Properties of Metallodielectric Nanostructures. Tese de doutorado, Rice University,
[22] OUBRE, C.
2005. [23] TENENBAUM, M., POLLARD, H.
Ordinary Dierential Equations.
Publications, 1985. [24] ZILL, D.
Equações Diferenciais.
Pearson, 2001.
57
Dover
[25] CHEN, W., SHAO, W., LI, J., et al. Numerical Dispersion Analysis and Key Parameter Selection in Laguerre-FDTD Method,
IEEE Microwave and
Wireless Components Letters VOL. 23. NO 12., 2013. [26] PRESS, W., TEUKOLSKY, S., VETTERLING, W. Cambridge University Press, 1992.
58
Numerical Recipes in C.
Apêndice A Critério de Von Neumann para Estabilidade Numérica O chamado Critério de Von Neumann consiste em escrever a distribuição espacial de uma grandeza física qualquer
U
como Série de Fourier complexa, como mostrado
na Equação (A.1) [18]:
U (x, t) =
∞ X
Um (x, t) =
Sendo Fourier,
λm
Cm (t) e−jkm x
(A.1)
m=−∞
m=−∞
√ j = −1, km =
∞ X
2π o número de onda referente ao termo da Série de λm
o comprimento de onda referente ao termo da Série de Fourier e
Cm (t)
uma função dependente do tempo, também referente ao termo da Série de Fourier.
Um único termo da Série de Fourier é, portanto:
Um (x, t) = Cm (t) ejkm x
(A.2)
Para análise de estabilidade, a Equação (A.2) é discretizada, fazendo o termo
m=1
para simpli car a notação:
U |ti = C|t ejk i∆x
(A.3)
Já que a grandeza foi descrita no espaço como uma Série de Fourier, um deslocamento no espaço equivale a uma defasagem nos termos da Série.
U |ti±1 = U |ti e±jk ∆x
(A.4)
A análise de estabilidade usando o Critério de Von Neumann consiste basicamente nos seguintes passos: 1. Substituir a versão discretizada da componente da Série de Fourier dada pela
59
Equação (A.3) na equação do método numérico a ser analisado. 2. De nir o chamado fator de ampli cação 3. Escrever os termos da forma
ejk∆x
q
como
q=
U |t+1 i U |ti
como senos e cossenos.
4. Analisar o valor assumido pelo fator de ampli cação de discretização espacial (por exemplo,
∆x)
q
para diferentes valores
e discretização temporal (∆t)
O método numérico analisado será estável se a aplicação destes passos resultar em
|q|≤ 1
para todo
k.
Como exemplo, aplicaremos o critério ao método
leapfrog.
A equação geral de
um método deste tipo é descrita abaixo:
U |t+1 i =
U |t−1 i −
v∆t ∆x
(U |ti ejk(i±1)∆x − U |ti−1 )
(A.5)
Aplicando a Equação (A.4) e a de nição do fator de ampli cação
U |t+1 i , U |ti
q =
obtemos a Equação (A.6):
U |t+1 i =
U |ti − q
v∆t ∆x
Escrevendo os termos da forma
(U |ti e+jk ∆x − U |ti e−jk ∆x )
ejk∆x
(A.6)
como senos e cossenos e manipulando o
lado direito da Equação (A.6) obtemos a Equação (A.7)
U |t+1 i =
1 − 2j q
v∆t ∆x
sin(k∆x) U |ti = q U |ti
(A.7)
Ou seja:
q= De nindo
A=
v∆t sin(k ∆x
1 − 2j q
∆x),
v∆t ∆x
sin(k∆x)
(A.8)
podemos reescrever a Equação (A.10) como:
q = −jA ± Lembrando, a exigência é
|q|≤ 1
√ 1 − A2
(A.9)
e portanto:
|q|2 = [Re(q)]2 + [Im(q)]2 = 1 − A2 + A2 = 1 Ou seja, a condição se traduz em de
A
|A|≤ 1,
como o termo
sin(k ∆x)
(A.10)
na de nição
é no máximo 1, a condição de estabilidade para um método do tipo
unidimensional se torna:
60
leapfrog
v ∆t ∆x ≤ 1 → ∆t ≤ ∆x v
(A.11)
Ou seja, utilizando o Critério de Von Neumann chega-se à condição CFL (Equação (2.57)).
61
Apêndice B Correlação Linear O coe ciente de correlação linear
RP earson
(também conhecido como coe ciente de
correlação de Pearson) é dado, para um par de grandezas
(xi , yi )
para
i = 1, 2 · · · N
pela fórmula:
RP earson Sendo
x¯
a média de
x
e
y¯ a
P (xi − x¯)(yi − y¯) pP = pP (xi − x¯)2 (yi − y¯)2 média de
y.
62
[26]
(B.1)
Apêndice C Código Desenvolvido para o Caso Bidimensional utilizando WLP-FDTD 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19 20
% vasculhar a malha de pra cima , esquerda pra direita . % Hz na extrema direita esta fora da malha % Hz na extremidade superior tambem esta fora % Ex na extremidade direita esta fora % Ey na extremidade superior esta fora sizet = 400; lim_t = 11.71*(10^ -9) ; % s t = linspace (0 , lim_t , sizet ) ; p = 151; % numero de polinomios de laguerre a serem utilizados , p = (( ordem do maior polinomio ) -1) deltax = 0.01; % m deltay = 0.001; % m l = 1; % m h = 0.1; % m ptslinha = round ( l / deltax ); ptscoluna = round ( h / deltay ) ; npontos = ptslinha * ptscoluna ; s = 6.07*(10^10) ; % fator de escala temporal mi0 = 4* pi *(10^ -7) ; % H /m eps0 = 8.854*(10^ -12) ; % F/ m c0 = 2.998*(10^8) ; % velocidade da luz
21 22
weight = exp ( - s * t *0.5) ;
63
23
24
Jq = zeros (1 , p ) ; % lista dos coeficientes de laguerre que descrevem a corrente nos pontos de entrada , a ser completada Jp = sparse ( zeros (2*( npontos ) ,p ) ) ; % lista dos coeficientes que descrevem correntes em todos os pontos da geometria
25 26 27 28 29 30
% dados do pulso , todos diretamente do Chung fc = 10^9; Td = 1/(2* fc ); Tc = 3* Td ; Jinput = exp ( -((t - Tc ) / Td ) .^2) .* sin (2* pi .* fc .*( t - Tc ) ) ;
31 32
33 34
laguerre = zeros (p , sizet ) ; % polinomios de laguerre no dominio do tempo laguerre (1 ,:) = 1; % ordem 0 laguerre (2 ,:) = 1 -( s * t ) ; % ordem 1
35 36
37 38 39 40
41 42
% determinando os polinomios de ordem maior que 2 , no dominio do tempo % lembrando , tamanho laguerre = (p , sizet ) if p >2 for q =3: p laguerre (q ,:) =(( laguerre (q -1 ,:) .*((2*( q -1) ) -1 - s *t ) ) -(q -2) .* laguerre (q -2 ,:) ) /( q -1) ; end end
43 44
45 46 47
48
% lista das bases de laguerre ( polinomios multiplicados pelo fator weight ) no dominio do tempo % lembrando tamanho phiq = (p , sizet ) for i =1: p disp ( sprintf ( ' Obtendo bases de laguerre : % d de % d ', i ,p)) phiq (i ,:) = ( weight ) .* laguerre (i ,:) ;
49 50
51
% obtendo termo a ser integrado para completar a lista de coeficientes de laguerre que definem a corrente Jq_pre_integ (i ,:) = Jinput .* phiq (i ,:) ; 64
52 53
end clear weight memoria
% limpando a variavel weight pra liberar
54 55
56 57
58 59 60
% integrando os termos obtidos anteriormente , verdadeiramente obtendo os coeficientes de laguerre que definem o sinal de corrente . for i =1: p disp ( sprintf ( ' Obtendo coeficientes de laguerre para a corrente : % d de % d ' , i , p ) ) % lembrando : Jq =(1 , p ) Jq (1 , i ) = trapz ( s *t , Jq_pre_integ (i ,:) ) ; end
61 62
clear Jq_pre_integ % limpando o termo pre - integracao , ele e inutil e so ocupa memoria
63 64
65
66
67 68 69
70 71
72 73
% ao criar os campos E para um sistema 2 D tenho que decidir qual e a ordem dos campos eletricos , ela sera a seguinte : % Primeiro : todos os campos Ex , percorrendo a geometria linha a linha de baixo pra cima % Depois : todos os campos Ey , percorrendo a geometria linha a linha de baixo pra cima % numeros uteis : % numero de linhas = npontos % todos os Ex = npontos , todos os Ey = npontos , todos os campos eletricos = 2*( npontos ) % todos os Hz = todos os campos magneticos = npontos % como cada campo vai ser a soma de p coeficientes de laguerre : % matriz E tera dimensoes p x 2*( npontos ) % matriz H tera dimensoes p x npontos
74 75 76 77 78 79
E = sparse ( zeros (p ,2*( npontos ) ) ) ; H = sparse ( zeros (p ,( npontos ) ) ) ; % lembrando : definicao de beta na eq 40 do Chung % beta tem o mesmo tamanho da matriz A beta = sparse ( zeros (1 ,2*( npontos ) ) ) ; 65
80 81
82
83
84 85
86 87 88 89 90
91 92 93 94 95 96
% para executar mais de uma vez com as mesmas matrizes e agilizar o processo disp ( ' Checando se existe alguma matriz A e de relacoes aproveitavel ') disp ( 'Se fez alguma mudanca recente em A lembre - se de limpar as variaveis ') flag_A = 1; if (( exist ( 'A ') & exist ( ' relat_betasomaH ') & exist ( ' relat_betasomaE ') ) == 0) | ( max ( size ( A ) ) ~= 2*( npontos ) ) % se matrizes A e de relacoes nao existem ainda ... disp ( ' Matriz A nao detectada ') disp ( ' criando matriz A ') A = sparse ( zeros (2*( npontos ) ) ) ; % A = ( zeros (2*( npontos ) ) ) ; % esses vetores de relacao sao um jeito rapido de implementar , por exemplo , a equacao 40 do Chung ( dentre outras ) disp ( ' criando matrizes de relacao ') relat_betasomaH = sparse ( zeros (2*( npontos ) , npontos ) ) ; relat_betasomaE = sparse ( zeros (2*( npontos ) ) ) ; relat_H_E = sparse ( zeros ( npontos ,2*( npontos ) )) ; flag_A = 0; end
97 98
99
100 101
% a cada iteracao do metodo sera necessario somar todos os coefs das iteracoes anteriores % como so trato Hz , somaH tem metade de somaE ( que inclui tanto Ex quanto Ey ) somaH = sparse ( zeros (1 ,( npontos ) )) ; somaE = sparse ( zeros (1 ,2*( npontos ) ) ) ;
102 103 104 105 106 107
% constantes definidas tambem no Chung CEx = 2/( deltax * eps0 * s ) ; CEy = 2/( deltay * eps0 * s ) ; CHx = 2/( deltax * mi0 * s ) ; CHy = 2/( deltay * mi0 * s ) ;
108
66
109
110 111
112
113
114
115
% existem duas componentes , Ex e Ey , cada uma delas responsavel por metade da matriz A ( npontos ) % preciso de dois loops , um pra cada metade da matriz % [ A ]* transpose ([ Ex (i ,1) Ex (i ,2) ... Ex (1 , j ) Ex (2 , j ) ... Ey (i ,1) Ey (i ,2) ... Ey (1 , j ) Ey (2 , j) ]) = transpose ([ Jx Jy ]) + Betaq -1 % criarei constantes pra lembrar as distancias ( os shifts ) dentro da matriz % sf = same field , ou seja , shift de Ex pra Ex ( termos da eq 39 dependendo de Ex e termos da eq 40 dependendo de Ey ) % Ey39 = distancia entre os termos da eq 39 e os termos Ey % Ex40 = distancia entre os termos da eq 40 e os termos Ex
116 117
118 119
120
% outra observacao e que a determinacao de J e muito mais complicada , como estou tentando o caso apresentado por Chung : % Jx = 0 sempre % e Jy = j( t ) , em uma faixa que corta a geometria verticalmente , o que faz com que os valores aparecam de forma intermitente no vetor Jy . % apareceriam continuamente no vetor Jy se a faixa cortasse a geometria horizontalmente
121 122
% importante notar que max ( size ( A ) ) /2 e o ultimo ponto dos Ex , o primeiro dos Ey e ( max ( size ( A ) ) /2) +1
123 124
125
sfvizinhox = 1; % exemplo , shift do Ex no ponto ( x =1 , y =1) pro Ex no ponto ( x =2 , y =1) sfvizinhoy = ptslinha ; % exemplo , shift do Ex no ponto ( x =1 , y =1) pro Ex no ponto ( x =1 , y =2)
126 127 128 129 130
sfdiagsupesq sfdiagsupdir sfdiaginfesq sfdiaginfdir
= = = =
sfvizinhoy - sfvizinhox ; sfvizinhoy + sfvizinhox ; - sfvizinhoy - sfvizinhox ; - sfvizinhoy + sfvizinhox ;
131
67
132 133
ptsquadrado = npontos ; shift_ExEy = ptsquadrado ;
134 135 136 137
138
139
140
pontos_geometria = 1: ptsquadrado ; pontos_inferiores = 2: ptslinha -1; pontos_superiores = ptsquadrado - ptslinha +2: ptsquadrado -1; % sem pontas pontos_esquerda = 1+ ptslinha : ptslinha : ptsquadrado -(2* ptslinha ) +1; % sem pontas pontos_direita = 2* ptslinha : ptslinha : ptsquadrado ptslinha ; % sem pontas pontos_pontas = [1 ptslinha ptsquadrado - ptslinha +1 ptsquadrado ];
141 142
143 144
pontos_limite = horzcat ( pontos_inferiores , pontos_superiores , pontos_esquerda , pontos_direita , pontos_pontas ) ; pontos_interior = setdiff ( pontos_geometria , pontos_limite ) ;
145 146
147
148
% Jy , todos os pontos com corrente estao na segunda metade do vetor pontos_com_corrente = ptsquadrado + ptslinha +1+(30) : ptslinha : 1+2* ptsquadrado - ptslinha -(70) ;
149 150 151 152 153
disp ( sprintf ( ' Construindo Jp ') ) for k = pontos_com_corrente Jp (k ,:) = - deltax * Jq ; end
154 155
disp ( sprintf ( ' Construindo matriz A ') )
156 157 158
159 160 161
for i = pontos_interior disp ( sprintf ( ' Construindo matriz A : % d de % d ' ,i , ptsquadrado ) ) % eq 39 do Chung Ey_i = i + shift_ExEy ; if flag_A == 1 68
162
163 164 165 166 167 168
disp ( ' Parece que ja existe uma matriz A no workspace , o programar tentara utiliza - la ') break end % Elementos Ex , normal A (i , i) = (1/ CEy ) + CHy + CHy ; A (i , i+ sfvizinhoy ) = - CHy ; A (i ,i - sfvizinhoy ) = - CHy ; % ok
169 170 171 172 173
A (i , Ey_i ) = - CHx ; A (i , Ey_i + sfvizinhox ) = CHx ; A (i , Ey_i + sfdiaginfdir ) = - CHx ; A (i , Ey_i - sfvizinhoy ) = CHx ;
174 175 176
177
178 179 180
% lembrando : relat_betasomaE e relat_betasomaH tem ndelinhas = ndelinhas de beta (2* npontos ) % relat_betasomaH tem menos colunas ( tantas colunas quanto H tem linhas ) relat_betasomaH (i ,i - sfvizinhoy ) = 2; relat_betasomaH (i , i ) = -2; relat_betasomaE (i , i ) = -2/ CEy ;
181 182 183 184 185
% Elementos Ey , normal A ( Ey_i , Ey_i ) = (1/ CEx ) + CHx + CHx ; A ( Ey_i , Ey_i - sfvizinhox ) = - CHx ; A ( Ey_i , Ey_i + sfvizinhox ) = - CHx ;
186 187 188 189 190
A ( Ey_i , i + sfvizinhoy ) = CHy ; A ( Ey_i , i + sfdiagsupesq ) = - CHy ; A ( Ey_i , i ) = - CHy ; A ( Ey_i ,i - sfvizinhox ) = CHy ;
191 192
193 194 195
% equivalente a equacao 40 , Chung ( existe um sinal errado no artigo ) relat_betasomaH ( Ey_i , i ) = 2; relat_betasomaH ( Ey_i ,i - sfvizinhox ) = -2; relat_betasomaE ( Ey_i , Ey_i ) = -2/ CEx ;
196
69
197 198 199 200 201 202 203 204 205 206 207 208 209 210
relat_H_E (i , i ) = - CHy ; relat_H_E (i , i + sfvizinhoy ) = CHy ; relat_H_E (i , Ey_i ) = CHx ; relat_H_E (i , Ey_i + sfvizinhox ) = - CHx ; end for i = pontos_esquerda Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Elementos Ex , normal A (i , i) = (1/ CEy ) + CHy + CHy ; A (i , i+ sfvizinhoy ) = - CHy ; A (i ,i - sfvizinhoy ) = - CHy ;
211 212 213 214 215 216
217
218 219 220
A (i , Ey_i ) = - CHx ; A (i , Ey_i + sfvizinhox ) = CHx ; A (i , Ey_i + sfdiaginfdir ) = - CHx ; A (i , Ey_i - sfvizinhoy ) = CHx ; % lembrando : relat_betasomaE e relat_betasomaH tem ndelinhas = ndelinhas de beta (2* npontos ) % relat_betasomaH tem menos colunas ( tantas colunas quanto H tem linhas ) relat_betasomaH (i ,i - sfvizinhoy ) = 2; relat_betasomaH (i , i ) = -2; relat_betasomaE (i , i ) = -2/ CEy ;
221 222 223 224 225 226 227
% Ey , % ABC A ( Ey_i , Ey_i ) = 0.25*( s / c0 ) +(1/ deltax ) ; A ( Ey_i , Ey_i + sfvizinhox ) = 0.25*( s / c0 ) -(1/ deltax ) ; relat_betasomaE ( Ey_i , Ey_i ) = -0.5*( s/ c0 ) ; relat_betasomaE ( Ey_i , Ey_i + sfvizinhox ) = -0.5*( s / c0 ) ; relat_betasomaH (i ,:) = 0;
228 229 230 231 232 233
% Relacao de H aos campos , normal relat_H_E (i , i ) = - CHy ; relat_H_E (i , i + sfvizinhoy ) = CHy ; relat_H_E (i , Ey_i ) = CHx ; relat_H_E (i , Ey_i + sfvizinhox ) = - CHx ; 70
234 235 236 237 238 239 240 241 242 243 244
end for i =1:1 % ponta inferior esquerda Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex , % PEC A (i ,:) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; relat_betasomaH (i ,:) = 0;
245 246 247 248 249 250 251
% Ey , % ABC A ( Ey_i , Ey_i ) = 0.25*( s / c0 ) +(1/ deltax ) ; A ( Ey_i , Ey_i + sfvizinhox ) = 0.25*( s / c0 ) -(1/ deltax ) ; relat_betasomaE ( Ey_i , Ey_i ) = -0.5*( s/ c0 ) ; relat_betasomaE ( Ey_i , Ey_i + sfvizinhox ) = -0.5*( s / c0 ) ; relat_betasomaH ( Ey_i ,:) = 0;
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
% Relacao de H aos campos , normal relat_H_E (i , i ) = - CHy ; relat_H_E (i , i + sfvizinhoy ) = CHy ; relat_H_E (i , Ey_i ) = CHx ; relat_H_E (i , Ey_i + sfvizinhox ) = - CHx ; end for i = pontos_inferiores Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex , % PEC A (i ,:) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; relat_betasomaH (i ,:) = 0;
269 270 271 272
% Elementos Ey , normal A ( Ey_i , Ey_i ) = (1/ CEx ) + CHx + CHx ; A ( Ey_i , Ey_i - sfvizinhox ) = - CHx ; 71
273
A ( Ey_i , Ey_i + sfvizinhox ) = - CHx ;
274 275 276 277 278
A ( Ey_i , i + sfvizinhoy ) = CHy ; A ( Ey_i , i + sfdiagsupesq ) = - CHy ; A ( Ey_i , i ) = - CHy ; A ( Ey_i ,i - sfvizinhox ) = CHy ;
279 280
281 282 283
% equivalente a equacao 40 , Chung ( existe um sinal errado no artigo ) relat_betasomaH ( Ey_i , i ) = 2; relat_betasomaH ( Ey_i ,i - sfvizinhox ) = -2; relat_betasomaE ( Ey_i , Ey_i ) = -2/ CEx ;
284 285 286 287 288 289
% Relacao entre H e os campos E relat_H_E (i , i ) = - CHy ; relat_H_E (i , i + sfvizinhoy ) = CHy ; relat_H_E (i , Ey_i ) = CHx ; relat_H_E (i , Ey_i + sfvizinhox ) = - CHx ;
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
%J Jp (i ,:) = 0; Jp ( Ey_i ,:) = 0; end for i = ptslinha : ptslinha % ponta inferior direita Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex , % Nao existe Ex aqui A (: , i) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i , i ) = 0; relat_betasomaE (i , i + sfvizinhoy ) = 0; relat_betasomaH (i ,:) = 0;
306 307 308 309 310
% Ey , % ABC A ( Ey_i , Ey_i ) = 0.25*( s / c0 ) +(1/ deltax ) ; A ( Ey_i , Ey_i - sfvizinhox ) = 0.25*( s / c0 ) -(1/ deltax ) ; relat_betasomaE ( Ey_i , Ey_i ) = -0.5*( s/ c0 ) ; 72
311 312
relat_betasomaE ( Ey_i , Ey_i - sfvizinhox ) = -0.5*( s / c0 ) ; relat_betasomaH ( Ey_i ,:) = 0;
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
% Relacao de H aos campos , nao existe Hz relat_H_E (i ,:) = 0; end for i = pontos_direita Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex , nao existe A (: , i) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; % relat_betasomaE (i , i) = 0; % relat_betasomaE (i , i+ sfvizinhoy ) = 0;
328 329 330 331 332 333 334
% Ey , ABC A ( Ey_i , Ey_i ) = 0.25*( s / c0 ) +(1/ deltax ) ; A ( Ey_i , Ey_i - sfvizinhox ) = 0.25*( s / c0 ) -(1/ deltax ) ; relat_betasomaE ( Ey_i , Ey_i ) = -0.5*( s/ c0 ) ; relat_betasomaE ( Ey_i , Ey_i - sfvizinhox ) = -0.5*( s / c0 ) ; relat_betasomaH ( Ey_i ,:) = 0;
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
% Relacao de H aos campos relat_H_E (i ,:) = 0; relat_H_E (i , i ) = 0; % nao existe Ex nem Hz relat_H_E (i , i + sfvizinhoy ) = 0; relat_H_E (i , Ey_i + sfvizinhox ) = 0; relat_H_E (i , Ey_i ) = 0; end for i = ptsquadrado : ptsquadrado % ponta superior direita Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex Nao existe A (: , i) = 0; 73
350 351 352
A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; relat_betasomaH (i ,:) = 0;
353 354 355 356 357 358
% Ey Nao existe A (: , Ey_i ) = 0; A ( Ey_i , Ey_i ) = (1/ CEx ) + CHx + CHx ; relat_betasomaE ( Ey_i ,:) = 0; relat_betasomaH ( Ey_i ,:) = 0;
359 360 361 362 363 364 365 366 367 368 369 370 371 372
% Relacao de H aos campos , H nao existe relat_H_E (i ,:) = 0; end for i = pontos_superiores Ey_i = i + shift_ExEy ; if flag_A == 1 break end % Ex , % PEC A (i ,:) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; relat_betasomaH (i ,:) = 0;
373 374 375 376 377 378
% Ey Nao existe A (: , Ey_i ) = 0; A ( Ey_i , Ey_i ) = (1/ CEx ) + CHx + CHx ; relat_betasomaE ( Ey_i ,:) = 0; relat_betasomaH ( Ey_i ,:) = 0;
379 380 381
% Relacao de H aos campos , H nao existe relat_H_E (i ,:) = 0;
382 383 384 385 386 387
%J Jp (i ,:) = 0; Jp ( Ey_i ,:) = 0; end for i = ptsquadrado - ptslinha +1: ptsquadrado - ptslinha +1 % ponta superior esquerda 74
388 389 390 391 392 393 394 395 396
Ey_i = i + shift_ExEy ; if flag_A == 1 break end % % Ex , % PEC A (i ,:) = 0; A (i , i) = (1/ CEy ) + CHy + CHy ; relat_betasomaE (i ,:) = 0; relat_betasomaH (i ,:) = 0;
397 398 399 400 401 402
% Ey Nao existe A (: , Ey_i ) = 0; A ( Ey_i , Ey_i ) = (1/ CEx ) + CHx + CHx ; relat_betasomaE ( Ey_i ,:) = 0; relat_betasomaH ( Ey_i ,:) = 0;
403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422
% Relacao de H aos campos , H nao existe relat_H_E (i ,:) = 0; end for j = 1: p disp ( sprintf ( ' Loop principal : % d de % d ' , j , p ) ) j; % j = 1 , ordem = 0 , j = 2 , ordem = 1... if j ==1 % lembrando H =( p ,2*( npontos ) ) % logo , sum ( H ) = (1 ,( npontos ) ) % lembrando beta = (1 ,2*( npontos ) ) rightside = Jp (: , j ) end if j > 1 % lembrando H =( p , npontos ) % logo , sum ( H ) = (1 , npontos ) % lembrando beta = (1 , npontos ) somaH = sum (H (1: j ,:) ) ; somaE = sum (E (1: j ,:) ) ;
423 424
425
beta (1 ,:) = transpose ( relat_betasomaH * transpose ( somaH ) + relat_betasomaE * transpose ( somaE ) ) ; rightside = Jp (: , j ) + transpose ( beta ) ; 75
426 427 428
429
end E (j ,:) = transpose ( mldivide (A , rightside ) ) ; H (j ,:) = transpose ( relat_H_E * transpose ( E (j ,:) ) - 2*( transpose ( somaH ) ) ) ; end
430 431 432 433 434 435 436 437 438
439
440 441
% Descobertos os coeficientes , voltando pro tempo Exrecuperado = zeros ( sizet ,( npontos ) ) ; Eyrecuperado = zeros ( sizet ,( npontos ) ) ; Hzrecuperado = zeros ( sizet ,( npontos ) ) ; disp ( sprintf ( ' Voltando pro dominio do tempo ... Ex ') ) for pos = 1:( npontos ) for laguerre = 1: p Exrecuperado (: , pos ) = Exrecuperado (: , pos ) + E ( laguerre , pos ) * transpose ( phiq ( laguerre ,:) ) ; Hzrecuperado (: , pos ) = Hzrecuperado (: , pos ) + H ( laguerre , pos ) * transpose ( phiq ( laguerre ,:) ) ; end end
442 443 444 445 446
447 448
disp ( sprintf ( ' Voltando pro dominio do tempo ... Ey ') ) for pos = (( npontos ) +1:2*( npontos ) ) for laguerre = 1: p Eyrecuperado (: , pos -( npontos ) ) = Eyrecuperado (: , pos -( npontos ) ) + E ( laguerre , pos ) * transpose ( phiq ( laguerre ,:) ) ; end end
76
Apêndice D Código Desenvolvido para o Caso Unidimensional utilizando WLP-FDTD 1 2 3 4 5 6 7 8 9
sizet = 8100; lim_t = 30 e -8; t = linspace (0 , lim_t , sizet ) ; p = 91; deltat = t (2) -t (1) deltax = 0.01; l = 10; npontos = ( l / deltax ) ; s = 1.1333 e +09;
10 11 12 13
mi0 = 4* pi *(10^ -7) ; % H /m eps0 = 8.854187817*(10^ -12) ; % F / c0 = 2.998*(10^8) ; % velocidade da luz
14 15 16
17
weight = exp ( - s * t *0.5) ; Jq = zeros (1 , p) ; % lista dos coeficientes de laguerre que descrevem a corrente no ponto de entrada Jp = zeros (1 , npontos ) ;
18 19 20 21
Tk = 1e -8; Jinput = exp ( -((t -5* Tk ) / Tk ) .^2) ; pt_entrada = npontos /2;
22
77
23 24 25 26 27 28
29 30 31 32 33
34
Eanalitico = zeros ( sizet , npontos ) ; Hanalitico = zeros ( sizet , npontos ) ; for i = 1: npontos /2 x = i* deltax ; tdeslocado = t +(( x -( deltax * npontos /2) ) / c0 ) ; Eanalitico (: , i ) = - exp ( -(( tdeslocado -5* Tk ) / Tk ) .^2) /(2* c0 * eps0 ) ; end for i = 1+ npontos /2: npontos x = i* deltax ; tdeslocado = t -(( x -( deltax * npontos /2) ) / c0 ) ; Eanalitico (: , i ) = - exp ( -(( tdeslocado -5* Tk ) / Tk ) .^2) /(2* c0 * eps0 ) ; end
35 36
37 38
laguerre = zeros (p , sizet ); % polinomios de laguerre no dominio do tempo laguerre (1 ,:) = 1; % ordem 0 laguerre (2 ,:) = 1 -( s * t ) ; % ordem 1
39 40 41
42 43 44
45 46
% lembrando , tamanho laguerre = (p , sizet ) % determinando os polinomios de ordem maior que 2 , no dominio do tempo if p >2 for q =3: p laguerre (q ,:) =(( laguerre (q -1 ,:) .*((2*( q -1) ) -1 - s* t ) ) -(q -2) .* laguerre (q -2 ,:) ) /( q -1) ; end end
47 48 49 50 51
52 53 54
% lista das bases de laguerre no dominio do tempo % lembrando tamanho phiq = (p , sizet ) for i =1: p disp ( sprintf ( ' Obtendo bases de laguerre : % d de % d ', i , p ) ) phiq (i ,:) = ( weight ) .* laguerre (i ,:) ; Jq_pre_integ (i ,:) = Jinput .* phiq (i ,:) ; end 78
55
clear weight
56 57 58
59 60 61
for i =1: p disp ( sprintf ( ' Obtendo coeficientes de laguerre para a corrente : % d de % d ' , i , p ) ) % lembrando Jq =(1 , p ) Jq (1 , i ) = trapz ( s *t , Jq_pre_integ (i ,:) ) ; end
62 63
clear Jq_pre_integ ;
64 65 66 67 68 69 70
A = sparse ( zeros ( npontos ) ); E = zeros (p , npontos ) ; H = zeros (p , npontos ) ; beta = zeros (1 , npontos ) ; somaH = zeros (1 , npontos ) ; somaE = zeros (1 , npontos ) ;
71 72
73 74 75
% os vetores de relacao sao um jeito rapido de implementar , por exemplo , equacoes da forma 40 do Chung relat_betasomaH = zeros ( npontos ) ; relat_betasomaE = zeros ( npontos ) ; relat_H_E = zeros ( npontos ) ;
76 77 78
CE = 2/( deltax * eps0 * s) ; CH = 2/( deltax * mi0 * s );
79 80 81
82 83 84 85
86 87 88
for i =2: npontos -1 % equivalente a equacao 39 , Chung . modificada pra TEMx A (i , i ) = ((1/ CE ) +2* CH ) ; A (i ,i -1) = - CH ; A (i , i +1) = - CH ; % equivalente a equacao 40 , Chung , modifcada pra TEMx relat_betasomaH (i , i ) = 2; relat_betasomaH (i ,i -1) = -2; relat_betasomaE (i , i ) = -2/ CE ; 79
89 90
91 92 93 94 95
% equivalente a equacao 34 , Chung . Modificada pra TEMx relat_H_E (i , i ) = CH ; relat_H_E (i , i +1) = - CH ; end relat_betasomaH = sparse ( relat_betasomaH ) ; relat_betasomaE = sparse ( relat_betasomaE ) ;
96 97 98 99
% adicionando fronteiras ABC nas duas extremidades A (1 ,1) = 0.25*( s / c0 ) +(1/ deltax ) ; A (1 ,2) = 0.25*( s/ c0 ) -(1/ deltax ) ;
100 101 102
A ( end , end -1) = 0.25*( s / c0 ) -(1/ deltax ); A ( end , end ) = 0.25*( s / c0 ) +(1/ deltax ) ;
103 104 105
relat_H_E (1 ,1) = CH ; relat_H_E (1 ,2) = - CH ;
106 107 108
relat_betasomaH (1 ,1) = 2; relat_betasomaH ( end , end ) = 2;
109 110 111 112 113
relat_betasomaE (1 ,1) = -0.5*( s / c0 ) ; relat_betasomaE (1 ,2) = -0.5*( s / c0 ) ; relat_betasomaE ( end , end ) = -0.5*( s / c0 ) ; relat_betasomaE ( end , end -1) = -0.5*( s/ c0 ) ;
114 115 116 117 118 119 120 121 122 123 124 125 126
% lembrando Jq = (1 , p) % %%%%%%% Loop Principal for j = 1: p j; disp ( sprintf ( ' Loop principal : % d de % d ', j , p ) ) % j = 1 , ordem = 0 , j = 2 , ordem = 1... if j ==1 % lembrando H =( p , npontos ) % logo , sum ( H ) = (1 , npontos ) % lembrando beta = (1 , npontos ) somaH = H (1 ,:) ; somaE = E (1 ,:) ; 80
Jp (1 , npontos ) = 0; % PEC ou ABC , as duas tem essa caracteristica relativa a Jp Jp (1 ,1) = 0; % PEC ou ABC Jp (1 , pt_entrada ) = - Jq (1 , j ) ; rightside = transpose ( Jp ) ;
127
128 129 130
end if j > 1
131 132
% lembrando H =( p , npontos ) % logo , sum ( H ) = (1 , npontos ) % lembrando beta = (1 , npontos ) somaH = sum (H (1: j ,:) ) ; somaE = sum (E (1: j ,:) ) ; beta (1 ,:) = relat_betasomaH * transpose ( somaH ) + relat_betasomaE * transpose ( somaE ) ; Jp (1 , npontos ) = 0; % PEC Jp (1 ,1) = 0; % PEC Jp (1 , pt_entrada ) = - Jq (1 , j ) ; rightside = transpose ( Jp + beta ) ;
133 134 135 136 137 138
139 140 141 142
end E (j ,:) = mldivide (A , rightside ) ; H (j ,:) = relat_H_E * transpose ( E (j ,:) ) - 2*( transpose ( somaH ) ) ;
143 144 145
146
end
147 148 149 150 151 152
153
154 155
Erecuperado = zeros ( sizet , npontos ) ; Hrecuperado = zeros ( sizet , npontos ) ; for pos = 1: npontos for laguerre = 1: p Erecuperado (: , pos ) = Erecuperado (: , pos ) + E ( laguerre , pos ) * transpose ( phiq ( laguerre ,:) ) ; Hrecuperado (: , pos ) = Hrecuperado (: , pos ) + H ( laguerre , pos ) * transpose ( phiq ( laguerre ,:) ) ; end end
81
