Estabilidad Y Exactitud De Una Extensión Del Método

Transcript

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 kz2∆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 ˜0x y E ˜0y del sistema homodo de Fourier. Eliminando E mogéneo resultante, llegamos a ˜ 0x (Z 2 − Z − 2 )2 H 1 1 ˜ 0y (Z 2 − Z − 2 )2 H 1 1 con ν 2z = 1 ¯ µ ˜ 0x + A(Z − Z −1 )H ˜ 0y , = −4ν 2z H ˜ 0y − A(Z − Z −1 )H ˜ 0x , = −4ν 2z H  ∆t ∆z 2 sin2  ˜z ∆z k 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 2 ω∆t z ∆z 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 =µ ¯±κ ˜, con 2 ∆t ˜± ± j2 ∆z     ω∆t ∆z sin−1 ¯ µ± sin , eff 2 ∆t 20 30 Frecuencia, [GHz] 40 50 de las constantes de propagación, utilizando para ello la relación: e−j kz ∆z = ± 10 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β j k˜z± ≡ α - 0 (9) ˜ son, respectivamente, las constantes de fase donde α ˜± y β 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.