Download Estabilidad y Exactitud de una Extensión del Método

1
Estabilidad y Exactitud de una Extensión del
Método FDTD Para la Incorporación de Ferritas
Parcialmente Magnetizadas
José A. Pereda, Ana Grande, Oscar González, y Ángel Vegas
Departamento de Ingeniería de Comunicaciones (DICom), Universidad de Cantabria.
Correo electrónico: [email protected]
Abstract – En esta comunicación se presenta una extensión del método de las diferencias finitas en el dominio del
tiempo (FDTD) para incorporar ferritas parcialmente magnetizadas. El comportamiento macroscópico de las ferritas
se describe a través de un tensor permeabilidad obtenido de
forma empírica. Con el ob jetivo de estudiar las características numéricas (estabilidad y excatitud) de la formulación
FDTD resultante, se considera el problema de la propagación de ondas planas en medios que contengan ferritas magnetizadas longitudinalmente. Utilizando el método de von
Neuman para estudiar la estabilidad numérica del algoritmo,
se demuestra analíticamente que la extensión propuesta conserva la condición de estabilidad del método FDTD convencional y que, además, el algoritmo no presenta el fenómeno
de disipación numérica. Por último, también se comprueba
que la extactitud de los resultados numéricos obtenidos no
se deteriora en comparación a la obtenida para el caso de
ferritas desmagnetizadas.
Keywords– Magnetized ferrite, FDTD methods, stability
I. I
incorporar ferritas parcialmente magnetizadas. Esta nueva
formulación se basa en las mismas ideas que la presentada
previamente en [7]. Sin embargo, el método se aplica ahora al caso de la propagación de ondas planas a lo largo de
ferritas longitudinalmente magnetizadas. Esta situación se
puede resolver con una formulación FDTD en 1D, lo que
nos permite estudiar las propiedades numéricas del algoritmo de una forma analítica.
Como veremos, a través de un estudio de la estabilidad
utilizando el método de von Neumann, se mostrará que
nuestra formulación mantiene la condición de estabilidad
del método FDTD convencional. Además, también obtendremos la ecuación de dispersión numérica, demostrando
que nuestra formulación no presenta disipación numérica.
Por último, para validar la exactitud del algoritmo, se calculan los errores absolutos de la constante de propagación
de las dos ondas propias en una ferrita magnetizada, es
decir, una onda circularmente polarizada a derechas (cuyas siglas en inglés son RHCP) y una onda circularmente
polarizada a izquierdas (LHCP). Los errores obtenidos se
comparan con el caso de la ferrita desmagnetizada.
N LOS últimos quince años, la extensión del método
de las diferencias finitas en el dominio del tiempo (cuyas siglas en inglés son FDTD) para incorporar ferritas
magnetizadas ha recibido mucha atención. La mayoría de
los trabajos presentados hasta ahora se centran en ferritas
saturadas, pudiéndose clasificar en dos grandes categorías:
los métodos que describen la ferrita en el dominio de la
frecuencia a través del tensor de Polder [1]-[3]; y aquellos
que plantean el problema directamente en el dominio del
tiempo a través de la ecuación de movimiento de Gilbert
[4]-[6]. En ambos casos y, posiblemente, debido a la complejidad de los algoritmos FDTD resultantes, el estudio de
las propiedades numéricas de las distintas formulaciones
nunca han sido estudiadas con la profundidad necesaria.
En muchas aplicaciones prácticas del campo de las microondas, las ferritas polarizadas no están saturadas, por
lo que las formulaciones anteriores dejan de ser válidas. En
[7] el método FDTD se extendió al caso de ferritas parcialmente magnetizadas y se aplicó al cálculo de diagramas de
dispersión en guías de onda cargadas con ferritas. En este
caso, la ferrita se describe, en el dominio de la frecuencia,
por un tensor permeabilidad empírico [8], [9].
El objetivo de este trabajo es estudiar la estabilidad y la
exactitud de una extensión de método FDTD que permite
donde es la permitividad y [µ(ω)] es el tensor permeabilidad. Asumiendo que la ferrita está, en promedio, magnetizada según la dirección z, los elementos del tensor [µ(ω)]
distintos de cero vienen dados por [9]
3/2 M
, (2a)
µxx = µyy = µ = µ0 µd + (1 − µd )
Ms
3/2
M
µzz = µ0 µd 1 −
,
(2b)
Ms
µ ωM
µxy = −µyx = jκ = j 0
,
(2c)
ω
Este trabajo ha sido financiado por la Dirección General de Investigación del MEC a través del proyecto TEC2006-13268-C03-03/TCM.
donde M es la magnetización, Ms es la magnetización de
saturación y ωM = γµ0 M, siendo γ la razón giromagnética.
E
II. T
Nuestro punto de partida son las ecuaciones de Maxwell del rotacional, en el dominio de la frecuencia, para un
medio lleno de ferrita magnetizada:
jω[µ(ω)]H(ω)
= −∇ × E(ω),
jωE(ω)
= ∇ × H(ω),
(1a)
(1b)
2
En (2), µd es la permeabilidad relativa de la ferrita para el
caso en la que ésta se encuentre desmagnetizada, y viene
dada por [8]
µd =
1 2
+
3 3
1−
ω
M
ω
2
.
Para desarrollar un modelo FDTD para la propagación
de ondas en ferritas magnetizadas, µd se aproxima por un
valor constante en la banda frecuencial de interés, es decir
µd µ̄d = cte. Como resultado de la aproximación, µ
en (2a) también se convierte en una constante, quedando
µ µ̄ = cte.
A partir de ahora consideramos el problema 1D de ondas
planas propagándose en la dirección z. En este caso, las
ecuaciones diferenciales que describen el problema son:
∂Hx
∂t
∂Hy
µ̄
∂t
∂Ex
∂t
∂Ey
∂t
µ̄
∂Ey
+ µ0 ω M Hy ,
∂z
∂Ex
= −
− µ0 ω M Hx ,
∂z
∂Hy
= −
,
∂z
∂Hx
=
.
∂z
=
(3a)
(3b)
(3c)
δt n
H (k+ 12 )
∆t x
δt
µ̄ Hyn (k+ 12 )
∆t
δt n+ 1
Ex 2 (k)
∆t
δt n+ 1
Ey 2 (k)
∆t
(3d)
δz n
E (k+ 12 ) + µ0 ω M µt Hyn (k+ 21 ), (4a)
∆z y
−δ z n
=
E (k+ 12 )−µ0 ωM µt Hxn (k+ 21 ), (4b)
∆z x
−δ z n+ 12
=
Hy (k),
(4c)
∆z
δ z n+ 12
=
Hx (k),
(4d)
∆z
=
donde µt es el operador promedio centrado en el tiempo,
y δt y δ z son, respectivamente, los operadores diferencia
centrada en el espacio y en el tiempo, definidos como
1
µ0 ωM ∆t
.
2µ̄
Estas ecuaciones dan lugar a un algoritmo FDTD que
implica los siguiente pasos en casa iteración temporal:
n+ 1
1) Hx 2 se calcula utilizando (5a).
n+ 1
2) Hy 2 se actualiza por medio de (5b).
3) Exn+1 se obtiene utilizando (5c).
4) Finalmente, Eyn+1 se calcula gracias a (5d).
Para mostrar la aplicación práctica de este algoritmo, hemos considerado la propagación de dos ondas planas linealmente polarizadas viajando en direcciones opuestas a través
de una lámina de ferrita. Para conseguirlo, hemos colocado dos fuentes, una a cada lado del dominio de cómputo.
Cada fuente excita una onda plana linealmente polarizada
según el eje x que se propaga hacia la lámina de ferrita,
tal y como se muestra en la figura 1. Se puede apreciar
como, si miramos en la dirección de propagación, el campo
eléctrico de la onda que se propaga en el sentido de la z
positiva gira en dirección contraria a las agujas del reloj.
Tal y como era de esperar, y de acuerdo al fenómeno de
rotación de Faraday, el campo eléctrico de la onda que se
propaga en la dirección de la z negativa gira en el sentido
de las agujas del reloj. Los parámetros del medio simulado
fueron M = 400 G, εr = 10, µ̄r = 0.9.
1
δt F n (k) = F n+ 2 (k) − F n− 2 (k),
δz F n (k) = F n (k + 21 ) − F n (k − 12 ),
1 n+ 12
1
µt F n (k) =
F
(k) + F n− 2 (k) .
2
Ahora cabe reseñar que (4a) está acoplada a (4b) ya que
n+ 1
donde
A=
Una versión en el dominio del tiempo discreto de las ecuaciones anteriores se obtiene reemplazando las derivadas por
diferencias centradas y utilizando promediado temporal para Hy en (3a) y para Hx en (3b). Las ecuaciones en diferencias resultantes, en forma operacional, son:
µ̄
desacoplarse antes de programarse, quedando
n− 1
1 n+ 1
Hx 2 (k + 12 ) =
1 − A2 Hx 2 (k + 12 )
2
1+A
∆t n
+
E (k + 1) − Eyn (k)
µ̄∆z y
∆t
[E n (k + 1) − Exn (k)]
−A
µ̄∆z x
n− 1
+ 2AHy 2 (k + 12 ) ,
(5a)
n− 1
1 n+ 1
Hy 2 (k + 12 ) =
1 + A2 Hy 2 (k + 12 )
2
1−A
∆t
−
[E n (k + 1) − Exn (k)]
µ̄∆z x
∆t n
Ey (k + 1) − Eyn (k)
+A
µ̄∆z
n+ 1
−2AHx 2 (k + 12 ) ,
(5b)
∆t n+ 12
n+ 1
Hy (k+ 12 )−Hy 2(k− 12 ) , (5c)
Exn+1 (k) = Exn (k)−
∆z
∆t n+ 12
n+ 1
Hx (k+ 12 )−Hx 2(k− 12 ) , (5d)
Eyn+1 (k) = Eyn (k)+
∆z
n+ 1
los términos Hx 2 (k + 12 ) y Hy 2 (k + 21 ) aparecen en las
dos ecuaciones. Afortunadamente, estas ecuaciones pueden
III. C
N
A. Estabilidad
Para realizar el estudio de la estabilidad del esquema
FDTD para ferritas que acabamos de presentar, utilizamos
el método de von Neumann. Tomando como punto de partida las ecuaciones en diferencias obtenidas anteriormente,
la técnica de von Neumann nos permite obtener un polinomio de estabilidad S(Z) [10]. La condición de estabilidad
3
donde
S ± (Z) = (1 ∓ jA) Z 2 + 2 2ν 2z − 1 Z + 1 ± jA.
(7)
Las raíces de S + (Z) son
+
Z1,2
=
1/2
2
− 2ν 2z − 1 ± j − 2ν 2z − 1 + 1 + A2
1 − jA
,
y las raíces de S − (Z)
−
Z1,2
=
Fig. 1. Rotación de Faraday
es que todas las raíces Zi de S(Z) queden dentro, o sobre,
el círculo unidad en el dominio de la transformada Z, es
decir |Zi | ≤ 1 ∀i.
De acuerdo con [10], para obtener S(Z) se comienza reemplazando en la ecuaciones del método (5) los campos por
→F
0 y los
sus correspondientes amplitudes complejas F
operadores en diferencias por sus correspondientes autovalores:
δ z → j2 sin k̃z2∆z ,
1
1
→ Z 2 − Z− 2 ,
1
1 12
→
Z + Z− 2 .
2
δt
µt
donde k̃z ∈ [−π/∆z , π/∆z ] es el número de onda de un
modo de Fourier. Eliminando Ẽ0x y Ẽ0y del sistema homogéneo resultante, llegamos a
1
1
1
1
(Z 2 − Z − 2 )2 H̃0x
(Z 2 − Z − 2 )2 H̃0y
con
ν 2z =
1
µ̄
= −4ν 2z H̃0x + A(Z − Z −1 )H̃0y ,
= −4ν 2z H̃0y − A(Z − Z −1 )H̃0x ,
∆t
∆z
2
sin2
k̃z ∆z
2
.
(6)
Igualando el determinante de la matriz de los coeficientes
a cero:
(Z 21 − Z − 12 )2 + 4ν 2z
−A(Z − Z −1 )
= 0,
S(Z) = 1
1
−1
−
2
2
A(Z − Z )
(Z 2 − Z 2 ) + 4ν z obtenemos el siguiente polinomio de estabilidad:
2
2
S(Z) = Z 2 + 4ν 2z − 2 Z + 1 + A2 Z 2 − 1 = 0.
Este polinomio se puede factorizar como
S(Z) = S + (Z)S − (Z),
1/2
2
− 2ν 2z − 1 ± j − 2ν 2z − 1 + 1 + A2
1 + jA
.
Se puede ver fácilmente que, para ν 2z ≤ 1 e independientemente del valor de A, el radicando es siempre positivo y,
además
±
Z = 1.
1,2
Por lo tanto, la formulación propuesta no introduce disipación numérica. Por otra parte, considerando la definición de ν z dada en (6) y tomando el peor caso para
sin2 (k̃z ∆z /2), es decir, k̃z ∆z = ±π, la condición ν 2z ≤ 1 se
puede expresar como:
1 ∆t
√
= s ≤ 1,
µ̄ ∆z
que inmediatamente se reconoce como el límite de estabilidad del método FDTD convencional, siendo s el factor de
estabilidad.
Para visualizar el resultado teórico que acabamos de obtener, hemos considerado una ferrita con M = 3000 G,
εr = 10 y µ̄r = 0.9. Para este medio, hemos calculado
numéricamente las raíces de S ± (Z). Los cálculos se han
realizado con ∆t = 0.89 ps y ν z variando desde 0.2 a 1.01.
La figura 2 muestra el lugar de las raíces de los polinomios S + (Z) y S − (Z). Se puede observar como todas las
raíces se sitúan sobre el círculo unidad cuando ν z ≤ 1, por
lo que el esquema FDTD correspondiente se dice que es
no disipativo y que mantiene el límite de estabilidad del
método FDTD convencional. Para ν z > 1, existe una raíz
de S + (Z) y otra de S − (Z) que salen del círculo unidad,
haciendo inestable el algoritmo FDTD.
B. Dispersión numérica
La ecuación de dispersión numérica puede obtenerse facilmente, sin mas que evaluar S(Z) en el círculo unidad
en el dominio de la transformada Z e igualar el resultado
a cero. En nuestro caso, haciendo S ± (ejω∆t ) = 0 en (7),
obtenemos la siguiente relación de dispersión:
± 1
1
2 k̃z ∆z
2 ω∆t
sin
= µ̃±
,
eff 2 sin
2
2
∆2z
∆t
(8)
donde k̃z± son los números de onda numéricos de las ondas
propias RHCP y LHCP. En el análisis de dispersión k̃z±
pueden ser, en general, números complejos.
4
90
4000
]
1.5
S+(Z)
-1
60
Constante de Propagación, [m
120
1
150
30
0.5
180
νz = 1
0
210
330
240
300
3000
90
β
+
β
-
β
2000
1000
0
270
1.5
120
β
+
β
β
α
α
0
60
-
10
20
30
40
50
Frecuencia [GHz]
1
S-(Z)
Fig. 3. Constantes de propagación exactas.
30
150
0.5
180
0
Error Abs β
15
300
270
Fig. 2. Lugar de las raíces de S ± (Z) para ν z variando de 0.2 a 1.01.
Error Absoluto, [m
240
Error Abs β
]
330
-1
νz = 1
210
Error Abs ββ
+
Error Abs α
+
-
10
β
β
-
5
α
En (8), la expresión de la permeabilidad numérica efectiva para las ondas propias viene dada por:
µ̃±
eff
κ̃ =
0
10
= µ̄ ± κ̃,
con
2
∆t
±
±
j2
∆z
ω∆t ∆z
sin−1
µ̄±
sin
,
eff
2
∆t
30
40
50
de las constantes de propagación, utilizando para ello la
relación:
e−j k̃z ∆z =
±
20
Frecuencia, [GHz]
Fig. 4. Error absoluto en el cálculo de la constante de propagación.
µ0 ω M
t .
tan ω∆
2
De (8) también obtenemos la siguiente expresión para la
constante de propagación numérica:
j k̃z± ≡ α̃± + j β̃ =
-
0
(9)
donde α̃± y β̃ son, respectivamente, las constantes de fase
y atenuación de las ondas propias.
Con el objeto de validar (8) hemos realizado simulaciones FDTD y hemos calculado numéricamente k̃z± para una
onda viajando en una ferrita con M = 3000 G, εr = 10 y
µ̄r = 0.9. La discretización espacial utilizada fue ∆z = 90
µm y el factor de estabilidad s = 0.99. Durante la simulación FDTD, se almacenaron las secuencias temporales de
las componentes de campo eléctrico Ex y Ey en dos nodos
espaciales consecutivos de la malla FDTD. Posteriormente, en la etapa de postprocesado, se obtuvieron los valores
DFT [Ex (k)] ± jDFT [Ey (k)]
,
DFT [Ex (k + 1)] ± jDFT [Ey (k + 1)]
donde DFT representa la transformada discreta de Fourier.
La figura 3 muestra los valores físicos exactos de α± y β ±
en función de la frecuencia. Además, como referencia, también se ha dibujado el valor de la constante de fase β para el
caso en el que la ferrita está desmagnetizada. Por último,
la figura 4 recoge los errores absolutos de las constantes
de propagación cuando éstas se calculan mediante FDTD.
Estos errores concuerdan con los calculados analíticamente
a partir de la expresión (9).
IV. C
En este trabajo se ha introducido una extensión del
método FDTD convencional que permite la simulación de
5
la propagación de ondas planas en medios que contienen ferritas parcialmente magnetizadas. Se ha demostrado que,
al menos para la propagación 1D, y para ferritas polarizadas longitudinalmente, la nueva formulación FDTD conserva las mismas propiedades numéricas del método FDTD
convencional: la misma condición de estabilidad, no presenta el fenómeno de disipación numérica y mantiene, en
esencia, el mismo grado de exactitud.
R
[1]
J. W. Schuster and R. J. Luebbers, “Finite difference time domain analysis of arbitrary biased magnetized ferrites,” Radio
Sci., vol. 31, pp. 923-930, Jul-Aug 1996.
[2] H. Kruger, H. Spachmann, T. Weiland, “Time domain modeling of gyromagnetic materials using the finite integration technique,” IEEE Trans. Magnetics, vol.37, pp. 3269-3272, Sep.
2001.
[3] W.K. Gwarek, A. Moryc, “An alternative approach to FD-TD
analysis of magnetized ferrites,” IEEE Microwave and Wireless
Components Lett.,vol.14, pp. 331- 333, July 2004.
[4] A. Reineix, T. Monediere, F. Jecko, “Ferrite analysis using the
finite-difference time-domain (FDTD) method” Microwave and
Optical Tech. Lett., vol. 5, no. 13, pp. 685-686, Dec. 1992,
[5] J. A. Pereda, L. A. Vielva, A. Vegas, and A. Prieto, “A treatment of magnetized ferrites using the FDTD method,” IEEE
Microwave Guided Wave Lett., vol. 3, pp. 136-138, May 1993.
[6] M. Okoniewski and E. Okoniewska, “FDTD analysis of magnetized ferrites: a more efficient approach,” IEEE Microwave
Guided Wave Lett., vol. 4, pp. 169-171, June 1994.
[7] J. A. Pereda, L. A. Vielva, A. Vegas and A. Prieto,“An extended FDTD method for the treatment of partially magnetized
ferrites,” IEEE Trans. Magnetics, vol. 31, pp. 1666-1669, May
1995.
[8] E. Schlomann, “Microwave behavior of partially magnetized ferrites,” J. Applied Physics, vol. 41, pp. 204-214, Jan. 1970.
[9] J. J. Green and F. Sandy, “Microwave characterization of partially magnetized ferrites,” IEEE Trans. Microwave Theory Tech.,
vol. 22, pp. 641-645, June 1974.
[10] J. A. Pereda, L. A. Vielva, A. Vegas and A. Prieto,“Analyzing
the stability of the FDTD technique by combining the von Neumann method with the Routh-Hurwitz criterion,” IEEE Trans.
Microwave Theory Tech., vol. 49, no. 2, pp. 377-381, Feb. 2001.