Embed Size (px)
Citation preview
Universidade Federal do Paraná
Setor de Tecnologia
Departamento de Engenharia Elétrica - DELT
IMPLEMENTAÇÃO DE BIPOLOS DISCRETOS NO MÉTODO DE
ANÁLISE ELETROMAGNÉTICA POR ELEMENTOS FINITOS NO
DOMÍNIO DO TEMPO USANDO A TÉCNICA DE NEWMARK
CURITIBA DEZEMBRO 2005
ISMAEL CHIAMENTI
IMPLEMENTAÇÃO DE BIPOLOS DISCRETOS NO MÉTODO DE
ANÁLISE ELETROMAGNÉTICA POR ELEMENTOS FINITOS NO
DOMÍNIO DO TEMPO USANDO A TÉCNICA DE NEWMARK
CURITIBA DEZEMBRO 2005
Projeto Final de Curso apresentado à disciplina Projeto de Graduação como requisito parcial à conclusão do Curso de Engenharia Elétrica, Setor de Tecnologia, Departamento de Engenharia Elétrica, Universidade Federal do Paraná.
Orientador: Prof. Dr. Wilson Arnaldo Artuzi Junior
ii
Pois disseste: “O Eterno é meu refúgio”, e
fizeste tua a morada do Altíssimo.
Salmo 91,9
iii
AGRADECIMENTOS
Bendito seja o D’US de Avraham, Yitschac e Yaacov, pelo livre arbítrio!
Ao Prof. Dr. Wilson A. Artuzi Jr. pela orientação e paciência.
À minha mãe Marta e o meu pai Vilson pelo apoio, amor e educação.
À minha tia Edite por acolher-me em sua casa e ser minha segunda mãe.
Ao meu irmão de sangue Moisés pela camaradagem e amizade.
À todos os Professores e Funcionários do Departamento de Engenharia Elétrica da
Universidade Federal do Paraná pelo ensino e serviços prestados.
Ao LACTEC, especialmente as unidades de Metrologia Elétrica – UTME e de
Compatibilidade Eletromagnética – UTCE, pelo apoio e infra-estrutura
disponibilizados.
iv
SUMÁRIO
LISTA DE FIGURAS..................................................................................................vi
LISTA DE TABELAS...............................................................................................viii
LISTA DE SÍMBOLOS E SIGLAS............................................................................ix
RESUMO......................................................................................................................xi
ABSTRACT.................................................................................................................xii
1 INTRODUÇÃO.........................................................................................................1
2 PARÂMETROS DE ADMITÂNCIA E DE ESPALHAMENTO.........................3
2.1 INTRODUÇÃO......................................................................................................3
2.2 PARÂMETROS Y: ADMITÂNCIAS...................................................................4
2.3 PARÂMETROS DE ESPALHAMENTO OU PARÂMETROS S....................6
2.3.1 ADMITÂNCIAS A PARTIR DOS PARÂMETROS S........................................9
3 APROXIMAÇÃO DA ADMITÂNCIA POR UMA SÉRIE DE
FRAÇÕES PARCIAS...........................................................................................11
3.1 INTRODUÇÃO....................................................................................................11
3.2 APROXIMAÇÃO DA RESPOSTA EM FREQÜÊNCIA DA
ADMITÂNCIA.....................................................................................................11
3.2.1 VERIFICAÇÃO DA APROXIMAÇÃO ............................................................14
4 MÉTODO FDTD e FETD.......................................................................................19
4.1 INTRODUÇÃO....................................................................................................19
4.2 DEFINIÇÕES GEOMÉTRICAS E DISCRETIZAÇÃO ESPACIAL............20
4.3 MATRIZES DA EQUAÇÃO DE ONDA...........................................................23
4.4 MODELAGEM DE SUPERFÍCIES CONDUTORAS.....................................26
4.5 MODELAGEM DE FIOS CONDUTORES......................................................27
4.6 DISCRETIZAÇÃO TEMPORAL......................................................................27
4.6.1 PARÂMETROS DE DISCRETIZAÇÃO...........................................................29
4.7 INCLUSÃO DE PÓLOS E RESÍDUOS.............................................................30
4.7.1 MÉTODO A – USANDO CONVOLUÇÃO NO TEMPO.................................30
v
4.7.2 MÉTODO B – USANDO A TÉCNICA DE NEWMARK.................................33
5 MODELO COMPUTACIONAL...........................................................................35
5.1 INTRODUÇÃO....................................................................................................35
5.2 LINHA DE TRANSMISSÃO STRIPLINE........................................................36
5.2.1 DIMENSÕES......................................................................................................37
5.3 FERRAMENTAS COMPUTACIONAIS..........................................................38
5.3.1 SIMULAÇÃO COMPUTACIONAL.................................................................39
5.3.1.1 Tipo de Problema.............................................................................................40
5.3.1.2 Geometria.........................................................................................................40
5.3.1.3 Materiais...........................................................................................................44
5.3.1.4 Condições de Contorno....................................................................................46
5.3.1.5 Geração da Malha Tetraédrica e do Cálculo Estrutural....................................51
5.3.1.6 Processamento e Pós-processamento Matemático no MATLAB......................52
5.4 OBTENÇÃO DOS PARÂMETROS A PARTIR DA SIMULAÇÃO.............53
5.4.1 DA ESTRUTURA DE SIMULAÇÃO DO BIPLO............................................53
5.4.2 DA ESTRUTURA DE SIMULAÇÃO DA CHAVE..........................................55
6 RESULTADOS........................................................................................................59
6.1 COMPARAÇÃO DOS MÉTODOS ATRAVÉS DA SIMULAÇÃO DO
BIPOLO...............................................................................................................59
6.2 COMPARAÇÃO DOS MÉTODOS ATRAVÉS DA SIMULAÇÃO DA
CHAVE................................................................................................................61
7 CONCLUSÃO..........................................................................................................66
7.1 TRABALHOS FUTUROS .................................................................................66
REFERÊNCIAS BIBLIOGRÁFICAS......................................................................67
ANEXO A: Especificações do Transistor MGF4400A Series da Mitsubishi.........69
ANEXO B: Algoritmo da Família de Newmark.......................................................72
ANEXO C: Diferenças Finitas Centradas.................................................................75
ANEXO D: Memorial de Trabalho e Cronograma..................................................78
vi
LISTA DE FIGURAS
FIGURA 2.1 – Analogia com a luz.................................................................................3
FIGURA 2.2 – Dispositivo de duas portas ou quadripólo...............................................4
FIGURA 2.3 – Quadripólo caracterizado por uma matriz de parâmetros-S...................6
FIGURA 3.1 – Magnitudes medida e da aproximação da Yin......................................15
FIGURA 3.2 – Fases medida e da aproximação da Yin................................................15
FIGURA 3.3 – Magnitudes medida e da aproximação da Yout....................................16
FIGURA 3.4 – Fases medida e da aproximação da Yout..............................................16
FIGURA 4.1 – Uma célula tetraédrica..........................................................................20
FIGURA 5.1 – Secção transversal da linha de transmissão stripline............................37
FIGURA 5.2 – Tela principal do software GID............................................................39
FIGURA 5.3 – Comandos de definição do tipo de problema.......................................40
FIGURA 5.4 – Comandos para a criação das linhas da estrutura.................................41
FIGURA 5.5 – Linhas de contorno de um quarto da estrutura simulada......................41
FIGURA 5.6 – Comandos de criação das superfícies da estrutura...............................42
FIGURA 5.7 – Estrutura com as superfícies definidas.................................................42
FIGURA 5.8 – Comandos de criação do volume da estrutura......................................43
FIGURA 5.9 – Estrutura com o volume definido.........................................................43
FIGURA 5.10 – Comandos para inclusão de novos materiais......................................44
FIGURA 5.11 – Preenchimento do volume com teflon................................................44
FIGURA 5.12 – Características do teflon......................................................................45
FIGURA 5.13 – Definição das superfícies condutoras com PEC.................................45
FIGURA 5.14 – Características do PEC.......................................................................46
FIGURA 5.15 – Local da fonte de excitação e da resistência de casamento................47
FIGURA 5.16 – Características da resistência da fonte................................................47
FIGURA 5.17 – Pulso de corrente de excitação do circuito.........................................48
FIGURA 5.18 – Parâmetros gerais da simulação..........................................................49
FIGURA 5.19 – Características de modelamento a) da chave e, b) do bipólo..............50
FIGURA 5.20 – Estrutura de modelamento a) da chave e, b) do bipólo......................51
FIGURA 5.21 – Malha tetraédrica da estrutura de simulação da chave.......................52
vii
FIGURA 5.22 – Amostras de tensão a) na chave e, b) no bipólo.................................53
FIGURA 5.23 – Ondas incidente e refletida segundo a fonte e o bipólo......................54
FIGURA 5.24 – Ondas incidente a1 e a2 segundo a posição da chave.........................56
FIGURA 5.25 – Circuito equivalente da chave modelada............................................57
FIGURA 5.26 – Detalhe das indutâncias e capacitâncias parasitas..............................58
FIGURA 6.1 – Magnitude do parâmetro S11 medido e dos métodos A e B..................59
FIGURA 6.2 – Fase do parâmetro S11 medido e dos métodos A e B............................60
FIGURA 6.3 – Magnitude do parâmetro S22 medido e dos métodos A e B..................60
FIGURA 6.4 – Fase do parâmetro S22 medido e dos métodos A e B............................61
FIGURA 6.5 – Coeficientes de Transmissão S21 da chave com LC =10nH, usando
o método A....................................................................................................................62
FIGURA 6.6 – Coeficientes de Transmissão S21 da chave com LC =10nH, usando
o método B....................................................................................................................62
FIGURA 6.7 – Coeficientes de Transmissão S21 da chave com LC =1nH, usando
o método A....................................................................................................................63
FIGURA 6.8 – Coeficientes de Transmissão S21 da chave com LC =1nH, usando
o método B....................................................................................................................63
FIGURA 6.9 – Coeficientes de Transmissão S21 da chave com LC =0,1nH, usando
o método A....................................................................................................................64
FIGURA 6.10 – Coeficientes de Transmissão S21 da chave com LC =0,1nH, usando
o método B....................................................................................................................64
viii
LISTA DE TABELAS
TABELA 3.1 – Erros das aproximações das magnitudes para Yin e Yout....................17
TABELA 3.2 – Erros das aproximações das fases para Yin e Yout..............................17
TABELA 3.3 – Pólos das aproximações das magnitudes das Yin e Yout.....................18
TABELA 3.4 – Resíduos das aproximações das magnitudes das Yin e Yout................18
TABELA B.1 – Métodos da família de Newmark........................................................74
ix
LISTA DE SÍMBOLOS E SIGLAS
e Vetor campo elétrico
h Vetor campo magnético
d Vetor densidade de fluxo elétrico
b Vetor densidade de fluxo magnético
j Vetor densidade de corrente elétrica
ρ Densidade de carga
σ Condutividade elétrica
mn,δ Delta de Kronecker
µ Permeabilidade magnética
ε Permissividade elétrica
c Velocidade da luz no vácuo
rε Permissividade elétrica relativa
ν Relutividade magnética
fmax Freqüência máxima de operação
fc Freqüência de corte
oµ Permeabilidade elétrica no vácuo
oε Permissividade elétrica no vácuo
s Passo do tempo normalizado
γ Constante de propagação da onda
i Corrente de excitação
v Vetor de tensão desconhecido
T Duração do pulso de excitação
t∆ Passo da simulação
n Número do passo do tempo
Rs Resistência superficial
RL Resistência linear
x
Zo Impedância característica
ZS Impedância da fonte de excitação
ZL Impedância da carga
ZC Impedância da chave eletrônica
ZC Impedância da chave eletrônica
LC Indutância da chave eletrônica
LP Indutância parasita da chave eletrônica
CP Capacitância parasita da chave eletrônica
pR Vetor de posicionamento vértice p
pR Vetor de posicionamento vértice q
pqL Comprimento vetorial da aresta primária
mL Comprimento médio das arestas
L Matriz de comprimentos de arestas primárias
SS Área vetorial de face primária
S Matriz de áreas primárias
pqS Área vetorial de face secundária
S’ Matriz de áreas secundárias
C Matriz de capacitância de uma célula
G Matriz de condutância de uma célula
K Matriz de relutância de uma célula
ppm Partes por milhão
S Condutância
PEC Condutor elétrico perfeito (Perfect electric conductor)
PMC Condutor magnético perfeito (Perfect magnetic conductor)
FETD Elementos Finitos no Domínio do Tempo (Finite Element Time Domain)
FDTD Diferenças Finitas no Domínio do Tempo (Difference Element Time
Domain)
TEFLON Politetrafluoroetileno (PTFE)
MESFET Transistor de Efeito de Campo Metal - Semicondutor
xi
RESUMO
Este projeto apresenta um estudo sobre diferentes formas de inclusão de bipólos
discretos no método de análise eletromagnética FETD (Finite Element Time Domain).
São apresentadas duas formas de implementação numérica: uma é feita a partir de uma
convolução no domínio do tempo [1] e a outra é feita na forma de inclusão de frações
parciais no domínio da freqüência usando a técnica de estabilidade de Newmark [2].
Foram criados modelos para a linha de transmissão stripline, usadas nas simulações,
usando um software de modelagem computacional GID (Geometry and Data).
Inicialmente foi verificada a precisão dos métodos estudados através dos parâmetros
de admitância de entrada e de saída de um transistor MESFET comercial caracterizado
por pólos e resíduos obtidos a partir de uma aproximação [1] por funções parciais da
resposta ao impulso. Depois foi modelada uma chave eletrônica para verificar a
estabilidade das duas formas de implementação ao operarem em sistemas com baixa
indutância. As duas formas apresentadas tiveram uma boa precisão, porém, a primeira
forma apresentou instabilidade em condições de operação com baixa indutância.
Palavras-chave: bipólos discretos, FETD, técnica de estabilidade de Newmark.
xii
ABSTRACT
This project shows a study about differents kinds of inclusion of discreet dipoles in the
method for electromagnetic analysis FETD (Finite Element Time Domain). It shows
two kinds of numerical implementions: one is made through a convolution in the time
domain [1] and another one is made in the form of inclusion of partial fractions in the
frequency domain using the Newmark stability technique [2]. It has been created a
model for a stripline transmission line to be used in simulations using a computational
modeling software GID (Geometry and Data). Initially it was verified the accuracy of
the studied methods using the input and output admittance parameters of a commercial
MESFET transistor which are characterized with poles and residues gotten from an
approach [1] for partial functions of the impulse response. Later an electronic switch
was modeled to verify the stability of the two forms of implementation when operating
in systems with low inductance. The two presented forms have shown good precision,
however, the first form presented instabilities in conditions of operation with low
inductance.
Keywords: discrete dipoles, FETD, Newmark stability technique.
1 INTRODUÇÃO
Neste trabalho são comparados os resultados entre duas formas de implementar
bipólos discretos operando em altas freqüências no método numérico de análise
eletromagnética FETD (Finite Element Time Domain).
A primeira forma de inclusão foi feita por uma convolução da transformada de
Laplace [1] das aproximações por frações parciais das admitâncias dos bipólos
analisados com as tensões amostradas.
A segunda forma foi feita pela inclusão dessas frações parciais utilizando a
técnica de estabilidade de Newmark [2].
As comparações foram realizadas de duas maneiras: usando as impedâncias de
entrada e de saída de um transistor MESFET comercial e usando diferentes estados de
operação de uma chave eletrônica.
O transistor utilizado é caracterizado por seus parâmetros de espalhamento ou
simplesmente parâmetros S, medidos no domínio da freqüência.
As simulações realizadas neste trabalho foram feitas no domínio do tempo e por
isso foi necessário obter parâmetros neste domínio que descrevessem o transistor
usado, para isso foram extraídos os parâmetros de admitância ou parâmetros Y a partir
dos parâmetros S [1] que caracterizam tal transistor.
Os parâmetros S foram desenvolvidos devido a necessidade de técnicas de
análise de redes elétricas que operam em altas freqüências e teve início em 1967,
quando o Hewlett-Packard Journal publicou um artigo de Dick Anderson, com as
colaborações de Lee Smith e Jeff Gruszynski. Nele eram definidos parâmetros de
espalhamento (scattering parameters) [3], que apresentam como principal vantagem o
fato de não necessitarem conexões em curto-circuito ou em circuito aberto para serem
determinados, possibilitando medidas precisas em altas freqüências.
Na construção das estruturas onde operaram os bipólos foram simuladas linhas
de transmissão do tipo stripline utilizando o software GID (Geometry and Data). Este
2
software usou como rotina de cálculo numérico o método FETD para a análise
eletromagnética das estruturas construídas.
O diferencial deste método numérico de análise eletromagnética, em contraste
com os métodos de análise de circuitos elétricos, é que ele considera a distribuição de
campo eletromagnético que gera elementos parasitas, não contemplados na análise de
circuitos elétricos.
Os dados fornecidos pelo GID foram processados pelo software MATLAB
(através do método FETD) onde foram gerados os gráficos usados nas comparações
entre os métodos A e B.
O capítulo dois apresenta uma descrição das características e formulações dos
parâmetros Y e S e a relação entre eles é obtida, pois ela será necessária nos capítulos
seguintes.
No capítulo três descreve-se a metodologia matemática usada na aproximação
da resposta em freqüência da admitância por uma série de frações parciais [1] e
verifica-se a melhor aproximação obtida.
O capítulo quatro descreve o método numérico de simulação eletromagnética
FETD e as formas de inclusão de pólos e resíduos que caracterizam os bipólos
simulados.
No capítulo cinco é descrita a metodologia de construção, as etapas de
processamento das estruturas simuladas e a obtenção dos parâmetros necessários para
as análises a partir das simulações.
O capítulo seis contém os resultados. Primeiramente é verificada a precisão dos
resultados obtidos das simulações realizadas com as duas diferentes formas de
implementação dos bipolos através de uma comparação com os dados retirados de um
catálogo comercial (Anexo A). Posteriormente é analisada a estabilidade das formas
de implementação fazendo a comparação entre os resultados da simulação de uma
chave com os resultados de um circuito equivalente que modela a estrutura simulada
da chave.
No capítulo sete são apresentadas as conclusões gerais do projeto bem como
perspectivas de trabalhos futuros.
3
2 PARÂMETROS DE ADMITÂNCIA E DE ESPALHAMENTO
2.1 INTRODUÇÃO
A caracterização de sistemas elétricos em altas freqüências envolve a medida de
ondas de sinal incidente, refletida e transmitida, em analogia com a onda da luz, como
mostra a FIGURA 2.1.
FIGURA 2.1 – Analogia com a luz.
Em altas freqüências o comprimento de onda do sinal é comparável ou menor
que o comprimento dos condutores num circuito elétrico, surge então a importância de
se estudar referências que caracterizem o componente de acordo com as ondas
existentes no circuito.
Neste capitulo é realizado uma abordagem sobre parâmetros de admitância e de
espalhamento, suas formulações e algumas de suas características.
4
2.2 PARÂMETROS Y: ADMITÂNCIAS
A FIGURA 2.2 mostra um dispositivo de duas portas, ou quadripólo, com as
tensões aplicadas e correntes que entram nas portas 1 e 2.
FIGURA 2.2 – Dispositivo de duas portas ou quadripólo.
Este dispositivo pode ser caracterizado pela matriz de admitâncias:
=
2221
1211
YY
YYY (2.1)
Pode-se excitar o quadripólo por uma tensão V1 na porta 1 e uma tensão V2 na
porta 2, tendo as duas correntes I1 e I2 como respostas da rede. Neste caso V1 e V2 são
as variáveis independentes e I1 e I2 são as variáveis dependentes, este funcionamento é
descrito pelas equações:
2121111 VYVYI += (2.2)
2221212 VYVYI += (2.3)
5
ou na forma matricial:
=
2
1
2221
1211
2
1
V
V
YY
YY
I
I (2.4)
sendo as admitâncias definidas por:
01
111
2 =
=V
V
IY (2.5)
02
112
1 =
=V
V
IY (2.6)
01
221
2 =
=V
V
IY (2.7)
02
222
1 =
=V
V
IY (2.8)
As Ys são conhecidas como admitâncias em curto-circuito (V1 ou V2 = 0), ou
parâmetros em curto-circuito, ou simplesmente parâmetros Y. Os parâmetros Y11 e Y22
são admitâncias vistas a partir de uma entrada, com a outra curto-circuitada, e Y12 e Y21
são admitâncias de transferência ou relações entre a corrente de uma entrada pela
tensão na outra, nas condições apropriadas de curto-circuito.
A obtenção de curto-circuito e as medidas de correntes em altas freqüências são
muito difíceis de se fazerem, tornando a obtenção dos parâmetros Y complicada.
6
2.3 PARÂMETROS DE ESPALHAMENTO OU PARÂMETROS S
Os parâmetros de espalhamento foram desenvolvidos para facilitar a
caracterização das redes elétricas em alta freqüência. Não se pode simplesmente
conectar um voltímetro ou uma ponta de prova para corrente e obter medidas precisas
devido à impedância das ponteiras e à dificuldade de colocar as ponteiras em posições
desejadas. No mais, dispositivos ativos podem oscilar ou se destruir com as conexões
de curto-circuito e circuito aberto. Tais parâmetros foram primeiramente mencionados
no artigo S-Parameters Theory and Applications publicado em fevereiro de 1967 pelo
Hewlett-Packard Journal [3].
Caracterizam-se dispositivos de duas ou mais portas, projetados para
trabalharem em altas freqüências, utilizando a matriz de espalhamento ou de
parâmetros S [4]:
=
2221
1211
SS
SSS (2.9)
A partir desta matriz são relacionas as ondas de sinal incidentes e refletidas nas
portas dos dispositivos, para a matriz (2.9) um quadripólo. A FIGURA 2.3 mostra um
quadripólo com a representação dessas ondas.
FIGURA 2.3 – Quadripólo caracterizado por uma matriz de parâmetros-S.
7
onde:
- a1 é a amplitude da onda de sinal incidente na porta 1;
- b1 é a amplitude da onda de sinal refletida na porta 1;
- a2 é a amplitude da onda de sinal incidente na porta 2;
- b2 é a amplitude da onda de sinal refletida na porta 1;
- Zo a impedância característica das linhas de transmissão;
- ZS a impedância da fonte de excitação;
- ZL a impedância da carga e;
- VS a fonte de excitação do circuito.
Descreve-se o funcionamento do quadripólo entre seus terminais pelas
equações:
2121111 aSaSb += (2.10)
2221212 aSaSb += (2.11)
Escrevendo as equações (2.10) e (2.11) na forma matricial tem-se:
=
2
1
2221
1211
2
1
a
a
SS
SS
b
b (2.12)
8
Sendo os parâmetros de espalhamento definidos como:
01
111
2 =
=a
a
bS (2.13)
02
112
1 =
=a
a
bS (2.14)
01
221
2 =
=a
a
bS (2.15)
02
222
1 =
=a
a
bS (2.16)
Eles representam os:
- Coeficiente de reflexão de entrada: S11;
- Coeficiente de transmissão inverso (isolamento): S12;
- Coeficiente de transmissão direto (ganho): S21;
- Coeficiente de reflexão de saída: S22.
Os parâmetros S relacionam ondas incidentes e refletidas em todas as portas e a
partir destas relações podem-se obter medidas importantes como ganho, atenuação,
impedância de entrada, etc. Podendo também ser obtidos os parâmetros Y, Z e H a
partir dos parâmetros S. Todas as portas estão terminadas nas impedâncias de
funcionamento, habitualmente com 50Ω e com menos freqüência 75Ω.
9
2.3.1 ADMITÂNCIAS A PARTIR DOS PARÂMETROS S
As admitâncias de entrada inY e de saída outY do quadripólo são obtidas a partir
dos parâmetros S [4] da seguinte forma:
0
1
1
1
ZY
in
in
inΓ+
Γ−= (2.17)
0
1
1
1
ZY
out
out
outΓ+
Γ−= (2.18)
sendo:
L
L
inS
SSS
Γ−
Γ+=Γ
22
122111 1
(2.19)
S
S
outS
SSS
Γ−
Γ+=Γ
11
122122 1
(2.20)
e
0
0
ZZ
ZZ
L
L
L+
−=Γ (2.21)
0
0
ZZ
ZZ
S
S
S+
−=Γ (2.22)
onde inΓ , outΓ , SΓ e LΓ são os coeficientes de reflexão da entrada, da saída, da fonte e da
carga, respectivamente.
10
Neste trabalho considerou-se a linha de transmissão casada com a fonte e com a
carga, ou seja, Ω=== 500ZZZ LS e portanto 0=Γ=Γ LS , tendo 11Sin =Γ e 22Sout =Γ
e:
011
11 1
1
1
ZS
SYin
+
−= (2.23)
022
22 1
1
1
ZS
SYout
+
−= (2.24)
11
3 APROXIMAÇÃO DA ADMITÂNCIA POR UMA SÉRIE DE
FRAÇÕES PARCIAS [1]
3.1 INTRODUÇÃO
Neste capítulo é feita uma abordagem sobre o modelo matemático para a
obtenção da resposta em freqüência da admitância.
A partir de funções racionais, observando-se os vários pontos de freqüência
que seriam estudados, pode-se reduzir a função racional à uma série de frações
parciais, obtendo-se seus pólos e resíduos.
3.2 APROXIMAÇÃO DA RESPOSTA EM FREQÜÊNCIA DA ADMITÂNCIA
Fez-se o modelamento das inY e outY , representadas por Y neste capítulo, a partir
de uma função definida por um quociente de dois polinômios como:
N
N
N
N
sasasasa
sbsbsbsbbsY
+++++
+++++=
L
L
33
221
33
2210
1)( (3.1)
sendo fjs π2= a freqüência complexa.
12
Considerou-se um componente cujo fabricante realizou as medidas de seus
parâmetros em função de vários pontos de freqüência (ver Anexo A) logo:
Mmsasasasa
sbsbsbsbbsY
N
mNmmm
N
mNmmm
m ,,1,1
)(3
32
21
33
2210
KL
L=
+++++
+++++= (3.2)
onde o índice m representa o número de pontos de freqüência usados pelo fabricante,
que vai de 1 até M.
A partir da equação (3.2) montou-se o seguinte sistema de equações:
)()()(
)()()(
)()()(
110
2222122102
1111111101
M
N
MNMM
N
MNMM
N
N
N
N
N
N
N
N
sYsasYsasbsbbsY
sYsasYsasbsbbsY
sYsasYsasbsbbsY
−−−+++=
−−−+++=
−−−+++=
LL
MM
LL
LL
(3.3)
Representando o sistema anterior na seguinte forma matricial:
−−−
−−−
−−−
=
NN
MMMMMM
N
MMM
NN
NN
M a
b
b
ssYssYssYsss
ssYssYssYsss
ssYssYssYsss
sY
sY
sY
1
0
22
22222222
222
11211111
211
2
1
)()()(1
)()()(1
)()()(1
)(
)(
)(
LL
MMM
LL
LL
M (3.4)
13
ou, resumidamente:
Ax =B (3.5)
sendo:
−−−
−−−
−−−
=
N
MMMMMM
N
MMM
NN
NN
ssYssYssYsss
ssYssYssYsss
ssYssYssYsss
A
)()()(1
)()()(1
)()()(1
22
22222222
222
11211111
211
LL
MMM
LL
LL
(3.6)
[ ]T
NN aabbx LL 10= (3.7)
[ ]T
MsYsYsYB )()()( 21 L= (3.8)
Na primeira linha de A, estão os elementos que vão de 1 até Ns1 , podendo ser
reais ou imaginários. Logo depois deste último, estão os elementos que vão de
11 )( ssY− até NssY 11 )(− , compostos de parte real e parte imaginária.
Foram colocados nas primeiras linhas a parte real da matriz e em seguida as
partes imaginárias, ficando a matriz A como:
=
ImIm00Im0
ImIm00Im0
ReReReRe01
ReReReRe01
LL
MMM
LL
LL
MMM
LL
A (3.9)
Desta forma trabalhou-se apenas com os coeficientes reais.
14
Foram obtidos dos coeficientes as e bs da função (3.2) resolvendo-se o sistema
de equações lineares supradeterminado (3.5) utilizando-se a rotina QR do software
MATLAB na forma:
x= B / A (3.10)
Tendo calculado os coeficientes as e bs fez-se uma aproximação da função (3.2)
por uma série de frações parciais e obteve-se a função:
Gps
r
ps
rsY
N
N +−
++−
= L
1
1)( (3.11)
onde ri e pi são os resíduos e os pólos, respectivamente e G um termo constante,
equivalente a uma condutância pura, ou seja, independente da freqüência complexa s.
Numericamente, estes valores também foram obtidos através do software MATLAB,
pelo comando residue.
3.2.1 VERIFICAÇÃO DA APROXIMAÇÃO
Verificou-se a validade do método matemático exposto comparando-o com as
impedâncias de entrada e de saída obtidas a partir de uma tabela dos parâmetros S, de
um catálogo da Mitsubishi, medidos no domínio da freqüência em módulo e fase para
o transistor GaAs FET MGF4403A (Anexo A). O fabricante usou uma linha
microstrip com Ω= 500Z e fez 40 amostras em diferentes freqüências que vão de 0,5 a
20GHz. Nas FIGURAS 3.1 e 3.2 estão plotadas as magnitudes e fases,
respectivamente, das aproximações da Yin segundo a equação (3.11) com 2, 3, 4 e 5
pólos e a magnitude e fase da Yin a partir da equação (2.23) usando os valores de S11
medidos do catálogo, tendo os mesmos 40 valores de freqüência.
15
FIGURA 3.1 – Magnitudes medida e da aproximação da Yin.
FIGURA 3.2 – Fases medida e da aproximação da Yin.
16
As FIGURAS 3.3 e 3.4 mostram as aproximações da magnitude e da fase, respectivamente, para a admitância de saída Yout.
FIGURA 3.3 – Magnitudes medida e da aproximação da Yout.
FIGURA 3.4 – Fases medida e da aproximação da Yout.
17
Para saber com qual número de pólos houve a melhor aproximação calculou-se
o erro médio quadrático em cada situação [1].
O erro médio quadrático é uma medida dos valores calculados em relação aos
valores originais. Este erro é estimado tomando-se uma amostra dos valores calculados
e comparando-a com seus valores reais. As diferenças entre elas são elevadas ao
quadrado e somadas. A soma é então dividida pelo número de medidas, para obter-se a
média cuja raiz quadrada fornece uma medida característica de erro na mesma unidade
das medidas originais e é dado por:
( )∑=
−=M
m
ii aproximadomedidoM
Erro1
21 (3.12)
As TABELAS 3.1 e 3.2 mostram os valores dos erros das aproximações das
magnitudes e fases para Yin e Yout.
TABELA 3.1 – Erros das aproximações das magnitudes para Yin e Yout.
Erros com: Admitância
2 pólos 3 pólos 4 pólos 5 pólos
Yin (S) 0,017188355352 0,001254200564 0,000122312797 0,000015496329
Yout (S) 0,000681394423 0,000076130879 0,000197219103 0,000296192659
TABELA 3.2 – Erros das aproximações das fases para Yin e Yout.
Erros com: Fase
2 pólos 3 pólos 4 pólos 5 pólos
Yin (o) 3,524463 3,747531 3,512576 3,727962
Yout (o) 9,553926 9,694467 9,694666 9,688485
18
As diferenças entre os erros das aproximações feitas para as fases de Yin e Yout
foram pequenas sendo, assim, irrelevantes na escolha da aproximação usada.
Porém, os menores erros para as magnitudes de Yin e Yout foram os das
aproximações feitas com 4 e com 5 pólos. Com 4 pólos foram de 122 ppm e 196 ppm
e com 5 pólos foram de 15 ppm e 296 ppm, respectivamente.
Escolheu-se a aproximação com 4 pólos porque os erros entre as aproximações
das magnitudes das Yin e Yout tiveram uma diferença de 74 ppm, enquanto que com 5
pólos a diferença foi de 285 ppm entre eles. Manteve-se assim uma simetria entre os
erros das aproximações.
As TABELAS 3.3 e 3.4 mostram os 4 pólos e os 4 resíduos obtidos para as
aproximações usadas na caracterização da resposta em freqüência da admitância
TABELA 3.3 – Pólos das aproximações das magnitudes das Yin e Yout.
Admitância Pólos
Yin Yout
p1 [rad/ns] -56,338 + j 153,78 -162,86
p2 [rad/ns] -56,338 - j 153,78 -23,158 + j 117,31
p3 [rad/ns] -14,092 + j 67,657 -23,158 - j 117,31
p4 [rad/ns] -14,092 - j 67,657 -66,953
TABELA 3.4 – Resíduos das aproximações das magnitudes das Yin e Yout.
Admitância Resíduos
Yin Yout
r1 [S/ns] 0,11869 + j 0,8255 1,7272
r2 [S/ns] 0,11869 - j 0,8255 2,1552 + j 0,29422
r3 [S/ns] 1,226 + j 0,30889 2,1552 - j 0,29422
r4 [S/ns] 1,226 - j 0,30889 -0,80174
Obteve-se também, conforme a equação (3.11), uma condutância pura de
0,010531 [S] para Yin e 0,003034 [S] para Yout.
19
4 MÉTODO FDTD e FETD
4.1 INTRODUÇÃO
Soluções numéricas no domínio do tempo das equações de Maxwell
demonstram ser poderosas ferramentas para prever fenômenos e analisar dispositivos
que envolvem a propagação de ondas eletromagnéticas. Dentre estes, destaca-se o
método FDTD (Finite Difference Time Domain) que teve origem em 1966. Neste ano,
Kane S. Yee desenvolveu um trabalho que chamou de Numerical Solution of Initial
Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media [5]. O
trabalho de Yee marcou o início do estudo do método FDTD.
O termo FDTD foi utilizado pela primeira vez por Allen Taflove em 1980 [6].
O método FDTD é baseado num arranjo de grades ortogonais entrelaçadas para
discretizar o espaço contínuo e transformá-lo num problema numérico. O uso do
método FDTD difundiu-se rapidamente com o avanço da tecnologia de computadores
com alta velocidade de processamento e grande capacidade de armazenamento de
dados. O maior entrave, entretanto, reside no fato de que estruturas de forma arbitrária
não casam com a grade de discretização ortogonal. Diversos algoritmos surgiram com
o objetivo de amenizar tal problema. Alguns mantiveram a base de discretização
ortogonal [7]. Outros utilizaram malhas não estruturadas, porém seu uso acabou não
sendo difundido devido à ocorrência de instabilidades numéricas associadas a
determinados tipos de malhas [8],[9].
Mais recentemente, métodos baseados em elementos finitos FETD [10] com
funções de aresta têm atraído a atenção dos pesquisadores devido a sua boa precisão
em malhas não estruturadas apesar de sua formulação ser aparentemente mais
complexa que a do método FDTD.
Usou-se neste trabalho o método FETD com malha tetraédrica de discretização
espacial [2].
20
A análise por elementos finitos de um problema qualquer envolve, basicamente,
quatro etapas:
- discretização do domínio em um número finito de sub-regiões ou elementos;
- obtenção das equações que regem um elemento típico;
- conexão de todos os elementos no domínio e;
- resolução do sistema de equações obtido.
4.2 DEFINIÇÕES GEOMÉTRICAS E DISCRETIZAÇÃO ESPACIAL
Usou-se no desenvolvimento do método FETD foi utilizada uma formulação
diferente daquela que é convencionalmente adotada por Taflove [6]. Tal formulação
foi inicialmente apresentada por Artuzi [2] e é detalhada neste capítulo.
A região de meio homogêneo compreendida entre determinado número de faces
chama-se de célula, contendo uma malha elementar primária (arestas que compõem as
faces) e outra secundária (arestas interiores que partem de um ponto no interior da
célula até um ponto localizado em cada uma das suas faces). O método foi baseado em
uma célula tetraédrica, como a da FIGURA 4.1, podendo ser ampliado a qualquer
célula poliédrica convexa.
FIGURA 4.1 - Uma célula tetraédrica.
21
Sendo que os pontos 1, 2, 3 e 4 compõem a malha primária (linhas sólidas) e os
pontos 0, 1`, 2`, 3` e 4` referem-se à malha secundária (linhas pontilhadas). As
superfícies hachuradas indicam a face primária 3` (azul) e face secundária 23 (laranja).
O ponto P` é o baricentro da aresta.
Os pontos 1, 2, 3 e 4 são fornecidos por meio de softwares geradores de malha e
a reta formada por dois destes pontos é chamada de aresta primária e é identificada
pelos dois índices dos pontos que a formam, sendo seu comprimento vetorial:
qppq RRL −= (4.1)
sendo pR e qR os vetores de posicionamento dos vértices p e q. Estes vetores foram
arranjados em uma matriz de comprimentos primários como:
=
zyx
zyx
zyx
zyx
zyx
zyx
LLL
LLL
LLL
LLL
LLL
LLL
343434
242424
232323
141414
131313
121212
L (4.2)
Três destas arestas formam uma face primária, que é denominada pelo índice do
ponto que não faz parte de nenhuma destas arestas. A área vetorial de uma face
primária, direcionada para fora da célula é:
( )qrpq
S
S LLS ×−
=2
1 (4.3)
onde p, q, r e s representam permutações cíclicas dos índices 1, 2, 3 e 4. Estas áreas
são reagrupadas em uma matriz de áreas primárias:
22
=
zyx
zyx
zyx
zyx
SSS
SSS
SSS
SSS
444
333
222
111
S (4.4)
A malha secundária é formada pelo ponto 0 no interior da célula (baricentro do
poliedro) e por um ponto no baricentro de cada face 1`, 2` ,3` e 4`. A reta formada pelo
ponto 0 e pelo baricentro da face primária é denominada aresta secundária, cujo índice
é dado pelo índice do ponto da face na qual está ligada, sendo seu comprimento
vetorial definido como:
12rsqspq
s
LLLL
++= (4.5)
Estes comprimentos formam a matriz de comprimentos secundários:
=
zyx
zyx
zyx
zyx
LLL
LLL
LLL
LLL
444
333
222
111
L` (4.6)
Duas destas arestas secundárias definem uma face secundária que tem como
índice os pontos que não fazem parte desta face, sendo a área vetorial:
12qp
pq
SSS
−= (4.7)
23
Que formam a matriz de áreas secundárias:
=
zyx
zyx
zyx
zyx
zyx
zyx
SSS
SSS
SSS
SSS
SSS
SSS
343434
242424
232323
141414
131313
121212
S` (4.8)
Os comprimentos e áreas vetoriais são úteis na discretização espacial das
equações de Maxwell descrita na próxima secção.
4.3 MATRIZES DA EQUAÇÃO DE ONDA
Como no método dos elementos finitos [11], a formulação matricial revela-se
como compacta e permite a implementação direta no algoritmo computacional,
portanto, tal notação foi aqui usada em lugar da notação usual das formulações em
FDTD.
As matrizes da equação de onda são derivadas das equações de Maxwell na
forma integral:
∫ ∫∂
∂−=
L SSdb
tLde .. (4.9)
SdjSddt
LdhSSL
... ∫∫∫ +∂
∂= (4.10)
dVSddVS
.. ∫∫ = ρ (4.11)
0. =∫ SdbS
(4.12)
24
sendo e e h as intensidades de campos elétrico e magnético, d e b as densidades de
fluxo elétrico e magnético, ρ e j as densidades de carga e corrente, respectivamente.
A discretização espacial é baseada em uma aproximação por partes constantes
das quantidades eletromagnéticas, ou seja, é assumido que elas sejam constantes
dentro de uma célula e, portanto, as equações (4.11) e (4.12) são automaticamente
satisfeitas uma vez que a acumulação de carga é desconsiderada.
A aplicação da equação (4.9) em cada uma das faces triangulares primárias
resulta em:
bSt
eL∂
∂−=Σ (4.13)
onde as letras minúsculas representam as matrizes dos vetores da célula dependentes
do tempo:
e =[ex ey ez ]T (4.14)
b =[bx by bz ]T (4.15)
e as letras maiúsculas e a letra grega Σ representam matrizes independentes do tempo.
A matriz Σ informa a relação entre as faces e as arestas:
−−
−
−−
−
=Σ
001011
010101
100110
111000
(4.16)
25
A primeira linha da matriz efetua a circulação sobre a face triangular 1` da
FIGURA 4.1, a qual é definida pelos vértices 2, 3, e 4, logo, as colunas com elementos
não nulos são aquelas cujos índices p e q contém dois destes vértices.
A orientação positiva é dada pela regra da mão direita em relação às faces que
são orientadas para fora da célula em conformidade com a equação (4.9). De maneira
análoga formam-se as demais linhas da matriz.
A aplicação de (4.10) em cada uma das faces triangulares secundárias resulta
em:
ΣT L` ν b =
t∂
∂S`ε e + S`σ e + S`j (4.17)
onde ν, ε e σ são a relutividade magnética (inverso da permeabilidade magnética), a
permissividade elétrica e a condutividade elétrica, respectivamente, e:
j = [ jx jy jz ]T (4.18)
é o vetor da densidade de corrente de excitação.
As matrizes de equações (4.13) e (4.17) são combinadas, sendo organizadas em
uma equação de onda matricial:
Cdt
dv + G v + K ∫
t
dt0v = i (4.19)
onde:
i = S`j e v = -Le (4.20)
são a corrente de excitação i e v o vetor de tensão desconhecido ao longo das arestas
primárias e:
26
C = S`ε (LTL)-1LT (4.21)
G = S`σ (LTL)-1LT (4.22)
K = ΣT L`ν (STS)-1ST Σ (4.23)
são as matrizes de capacitância, condutância e relutância de cada uma das células,
respectivamente. Não sendo L e S matrizes quadradas, suas inversas não são possíveis
de se realizar. Contudo, a solução pode ser obtida através da aplicação da inversão
generalizada de Moore-Penrose .
A equação global da onda é formada pela superposição dos elementos das
matrizes C, G e K que tem todas as posições das arestas submetidas a um
procedimento de numeração global [12].
4.4 MODELAGEM DE SUPERFÍCIES CONDUTORAS
Para se fazer a modelagem das superfícies, foi utilizada a equação:
pqrs
pqS
pqrsLR
SG ,2,
3δ= (4.24)
onde RS é a resistência superficial do material e S é a área da face em contato com a
superfície condutora. Os valores resultantes desta matriz são superpostos na matriz
global G. Valores de RS na ordem de 10-6 Ω mostram ser adequados na prática para a
implementação de condutores elétricos perfeitos PEC (Perfect Electric Conductor).
27
Para o condutor magnético perfeito PMC (Perfect Magnetic Conductor) o valor
utilizado para a RS é ∞ .
4.5 MODELAGEM DE FIOS CONDUTORES
No caso dos fios, é feita a modelagem através de:
L
pq
pqrsR
LG =, (4.25)
em que LR é a resistência linear do material [1].
4.6 DISCRETIZAÇÃO TEMPORAL
A discretização temporal baseia-se na amostragem de v a intervalos regulares
t∆ , supondo que v varia linearmente dentro de cada intervalo. Usando o método de
estabilidade de Newmark [2] e reescrevendo a equação de onda global (4.19) tem-se:
C2
2
dt
d w + G
dt
d w + K w = i (4.26)
28
com:
w ∫=t
0v dt (4.27)
A discretização temporal é realizada através de diferenças finitas centradas (ver
Apêndice C) e uma média ponderada de w (Técnica de Newmark) conforme:
n
nnnnnnnn iwww
Kt
wwG
t
wwwC =
+++
∆
−+
∆
+− −+−+−+
4
2
2
2 11112
11 (4.28)
onde t∆ é o passo do tempo e o índice n representa o instante de tempo tnt ∆= . A
equação diferencial (4.28) foi reorganizada como o seguinte conjunto de equações
recursivas:
( )nnnn KwGviKt
Gt
Cu −−
∆+
∆+= −
−
2/1
12
.42
(4.29)
nnn tuvv ∆+= −+ 2/12/1 (4.30)
2/11 ++ ∆+= nnn tvww (4.31)
Como as matrizes C e G são diagonais, suas inversas são triviais.
29
4.6.1 PARÂMETROS DE DISCRETIZAÇÃO
Na divisão do espaço computacional em poliedros, neste caso em tetraedros,
são respeitados alguns critérios de dispersão e estabilidade numérica do método [13].
Para estabelecer esses critérios foi determinado o limite máximo de freqüência da
excitação da simulação. O princípio básico para o dimensionamento da célula é que a
célula na qual a estrutura foi discretizada deve ser bem menor que o menor
comprimento de onda excitado na simulação. Quanto menor for à célula maior será a
precisão da resposta. O que se deve ter em mente é que deve-se respeitar o limite de
amostragem de Nyquist, onde, mL∆= 2λ , com mL∆ sendo o comprimento médio das
arestas. Porém, adotam-se dimensões de células bem menores do que o comprimento
de onda, fazendo 10/λ≈∆ mL . Isso foi verificado a partir de resultados numéricos
obtidos segundo Karl [13] que, para se respeitar o limite de dispersão numérica, deve
ser observada a seguinte equação:
rmáx
m
f
cL
ε.10≅ (4.32)
onde, através da equação (4.32), definindo o limite máximo de freqüência de excitação
máxf , obtém-se o comprimento médio das arestasmL . Onde c é a velocidade da luz no
vácuo e rε é a permissividade elétrica relativa do material dielétrico de maior
permissividade aplicado na estrutura a ser simulada.
O valor do passo do tempo foi normalizado de acordo com [14]:
µεminL
ts
∆= (4.33)
30
onde Lmin é o comprimento da menor aresta primária. Esta fórmula está de acordo com
a definição de fator de estabilidade para métodos condicionalmente estáveis [6] Neste
trabalho o passo de tempo usado foi s = 1.
4.7 INCLUSÃO DE PÓLOS E RESÍDUOS
Os métodos de inclusão dos pólos e resíduos descritos a seguir foram realizados
com o objetivo de verificar estabilidade do método FETD usando a técnica de
Newmark através da comparação com outro método proposto [1].
A inclusão dos pólos e resíduos realizou-se de forma a não interferir
eletromagneticamente no restante da estrutura em análise. Assim, os dispositivos
puderam ser modelados através de sua admitância no domínio da freqüência.
No método de implementação A admitância foi convertida [1] de forma a
realizar uma eficiente convolução no domínio do tempo.
No método de inclusão B a equação (4.19) foi passada para o domínio da
freqüência e os pólos e resíduos inseridos usando a técnica de Newmark.
Estes métodos de inclusão permitem que dispositivos cujo comportamento seja
obtido através de medições também possam ser parte integrante da simulação.
4.7.1 MÉTODO A - USANDO CONVOLUÇÃO NO TEMPO
A modelagem seria igual à modelagem do fio condutor caso sua admitância não
apresentasse variação com a freqüência. A fim de incluir tal comportamento, a
31
admitância deve ser transformada para o domínio do tempo e a convolução da mesma
deve ser realizada com v. Assim a equação de onda (4.19) assume a forma:
ivdtKvtyGvvdt
dC
t
=+∗++ ∫0
)( (4.34)
onde:
)]([£)( 1 sYty −= (4.35)
sendo Y(s) a admitância do bipolo no domínio da freqüência complexa e £-1 representa
a transformada inversa de Laplace. Utilizando a aproximação de Y(s) na forma de uma
série de frações parciais no domínio s, sua transformada inversa de Laplace passa a
ser:
]...[)( 2121
tp
n
tptp nerererty +++= (4.36)
Para demonstrar o procedimento da convolução, apenas o primeiro termo de y(t)
será considerado, ou seja, com apenas um pólo e resíduo, valendo o princípio da
superposição para os demais termos. Utilizando a definição de convolução tem-se:
∫ −=∗=t
dtytvtvtytj0
)()()()()( ττ (4.37)
Assumindo que a discretização temporal se dá em intervalos regulares ∆t e que
a variação de v é linear dentro de um intervalo, a convolução pode ser aproximada pelo
somatório [1]:
32
−
−−−−
−∆−+
∆+
∆+= ∆
−−
−
n
tp
n
tp
nnn
tp
n jevrp
eKwGvir
p
tpeK
tG
tCu 1
11
2/111
2/1
1
121
12 11
42
m
tp
n
tp
n
tp
n jetutp
tpev
tp
etrj
∆
−+ +
∆
∆
−∆−+
∆
−∆= 1
11
21
12/1
111 )(
11
nnn tuvv ∆+= −+ 2/12/1
2/11 ++ ∆+= nnn tvww
[ ]∑∫=
∆+
∆−−−+ −∆++=
n
m
tm
tm
p
mnmnn dertmuvj0
)1(
12/111)1( ττ τ (4.38)
Na expressão (4.38) o termo )( τ−tv foi aproximado pela equação de uma reta
com tnt ∆= e tmtm ∆+≤≤∆ )1(τ .Resolvendo-se a integral tem-se:
tmpn
m
mn
tp
mn
tp
n etutp
tpev
tp
etrj
∆
=
−−−+ ∑
∆
∆
−∆−+
∆
−∆= 1
11
02
1
12/1
111 )(
11 (4.39)
Cujo somatório pode ser eliminado através da forma recursiva, ficando:
m
tp
n
tp
n
tp
n jetutp
tpev
tp
etrj
∆−+ +
∆
∆
−∆−+
∆
−∆= 1
11
21
12/1
111 )(
11 (4.40)
Juntando esta expressão ao conjunto de equações recursivas (4.29), (4.30) e
(4.31) obtém-se o novo arranjo:
(4.41)
33
Nota-se que um aumento do esforço computacional é esperado pela inclusão do
cálculo de 1+nj a cada interação, entretanto, caso o somatório presente na equação
(4.39) não fosse eliminado, o aumento do esforço computacional inviabilizaria sua
aplicação na prática [1].
4.7.2 MÉTODO B – USANDO A TÉCNICA DE NEWMARK
Passando a equação de onda (4.19) para o domínio da freqüência, tem-se:
ivKs
GsC =
++
1 (4.42)
As características do elemento a ser implementado são inseridas na forma de
frações parciais. Para demonstrar o procedimento será considerada uma modelagem
feita por apenas um pólo e um resíduo, valendo o princípio da superposição para os
demais casos. Assim, a equação (4.42) será:
ivps
rK
sGsC =
−+++
1
11 (4.43)
Considerando a fração parcial inserida na equação (4.43) como uma nova
variável e organizando-a em termos independentes tem-se:
34
vps
rj
1
1
−= →
s
jp
s
vrj 11 += (4.44)
Organizando o sistema de duas equações e duas incógnitas formado a partir das
equações (4.43) e (4.44) tem-se:
( )s
jpivrK
sGsC 11
1−=
+++
(4.45)
s
vr
s
jp
p
s11
1
1 =
−
As seguintes considerações são feitas para a discretização temporal da equação
(4.45):
svu =
s
vw =
(4.46)
s
jpk 1=
Por analogia ao método de estabilidade de Newmark [2] explicado
anteriormente e reescrevendo o sistema de equações (4.45) na forma recursiva tem-se:
35
+
∆−
∆+
∆−
∆+= +
+ 22
2
2
2 11
1
1
1
11
nn
nn
wwr
tp
tpk
tp
tpk
(4.47)
nnn tuvv ∆+= −+ 2/12/1
2/11 ++ ∆+= nnn tvww
A segunda equação do sistema de equações (4.47) representa a inclusão do
bipólo pela técnica de Newmark aplicada na segunda equação do sistema de equações
(4.45).
5 MODELO COMPUTACIONAL
5.1 INTRODUÇÃO
Para a análise dos métodos de inclusão a partir das características das Yin e Yout
do quadripólo construiu-se uma linha de transmissão do tipo stripline [4] com a
mesma impedância característica da linha usada pelo fabricante para realizar as
medidas dos parâmetros S. Esses parâmetros independem do tipo da linha de
transmissão utilizada, logo, não há problema de incompatibilidade.
Também foi modelada uma chave em diferentes estados de operação na linha de
transmissão stripline a fim de verificar as respostas dos métodos em diferentes
( ) ( )nnnnn kwrKGvirKt
Gt
Cu −+−−
+
∆+
∆+= −
−
)(42 12/1
1
1
2
36
situações de operação, principalmente naquelas que possuem uma baixa impedância
agregada.
5.2 LINHA DE TRANSMISSÃO STRIPLINE
Escolheu-se a linha de transmissão stripline pelo fato de sua simetria auxiliar a
implementação e por não apresentar dispersão já que é construída com material
homogêneo. É uma das estruturas comumente utilizadas como linha de transmissão em
freqüências de microondas. Em geral ela é leve, pequena e tem a grande vantagem na
fabricação e no seu custo. Ela consiste de uma fita condutora situada simetricamente
entre dois planos de terra, e o espaço entre os mesmos é preenchido com um dielétrico
homogêneo.
A FIGURA 5.1 mostra a secção transversal de uma stripline.
FIGURA 5.1 – Secção transversal da linha de transmissão stripline.
Um parâmetro importante verificado foi a impedância característica Z0 [14],
como:
37
+
−+
−−+= 27,6
8841ln
4
12
0w
tb
w
tb
w
tbZ
πππε
µ
π (5.1)
onde 0µµ = e 0εεε r= , b é a distância entre os planos condutores e w e t são as
larguras e a espessura da fita condutora, respectivamente. A fórmula utilizada para o
cálculo de Z0 , possui um erro menor que 0,5%, quando 10)( <− tbw .
A freqüência de corte, a qual é a máxima freqüência para operação da stripline,
obedece a seguinte equação:
+
=
4
15
πε
b
wb
f
r
c (5.2)
sendo b e w em centímetros e f em GHz.
5.2.1 DIMENSÕES
Escolheu-se as dimensões de acordo com a equação (5.1) para que as estruturas
tivessem uma impedância característica de Ω= 500Z . As dimensões usadas para as
linhas de transmissão foram w = 0,1 cm, b = 0,14 cm e t = 0. A freqüência de corte foi
verificada para garantir o comportamento esperado nas estruturas e de acordo com a
equação (5.2) fc = 49,3 GHz.
Utilizou-se o teflon como material dielétrico nas linhas. Ele possui
permissividade elétrica relativa 2,2=rε e a permeabilidade magnética relativa 1=rµ .
38
5.3 FERRAMENTAS COMPUTACIONAIS
Neste trabalho foram utilizados dois softwares. Para a fase de pré-
processamento foi usado o software GID (Geometry and Data) e para o processamento
e o pós-processamento foi utilizado o software MATLAB. O software GID fez a
estruturação física e o software MATLAB os cálculos
O software GID [15] é uma interface gráfica interativa usada para definição,
preparação e visualização de todos os dados relacionados à simulação numérica. Estes
dados incluem a definição da geometria, materiais, condições de contorno e outros
parâmetros. Apesar de ser um software de pré e pós-processamento, neste projeto ele
foi utilizado apenas na fase de pré-processamento. O GID é um recurso computacional
muito útil e versátil, pois através dele a geometria das estruturas a serem simuladas
podem ser esquematizadas com precisão. Além disto, através do GID, os materiais que
foram atribuídos aos componentes das estruturas podem ter suas características
eletromagnéticas completamente definidas, assim como as condições de contorno da
situação simulada. O MATLAB é um software matemático interativo, de alta
performance, utilizado em cálculos numéricos de cientistas, engenheiros,
pesquisadores, estudantes, etc. Ele é composto de um módulo matemático básico ao
qual podem-se agregar as mais variadas caixas de ferramentas como: estatística,
matemática financeira, matemática simbólica, otimização, etc.
Concluindo, com o software GID, usando o método FETD, esquematizou-se e
atribuíram-se materiais e condições de contorno às estruturas simuladas (fase de pré-
processamento), enquanto que o MATLAB realizou o processamento [2] dos dados
gerados no pré-processamento. Com os dados do processamento foram realizados os
cálculos necessários na obtenção dos parâmetros estudados, possibilitando a
visualização gráfica dos resultados.
39
5.3.1 SIMULAÇÃO COMPUTACIONAL
A seguir são descritos os principais passos da construção da estrutura usada
para a análise com os bipólos, sendo citadas as diferenças com a estrutura de analise da
chave quando necessário. A FIGURA 5.2 mostra a tela principal do software GID com
suas principais ferramentas.
FIGURA 5.2 – Tela principal do software GID.
As cinco primeiras fases descritas a seguir são de pré-processamento e a última
a de processamento e pós-processamento:
- definição do tipo de problema;
- esquematização da geometria da estrutura;
- atribuição dos materiais de cada parte da estrutura;
- atribuição das condições de contorno às entidades geométricas;
- geração da malha e calculo estrutural da simulação;
- processamento e pós-processamento matemático.
40
5.3.1.1 Tipo de Problema
Neste trabalho usaram-se dois tipos de problemas: a forma de implementação A
e a forma de implementação B.
A FIGURA 5.3 mostra a localização destas opções no GID a partir dos
comandos: Data→Problem type.
FIGURA 5.3 – Comandos de definição do tipo de problema.
5.3.1.2 Geometria
Consideraram-se as simetrias geométricas, elétricas e magnéticas das estruturas
para economizar substancialmente os recursos computacionais sendo utilizado apenas
um quarto da estrutura.
O GID obedece à hierarquia de identidades geométricas: pontos, linhas,
superfícies, volumes – nesta ordem, ou seja, primeiro construiu-se os pontos, em
seguida as linhas, depois as superfícies e então os volumes.
41
A construção das linhas de contorno da estrutura foi feita a partir do menu de
comandos: Geometry→Create →Line, como mostra a FIGURA 5.4. As linhas foram
formadas pela ligação entre dois pontos definidos pelas coordenadas x, y e z digitadas
na linha de comando, mostrada na FIGURA 5.2.
FIGURA 5.4 – Comandos para a criação das linhas da estrutura.
A FIGURA 5.5 mostra o inicio da construção da estrutura, sendo somente as
linhas de contorno e os pontos definidos.
FIGURA 5.5 – Linhas de contorno de um quarto da estrutura simulada.
Definiram-se as superfícies da estrutura a partir das linhas de contorno. Para
criá-las, utilizou-se os comandos: Geometry →Create →NURBSsurface→By
contour, mostrados na FIGURA 5.6.
42
FIGURA 5.6 – Comandos de criação das superfícies da estrutura.
No GID as superfícies são identificadas pela cor rosa. A estrutura com as
superfícies definidas é mostrada na FIGURA 5.7.
FIGURA 5.7 – Estrutura com as superfícies definidas.
43
Criou-se o volume da estrutura segundo os seguintes comandos:
Geometry →Create →Volume→By contour, mostrada na FIGURA 5.8.
FIGURA 5.8 – Comandos de criação do volume da estrutura.
No GID o volume é identificado pela cor azul clara. A estrutura com o volume
definido é mostrada na FIGURA 5.9.
FIGURA 5.9 – Estrutura com o volume definido.
44
5.3.1.3 Materiais
Precisaram-se criar os materiais necessários, pois inicialmente existe apenas o
AR definido como material no GID. Os materiais foram adicionados seguindo os
comandos Data → Materials e a partir do ícone New material abriu-se uma janela
como a mostrada na FIGURA 5.10. Nela foram inseridos os materiais necessários e
suas características.
FIGURA 5.10 – Comandos para inclusão de novos materiais.
A FIGURA 5.11 mostra o volume da estrutura preenchida com o material
teflon.
FIGURA 5.11 – Preenchimento do volume com teflon.
45
A FIGURA 5.12 mostra as características do teflon inseridas na ferramenta
GID.
FIGURA 5.12 – Características do teflon.
A FIGURA 5.13 mostra as superfícies condutoras definidas com o material
PEC.
FIGURA 5.13 – Definição das superfícies condutoras com PEC.
46
A FIGURA 5.14 mostra as características do PEC inseridas na ferramenta GID.
FIGURA 5.14 – Características do PEC.
5.3.1.4 Condições de Contorno
Utilizou-se um pulso de corrente como fonte de excitação dos circuitos
colocada na extremidade da stripline como mostra a FIGURA 5.15, com duração de
0,15ns.
Foi conectada uma resistência de Ω50 em paralelo com a fonte de corrente de
excitação com o objetivo de manter a rede casada.
47
FIGURA 5.15 – Local da fonte excitação e da resistência de casamento.
A FIGURA 5.16 mostra as características da resistência inserida.
FIGURA 5.16 – Características da resistência da fonte.
48
O pulso comporta-se como mostra a FIGURA 5.17.
FIGURA 5.17 – Pulso de corrente de excitação do circuito.
O pulso de corrente aplicado é descrito pela equação:
( )( )
−=
2cos
1)(
2
ππ
π
u
uusenti (5.3)
sendo ( )[ ]3/6 −= Ttu , t é o tempo e T é a duração do pulso.
Na análise feita a partir da chave, a fonte de excitação utilizada teve um tempo
de duração de 0,15 ns e conforme a equação (5.4) uma freqüência máxima de
aproximadamente 20 GHz.
49
max
3
fT = (5.4)
Na análise feita a partir dos parâmetros S a freqüência máxima de operação do
circuito também foi de 20 GHz. Esse valor de freqüência, de acordo com a equação
(5.2), está abaixo da freqüência de corte da linha de transmissão utilizada que é de 49,3
GHz.
O tempo de duração da fonte e outros parâmetros para simular os modelos
propostos foram especificados na janela Problem Data, a partir do menu de comandos
Data→Problem Data, como mostra a FIGURA 5.18.
FIGURA 5.18 – Parâmetros gerais da simulação.
Os pólos e resíduos foram inseridos a partir dos comandos Data→Materials
para o método de implementação A e para o método de implementação B a partir dos
comandos Data → Conditions. A FIGURA 5.19 a) mostra as características e onde
50
foram inseridos o pólo e o resíduo que simularam a chave e a FIGURA 5.19 b) mostra
onde foram inseridos os pólos e resíduos usados na aproximação do bipolo.
FIGURA 5.19 – Características de modelamento a) da chave e, b) do bipolo.
A FIGURA 5.20 a) mostra a estrutura usada na simulação da chave e a
FIGURA 20 b) mostra a estrutura usada na simulação do bipolo.
51
FIGURA 5.20 – Estrutura de modelamento a) da chave e, b) do bipolo.
5.3.1.5 Geração da Malha Tetraédrica e do Cálculo Estrutural
Fez-se a discretização das estruturas por uma malha tetraédrica gerada após as
atribuições dos materiais e das condições de contorno, explicadas anteriormente.
Inicialmente estabeleceram-se as entidades (linhas e superfícies) para as quais o
software deve gerar elementos (linhas e tetraedros) de malha de discretização da
estrutura. Gerou-se a malha a partir dos comandos Meshing →Generate.
A FIGURA 5.21 mostra o resultado da geração na estrutura de simulação da
chave, sendo a geração da malha na estrutura de simulação do bipolo semelhante.
52
FIGURA 5.21 – Malha tetraédrica da estrutura de simulação da chave.
Após a geração da malha de discretização da estrutura, realizou-se o cálculo
estrutural a partir dos comandos Calculate →Calculate. O programa criou três
arquivos relacionados à fonte de excitação, o tipo de problema e à estrutura simulada.
Estes arquivos são fundamentais para o processamento e pós-processamento dos
resultados no MATLAB.
5.3.1.6 Processamento e Pós-processamento Matemático no MATLAB
No software MATLAB foram feitos o processamento e o pós-processamento da
simulação digitando-se os comandos, FETD e POSTP, respectivamente, aberto no
diretório com o nome do simulação (arquivo que tem a extensão .GID).
53
5.4 OBTENÇÃO DOS PARÂMETROS A PARTIR DAS SIMULAÇÕES
Foram coletadas três amostras de tensão para a obtenção das constantes de
propagação e das amplitudes das ondas de sinal incidentes e refletidas. A FIGURA
5.22 mostra o local dessas amostras nas duas estruturas simuladas.
FIGURA 5.22 – Amostras de tensão a) na chave e, b) no bipolo.
5.4.1 DA ESTRUTURA DE SIMULAÇÃO DO BIPOLO
No cálculo da constante de propagação da onda na estrutura de simulação do
bipolo, considerou-se uma linha de transmissão com uma fonte alocada em uma de
suas extremidades e uma carga na outra extremidade; e que a direção de propagação da
onda incidente é + z e da onda refletida - z. A origem da coordenada z foi estabelecida
em um ponto no centro da linha de transmissão que coincide com a amostra de tensão
V2, como mostra a FIGURA 5.23.
54
FIGURA 5.23 – Ondas incidente e refletida segundo a fonte e o bipólo.
Gerou-se um sinal a partir da fonte (pulso Gaussiano) representado na FIGURA
5.23 pela onda incidente - a e teve-se conseqüentemente uma onda refletida – b. As
tensões Vs amostradas no domínio do tempo na origem, em um ponto a uma distância
–d e em um ponto a uma distância d são dadas por [4]:
dd beZaeZV γγ001 += − (5.5)
bZaZV 002 += (5.6)
ddbeZaeZV
γγ −+= 003 (5.7)
onde a e b são as amplitudes da ondas incidente e refletida, respectivamente, no ponto
zero e γ é a constante de propagação, calculada por:
2
13
2cosh
1
V
VVar
d
+−=γ (5.8)
sendo d a distância entre as amostras, no caso 0,2 cm.
55
Isolaram-se os termos de interesse das equações (5.5) e (5.7) tendo:
[ ]4
231
0)(1
)(d
dd
e
eeVVZa
γ
γγ
−
−= (5.9)
[ ]4
213
0)(1
)(d
dd
e
eeVVZb
γ
γγ
−
−= (5.10)
Aplicou-se transformada de Fourier nas amostras Vs obtidas no domínio do
tempo para que pudessem ser realizadas as análises no domínio da freqüência.
A partir das equações (5.9) e (5.10) obteve-se, de acordo com as equações
(2.13) e (2.16) os parâmetros S11 e S22, respectivamente, que são aqui generalizados
pos S:
a
b
eVV
eVVS
d
d
=−
−=
231
213
)(
)(γ
γ
(5.11)
Com esta equação calculou-se os parâmetros necessários para a definição das
Yin e Yout, conforme as equações (2.23) e (2.24), respectivamente. Estes resultados
foram comparados com as Yin e Yout obtidas a partir dos valores de S11 e S22 medidos
pelo fabricante.
5.4.2 DA ESTRUTURA DE SIMULAÇÃO DA CHAVE
Na comparação dos métodos realizada através do modelamento de uma chave a
análise da onda refletida não foi necessária. Neste caso o interesse esteve na onda que
56
é transmitida depois da chave e na onda que incidente nela, ou seja, o ganho. Tal
situação é mostrada na FIGURA 5.24.
FIGURA 5.24 – Ondas incidente a1 e a2 segundo a posição da chave.
As duas ondas são descritas pela equação (5.9). Dividiu-se a onda incidente a2
pela onda incidente a1 para obter o parâmetro desejado S21, conforme a equação (2.15)
tendo:
[ ][ ] 1
2
231
264
21)(
)(
a
a
eVV
eVVS
d
d
=−
−=
γ
γ
(5.12)
Com a equação (5.12) compararam-se os coeficientes de transmissão obtidos
pelos métodos A e B com uma curva de referência formada a partir de um circuito
equivalente com parâmetros concentrados que simulou as condições de modelamento
da chave, conforme a FIGURA 5.25.
57
FIGURA 5.25 – Circuito equivalente da chave modelada.
A seguinte equação descreve o coeficiente de transmissão (ganho) entre os
terminais da chave do circuito mostrado na FIGURA 5.25:
( )( ) 12
21
+++
+++
+=
PCPPC
PC
LS
LS
CsZCLLs
LLsZZ
ZZS (5.13)
onde Lp e Cp são as indutâncias e capacitâncias parasitas, respectivamente. SZ e LZ
são as impedâncias da fonte e da carga, respectivamente, e valem Ω50 . Zc e Lc são a
impedância e indutância que simularam a chave. Foram usados os seguintes valores
para Lc = 10 nH, 1 nH e 0,1 nH e para cada valor de indutância foram utilizados os
valores Zc = 0 Ω , 50 Ω , 150 Ω , 500 Ω e 10000 Ω . O primeiro valor da Zc simulou o
estado fechado da chave e o último o estado aberto. Os outros valores simularam
estados intermediários entre os dois citados.
Os valores das capacitâncias e indutâncias parasitas adicionadas no circuito
equivalente da chave foram extraídas a partir de ajustes para fazer coincidir a curva do
coeficiente de transmissão calculada pelo circuito equivalente com a curva obtida
através da simulação.
O valor da capacitância parasita foi ajustado considerando-se o caso simulado
com Zc = 10000Ω e Lc = 1 nH, fornecendo um valor de Cp = 58 fF. E o valor da
indutância parasita foi ajustado considerando-se o caso simulado com Zc = 0 e
58
Lc = 1 nH, em paralelo com a capacitância parasita obtida anteriormente, fornecendo
um valor de Lp = 0,2 nH.
A FIGURA 5.26 mostra o local de implementação da chave e as indutâncias e
capacitância parasitas existentes.
FIGURA 5.26 – Detalhe das as indutâncias e capacitâncias parasitas.
59
6 RESULTADOS
6.1 COMPARAÇÃO DOS MÉTODOS ATRAVÉS DA SIMULAÇÃO DO
BIPOLO
As FIGURAS 6.1 à 6.4 mostram os resultados obtidos da comparação dos
métodos A e B através da simulação do bipólo, sendo a legenda medido usada para os
coeficientes calculados a partir dos dados do catálogo do fabricante do transistor
(Anexo A).
FIGURA 6.1 – Magnitude do parâmetro S11 medido e dos métodos A e B.
60
FIGURA 6.2 – Fase do parâmetro S11 medido e dos métodos A e B.
FIGURA 6.3 – Magnitude do parâmetro S22 medido e dos métodos A e B.
61
FIGURA 6.4 – Fase do parâmetro S22 medido e dos métodos A e B.
As simulações realizadas mostram que os dois métodos utilizados tiveram uma
boa precisão, aproximando-se dos valores medidos pelo fabricante.
Nas FIGURAS 6.1 à 6.4 é observado que em baixas freqüências ocorre uma
distorção a qual é proveniente de uma imprecisão numérica no cálculo do coeficiente
de propagação γ através da equação (5.8). Nas FIGURAS 6.2 e 6.3 a simulação foi
feita sem levar em conta os elementos parasitas da estrutura que podem ser os
causadores do aumento do erro com a freqüência.
6.2 COMPARAÇÃO DOS MÉTODOS ATRAVÉS DA SIMULAÇÃO DA
CHAVE
As FIGURAs 6.5 à 6.10 mostram os resultados da comparação entre os métodos
A e B através da simulação da chave eletrônica, sendo a legenda circuito usada para os
dados obtidos a partir da equação (5.13).
62
FIGURA 6.5 – Coeficientes de Transmissão S21 da chave com LC =10nH, usando o método A.
FIGURA 6.6 – Coeficientes de Transmissão S21 da chave com LC =10nH, usando o método B.
63
FIGURA 6.7 – Coeficientes de Transmissão S21 da chave com LC =1nH, usando o método A.
FIGURA 6.8 – Coeficientes de Transmissão S21 da chave com LC =1nH, usando o método B.
64
FIGURA 6.9 – Coeficientes de Transmissão S21 da chave com LC =0,1nH, usando o método A.
FIGURA 6.10 – Coeficientes de Transmissão S21 da chave com LC =0,1nH, usando o método B.
65
Nas FIGURAS 6.5, 6.7 e 6.9 foi usado um valor de 1m Ω ao invés de zero,
porque a solução pelo método A não convergiu para este valor. Também na FIGURA
6.9 houve uma instabilidade deste método para o valor de 10k Ω , não sendo possível a
obtenção da curva.
Percebe-se que o método A apresenta um aumento gradual de divergência com
a diminuição dos valores das indutâncias usadas no circuito equivalente. Esse
comportamento é confirmado por outros autores e sua causa é ainda pesquisada [16].
66
7 CONCLUSÃO
Duas formas de inclusão de bipólos discretos no método numérico de análise
eletromagnética FETD foram estudados neste trabalho. O primeiro método estudado,
chamado de método A, realizou a inclusão dos bipólos no domínio do tempo [1] e o
segundo, chamado de método B, fez a inclusão no domínio da freqüência.
O método A realizou a inclusão através de uma convolução no tempo, enquanto
o método B apresentou uma forma mais simples de implementação, sendo incluído os
bipólos através de frações parciais, utilizando a técnica de estabilidade de Newmark
[14].
A precisão e a estabilidade desses dois métodos foram comprovadas através de
suas aplicações em simulações das impedâncias de entrada e de saída de um transistor
e nas simulações de diferentes estados de operação de uma chave eletrônica.
Os dois métodos tiveram uma boa precisão nas simulações que utilizaram as
impedâncias de entrada e de saída do transistor comparados com os valores medidos
pelo fabricante.
Nas simulações com a chave eletrônica o método A apresentou instabilidade à
medida que os valores de indutância e admitância resistiva diminuíram [16], enquanto
o método B manteve-se estável por utilizar a técnica de estabilidade de Newmark [14].
Com isso o trabalho alcançou os objetivos desejados, que eram verificar a
precisão e estabilidade dos dois métodos estudados.
7.1 TRABALHOS FUTUROS
Para trabalhos futuros é sugerido usar o método B para analises dos parâmetros
S21 e S12 de um quadripolo e da operação dinâmica de uma chave eletrônica.
67
REFERÊNCIAS BIBLIOGRÁFICAS [1] OLIVEIRA, C. H.: “Modelagem de Bipólos Lineares no Domínio do Tempo Baseada em Medidas no Domínio da Freqüência”, Dissertação de Mestrado, Programa de Pós-Graduação em Engenharia Elétrica (PPGEE) UFPR – 2004, pp.6-18; 28-31. [2] W. A. ARTUZI JR: “An Unconditionally Stable FETD Method Using Tetrahedral Cells”, Departamento de Engenharia Elétrica UFPR, pp. 1-3. [3] ANDERSON, D. L et al: “S-Parameter Theory and Applications”, Hewlett-Packard Journal, 1966. [4] LUDWIG, R. and BRETCHKO, P.: RF Cricuit Design: Theory and Applications. Prentice Hal, 2000, pp. 168-175. [5] YEE, K. S.: Numerical Solution of Initial Boundary Value Problems Involving Maxwell s Equations in Isotropic Media. IEEE Trans. Antennas Propagat., vol. 14, pp. 302-307, May 1966. [6] TAFLOVE, A.: Computational Electrodynamics: The Finite Difference Time-Domain Method, Artech House, Boston, 1995. [7] RAILTON, C. J. and SCHNEIDER, J. B.: "An analytical and numerical analysis of several locally conformal FDTD schemes," IEEE Transactions on Microwave Theory and Tech., vol.47, pp.56-66, Jan. 1999. [8] MADSEN, N. K.: "Divergence Preserving Discrete Surface Integral Methods for Maxwell's curl Equations using non-orthogonal Unstructured Grids " Journal of Computational Physics, vol.119, pp.34-45, 1995. [9] SANGANI, H. et al: "An Explicit Time-Domain Method Using Three-Dimensional Whitney Elemants" Microwave Theory and Opt. Technol. Lett., vol. 7, no.13, pp. 607-609. [10] S. D. GEDNE, F. S. LANSING and D. L. RASCOE: "Full wave analysis of microwave monolithic circuit devices using a generalized yee-algorithm based on unstructured grid", IEEE Trans. Microwave Theory and Tech., vol.44, pp.1393-1400, Aug. 1996. [11] SADIKU, M. N. O.: Elements of Electromagnetics, Oxford University Press, 2000, pp. 598-636.
68
[12] TAFLOVE, A. and BRODWIN, M. E.: “Numerical Solutions of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations” IEEE Trans. Microwave Theory Tech., vol. MTT-23, pp. 623-630, agosto de 1975. [13] KARL S. KUNZ; RAYMOND J. LUEBBERS, The Finite Difference Time Domain Method for Eletromagnetics, New York: CRC Press LLC, 1993. [14] W. A. ARTUZI JR: “Improving the Newmark Time Integration Scheme in Finite Element Time Domain Methods”, IEEE Transactions on Microwave Theory and Tech., pp.898-900, Nov 2005. [15] SCHMIDKE, W.: Simulação Eletromagnética: Curso de Simulação de Interferência Eletromagnética pelo Método FDTD - Esquematização da Geometria, Dos Materiais e Condições de Contorno Da Estrutura Simulada, 2004, pp. 1-30. [16] RILEY, D. J. and RILEY, N. W.: "Firs Order Models for Thin-Material Sheets and Coatings in the Finite-Element Time Domain Method" Northrop Grumman Mission Systems
69
ANEXO A
70
71
72
ANEXO B
73
ALGORITMOS DA FAMÍLIA NEWMARK
Dentre as famílias de algoritmos de integração, a família NEWMARK (1959)
possui os mais populares algoritmos de integração. Consiste basicamente em
expressar as velocidades e deslocamentos segundo aproximações por diferenças
finitas no domínio do tempo, dadas por:
( )
+−
∆+∆+= +
•••••
+ 1
2
1221
2. nnnn uu
tutuu ββ (B.1)
( )
+−∆+= +
•••••
+
•
11 1 nnnn uutuu γγ (B.2)
Os parâmetros β e γ determinam propriedades de estabilidade e precisão dos
métodos. HUGHES (1987) apresenta uma análise de estabilidade para os métodos
da família Newmark, para problemas lineares, colocando que se:
2
12 ≥≥ γβ (A.3)
os algoritmos são incondicionalmente estáveis.
Muitos métodos conhecidos podem ser gerados a partir de particularizações
das expressões (B.1) e (B.2). A TABELA B.1 apresenta alguns desses métodos e
suas características.
74
TABELA B.1 – Métodos da família Newmark
Método TIPO β γ CONDIÇÃO DE
ESTABILIDADE
ORDEM DE
PRECISÃO
Aceleração Média
(Regra Trapezoidal) Implicito 1/4 1/2 Incondicional 2
Aceleração Linear Implicito 1/6 1/2 Condicional 2
Diferença Central Explícito 0 1/2 Condicional 2
Outras informações sobre cada método em particular podem ser encontradas
em literaturas tradicionais como HUGHES (1987), COOK e outros (1989), BATHE
(1996) entre outros.
BATHE, K.J. (1996) Finite Element Procedures. Prentice-Hall International
Editions.
COOK, R.D.; MALKUS, D.S. e PLESHA, M.E. (1989) Concepts and
Applications of
Finite Element Analysis, Third Edition, Jonh Wiley & Sons.
HUGHES, T.J.R. (1987) The Finite Element Method – Linear static and
dynamic finite
element analysis. Prentice-Hall International Editions.
NEWMARK, N. M. (1959) A Method of Computation for Structural Dynamics,
ASCE
Journal of Engineering Mechanics Division, Vol. 85, 67-94.
75
ANEXO C
76
DIFERENÇAS FINITAS CENTRADAS
Problemas de base matemática, envolvendo a solução de equações
diferenciais, são o núcleo da maioria dos problemas de exatas, tecnologia e
indústria, ou seja, daquela classe de problemas associados a meios contínuos, à
evolução no tempo de corpos imersos em meios contínuos, ou associados a
evolução temporal do próprio meio contínuo. Os fenômenos associados com
contínuos resultam, quando modelados matematicamente, em equações ou
sistemas de equações diferenciais parciais, exigindo a solução destas em meios
reais com infinitas variáveis, de coordenadas infinitas. Tal modelo requer
discretização do contínuo em um sistema discreto composto por uma malha
geométrica (1D, 2D ou 3D), que reduz o número de variáveis do problema a uma
quantidade finita tornando tratável um sistema contínuo a partir de um modelo
discreto equivalente. As equações diferenciais representam assim o fenômeno
associado ao meio contínuo discretizado. A solução numérica de equações
diferenciais está relacionada a métodos de discretização e integração de variáveis
sobre um contínuo, como o Método da Diferenças Finitas e o Método dos Elementos
Finitos. Outros Métodos, além desses historicamente mais antigos, surgiram nas
últimas décadas, podendo-se referenciá-los como Diferenças Finitas Energéticas,
Elementos de Contorno, Tiras Finitas, dentre outros.
O Método das Diferenças Finitas busca efetuar diretamente a substituição das
formas derivadas da função / variável dependente, por suas formas de diferenças
finitas, obtidas pela expansão em Série de Taylor e truncamento a nível da ordem de
erro desejada:
)()(')()(2
hOxxfxfxxf +∆+=∆+
(C.1)
)()(')()(2
hOxxfxfxxf +∆−=∆−
As duas equações em (C.1) são o desenvolvimento em Taylor,
respectivamente a vante e a ré, em torno do ponto x da função f(x), assim,
considerando a soma da 1ª com o simétrico da 2ª,
)()('2)()(2
hOxxfxxfxxf +∆=∆−−∆+ (C.2)
77
resultando a forma em diferença finita central, da derivada primeira da função, a
menos de um erro de 2ª ordem,
)(2
)()()('
2hO
x
xxfxxfxf +
∆
∆−−∆+= (C.3)
analogamente, tomando as equações em (C.1), resultam, respectivamente, as
formas em diferenças finitas a vante e a ré da primeira derivada,
)(2
)()()(' hO
x
xfxxfxf +
∆
−∆+= (C.4)
)(2
)()()(' hO
x
xxfxfxf +
∆
∆+−= (C.5)
Formulações análogas podem ser desenvolvidas para formas em diferenças
finitas centradas, a vante e a ré, para derivadas de ordem superior, assim para a 2ª
derivada resulta a forma central:
2
),(),(2),(),(''
x
yxxfyxfyxxfyxf xx
∆
∆−+−∆+= (C.5)
2
),(),(2),(),(''
y
yyxfyxfyyxfyxf yy
∆
∆−+−∆+= (C.6)
A aplicação dessas formas em diferenças finitas em equações diferenciais resulta
em obter uma formulação algébrica dessas equações, que integram contribuições
pontuais do valor da função no domínio contínuo agora discretizado em uma malha
de pontos, em cujos vértices são tomados os valores da função (ver FIGURA 5.21)
http://www.uesc.br/arbelos/arquivo/sm/2002/mc.03.doc, em 12/12/2005 às 7:40 hrs.
MEMORIAL DE TRABALHO
UNIVERSIDADE FEDERAL
DO PARANÁ DATA
23/11/2005 PÁGINA
79
ASSUNTO:
Histórico de desenvolvimento do projeto de final. RELATOR:
Ismael Chiamenti
PARTICIPANTES COMO RUBRICA PARTICIPANTES COMO RUBRICA
Wilson Arnaldo Artuzi Junior Orientador Ismael Chiamenti Orientado
ASSUNTOS TRATADOS E DECISÕES TOMADAS
02/06 Confirmação de orientação para Projeto de Final de Curso para o segundo semestre de 2005; 27/06 Decisão sobre qual as configurações físicas da stripline e o quadripólo usado; 02/08 Tirado dúvidas sobre o método de extração da admitância através de frações parciais e análise do método matemático a ser usado para simular o quadripólo; 10/08 Feito o cronograma de trabalho, datas das avaliações e prazo de entrega; 17/08 Decidido refazer o programa de extração das admitâncias; 24/08 Entrega do resumo da proposta de trabalho para o Prof. Horácio; 25/08 Confirmação das extrações ; 30/08 Verificação das aproximações das admitâncias, decidido usar 110,7 graus. Decidido usar as dimensões de estruturas conforme trabalho anteriormente realizado [1]. Início da modulação computacional; 20/09 Esclarecido dúvidas sobre dimensões usadas. Definição dos materiais e suas características, tempo de duração do pulso e construção da malha estrutural; 26/09 PRIMEIRA AVALIAÇÃO PARCIAL (1/2) 27/09 Término da etapa anterior. Início da extração dos parâmetros a partir das simulações; 06/10 Verificação do calculo do coeficiente de reflexão para carga de 50 ohms. Corrigidos erros de fórmulas; 11/10 Corrigido erro na altura da stripline; 14/10 Verificação da aproximação da fase e tentativa de correção da mesma. Análise sobre o andamento do projeto. Comparação dos métodos e diferença do grau do numerador, aparecendo os termos sA+B, sendo B a constante já existente e A a capacitância. Neste método não é preciso multiplicar os parâmetros por 0,1 (gap de inserção) somente pela ordem da fração da estrutura simulada; 18/10 Início da análise do segundo método. Feita a análise teórica e tentativa de correção das curvas obtidas; 20/10 Verificação do método. Revisão do algoritmo pelo professor. Início da extração dos pólos e resíduos com o grau do numerador maior que o do denominador, tendo o objetivo de considerar as capacitâncias existentes; 25/10 Constatado a não solução do sistema proposto. Feita correção da fase. Decisão de realizar a análise de uma chave eletrônica com estudo do parâmetro S21 (coeficiente de transferência); 27/10 Sanada dúvidas sobre a extração do parâmetro S21 e definido método de análise das capacitâncias e indutâncias parasita da estrutura; 01/11 Verificado o porque da ondulação inesperada na magnitude da capacitância parasita: largura do pulso
MEMORIAL DE TRABALHO
UNIVERSIDADE FEDERAL
DO PARANÁ DATA
23/11/2005 PÁGINA
80
de corrente. Foi diminuído para aumentar a largura de banda. Realizado tentativas de obter os resultados diretamente das amostras de tensão e não dos coeficientes de propagação. Permaneceram os coeficientes; 03/11 Realizados os ajustes entre os resultados simulados e o circuito equivalente, tendo extraído os valores dos elementos parasitas; 07/11 Decidido à forma de exposição dos resultados obtidos; 08/11 Apresentação dos gráficos finais e início da elaboração do texto; 11/11 Entregado a primeira parte do texto para revisão; 18/11 Entregado a segunda parte do texto para revisão; 21/11 SEGUNDA AVALIAÇÃO PARCIAL (2/2); 23/11 Agendada a data de apresentação para o dia 05/12/2005 ás 9:00 horas; 24/11 Entregada a versão final; 25/11 FIM DO PRAZO PARA ENTREGAR VERSÃO FINAL.
