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 DE GRADUAÇÃO SUBMETIDO AO CORPO DOCENTE DO

CURSO DE ENGENHARIA ELÉTRICA DA ESCOLA POLITÉCNICA

DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE

DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO 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. 2. FDTD. 3. Eletromagnetismo

Computacional. I. Tiago Batista de Sousa, Wescley

et al. 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 mo-

mentaneamente o seu computador, muito obrigado! Também devo agradecer

ao Felipe Costa por ter me recomendado muito material sobre diferenças �ni-

tas 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 escre-

ver 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 pri-

meiro período. Desejo sucesso e felicidades a todos vocês!

� Ao meu amigo, Fred Guimarães, pelas dicas que me deu relativas a este tra-

balho.

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: Finite

Di�erence Time Domain), 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 ver-

sõ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 bidi-

mensional. 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-Di�erence Time-Domain method based on the Yee Algorithm is a

traditional and relatively simple numerical technique to solve the Maxwell Equa-

tions. 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 Courant-

Friedrichs-Lewy Condition, that introduces a coupling requirement between the dis-

cretization 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 eletro-

dynamics 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) . . . . . . . . . . . 11

2.4 Algoritmo de Yee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.4.1 Condições de Fronteira . . . . . . . . . . . . . . . . . . . . . . 15

2.5 Estabilidade e a Condição Courant-Friedrichs-Lewy . . . . . . . . . . 17

3 O Método FDTD baseado em Polinômios de Laguerre 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

3.4.1 Modo TE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

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 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.4 Simulação Unidimensional: Propagação de Onda Transversal Eletro-

magnética . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.4.1 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 WLP-

FDTD 63

D Código Desenvolvido para o Caso Unidimensional utilizando WLP-

FDTD 77

ix

Lista de Figuras

2.1 Aumento no número de artigos publicados em IEEE Transactions on

Microwaves utilizando FDTD a partir dos anos 90. . . . . . . . . . . 4

2.2 Diferenças centradas (DC), regressivas (DR) e progressivas (DP). . . 6

2.3 Malha de Yee para o modo TEz. . . . . . . . . . . . . . . . . . . . . . 12

2.4 Malha de Yee para o modo TEM em relação a x polarizado em y. . . 13

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). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

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

2.7 Simulação de Campo Elétrico Ey utilizando método FDTD com Sc =

1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.8 Simulação de Campo Elétrico Ey utilizando método FDTD com Sc =

1.005. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.8 Simulação de Campo Elétrico Ey utilizando método FDTD com Sc =

1.005. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.1 Polinômios de Laguerre de ordem 0 a 4. . . . . . . . . . . . . . . . . . 23

3.2 Bases de Laguerre até ordem 4, escala temporal s = 1. . . . . . . . . 26

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. . . . . . . . . . 27

4.1 Fluxograma do Algoritmo Implementado para o Método WLP-FDTD. 40

4.2 Simulação bidimensional para validação do método. . . . . . . . . . . 41

4.3 Forma da densidade de corrente elétrica inserida no problema bidi-

mensional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.4 Resultado da simulação bidimensional mostrada na Figura 4.2. . . . . 43

x

4.5 Resultados para a simulação bidimensional utilizando o métodoWLP-

FDTD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.5 Resultados para a simulação bidimensional utilizando o métodoWLP-

FDTD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.6 Estrutura do problema unidimensional. . . . . . . . . . . . . . . . . . 48

4.7 Forma de Onda da Densidade Super�cial de Corrente Ky. . . . . . . . 49

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

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

4.8 Sobreposição dos resultados da simulação unidimensional utilizando

o método WLP-FDTD e dos apresentados por Feynman et al.[1]. . . . 52

xi

Lista de Tabelas

2.1 Algumas das Grandezas Físicas Eletromagnéticas. . . . . . . . . . . . 7

4.1 Comparação entre os valores obtidos por Chung et al. e os obtidos

por algoritmo próprio, para os pontos a, b, c e d. . . . . . . . . . . . . 44

4.2 Comparação entre os valores analíticos e os obtidos pelo método

WLP-FDTD para campo elétrico no sistema físico unidimensional. . . 52

xii

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

CHx,y, C

Ex,y Constantes utilizadas no método WLP-FDTD, p. 32

Ln Polinômio de Laguerre de ordem n, p. 22

NL Número necessário de bases de Laguerre para sintetizar uma

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

RPearson Coe�ciente de Correlação de Pearson, p. 38

Tf Duração, em segundos, de uma forma de onda., p. 28

U Função qualquer descrita como uma série., p. 29

∆t Discretização Temporal, p. 13

∆x Discretização Espacial, coordenada x, p. 13

∆y Discretização Espacial, coordenada y, p. 13

Φ0,1,2... Conjunto ortogonal in�nito de funções, p. 24

βq−1 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: F/m, p. 7

µ Permeabilidade Magnética, Unidade SI: H/m, p. 7

xiii

φ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: Wb/m2, p. 7

~D Densidade de Fluxo Elétrico, Unidade SI: C/m2, p. 7

~E Intensidade de Campo Elétrico, Unidade SI: V/m, p. 7

~H Intensidade de Campo Magnético, Unidade SI: A/m, p. 7

~J Densidade de Corrente, Unidade SI: A/m2, p. 7

~K Densidade Super�cial de Corrente, p. 46

c0,1,2... 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 Vetores quaisquer., p. 23

v Velocidade qualquer de propagação de uma onda eletromagné-

tica, dada em m/s, p. 16

w(t) Função peso., p. 25

x0 Ponto central de um intervalo, ponto este no qual se deseja

avaliar o operador diferencial., p. 4

xNLmax 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. 15

CAE Computer Assisted Engineering, p. 1

CFL Courant-Friedrichs-Lewy, p. 1

FDTD Finite Di�erence Time Domain, p. vi

IEEE Institute of Electrical and Electronic Engineers, p. 3

LASUP Laboratório de Aplicações de Supercondutores, p. 2

PEC Perfect Electric Conductor, p. 15

PML Perfectly Matched Layer, p. 16

TE Transversal Elétrica(o), p. 9

TM Transversal Magnética(o), p. 9

WLP-FDTD Weighted Laguerre Polynomials Finite Di�erence Time Do-

main, 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, utili-

zando 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

FDTD, do inglês Finite-Di�erence Time Domain) é largamente empregado e possui

diversas implementações na forma de softwares, comerciais e gratuitos. Este método

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 multi-

escala, 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 La-

guerre (também conhecido como Laguerre-FDTD ou WLP-FDTD, do inglês Weigh-

ted Laguerre Polynomials - Finite-Di�erence 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 par-

cela 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 impor-

tâ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 normal-

mente encontradas analiticamente, utilizando séries in�nitas ou formas fechadas.

Tais abordagens possuem algumas limitações no tratamento de materiais com geo-

metrias 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 dire-

tamente 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 tratando do assunto entre os anos

80 e 90 [13].

3

1965 1970 1975 1980 1985 1990 19950

100

200

300

400

500

Publicações Utilizando FDTD

Nú

me

ro d

e P

ub

lica

çõ

es

Ano

Figura 2.1: Aumento no número de artigos publicados em IEEE Transactions onMicrowaves utilizando FDTD a partir dos anos 90.

2.2 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 leapfrog [14]. Tomamos uma função f e um intervalo de

raio h ao redor do ponto x0, 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]:

f(x0 + h) = f(x0) + h∂f(x)

∂x

∣∣∣∣x=x0

+h2

2!

∂2f(x)

∂x2

∣∣∣∣x=x0

+h3

3!

∂3f(x)

∂x3

∣∣∣∣x=x0

+ · · · (2.1)

Rearranjando termos e dividindo por h os dois lados da Equação (2.1) obtemos:

f(x0 + h)− f(x0)

h=

∂f(x)

∂x

∣∣∣∣x=x0

+h

2!

∂2f(x)

∂x2

∣∣∣∣x=x0

+h2

3!

∂3f(x)

∂x3

∣∣∣∣x=x0

+ · · · (2.2)

Truncando o lado direito desta equação no segundo termo conseguimos uma aproxi-

mação pra derivada conhecida como aproximação por diferenças progressivas (Equa-

ção (2.3)):

4

∂f(x)

∂x

∣∣∣∣x=x0

=f(x0 + h)− f(x0)

h+O(h) (2.3)

O termo O(h) está representado na notação Big Oh e simboliza o erro devido ao

truncamento. Como neste caso o primeiro termo do truncamento é proporcional

a h, 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):

f(x0 − h) = f(x0)− h∂f(x)

∂x

∣∣∣∣x=x0

+h2

2!

∂2f(x)

∂x2

∣∣∣∣x=x0

− h3

3!

∂3f(x)

∂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)

∂x

∣∣∣∣x=x0

=f(x)− f(x0 − h)

h+O(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(x0 + h)− f(x0 − h) = h∂f(x)

∂x

∣∣∣∣x=x0

+2 h3

3!

∂3f(x)

∂x3

∣∣∣∣x=x0

+ · · · (2.6)

Dividindo por h nos dois lados da equação obtemos:

f(x0 + h)− f(x0 − h)

h=∂f(x)

∂x

∣∣∣∣x=x0

+2 h2

3!

∂3f(x)

∂x3

∣∣∣∣x=x0

+ · · · (2.7)

Importante notar que, neste caso, o termo a ser truncado depende de 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)

h=∂f(x)

∂x

∣∣∣∣x=x0

+O(h2) (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((x0 + h), y)− f((x0 − h), y)

h=∂f(x, y)

∂x

∣∣∣∣x=x0

+O(h2) (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 centra-

das.

DC

DR

DP

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] :

∇× ~H = ~J +∂ ~D

∂t(2.10)

∇× ~E = −∂~B

∂t(2.11)

∇ · ~D = 0 (2.12)

∇ · ~B = 0 (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).

~D = ε ~E (2.14)

~B = µ ~H (2.15)

Na presença de materiais condutores é preciso adicionar também a Equação (2.16).

~J = ~Jindependente + σ ~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 Campo Elétrico Vm

~D Densidade de Fluxo Elétrico Cm2

~H Campo Magnético Am

~B Densidade de Fluxo Magnético Wbm2

~J Densidade de Corrente Elétrica Am2

ε Permissividade Elétrica Fm

µ Permeabilidade Magnética Hm

σ Condutividade Elétrica S

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

∂t=

1

ε(∇× ~H − ~J) (2.17)

∂ ~H

∂t= − 1

µ(∇× ~E) (2.18)

∇ · ε ~E = 0 (2.19)

∇ · µ~B = 0 (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

∂t=

1

µ

(∂Ey∂z− ∂Ez

∂y

)(2.21)

∂Hy

∂t=

1

µ

(∂Ez∂x− ∂Ex

∂z

)(2.22)

∂Hz

∂t=

1

µ

(∂Ex∂y− ∂Ey

∂x

)(2.23)

∂Ex∂t

=1

ε

(∂Hz

∂y− ∂Hy

∂z− Jx

)(2.24)

∂Ey∂t

=1

ε

(∂Hx

∂z− ∂Hz

∂x− Jy

)(2.25)

∂Ez∂t

=1

ε

(∂Hy

∂x− ∂Hx

∂y− Jz

)(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 cor-

rente 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 pre-

sentes 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 y e z. Este procedimento

8

resulta nas Equações (2.31), (2.32), (2.33) e (2.34):

∂Ez∂t

=1

ε

(∂Hy

∂x− ~Jz

)(2.27)

∂Hy

∂t=

1

µ

(∂Ez∂t

)(2.28)

∂Ey∂t

=1

ε

(− ∂Hz

∂x− Jy

)(2.29)

∂Hz

∂t=

1

µ

(− ∂Ey

∂x

)(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.

∂Ez∂t

=1

ε

(∂Hy

∂x− ~Jz

)(2.31)

∂Hy

∂t=

1

µ

(∂Ez∂t

)(2.32)

O conjunto formado pelas Equações (2.33) e (2.34) é chamado de modo trans-

versal 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 TEMx polar. y.

∂Ey∂t

=1

ε

(− ∂Hz

∂x− Jy

)(2.33)

∂Hz

∂t=

1

µ

(− ∂Ey

∂x

)(2.34)

Os dois modos apresentados, TEMx polar. z e TEMx polar. y, 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 cor-

rente 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

∂t=

1

µ

(− ∂Ez

∂y

)(2.35)

∂Hy

∂t=

1

µ

(∂Ez∂x

)(2.36)

∂Hz

∂t=

1

µ

(∂Ex∂y− ∂Ey

∂x

)(2.37)

∂Ex∂t

=1

ε

(∂Hz

∂y− Jx

)(2.38)

∂Ey∂t

=1

ε

(− ∂Hz

∂x− Jy

)(2.39)

∂Ez∂t

=1

ε

(∂Hy

∂x− ∂Hx

∂y− Jz

)(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).

Como é na direção z que não há qualquer variação, este modo é ditomodo transversal

elétrico em relação a z, ou modo TEz [7].

∂Ex∂t

=1

ε

(∂Hz

∂y− Jx

)(2.41)

∂Ey∂t

=1

ε

(− ∂Hz

∂x− Jy

)(2.42)

10

∂Hz

∂t=

1

µ

(∂Ex∂y− ∂Ey

∂x

)(2.43)

2.3.3 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 transversal magnético em relação a z, ou

Modo TMz:

∂Hx

∂t=

1

µ

(− ∂Ez

∂y

)(2.44)

∂Hy

∂t=

1

µ

(∂Ez∂x

)(2.45)

∂Ez∂t

=1

ε

(∂Hy

∂x− ∂Hx

∂y− Jz

)(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

bidimensionais (modo TEz) e a Figura 2.4 para estudos unidimensionais, modo TEM

em relação a x, polarizado em y.

Figura 2.3: Malha de Yee para o modo TEz.

12

Figura 2.4: Malha de Yee para o modo TEM em relação a x polarizado em y.

Os índices i e j 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 (x, y) possa ter sua posição descrita com

uma expressão da forma (i∆x, j∆y) [7].

É importante observar que campos elétricos e magnéticos estão deslocados tam-

bé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 Ez e a outra (Equa-

ção (2.48)) atualizará os valores para os campos magnéticos Hy. 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.

Ey|t+1i = Ey|ti −

∆t

ε∆x

(Hz|ti−Hz|ti−1

)(2.47)

Hz|t+1i = Hz|ti −

∆t

µ∆x

(Ey|ti+1−Ey|ti

)(2.48)

13

Figura 2.5: Visualização da Malha de Yee no espaço-tempo, o círculo vermelhomostra o ponto onde o operador diferencial é discretizado para que se obtenha aEquaçã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 vermelhomostra o ponto onde o operador diferencial é discretizado para que se obtenha aEquação de Atualização para Campos Magnéticos (Equação 2.44).

14

Adicionando Jy(t) como função forçante obtemos:

Ey|t+1i = Ey|ti −

∆t

ε∆x

(Hz|ti −Hz|ti−1

)+

∆t

εJy|t+1

i (2.49)

Hz|t+1i = Hz|ti −

∆t

µ∆x

(Ey|ti+1 − Ey|ti

)(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 leapfrog.

2.4.1 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 fron-

teiras da simulação.

Condições do tipo PEC

Em alguns casos deseja-se determinar um espaço de simulação cercado por conduto-

res 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: Perfect Elec-

tric Conductor). 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�nitiva-

mente, 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 ar-

ti�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 ab-

sorvem 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]. Este tipo de condição de fronteira é conhecida como

condição ABC (do inglês: Absorbing Boundary Conditions).

15

(∂

∂x± 1

v

∂

∂t

)Ey(x, y, t) = 0 (2.51)

sendo v a velocidade de propagação da onda eletromagnética no meio estudado.

A Equação (2.51) pode ter duas soluções, a saber:(∂

∂x+

1

v

∂

∂t

)Ey(x, t) = 0 → Ey(x, t) = f(x − vt) (2.52)

(∂

∂x− 1

v

∂

∂t

)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 di-

reçã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-

plo, para uma fronteira vertical presente em x = 0, equivalente a i = 1. O ponto de

discretização espacial utilizado é normalmente i = 12, de modo a caracterizar o início

da fronteira no ponto médio entre dois campos elétricos, o ponto de discretização

temporal é t = 12, ou seja, também no ponto médio entre dois campos elétricos,

mas no tempo. Estas escolhas para realizar a aproximação por diferenças centrais

resultam na Equação (2.54) [18]:

Ey|t+1i=1+1/2 − Ey|ti=1+1/2

v∆t=Ey|t+1/2

i=2 − Ey|t+1/2i=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:

(Ey|t+1

i=1+Ey|t+1i=2

)− (Ey|ti=1+Ey|ti=2)

2v∆t=

(Ey|t+1

i=2+Ey|ti=2

)−(Ey|t+1

i=1+Ey|ti=1

)2∆x

(2.55)

Finalmente, rearranjando para o campo elétrico na fronteira, conseguimos chegar à

Equação de Atualização (2.56):

Ey|t+1i=1= Ey|ti=2 +

(v∆t−∆x

v∆t+ ∆x

)(Ey|t+1

i=2−Ey|ti=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 per-

pendicular às fronteiras da simulação. Caso este não seja o caso, outros métodos

posteriores tais como o Perfectly-Matched Layer (PML) são recomendados [8].

16

2.5 Estabilidade e a Condição Courant-Friedrichs-

Lewy

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 Courant-

Friedrichs-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].

∆t ≤ 1

v√

1∆x2

+ 1∆y2

(2.57)

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 discreti-

zaçã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 on-

das 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 tem-

poral, 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(2.58)

A Figura 2.7 mostra resultados para campo elétrico Ey de uma simulação uni-

dimensional da propagação de uma onda transversal eletromagnética utilizando o

17

método FDTD com Sc = 1.0 (valor considerado ideal para este tipo de simulação

[9]). Note que, utilizando Número de Courant unitário, a onda eletromagnética

percorre a geometria sem sofrer distorções.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Distribuição de Ey ao longo da geometria em diferentes instantes

Eixo x (cm)

Ey (

V/m

)

t = 0.93 ns

t = 1.06 ns

t = 1.19 ns

t = 1.33 ns

Figura 2.7: Simulação de Campo Elétrico 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: Sc = 1.005. A Figura

2.8 mostra a intensidade do campo elétrico Ey 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

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(a) t = 0.134 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(b) t = 0.335 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(c) t = 0.469 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(d) t = 0.536 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(e) t = 0.603 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(f) t = 0.670 ns

Figura 2.8: Simulação de Campo Elétrico Ey utilizando método FDTD com Sc =1.005.

19

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−10

−8

−6

−4

−2

0

2

4

6

8

10

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(g) t = 0.804 ns

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4−50

−40

−30

−20

−10

0

10

20

30

40

50

Distribuição de Ey ao longo da geometria em um certo instante.

Eixo x (cm)

Ey (

V/m

)

(h) t = 0.938 ns

Figura 2.8: Simulação de Campo Elétrico Ey utilizando método FDTD com Sc =1.005.

20

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-

FDTD (do inglês:Alternate-Direct Implicit FDTD). 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 Laguerre-

FDTD 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∂2y

∂x2+ (1− x)

∂y

∂x+ ky = 0 (3.1)

sendo x ≥ 0 e k um número natural, a solução em termos de séries in�nitas desta

equação diferencial é descrita pela Equação (3.2).

yk(x) = a0

(1− kx+

k(k − 1)

22x2 − k(k − 1)(k − 2)

22 · · · 32x3

+k(k − 1)(k − 2)(k − 3)

22 · · · 32 · · · 42x4 + . . .

(−1)n k!

(n! )2(k − n)!xn) (3.2)

Aplicar k = 0, 1, 2, 3... e a0 = k! 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

. . .

Lk(x) = (k! )2

k∑n=0

(−1)n

(n! )2(k − n)!xn

(3.3)

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) =(2k − 1− x) Lk−1(x)− (k − 1) Lk−2(x)

kpara k ≥ 2

(3.4)

22

A Figura 3.1 mostra os cinco primeiros polinômios de Laguerre.

0 5 10 15 20−10

−5

0

5

10

15

20

x

L(x

)

Ordem 0

Ordem 1

Ordem 2

Ordem 3

Ordem 4

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, t.

3.3 Ortogonalidade

Os conceitos vetoriais de produto interno escalar e ortogonalidade podem ser gene-

ralizados para funções. Consideremos dois vetores u e v, de�nidos, por exemplo, no

espaço vetorial IR2. O produto interno entre esses dois vetores, representado u · v ,é de�nido pelas seguintes propriedades [24]:

i. u · v = v · u

ii. ku · v = k(u · v), sendo k escalar.

iii. u · u = 0 se u = 0 e u · u > 0 se u 6= 0

iv. u+ v · w = u · w + v · w, sendo w 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

f1 · f2 =

∫ b

a

f1(t)f2(t)dt (3.5)

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 u e v formam base no espaço bidi-

mensional IR2, qualquer vetor 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 Φn de�nido no intervalo t = [a, b], ou seja, um conjunto da forma

[Φ0,Φ1,Φ2...] que respeita a condição (3.7).∫ b

a

Φm(t) Φn(t) dt = δmn (3.7)

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 (3.9)

onde cn 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

∫ b

a

w(t) Φm(t) Φn(t) dt = δmn (3.10)

O conjunto de polinômios de Laguerre se encaixa no caso (3.10). Estas tem uma

relação de ortogonalidade para uma função peso igual a e−t. Observando também

que o intervalo anteriormente representado por t = [a, b] está de�nido para este

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. ∫ ∞0

e−t Lp Lq dt = δpq (3.11)

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−s.t2 Lp(s.t) (3.12)

25

0 5 10 15 20−0.5

0

0.5

1

t

φp(t)

Ordem 0

Ordem 1

Ordem 2

Ordem 3

Ordem 4

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 utili-

zar 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−155

)2 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

0 5 10 15 20 25 30 35 40−0.5

0

0.5

1

t

Sin

al, f

(t)

Sinal Original

Sintese com 5 bases de Laguerre

Sintese com 8 bases de Laguerre

Sintese com 12 bases de Laguerre

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 bases

de Laguerre.

0 5 10 15 20 25 30 35 40−0.5

0

0.5

1

t

Sin

al, f(t

)

Sinal Original

Sintese com 100 bases de Laguerre

Figura 3.4: Síntese de um sinal utilizando 100 bases de Laguerre.

27

As escolhas para o número NL de bases de Laguerre e para o valor da escala

temporal s 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 =xNLmax

Tf(3.14)

sendo:

xNLmax o maior zero dentre os polinômios de Laguerre utilizados.

Uma pergunta razoável é: Por que utilizar especi�camente polinômios de La-

guerre? 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 decres-

cente, 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) =∞∑p=0

Epx(r) φp(t) (3.15)

sendo:

Ex(r, t) o campo elétrico na direção x em um determinado ponto r, função da

variável tempo t.

t = s.t simplesmente uma representação para o tempo escalado pelo fator s.

r a variável posição no espaço.

φp a base de Laguerre de ordem p.

Epx(r) o coe�ciente de ponderação da base Laguerre de ordem p para a repre-

sentação do campo elétrico Ex em uma posição especí�ca 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]:

∂U(r, t)

∂t= s

∞∑p=0

(0.5 Up(r) +

p−1∑k=0,p>0

Uk(r)

)φp (t) (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 TEM).

3.4.1 Modo TE

Primeiro, reapresentamos as equações que representam o Modo TEz:

∂Ex∂t

=1

ε

(∂Hz

∂y− Jx

)(3.17)

∂Ey∂t

=1

ε

(− ∂Hz

∂x− Jy

)(3.18)

29

∂Hz

∂t=

1

µ

(∂Ex∂y− ∂Ey

∂x

)(3.19)

O primeiro passo é utilizar as Equações (3.15) e (3.16) para reescrever as gran-

dezas 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∞∑p=0

(0.5Ep

x(x, y) +

p−1∑k=0

Ekx(x, y)

)φp(t)

=1

ε

∞∑p=0

∂Hpz

∂y(x, y) φp(t)−

Jx(x, y, t)

ε

(3.20)

s∞∑p=0

(0.5Ep

y(x, y) +

p−1∑k=0

Eky (x, y)

)φp(t)

= −1

ε

∞∑p=0

∂Hpz

∂y(x, y) φp(t)−

Jy(x, y, t)

ε

(3.21)

s∞∑p=0

(0.5Hp

z (x, y) +

p−1∑k=0

Hkz (x, y)

)φp(t)

= − 1

µ

∞∑p=0

(∂Ep

x

∂y(x, y)−

∂Epy

∂x(x, y)

)φp(t)

(3.22)

É 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 φq e integrar

em função de t no intervalo [0, Tf ), sendo Tf 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 q 6= p desaparecem, e quando q = p,

o termo resulta simplesmente em 1.

� Os termos relativos a densidade de corrente passam a ser:

Jqx(x, y) =

∫ TF

0

J(x)(x, y, t) φq(t) dt

Jqy (x, y) =

∫ TF

0

J(y)(x, y, t) φq(t) dt

(3.23)

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:

s

(0.5Eq

x(x, y) +

q−1∑k=0

Ekx(x, y)

)=

1

ε

∂Hqz

∂y(x, y) − Jqx(x, y)

ε(3.24)

s

(0.5Eq

y(x, y) +

q−1∑k=0

Eky (x, y)

)= −1

ε

∂Hqz

∂y(x, y)−

Jqy (x, y)

ε(3.25)

s

(0.5Hq

z (x, y) +

q−1∑k=0

Hkz (x, y)

)= − 1

µ

(∂Eq

x

∂y(x, y)−

∂Eqy

∂x(x, y)

)(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).

Eqx|i,j = CE

y [Hqz |i,j−Hq

z |i,j−1]− 2

sεJqx|i,j − 2

q−1∑k=0

Ekx |i,j (3.27)

Eqy |i,j = −CE

x [Hqz |i,j−Hq

z |i−1,j]−2

sεJqy |i,j − 2

q−1∑k=0

Eky |i,j (3.28)

31

Hqz |i,j = −CH

x

[Eqy |i+1,j−Eq

y |i,j]

+ CHy

[Eqx|i,j+1−Eq

y |i,j]− 2

q−1∑k=0

Hkz |i,j (3.29)

sendo:

CEx = 2

sε∆x.

CEy = 2

sε∆y.

CHx = 2

sµ∆x.

CHy = 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:

Eqx|i,j

[1

CEy

+ CHy + CH

y

]+ Eq

x|i,j−1

[−CH

y

]+ Eq

x|i,j+1

[−CH

y

]+

Eqy |i,j

[−CH

x

]+ Eq

y |i,j−1

[CHx

]+ Eq

y |i+1,j

[CHx

]+ Eq

y |i+1,j−1

[−CH

x

]= −∆y Jqx|i,j −

2

CEy

q−1∑k=0

Ekx |i,j + 2

q−1∑k=0

Hkz |i,j−1 − 2

q−1∑k=0

Hkz |i,j

(3.30)

Eqy |i,j

[1

CEx

+ CHx + CH

x

]+ Eq

y |i−1,j

[−CH

x

]+ Eq

y |i+1,j

[−CH

x

]+

Eqx|i,j

[−CH

y

]+ Eq

x|i−1,j

[CHy

]+ Eq

x|i−1,j+1

[−CH

y

]+ Eq

x|i,j+1

[CHy

]= −∆x Jqy |i,j −

2

CEx

q−1∑k=0

Eky |i,j − 2

q−1∑k=0

Hkz |i−1,j + 2

q−1∑k=0

Hkz |i,j

(3.31)

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. Isto dá

32

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 [Eq]2m× 1 = [Jq]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.

[Eq]2m× 1 o vetor coluna contendo todos os termos de Laguerre para os campos

elétricos.

[Jq]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 Jq, 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, [βq−1] é inexistente e a Equação (3.32) se torna simplesmente

[A] [E0] = [J0]. 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∂t

=1

ε

(∂Hz

∂x− Jy

)(3.33)

∂Hz

∂t=

1

µ

(∂Ey∂x

)(3.34)

Aplicando as Equações (3.15) e (3.16) obtemos as Equações (3.35) e (3.36):

s

∞∑p=0

(0.5Ep

y(x) +

q−1∑k=0

Eky (x)

)φp(t) = −1

ε

∞∑p=0

∂Hpz

∂x(x) φp(t)−

Jy(x, t)

ε(3.35)

s∞∑p=0

(0.5Hp

z (x) +

q−1∑k=0

Hkz (x)

)φp(t) = − 1

µ

∞∑p=0

∂Epy

∂x(x) φp(t) (3.36)

Aplicando o Procedimento de Teste Temporal de Galerkin:

s

(0.5Eq

y(x) +

q−1∑k=0

Eky (x)

)= −1

ε

∂Hqz

∂x(x)−

Jqy (x)

ε(3.37)

s

(0.5Hq

z (x) +

q−1∑k=0

Hkz (x)

)= − 1

µ

∂Eqy

∂x(x) (3.38)

Aplicando a discretização espacial por diferenças centradas e rearranjando termos

�camos com as Equações (3.39) e (3.40):

Eqy |i = −CE

x (Hqz |i −Hq

z |i−1)− 2

sεJqy |i − 2

q−1∑k=0

Eky |i (3.39)

Hqz |i = −CH

x

(Eqy |i+1 − Eq

y |i)− 2

q−1∑k=0

Hkz |i (3.40)

sendo:

CEx = 2

sε∆x.

CHx = 2

sµ∆x.

34

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:

Eqy |i[

1

CEx

+ CHx + CH

x

]+ Eq

y |i−1

[−CH

x

]+ Eq

y |i+1

[−CH

x

]= −∆x Jqy |i −

2

CEx

q−1∑k=0

Eky |i − 2

q−1∑k=0

Hkz |i−1 + 2

q−1∑k=0

Hkz |i,j

(3.41)

Assim como no caso anterior, chega-se a um sistema linear. O sistema para o

modo TEMx polar. y é da forma mostrada na Equação (3.42).

[A]m×m [Eq]m× 1 = [Jq]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.

[Eq]m× 1 o vetor coluna contendo todos os termos de Laguerre para os campos

elétricos.

[Jq]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 imple-

mentadas no método WLP-FDTD. Estas condições são adicionadas exatamente da

mesma forma tanto para problemas bidimensionais (modos TE e TM) quanto para

unidimensionais (modo TEM).

35

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 [Jq]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 = 0 ou

x = L para uma geometria de�nida no intervalo [0, L]) é 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):

∂Eqy(x, y)

∂x+s

v

(Eqy(x, y)

2+

q−1∑k=0

Eky (x, y)

)= 0 (3.43)

Aplicando a discretização espacial por diferenças centrais no ponto 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):

∂Ey∂x|1+ 1

2,j=

Eqy |2,j−Eq

y |1,j2

(3.44)

Eqy |1,j

[s

4v+

1

∆x

]+ Eq

y |2,j[s

4v− 1

∆x

]= − s

2v

q−1∑k=0

(Eky |2,j + Ek

y |1,j)

(3.45)

Aplicando o mesmo procedimento no outro extremo da geometria (x = L), ob-

temos a relação dada pela Equação (3.46), semelhante à Equação (3.45).

Eqy |L,j

[s

4v+

1

∆x

]+ Eq

y |L−1,j

[s

4v− 1

∆x

]= − s

2v

q−1∑k=0

(Eky |L,j + Ek

y |L−1,j

)(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 A e

seu vetor coluna β 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ís-

ticas de estabilidade incondicional e precisão do método WLP-FDTD. Para isso,

foram desenvolvidos algoritmos próprios utilizando linguagem MATLAB para os ca-

sos 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. Além

disso, o coe�ciente de correção linear RPearson entre as duas respostas é obtido. Este

coe�ciente está limitado entre -1 e 1, inclusive. RPearson = 1 denota correlação

positiva completa entre as grandezas, RPearson = −1 denota correlação negativa

completa e RPearson = 0 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, fundamental-

mente, em quatro loops :

� Loop de de�nição das bases de Laguerre (NL iterações): Neste loop, as equações

3.4 e 3.12 são utilizadas para gerar as bases de Laguerre. Neste mesmo loop já

são obtidos os coe�cientes de Laguerre Jqx e Jqy para a densidade de corrente

escolhida, utilizando a Equação (3.23).

38

� Loop de de�nição da matriz A (2m 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 TEMx polar. y) são

utilizadas para construir a matriz A do sistema linear apresentado na Equação

(3.32) (para o caso TEz) ou na Equação (3.42) (para o caso TEMx polar. y).

Existem 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 loop, três processos ocorrem; Primeiro,

o vetor β, dependente do somatório dos termos já descobertos, é atualizado

com os coe�cientes da última iteração, se esta é a primeira iteração, este ve-

tor simplesmente não existe. Em seguida, utilizando o vetor β, os termos

Jqx e Jqy descobertos no primeiro Loop, e a matriz A construída no segundo

Loop, o sistema linear ((3.32) ou (3.42)) é resolvido. Em seguida, utilizando os

coe�cientes Eq 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 descober-

tos 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 definição da Matriz A

Loop principal

Loop de retorno ao domínio do tempo

Dados deGeometria e Discretização

Geradorde

Gráficos

Figura 4.1: Fluxograma do Algoritmo Implementado para o Método WLP-FDTD.

4.3 Simulação Bidimensional: Sistema Físico Apre-

sentado 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. A componente de campo elétrico na direção y (Ey) é medida exatamente no

centro da geometria, denominada aqui como ponto p1. Esta con�guração exige a

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 1m na direção da coordenada x

e 0.1 m na direção da coordenada y. 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.001m 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 s e NL foram replicados diretamente

do artigo e são s = 6.07 × 1010 e NL = 150. A corrente elétrica é estabelecida em

x = 0.3m.

41

0 0.2 0.4 0.6 0.8 1 1.2

x 10−8

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

Forma de onda de Jy(t)

Tempo (segundos)

A/m

2

Figura 4.3: Forma da densidade de corrente elétrica inserida no problema bidimen-sional.

J(t) = e

(−(

t−TcTd

)2)

sin(2πfc(t− Tc)) (4.1)

Sendo fc = 1GHz, Td = 12fc

= 50 ns e Tc = 3Td = 150 ns.

O código do algoritmo próprio escrito em linguagem MATLAB para esta simu-

laçã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

0 0.2 0.4 0.6 0.8 1 1.2

x 10−8

−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Ey(t) no ponto p1

Tempo (segundos)

V/m

b

a

c

d

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 Ey nos pontos a, b, c e d obtidos tanto neste trabalho quanto na

referência. 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 poralgoritmo próprio, para os pontos a, b, c e d.

Ponto Grandeza Chung et al. Algoritmo Próprio

aEy −0.3± 0.1 V/m −0.28 V/mt 1.5± 0.2 ns 1.5 ns

bEy 1.5± 0.1 V/m 1.53 V/m

t 2.0± 0.2 ns 1.93 ns

cEy −1.5± 0.1 V/m −1.54 V/mt 2.4± 0.2 ns 2.35 ns

dEy 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

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(a) t = 0.878 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(b) t = 1.17 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(c) t = 1.46 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(d) t = 1.75 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(e) t = 2.05 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(f) t = 2.34 ns

Figura 4.5: Resultados para a simulação bidimensional utilizando o método WLP-FDTD.

45

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(g) t = 2.63 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(h) t = 3.22 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(i) t = 3.51 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(j) t = 3.80 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(k) t = 4.10 ns

Eixo x (cm)

Eix

o y

(cm

)

Ey(t)

10 20 30 40 50 60 70 80 90

9

8

7

6

5

4

3

2

1

0

V/m

−3

−2

−1

0

1

2

3

(l) t = 4.39 ns

Figura 4.5: Resultados para a simulação bidimensional utilizando o método WLP-FDTD.

46

Dos resultados apresentados é possível concluir que, primeiro, o método se man-

teve 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 y, o que signi�ca

dizer que esta densidade super�cial estará na direção y, como mostrado na Figura

4.6.

47

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 ele-

tromagné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) = −Ky(t± x

c)

2ε0 c(4.2)

sendo K a densidade super�cial de corrente, dada no SI por Am.

Os termos ± 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

Ey(t) = c µHz(t) = −Ky(t− x

c)

2ε0 c(4.3)

Caso a propagação seja no sentido negativo:

Ey(t) = −c µHz(t) = −Ky(t+ x

c)

2ε0 c(4.4)

O desenvolvimento do método WLP-FDTD para o caso TEM polarizado na dire-

çã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, con-

sistindo 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 cor-

rente ~K na direção y é colocado em x = 5 cm. A forma escolhida para a densidade

de corrente super�cial ~K é a descrita pela Equação (4.5) e mostrada na Figura 4.7.

0 0.5 1 1.5 2 2.5 3

x 10−7

−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Forma de onda de Ky(t)

Tempo (segundos)

A/m

Figura 4.7: Forma de Onda da Densidade Super�cial de Corrente Ky.

49

J = e−(

t−5TkTk

)2

(4.5)

sendo Tk = 10 ns.

O número de bases de Laguerre NL necessário para sintetizar esta forma de onda

foi obtido utilizando a Equação (3.13), e o fator de escala s foi obtido utilizando a

Equação (3.14). Os seguintes resultados foram obtidos: s = 1.13 × 109 e NL =

91. A discretização foi feita da seguinte maneira: ∆t = 0.037 ns e ∆x = 0.01 m.

Discretizações estas que também desrespeitam a condição de Courant, o que pode

ser visto facilmente pela Equação (2.58). Utilizando estes valores, obtemos Sc =

1.11 > 1.0. 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 WLP-

FDTD.

0 2 4 6 8 10

−4

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(a) t = 30.0 ns

0 2 4 6 8 10

−30

−25

−20

−15

−10

−5

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(b) t = 36.0 ns

Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando ométodo WLP-FDTD e dos apresentados por Feynman et al. [1].

50

0 2 4 6 8 10

−110

−100

−90

−80

−70

−60

−50

−40

−30

−20

−10

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(c) t = 42.0 ns

0 2 4 6 8 10

−200

−180

−160

−140

−120

−100

−80

−60

−40

−20

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(d) t = 48.0 ns

0 2 4 6 8 10

−220

−200

−180

−160

−140

−120

−100

−80

−60

−40

−20

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(e) t = 54.0 ns

0 2 4 6 8 10

−220

−200

−180

−160

−140

−120

−100

−80

−60

−40

−20

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(f) t = 60.0 ns

0 2 4 6 8 10

−220

−200

−180

−160

−140

−120

−100

−80

−60

−40

−20

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(g) t = 66.0 ns

0 2 4 6 8 10

−160

−140

−120

−100

−80

−60

−40

−20

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(h) t = 72.0 ns

Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando ométodo WLP-FDTD e dos apresentados por Feynman et al. [1].

51

0 2 4 6 8 10

−60

−50

−40

−30

−20

−10

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(i) t = 78.0 ns

0 2 4 6 8 10

−11

−10

−9

−8

−7

−6

−5

−4

−3

−2

−1

0

Distribuição do campo eletrico Ey ao longo da geometria

Ey (

V/m

)

Eixo x (cm)

Feynman et. al

WLP−FDTD

(j) t = 84.0 ns

Figura 4.8: Sobreposição dos resultados da simulação unidimensional utilizando omé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 |erro|max) encontrado dentre todos os

pontos da geometria, o ponto xerromax no qual este erro ocorre, o valor de pico

(em módulo) |Ey 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 RPearson entre as

duas respostas.

Tabela 4.2: Comparação entre os valores analíticos e os obtidos pelo método WLP-FDTD para campo elétrico no sistema físico unidimensional.

Instante |erro|max Ponto xerromax |Ey pico| RPearson

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

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 bidimen-

sional. Os resultados se revelaram próximos e, portanto, é possível concluir que o

algoritmo desenvolvido utilizando o método WLP-FDTD de fato apresenta estabi-

lidade 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 re-

centes 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. The Feynman Lectures on Phy-

sics, v. 2: Mainly Electromagnetics and Matter. Addison-Wesley, 1964.

[2] KASAL, R. Simulação de Supercondutores pelo Modelo do Estado Crítico. Dis-

sertação de mestrado, COPPE/UFRJ, 2006.

[3] SOTELO, G. Modelagem de Supercondutores Aplicada ao Projeto de Mancais

Magnéticos. Tese de doutorado, COPPE/UFRJ, 2007.

[4] SASS, F. Modelagem do Comportamento de Mancais Magnéticos Utilizando Fi-

tas e Blocos Maciços Supercondutores. Tese de doutorado, COPPE/UFRJ,

2015.

[5] DE SOUSA, W. Simulações de Regime Transitório de Limitadores de Corrente

de Curto-Circuito Supercondutores. Tese de doutorado, COPPE/UFRJ,

2015.

[6] OLIVEIRA SANTOS, B. Simulação de Sistemas Eletromagnéticos usando Mé-

todo de Elementos Finitos e Método de Diferenças Finitas no Tempo.

Trabalho de conclusão de curso, Universidade Federal do Rio de Janeiro,

2017.

[7] TAFLOVE, A., HAGNESS, S. Computational Electrodynamics: The Finite-

Di�erence Time-Domain Method. Artech House, 2000.

[8] VELOSO, L. Implementação do Método FDTD para as Equações de Maxwell

em Ambientes de Programação Paralela. Dissertação de mestrado, Uni-

versidade Federal do Rio de Janeiro, 2015.

[9] SCHNEIDER, J. Understanding the Finite-Di�erence Time-Domain Method.

2015.

[10] HA, M. EM Simulation Using the Laguerre-FDTD Scheme for Multiscale 3D

Interconnections. Tese de doutorado, Georgia Institute of Technology,

2011.

56

[11] CHUNG, Y., SARKAR, T., JUNG, B., et al. �An Unconditionally Stable

Scheme for the Finite-Di�erence Time-Domain Method�, IEEE Transac-

tions on Microwave Theory and Techniques, VOL 51. NO. 3, 2003.

[12] DE SOUSA, W., POLASEK, A., DIAS, R., et al. �Thermal�electrical ana-

logy for simulations of superconducting fault current limiters�, Cryogenics,

2014.

[13] SHLAGER, K., SCHNEIDER, J. �A Survey of the Finite-Di�erence Time

Domain Literature�. In: Ta�ove, A. (Ed.), Advances in Computational

Electrodynamics: The Finite-Di�erence Time-Domain, Artech House.

[14] CUMINATO, J., MENEGUETTE, M. Discretização de Equações Diferenci-

ais Parciais - Técnicas de Diferenças Finitas. Sociedade Brasileira de

Matematica, 2013.

[15] FARLOW, S. Partial Di�erential Equations for Scientists and Engineers. Dover

Publications, Inc., 1993.

[16] GRIFFITHS, D. Introduction to Electrodynamics. Prentice-Hall, 1999.

[17] KRAUS, J., CARVER, K. Electromagnetics. McGraw-Hill, 1973.

[18] INAN, U., MARSHALL, R. Numerical Electromagnetics: The FDTD Method.

Cambridge University Press, 2011.

[19] MUR, G. �Absorbing Boundary Conditions for the Finite-Di�erence Appro-

ximation of the Time-Domain Electromagnetic-Field Equations�, IEEE

Transactions on Electromagnetic Compatibility VOL. EMC 23. NO 4.,

1981.

[20] NAJM, F. Circuit Simulation. Wiley-IEEE, 2010.

[21] KUNG, F. Modeling of Electromagnetic Wave Propagation in Printed Circuit

Board and Related Structures. Tese de doutorado, Multimedia University,

Malasia, 2003.

[22] OUBRE, C. Finite-Di�erence Time-Domain Studies of the Optical Properties

of Metallodielectric Nanostructures. Tese de doutorado, Rice University,

2005.

[23] TENENBAUM, M., POLLARD, H. Ordinary Di�erential Equations. Dover

Publications, 1985.

[24] ZILL, D. Equações Diferenciais. Pearson, 2001.

57

[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. Numerical Recipes in C.

Cambridge University Press, 1992.

58

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) =∞∑

m=−∞

Um(x, t) =∞∑

m=−∞

Cm(t) e−jkmx (A.1)

Sendo j =√−1, km = 2π

λmo número de onda referente ao termo da Série de

Fourier, λ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) ejkmx (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 deslo-

camento 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 basica-

mente 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 q como q =U |t+1

i

U |ti

3. Escrever os termos da forma ejk∆x como senos e cossenos.

4. Analisar o valor assumido pelo fator de ampli�cação q para diferentes valores

de discretização espacial (por exemplo, ∆x) 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+1i = U |t−1

i −(v∆t

∆x

)(U |tiejk(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 q =U |t+1

i

U |ti,

obtemos a Equação (A.6):

U |t+1i =

U |tiq−(v∆t

∆x

)(U |ti e+jk ∆x − U |ti e−jk ∆x) (A.6)

Escrevendo os termos da forma ejk∆x como senos e cossenos e manipulando o

lado direito da Equação (A.6) obtemos a Equação (A.7)

U |t+1i =

(1

q− 2j

(v∆t

∆x

)sin(k∆x)

)U |ti = q U |ti (A.7)

Ou seja:

q =

(1

q− 2j

(v∆t

∆x

)sin(k∆x)

)(A.8)

De�nindo A = v∆t∆xsin(k ∆x), podemos reescrever a Equação (A.10) como:

q = −jA±√

1− A2 (A.9)

Lembrando, a exigência é |q|≤ 1 e portanto:

|q|2= [Re(q)]2 + [Im(q)]2 = 1− A2 + A2 = 1 (A.10)

Ou seja, a condição se traduz em |A|≤ 1, como o termo sin(k ∆x) na de�nição

de A é no máximo 1, a condição de estabilidade para um método do tipo leapfrog

unidimensional se torna:

60

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 RPearson (também conhecido como coe�ciente de

correlação de Pearson) é dado, para um par de grandezas (xi, yi) para i = 1, 2 · · ·Npela fórmula:

RPearson =

∑(xi − x)(yi − y)√∑

(xi − x)2√∑

(yi − y)2(B.1)

Sendo x a média de x e y a média de y. [26]

62

Apêndice C

Código Desenvolvido para o Caso

Bidimensional utilizando

WLP-FDTD

1 %vasculhar a malha de pra cima , esquerda pra direita.

2 %Hz na extrema direita esta fora da malha

3 %Hz na extremidade superior tambem esta fora

4 %Ex na extremidade direita esta fora

5 %Ey na extremidade superior esta fora

6 sizet = 400;

7 lim_t = 11.71*(10^ -9); %s

8 t = linspace(0,lim_t ,sizet);

9 p = 151; % numero de polinomios de laguerre a serem

utilizados , p = (( ordem do maior polinomio) -1)

10 deltax = 0.01; %m

11 deltay = 0.001; %m

12 l = 1; %m

13 h = 0.1; %m

14 ptslinha = round(l/deltax);

15 ptscoluna = round(h/deltay);

16 npontos = ptslinha*ptscoluna;

17 s = 6.07*(10^10); %fator de escala temporal

18 mi0 = 4*pi*(10^ -7); %H/m

19 eps0 = 8.854*(10^ -12); % F/m

20 c0 = 2.998*(10^8); % velocidade da luz

21

22 weight = exp(-s*t*0.5);

63

23 Jq = zeros(1,p); %lista dos coeficientes de laguerre

que descrevem a corrente nos pontos de entrada , a

ser completada

24 Jp = sparse(zeros (2*( npontos),p)); %lista dos

coeficientes que descrevem correntes em todos os

pontos da geometria

25

26 %dados do pulso , todos diretamente do Chung

27 fc = 10^9;

28 Td = 1/(2*fc);

29 Tc = 3*Td;

30 Jinput = exp(-((t-Tc)/Td).^2).*sin (2*pi.*fc.*(t-Tc));

31

32 laguerre = zeros(p,sizet); %polinomios de laguerre no

dominio do tempo

33 laguerre (1,:) = 1; %ordem 0

34 laguerre (2,:) = 1-(s*t); % ordem 1

35

36 %determinando os polinomios de ordem maior que 2, no

dominio do tempo

37 %lembrando , tamanho laguerre = (p,sizet)

38 if p>2

39 for q =3:p

40 laguerre(q,:)=(( laguerre(q-1,:) .*((2*(q-1))-1-s*t)) -(q

-2).* laguerre(q-2,:))/(q-1);

41 end

42 end

43

44 %lista das bases de laguerre (polinomios multiplicados

pelo fator weight) no dominio do tempo

45 %lembrando tamanho phiq = (p,sizet)

46 for i=1:p

47 disp( sprintf( 'Obtendo bases de laguerre: %d de %d', i

,p))

48 phiq(i,:) = (weight).* laguerre(i,:);

49

50 % obtendo termo a ser integrado para completar a lista

de coeficientes de laguerre que definem a corrente

51 Jq_pre_integ(i,:) = Jinput .*phiq(i,:);

64

52 end

53 clear weight %limpando a variavel weight pra liberar

memoria

54

55 %integrando os termos obtidos anteriormente ,

verdadeiramente obtendo os coeficientes de laguerre

que definem o sinal de corrente.

56 for i=1:p

57 disp( sprintf( 'Obtendo coeficientes de laguerre para a

corrente: %d de %d', i,p))

58 %lembrando: Jq =(1,p)

59 Jq(1,i) = trapz(s*t,Jq_pre_integ(i,:));

60 end

61

62 clear Jq_pre_integ %limpando o termo pre -integracao ,

ele e inutil e so ocupa memoria

63

64 %ao criar os campos E para um sistema 2D tenho que

decidir qual e a ordem dos campos eletricos , ela

sera a seguinte:

65 %Primeiro: todos os campos Ex, percorrendo a geometria

linha a linha de baixo pra cima

66 %Depois: todos os campos Ey, percorrendo a geometria

linha a linha de baixo pra cima

67 %numeros uteis:

68 %numero de linhas = npontos

69 %todos os Ex = npontos , todos os Ey = npontos , todos os

campos eletricos = 2*( npontos)

70 %todos os Hz = todos os campos magneticos = npontos

71 %como cada campo vai ser a soma de p coeficientes de

laguerre:

72 %matriz E tera dimensoes p x 2*( npontos)

73 %matriz H tera dimensoes p x npontos

74

75 E = sparse(zeros(p,2*( npontos)));

76 H = sparse(zeros(p,( npontos)));

77 %lembrando: definicao de beta na eq 40 do Chung

78 %beta tem o mesmo tamanho da matriz A

79 beta = sparse(zeros (1,2*( npontos)));

65

80

81 %para executar mais de uma vez com as mesmas matrizes e

agilizar o processo

82 disp(' Checando se existe alguma matriz A e de relacoes

aproveitavel ')

83 disp('Se fez alguma mudanca recente em A lembre -se de

limpar as variaveis ')

84 flag_A = 1;

85 if ((exist ('A') & exist ('relat_betasomaH ') & exist('

relat_betasomaE ')) == 0) | (max(size(A)) ~= 2*(

npontos)) %se matrizes A e de relacoes nao existem

ainda ...

86 disp('Matriz A nao detectada ')

87 disp('criando matriz A')

88 A = sparse(zeros (2*( npontos)));

89 % A = (zeros (2*( npontos)));

90 % esses vetores de relacao sao um jeito rapido de

implementar , por exemplo , a equacao 40 do Chung (

dentre outras)

91 disp('criando matrizes de relacao ')

92 relat_betasomaH = sparse(zeros (2*( npontos),npontos));

93 relat_betasomaE = sparse(zeros (2*( npontos)));

94 relat_H_E = sparse(zeros(npontos ,2*( npontos)));

95 flag_A = 0;

96 end

97

98 % a cada iteracao do metodo sera necessario somar todos

os coefs das iteracoes anteriores

99 %como so trato Hz , somaH tem metade de somaE (que

inclui tanto Ex quanto Ey)

100 somaH = sparse(zeros (1,( npontos)));

101 somaE = sparse(zeros (1,2*( npontos)));

102

103 %constantes definidas tambem no Chung

104 CEx= 2/( deltax*eps0*s);

105 CEy = 2/( deltay*eps0*s);

106 CHx = 2/( deltax*mi0*s);

107 CHy = 2/( deltay*mi0*s);

108

66

109 % existem duas componentes , Ex e Ey, cada uma delas

responsavel por metade da matriz A (npontos)

110 % preciso de dois loops , um pra cada metade da matriz

111 % [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

112 %criarei constantes pra lembrar as distancias (os

shifts) dentro da matriz

113 %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)

114 % Ey39 = distancia entre os termos da eq 39 e os termos

Ey

115 % Ex40 = distancia entre os termos da eq 40 e os termos

Ex

116

117 %outra observacao e que a determinacao de J e muito

mais complicada , como estou tentando o caso

apresentado por Chung:

118 %Jx = 0 sempre

119 %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.

120 %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 sfvizinhox = 1; %exemplo , shift do Ex no ponto (x=1,y

=1) pro Ex no ponto (x=2,y=1)

125 sfvizinhoy = ptslinha; %exemplo , shift do Ex no ponto (

x=1,y=1) pro Ex no ponto (x=1,y=2)

126

127 sfdiagsupesq = sfvizinhoy -sfvizinhox;

128 sfdiagsupdir = sfvizinhoy+sfvizinhox;

129 sfdiaginfesq = -sfvizinhoy -sfvizinhox;

130 sfdiaginfdir = -sfvizinhoy+sfvizinhox;

131

67

132 ptsquadrado = npontos;

133 shift_ExEy = ptsquadrado;

134

135 pontos_geometria = 1: ptsquadrado;

136 pontos_inferiores = 2:ptslinha -1;

137 pontos_superiores = ptsquadrado -ptslinha +2: ptsquadrado

-1; %sem pontas

138 pontos_esquerda = 1+ ptslinha:ptslinha:ptsquadrado -(2*

ptslinha)+1; %sem pontas

139 pontos_direita = 2* ptslinha:ptslinha:ptsquadrado -

ptslinha; %sem pontas

140 pontos_pontas = [1 ptslinha ptsquadrado -ptslinha +1

ptsquadrado ];

141

142 pontos_limite = horzcat(pontos_inferiores ,

pontos_superiores ,

143 pontos_esquerda ,pontos_direita ,pontos_pontas);

144 pontos_interior = setdiff(pontos_geometria ,

pontos_limite);

145

146 %Jy , todos os pontos com corrente estao na segunda

metade do vetor

147 pontos_com_corrente = ptsquadrado+ptslinha +1+(30):

ptslinha:

148 1+2* ptsquadrado -ptslinha -(70);

149

150 disp(sprintf( 'Construindo Jp'))

151 for k = pontos_com_corrente

152 Jp(k,:) = -deltax*Jq;

153 end

154

155 disp( sprintf( 'Construindo matriz A'))

156

157 for i=pontos_interior

158 disp( sprintf('Construindo matriz A: %d de %d',i,

ptsquadrado))

159 %eq 39 do Chung

160 Ey_i = i+shift_ExEy;

161 if flag_A == 1

68

162 disp('Parece que ja existe uma matriz A no workspace , o

programar tentara utiliza -la')

163 break

164 end

165 %Elementos Ex, normal

166 A(i,i) = (1/CEy)+CHy+CHy;

167 A(i,i+sfvizinhoy) = -CHy;

168 A(i,i-sfvizinhoy) = -CHy; % ok

169

170 A(i,Ey_i) = -CHx;

171 A(i,Ey_i+sfvizinhox) = CHx;

172 A(i,Ey_i+sfdiaginfdir) = -CHx;

173 A(i,Ey_i -sfvizinhoy) = CHx;

174

175

176 %lembrando: relat_betasomaE e relat_betasomaH tem

ndelinhas = ndelinhas de beta (2* npontos)

177 %relat_betasomaH tem menos colunas (tantas colunas

quanto H tem linhas)

178 relat_betasomaH(i,i-sfvizinhoy) = 2;

179 relat_betasomaH(i,i) = -2;

180 relat_betasomaE(i,i) = -2/CEy;

181

182 %Elementos Ey, normal

183 A(Ey_i ,Ey_i) = (1/CEx)+CHx+CHx;

184 A(Ey_i ,Ey_i -sfvizinhox) = -CHx;

185 A(Ey_i ,Ey_i+sfvizinhox) = -CHx;

186

187 A(Ey_i ,i+sfvizinhoy) = CHy;

188 A(Ey_i ,i+sfdiagsupesq) = -CHy;

189 A(Ey_i ,i) = -CHy;

190 A(Ey_i ,i-sfvizinhox) = CHy;

191

192 %equivalente a equacao 40, Chung (existe um sinal

errado no artigo)

193 relat_betasomaH(Ey_i ,i) = 2;

194 relat_betasomaH(Ey_i ,i-sfvizinhox) = -2;

195 relat_betasomaE(Ey_i ,Ey_i) = -2/CEx;

196

69

197 relat_H_E(i,i) = -CHy;

198 relat_H_E(i,i+sfvizinhoy) = CHy;

199 relat_H_E(i,Ey_i) = CHx;

200 relat_H_E(i,Ey_i+sfvizinhox) = -CHx;

201 end

202 for i =pontos_esquerda

203 Ey_i = i+shift_ExEy;

204 if flag_A == 1

205 break

206 end

207 %Elementos Ex, normal

208 A(i,i) = (1/CEy)+CHy+CHy;

209 A(i,i+sfvizinhoy) = -CHy;

210 A(i,i-sfvizinhoy) = -CHy;

211

212 A(i,Ey_i) = -CHx;

213 A(i,Ey_i+sfvizinhox) = CHx;

214 A(i,Ey_i+sfdiaginfdir) = -CHx;

215 A(i,Ey_i -sfvizinhoy) = CHx;

216 %lembrando: relat_betasomaE e relat_betasomaH tem

ndelinhas = ndelinhas de beta (2* npontos)

217 %relat_betasomaH tem menos colunas (tantas colunas

quanto H tem linhas)

218 relat_betasomaH(i,i-sfvizinhoy) = 2;

219 relat_betasomaH(i,i) = -2;

220 relat_betasomaE(i,i) = -2/CEy;

221

222 %Ey , %ABC

223 A(Ey_i ,Ey_i) = 0.25*(s/c0)+(1/ deltax);

224 A(Ey_i ,Ey_i+sfvizinhox) = 0.25*(s/c0) -(1/ deltax);

225 relat_betasomaE(Ey_i ,Ey_i) = -0.5*(s/c0);

226 relat_betasomaE(Ey_i ,Ey_i+sfvizinhox) = -0.5*(s/c0);

227 relat_betasomaH(i,:) = 0;

228

229 %Relacao de H aos campos , normal

230 relat_H_E(i,i) = -CHy;

231 relat_H_E(i,i+sfvizinhoy) = CHy;

232 relat_H_E(i,Ey_i) = CHx;

233 relat_H_E(i,Ey_i+sfvizinhox) = -CHx;

70

234 end

235 for i =1:1 %ponta inferior esquerda

236 Ey_i = i+shift_ExEy;

237 if flag_A == 1

238 break

239 end

240 %Ex , %PEC

241 A(i,:) = 0;

242 A(i,i) = (1/CEy)+CHy+CHy;

243 relat_betasomaE(i,:) = 0;

244 relat_betasomaH(i,:) = 0;

245

246 %Ey , %ABC

247 A(Ey_i ,Ey_i) = 0.25*(s/c0)+(1/ deltax);

248 A(Ey_i ,Ey_i+sfvizinhox) = 0.25*(s/c0) -(1/ deltax);

249 relat_betasomaE(Ey_i ,Ey_i) = -0.5*(s/c0);

250 relat_betasomaE(Ey_i ,Ey_i+sfvizinhox) = -0.5*(s/c0);

251 relat_betasomaH(Ey_i ,:) = 0;

252

253 %Relacao de H aos campos , normal

254 relat_H_E(i,i) = -CHy;

255 relat_H_E(i,i+sfvizinhoy) = CHy;

256 relat_H_E(i,Ey_i) = CHx;

257 relat_H_E(i,Ey_i+sfvizinhox) = -CHx;

258 end

259 for i =pontos_inferiores

260 Ey_i = i+shift_ExEy;

261 if flag_A == 1

262 break

263 end

264 %Ex , %PEC

265 A(i,:) = 0;

266 A(i,i) = (1/CEy)+CHy+CHy;

267 relat_betasomaE(i,:) = 0;

268 relat_betasomaH(i,:) = 0;

269

270 %Elementos Ey, normal

271 A(Ey_i ,Ey_i) = (1/CEx)+CHx+CHx;

272 A(Ey_i ,Ey_i -sfvizinhox) = -CHx;

71

273 A(Ey_i ,Ey_i+sfvizinhox) = -CHx;

274

275 A(Ey_i ,i+sfvizinhoy) = CHy;

276 A(Ey_i ,i+sfdiagsupesq) = -CHy;

277 A(Ey_i ,i) = -CHy;

278 A(Ey_i ,i-sfvizinhox) = CHy;

279

280 %equivalente a equacao 40, Chung (existe um sinal

errado no artigo)

281 relat_betasomaH(Ey_i ,i) = 2;

282 relat_betasomaH(Ey_i ,i-sfvizinhox) = -2;

283 relat_betasomaE(Ey_i ,Ey_i) = -2/CEx;

284

285 %Relacao entre H e os campos E

286 relat_H_E(i,i) = -CHy;

287 relat_H_E(i,i+sfvizinhoy) = CHy;

288 relat_H_E(i,Ey_i) = CHx;

289 relat_H_E(i,Ey_i+sfvizinhox) = -CHx;

290

291 %J

292 Jp(i,:) = 0;

293 Jp(Ey_i ,:) = 0;

294 end

295 for i =ptslinha:ptslinha %ponta inferior direita

296 Ey_i = i+shift_ExEy;

297 if flag_A == 1

298 break

299 end

300 %Ex , %Nao existe Ex aqui

301 A(:,i) = 0;

302 A(i,i) = (1/CEy)+CHy+CHy;

303 relat_betasomaE(i,i) = 0;

304 relat_betasomaE(i,i+sfvizinhoy) = 0;

305 relat_betasomaH(i,:) = 0;

306

307 %Ey , %ABC

308 A(Ey_i ,Ey_i) = 0.25*(s/c0)+(1/ deltax);

309 A(Ey_i ,Ey_i -sfvizinhox) = 0.25*(s/c0) -(1/ deltax);

310 relat_betasomaE(Ey_i ,Ey_i) = -0.5*(s/c0);

72

311 relat_betasomaE(Ey_i ,Ey_i -sfvizinhox) = -0.5*(s/c0);

312 relat_betasomaH(Ey_i ,:) = 0;

313

314 %Relacao de H aos campos , nao existe Hz

315 relat_H_E(i,:) = 0;

316 end

317 for i =pontos_direita

318 Ey_i = i+shift_ExEy;

319 if flag_A == 1

320 break

321 end

322 % Ex, nao existe

323 A(:,i) = 0;

324 A(i,i) = (1/CEy)+CHy+CHy;

325 relat_betasomaE(i,:) = 0;

326 % relat_betasomaE(i,i) = 0;

327 % relat_betasomaE(i,i+sfvizinhoy) = 0;

328

329 %Ey , ABC

330 A(Ey_i ,Ey_i) = 0.25*(s/c0)+(1/ deltax);

331 A(Ey_i ,Ey_i -sfvizinhox) = 0.25*(s/c0) -(1/ deltax);

332 relat_betasomaE(Ey_i ,Ey_i) = -0.5*(s/c0);

333 relat_betasomaE(Ey_i ,Ey_i -sfvizinhox) = -0.5*(s/c0);

334 relat_betasomaH(Ey_i ,:) = 0;

335

336 %Relacao de H aos campos

337 relat_H_E(i,:) = 0;

338 relat_H_E(i,i) = 0; %nao existe Ex nem Hz

339 relat_H_E(i,i+sfvizinhoy) = 0;

340 relat_H_E(i,Ey_i+sfvizinhox) = 0;

341 relat_H_E(i,Ey_i) = 0;

342 end

343 for i =ptsquadrado:ptsquadrado %ponta superior direita

344 Ey_i = i+shift_ExEy;

345 if flag_A == 1

346 break

347 end

348 %Ex Nao existe

349 A(:,i) = 0;

73

350 A(i,i) = (1/CEy)+CHy+CHy;

351 relat_betasomaE(i,:) = 0;

352 relat_betasomaH(i,:) = 0;

353

354 %Ey Nao existe

355 A(:,Ey_i) = 0;

356 A(Ey_i ,Ey_i) = (1/CEx)+CHx+CHx;

357 relat_betasomaE(Ey_i ,:) = 0;

358 relat_betasomaH(Ey_i ,:) = 0;

359

360 %Relacao de H aos campos , H nao existe

361 relat_H_E(i,:) = 0;

362 end

363 for i =pontos_superiores

364 Ey_i = i+shift_ExEy;

365 if flag_A == 1

366 break

367 end

368 %Ex , %PEC

369 A(i,:) = 0;

370 A(i,i) = (1/CEy)+CHy+CHy;

371 relat_betasomaE(i,:) = 0;

372 relat_betasomaH(i,:) = 0;

373

374 %Ey Nao existe

375 A(:,Ey_i) = 0;

376 A(Ey_i ,Ey_i) = (1/CEx)+CHx+CHx;

377 relat_betasomaE(Ey_i ,:) = 0;

378 relat_betasomaH(Ey_i ,:) = 0;

379

380 %Relacao de H aos campos , H nao existe

381 relat_H_E(i,:) = 0;

382

383 %J

384 Jp(i,:) = 0;

385 Jp(Ey_i ,:) = 0;

386 end

387 for i =ptsquadrado -ptslinha +1: ptsquadrado -ptslinha +1 %

ponta superior esquerda

74

388 Ey_i = i+shift_ExEy;

389 if flag_A == 1

390 break

391 end

392 % %Ex, %PEC

393 A(i,:) = 0;

394 A(i,i) = (1/CEy)+CHy+CHy;

395 relat_betasomaE(i,:) = 0;

396 relat_betasomaH(i,:) = 0;

397

398 %Ey Nao existe

399 A(:,Ey_i) = 0;

400 A(Ey_i ,Ey_i) = (1/CEx)+CHx+CHx;

401 relat_betasomaE(Ey_i ,:) = 0;

402 relat_betasomaH(Ey_i ,:) = 0;

403

404 %Relacao de H aos campos , H nao existe

405 relat_H_E(i,:) = 0;

406 end

407 for j = 1:p

408 disp( sprintf('Loop principal: %d de %d', j,p))

409 j;

410 % j = 1, ordem = 0, j = 2, ordem = 1...

411 if j==1

412 %lembrando H =(p,2*( npontos))

413 %logo , sum(H) = (1,( npontos))

414 %lembrando beta = (1,2*( npontos))

415 rightside = Jp(:,j)

416 end

417 if j > 1

418 %lembrando H =(p,npontos)

419 %logo , sum(H) = (1,npontos)

420 %lembrando beta = (1,npontos)

421 somaH = sum(H(1:j,:));

422 somaE = sum(E(1:j,:));

423

424 beta (1,:) = transpose(relat_betasomaH*transpose(somaH)

+ relat_betasomaE*transpose(somaE));

425 rightside = Jp(:,j)+transpose(beta);

75

426 end

427 E(j,:) = transpose(mldivide(A,rightside));

428 H(j,:) = transpose(relat_H_E*transpose(E(j,:)) - 2*(

transpose(somaH)));

429 end

430

431 % Descobertos os coeficientes , voltando pro tempo

432 Exrecuperado = zeros(sizet ,( npontos));

433 Eyrecuperado = zeros(sizet ,( npontos));

434 Hzrecuperado = zeros(sizet ,( npontos));

435 disp(sprintf('Voltando pro dominio do tempo ...Ex'))

436 for pos = 1:( npontos)

437 for laguerre = 1:p

438 Exrecuperado (:,pos) = Exrecuperado (:,pos) + E(laguerre ,

pos)*transpose(phiq(laguerre ,:));

439 Hzrecuperado (:,pos) = Hzrecuperado (:,pos) + H(laguerre ,

pos)*transpose(phiq(laguerre ,:));

440 end

441 end

442

443 disp(sprintf('Voltando pro dominio do tempo ...Ey'))

444 for pos = (( npontos)+1:2*( npontos))

445 for laguerre = 1:p

446 Eyrecuperado (:,pos -( npontos)) = Eyrecuperado (:,pos -(

npontos)) + E(laguerre ,pos)*transpose(phiq(laguerre

,:));

447 end

448 end

76

Apêndice D

Código Desenvolvido para o Caso

Unidimensional utilizando

WLP-FDTD

1 sizet = 8100;

2 lim_t = 30e-8;

3 t = linspace(0,lim_t ,sizet);

4 p = 91;

5 deltat = t(2)-t(1)

6 deltax = 0.01;

7 l = 10;

8 npontos = (l/deltax);

9 s = 1.1333e+09;

10

11 mi0 = 4*pi*(10^ -7); %H/m

12 eps0 = 8.854187817*(10^ -12); % F/

13 c0 = 2.998*(10^8); % velocidade da luz

14

15 weight = exp(-s*t*0.5);

16 Jq = zeros(1,p); %lista dos coeficientes de laguerre que

descrevem a corrente no ponto de entrada

17 Jp = zeros(1,npontos);

18

19 Tk = 1e-8;

20 Jinput = exp(-((t-5*Tk)/Tk).^2);

21 pt_entrada = npontos /2;

22

77

23 Eanalitico = zeros(sizet ,npontos);

24 Hanalitico = zeros(sizet ,npontos);

25 for i = 1: npontos /2

26 x = i*deltax;

27 tdeslocado = t+((x-( deltax*npontos /2))/c0);

28 Eanalitico (:,i) = -exp(-(( tdeslocado -5*Tk)/Tk)

.^2) /(2*c0*eps0);

29 end

30 for i = 1+ npontos /2: npontos

31 x = i*deltax;

32 tdeslocado = t-((x-( deltax*npontos /2))/c0);

33 Eanalitico (:,i) = -exp(-(( tdeslocado -5*Tk)/Tk)

.^2) /(2*c0*eps0);

34 end

35

36 laguerre = zeros(p,sizet); %polinomios de laguerre no

dominio do tempo

37 laguerre (1,:) = 1; %ordem 0

38 laguerre (2,:) = 1-(s*t); % ordem 1

39

40 %lembrando , tamanho laguerre = (p,sizet)

41 %determinando os polinomios de ordem maior que 2, no

dominio do tempo

42 if p>2

43 for q =3:p

44 laguerre(q,:)=(( laguerre(q-1,:) .*((2*(q

-1))-1-s*t)) -(q-2).* laguerre(q-2,:))/(

q-1);

45 end

46 end

47

48 %lista das bases de laguerre no dominio do tempo

49 %lembrando tamanho phiq = (p,sizet)

50 for i=1:p

51 disp( sprintf( 'Obtendo bases de laguerre: %d de

%d', i,p))

52 phiq(i,:) = (weight).* laguerre(i,:);

53 Jq_pre_integ(i,:) = Jinput .*phiq(i,:);

54 end

78

55 clear weight

56

57 for i=1:p

58 disp( sprintf( 'Obtendo coeficientes de laguerre

para a corrente: %d de %d', i,p))

59 %lembrando Jq =(1,p)

60 Jq(1,i) = trapz(s*t,Jq_pre_integ(i,:));

61 end

62

63 clear Jq_pre_integ;

64

65 A = sparse(zeros(npontos));

66 E =zeros(p,npontos);

67 H = zeros(p,npontos);

68 beta = zeros(1,npontos);

69 somaH = zeros(1,npontos);

70 somaE = zeros(1,npontos);

71

72 % os vetores de relacao sao um jeito rapido de

implementar , por exemplo , equacoes da forma 40 do

Chung

73 relat_betasomaH = zeros(npontos);

74 relat_betasomaE = zeros(npontos);

75 relat_H_E = zeros(npontos);

76

77 CE= 2/( deltax*eps0*s);

78 CH= 2/( deltax*mi0*s);

79

80 for i=2: npontos -1

81 %equivalente a equacao 39, Chung. modificada pra

TEMx

82 A(i,i) = ((1/CE)+2*CH);

83 A(i,i-1) = -CH;

84 A(i,i+1) = -CH;

85 %equivalente a equacao 40, Chung , modifcada pra

TEMx

86 relat_betasomaH(i,i) = 2;

87 relat_betasomaH(i,i-1) = -2;

88 relat_betasomaE(i,i) = -2/CE;

79

89

90 %equivalente a equacao 34, Chung. Modificada pra

TEMx

91 relat_H_E(i,i) = CH;

92 relat_H_E(i,i+1) = -CH;

93 end

94 relat_betasomaH = sparse(relat_betasomaH);

95 relat_betasomaE = sparse(relat_betasomaE);

96

97 %adicionando fronteiras ABC nas duas extremidades

98 A(1,1) = 0.25*(s/c0)+(1/ deltax);

99 A(1,2) = 0.25*(s/c0) -(1/ deltax);

100

101 A(end ,end -1) = 0.25*(s/c0) -(1/ deltax);

102 A(end ,end) = 0.25*(s/c0)+(1/ deltax);

103

104 relat_H_E (1,1) = CH;

105 relat_H_E (1,2) = -CH;

106

107 relat_betasomaH (1,1) = 2;

108 relat_betasomaH(end ,end) = 2;

109

110 relat_betasomaE (1,1) = -0.5*(s/c0);

111 relat_betasomaE (1,2) = -0.5*(s/c0);

112 relat_betasomaE(end ,end) = -0.5*(s/c0);

113 relat_betasomaE(end ,end -1) = -0.5*(s/c0);

114

115 %lembrando Jq = (1,p)

116 %%%%%%%% Loop Principal

117 for j = 1:p

118 j;

119 disp( sprintf('Loop principal: %d de %d', j,p))

120 % j = 1, ordem = 0, j = 2, ordem = 1...

121 if j==1

122 %lembrando H =(p,npontos)

123 %logo , sum(H) = (1,npontos)

124 %lembrando beta = (1,npontos)

125 somaH = H(1,:);

126 somaE = E(1,:);

80

127 Jp(1,npontos) = 0; %PEC ou ABC , as duas

tem essa caracteristica relativa a Jp

128 Jp(1,1) = 0; %PEC ou ABC

129 Jp(1, pt_entrada) = -Jq(1,j);

130 rightside = transpose(Jp);

131 end

132 if j > 1

133 %lembrando H =(p,npontos)

134 %logo , sum(H) = (1,npontos)

135 %lembrando beta = (1,npontos)

136 somaH = sum(H(1:j,:));

137 somaE = sum(E(1:j,:));

138 beta (1,:) = relat_betasomaH*transpose(

somaH) + relat_betasomaE*transpose(

somaE);

139 Jp(1,npontos) = 0; %PEC

140 Jp(1,1) = 0; %PEC

141 Jp(1, pt_entrada) = -Jq(1,j);

142 rightside = transpose(Jp+beta);

143 end

144 E(j,:) = mldivide(A,rightside);

145 H(j,:) = relat_H_E*transpose(E(j,:)) - 2*(

transpose(somaH));

146 end

147

148 Erecuperado = zeros(sizet ,npontos);

149 Hrecuperado = zeros(sizet ,npontos);

150 for pos = 1: npontos

151 for laguerre = 1:p

152 Erecuperado (:,pos) = Erecuperado (:,pos) + E(

laguerre ,pos)*transpose(phiq(laguerre ,:));

153 Hrecuperado (:,pos) = Hrecuperado (:,pos) + H(

laguerre ,pos)*transpose(phiq(laguerre ,:));

154 end

155 end

81