Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

4

1. Propagación en el espacio libre

En este capítulo hacemos una presentación sencilla del MDFDT discutiendo la

formulación para simular la propagación del campo electromagnético en un medio

homogéneo. Ciertamente, la propagación en un medio homogéneo o en el vacío no es un

problema que implique un gran interés físico. Sin embargo, el objetivo es la ilustración

del algoritmo.

1.1 Las ecuaciones de Maxwell en forma discreta

EL MDFDT analiza el problema de la propagación electromagnética en pequeñas

particiones espaciales. En estas celdas los campos eléctricos y magnéticos están

alternadamente distribuidos. Las celdas o nodos también estas intercaladas en el tiempo.

Para resolver este conjunto de ecuaciones, se desarrolla un proceso iterativo en el tiempo.

Comenzamos con las ecuaciones de Maxwell en el sistema (cgs).5

JxDxHct

tc

t π4),(1),( +∂

∂=×∇ (1.1-1)

t

tc

t∂

∂−=×∇

),(1),( xBxE (1.1-2)

πρ4),( =⋅∇ txD (1.1-3)

0),( =⋅∇ txB (1.1-4)

Donde B corresponde a la inducción magnética, E corresponde al campo

eléctrico, H corresponde al campo magnético, D es el desplazamiento eléctrico, ρ es la

densidad de carga, J es la densidad de corriente eléctrica y c corresponde a la velocidad

de la luz en el vacío. Debemos hacer notar que para medios isotrópicos y lineales las

relaciones constitutivas se escriben:

ED ε=

HB μ=

EJ .σ=

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

5

Siendo ε la constante dieléctrica del material, μ la permeabilidad magnética del material y

σ la conductividad.

Para medios dispersivos

)()()( ωωεω ED =

)()()( ωωμω HB =

)()()( ωωσω EJ =

Para ondas que se propagan en medios homogéneos y en dieléctricos y en donde

no existen fuentes de carga ni de corrientes (ρ = J = 0), las ecuaciones (1.1-1) y (1.1-2)

toman la forma:

tt

ct

∂∂

=×∇),(),( xExH ε

(1.2-1)

t

tc

t∂

∂−=×∇

),(),( xHxE μ (1.2-2)

Si escogemos la dirección x para la polarización del campo eléctrico y la

dirección z para la dirección de la propagación, entonces el campo magnético tiene la

polarización y. Las ecuaciones rotacionales vectoriales (1.2-1) y (1.2-2) son ahora las

ecuaciones escalares,

),(),( tzHz

ctzEt yx ∂

∂−=

∂∂

ε (1.3-1)

),(),( tzEz

ctzHt xy ∂

∂−=

∂∂

μ (1.3-2)

La primera de estas ecuaciones es también llamada de Maxwell-Ampere y la

segunda se conoce como la ecuación de Faraday. Estas ecuaciones son ecuaciones

puntuales, es decir, que son válidas para cada valor continuo de las coordenadas espacial

y temporal. Para lograr esto tenemos que trabajar en un espacio discreto en lugar de

hacerlo en un espacio continuo.

Tomemos la definición de derivada en la variante de diferencias centrales un

punto zo,6

z

zzfzzfdz

zdf oo

oz ΔΔ−−Δ+

=→Δ

)2/()2/(lim

)( 0 (1.3-3)

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

6

f(z - z )

z0

f(z )0

z 0

0 2

z

f(z + z )0 2

f (z )

(z - z )20 0(z + z )

2

f '(z )0

Esta ecuación es valida para todos los puntos del espacio z. En esta ecuación, para un

valor finito de Δz, la ecuación diferencial se transforma en una ecuación de diferencias

finitas. La derivada esta definida por el valor de la función en puntos discretos contiguos

a zo.

Figura 1.1. Representación gráfica de derivada en diferencias centrales.

Es necesario discretizar las ecuaciones de Maxwell. La estrategia consiste en

definir un punto “z” por medio de un valor discreto “z = kΔz” y de igual forma la

coordenada temporal “t” toma valores discretos de la forma “t = n Δt”. Para ilustrar estas

ideas obsérvese la figura en donde los valores accesibles de z y t están determinados por

los saltos de rana discretos que se muestran en la Fig. 1.2, cada salto de rana simula un

paso espacial o temporal según sea el caso.

Figura 1.2. Representación de la utilización de valores discretos de z y t.

kΔz (k+1)Δz z

nΔt (n+1)Δt t

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

7

Con estas consideraciones, la forma discreta de la derivada espacial para el campo

eléctrico (E) es:

ztnzkEtnzkE

tnzkEz

xxx Δ

ΔΔ−−ΔΔ+≅ΔΔ

∂∂ ],)2/1[(],)2/1[(

),(

Similarmente, los valores discretos de la coordenada temporal están indexados por el

entero ‘n’ por medio de la relación t = nΔt y la forma discreta temporal queda

ttnzkEtnzkE

tnzkEt

xxx Δ

Δ−Δ−Δ+Δ≅ΔΔ

∂∂ ])2/1(,[])2/1(,[

),(

Vamos a usar la convención de escribir en el superíndice el valor discreto de la

variable temporal y entre paréntesis el valor discreto de la variable espacial, es decir,

escribimos

)(),(),( kEtnzkEtzE nxxx =ΔΔ≅

Relaciones similares pueden obtenerse para el campo Hy. Al implementar las

ecuaciones en diferencias finitas para Ex(z,t) y Hy(z,t) en la ec. (1.1-7) se obtiene la forma

discreta de la ecuación de Maxwell-Ampere:

z

kHkHct

kEkE ny

ny

nx

nx

Δ

−−+−=

Δ− −+ )2/1()2/1()()( 2/12/1

ε (1.4-1)

Al igual que hemos hecho la consideración de los campos en un punto (z,t), ahora

vamos a considerar la siguiente ecuación rotacional (1.1-8) en el punto

)2/,2/( ttzz Δ+Δ+ , así la ecuación de Faraday nos queda:

z

kEkEct

kHkH nx

nx

ny

ny

Δ−+

−=Δ

+−+ +++ )()1()2/1()2/1( 2/12/11

μ (1.4-2)

Al considerar a los campos en coordenadas diferentes, es decir (z,t) y

)2/,2/( ttzz Δ+Δ+ , logramos describir una situación de evolución espacio temporal.

Estas ecuaciones pueden expresarse de la forma:

[ ])2/1()2/1()()( 2/12/1 −−+ΔΔ

−= −+ kHkHztckEkE n

yny

nx

nx ε

(1.5-1)

[ ])()1()2/1()2/1( 2/12/11 kEkEztckHkH n

xnx

ny

ny

+++ −+ΔΔ

−+=+μ

(1.5-2)

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

8

Estas ecuaciones están intercaladas en ambas coordenadas: espacio y tiempo. Esta forma

de escribir las ecuaciones de Maxwell tiene la ventaja que asemeja la forma usual en que

visualizamos la propagación de un campo electromagnético: un campo eléctrico induce

a un campo magnético, el cual induce un campo eléctrico, y así al infinito. Este

mecanismo es esquematizado en la Fig. 1.3. En esta figura, las flechas indican la forma

en como el nuevo valor del campo eléctrico proviene del campo eléctrico en ese punto un

paso temporal anterior más la contribución de los campos Hy en puntos contiguos en el

espacio un medio tiempo atrás. Este es el paradigma fundamental del MDFDT.

Hn(k+3/2)Hn(k-1/2)Hn(k-3/2) Hn(k+1/2)

En-1/2(k-2) En-1/2(k-1) En-1/2(k+2)En-1/2(k+1)

(n+2)Δt

(n+3/2)Δt

(n+1)Δt

(n+1/2)Δt

nΔt

kΔz-2Δz kΔz-Δz kΔz+2Δz kΔz+Δz

tiem

po (t

)

distancia (z)

E H

kΔz

(n-1/2)ΔtEn-1/2(k)

En+1/2(k-2) En+1/2(k-1) En+1/2(k+2)En+1/2(k+1)En+1/2(k)

Hn+1(k+3/2)Hn+1(k-1/2)Hn+1(k-3/2) Hn+1(k+1/2)

En+3/2(k-2) En+3/2(k-1) En+3/2(k+2)En+3/2(k+1)En+3/2(k)

Figura 1.3. Representación de la evolución temporal de un campo electromagnético.

En este punto, es necesario determinar las cantidades que afectan a los corchetes

en las ecs. (1.5-1) y (1.5-2). En el vacío, la función dieléctrica y la permeabilidad

magnética son iguales a 1; ε = μ =1. Una onda electromagnética necesita un mínimo de

tiempo Δt =Δz/c para propagarse dentro de la malla discreta espacial Δz. Si tomamos un

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

9

paso temporal mas grande, la onda ya habría pasado la distancia Δz y ya no estaría en

nuestra malla de simulación. Por esta razón es necesario escoger un paso temporal de la

forma: 1,7

czt Δ

≤Δ (1.6-1)

En este trabajo vamos a utilizar la convención de escoger un paso temporal

czt

2Δ

=Δ (1.7-1)

Para la propagación en el vacío, el paso espacial lo escogemos tomando 10

puntos de la longitud de onda de la luz en el vacío λ0,

100λ

=Δz . (1.8-1)

1.2 Implementación Computacional

El primer campo electromagnético que vamos a simular es un pulso gaussiano que

tiene la forma:

⎪⎪

⎩

⎪⎪

⎨

⎧

>

<=

⎟⎟⎠

⎞⎜⎜⎝

⎛ −−

0

00)(

20

21

te

ttE

t

ttσ

Esta es una función temporal que asemeja un pulso de luz y esta ilustrada en la Fig.

1.4. Tiene la característica de que su amplitud máxima es E0. Solo toma valores

diferentes de cero en la vecindad de t0, tiempo al cual esta centrada y tiene una duración

media σ t.

La implementación de las ecs. (1.5-1) y (1.5-2) en un algoritmo computacional

requiere dos vectores unidimensionales, uno para Ex y otro para Hy, para describir la

coordenada espacial. En cada tiempo ‘n’ los campos son calculados en Matlab mediante

las líneas de programa:

ex(k) = ex(k) + 0.5*( hy(k-1) – hy(k) ) (1.9-1)

hy(k) = hy(k) + 0.5*( ex(k) – ex(k+1) ) (1.9-2)

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

10

0 1 2 3 4 5 6 7 8 9 100,0

0,2

0,4

0,6

0,8

1,0

t0=5τ

σt

C

ampo

ele

ctric

o (1

/E0)

tiempo (t/τ)

Figura 1.4 Ilustración de un pulso gaussiano.

En estas relaciones, ya no existe el superíndice. El tiempo es una variable

implícita en el MDFDT. En la primera de estas ecuaciones, el nuevo valor de ex(k) a la

izquierda de la relación [para el tiempo (n+1/2)] proviene del antiguo valor ex(k) [al

tiempo(n-1/2)] que la máquina ya ha calculado y tenia en la memoria más los campos

magnéticos hy(k-1)y hy(k) [propios del tiempo n]. Para el vector hy los índices

se han redondeado de los índices k+1/2 y k-1/2 a k y k-1 para poder especificar la

posición de los campos en el arreglo computacional.

En el apéndice 1 mostramos el programa ‘simple.m’ que es la primera

implementación del algoritmo MDFDT en Matlab. Ya que el MDFDT es un proceso

iterativo, puede ser presentado en forma concisa en un diagrama de flujo.

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

11

Como se aprecia en el diagrama el programa comienza inicializando constantes y

parámetros. Los campos eléctricos y magnéticos son inicializados a cero. El cálculo

comienza obteniendo los campos eléctricos en cada nodo usando la relación (1.6-1). A

continuación, se introduce el valor de la fuente de campo eléctrico. El siguiente paso es

calcular los campos magnéticos. Finalmente, el programa grafica e itera el procedimiento

al siguiente valor temporal.

En la Fig. 1.5 mostramos el resultado del programa para 100 iteraciones

temporales. El comportamiento del pulso en el espacio y en el tiempo puede ser

visualizado en la Fig. 1.6. A medida que el pulso se propaga en el tiempo, se extiende

hacia fuera del centro de la malla espacial.

Definición de parámetros e inicialización

de los campos

Cálculo del campo eléctrico

Inserción de la fuente

Cálculo del campo magnético

Graficación

Fín

Inicio

¿Fin de pasos

temporales?

No

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

12

El valor de la longitud de onda utilizado en las simulaciones para este capítulo es

del orden de 632.8x10-9 m.

El programa detiene la simulación antes de tocar el límite de la malla de

simulación. En este método, para calcular el campo E en un punto, es necesario conocer

los valores H que rodean dicho punto. Sin embargo, en el borde de la malla de simulación

no conocemos el punto correspondiente que esta fuera de la malla. De esta forma, no es

posible calcular el campo E en los valores límites de la malla.

Si el programa continuase, es decir, si hacemos el ejercicio de modificar el

número de pasos temporales a N_t=200, veremos en la evolución temporal los valores

falsos que toma el campo en el borde de la malla. El efecto neto es una reflexión hacia el

interior de la malla. Esto es algo que es necesario evitar ya que esto no es físicamente

aceptable. Por ello introducimos las condiciones de frontera absorbentes.

Figura 1.5. Propagación de un pulso gaussiano en el espacio libre.

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

13

Figura 1.6. Variación de un pulso gaussiano en el espacio y en el tiempo

1.3 Condiciones de frontera absorbentes

Las condiciones de frontera absorbentes son necesarias para simular

adecuadamente la propagación de los campos electromagnéticos una vez que estos han

llegado a la frontera de nuestra malla de simulación. En el algoritmo del MDFDT, los

valores de los campos se determinan mediante un promedio de los campos en los puntos

vecinos. El problema consiste en que en la frontera de la malla de simulación, este

promedio no se puede dar porque no conocemos el valor del campo fuera de la malla. De

esta forma, si la simulación continua, los campos toman valores no válidos.

Físicamente, esperamos que los campos se propaguen hacia fuera de la malla de

simulación, ya que consideramos que fuera de nuestra malla no existen fuentes. La

distancia que la onda viaja en un paso temporal esta determinado por la relación (1.7-1).

distancia = c0Δt=c0(Δz/2c0)=Δz/2 (1.7-1)

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

14

Figura 1.7 Pulso gaussiano después de 140 pasos temporales incluyendo

condiciones absorbentes de frontera.

Esta relación implica que al campo le toma un par de pasos temporales viajar un

paso espacial. La forma de expresar esta condición es

)2()1( 2−= nx

nx EE

Con esta relación estamos asignando el valor del campo en la frontera en lugar de

calcularlo. Esta condición es fácil de implementar. Solamente es necesario guardar el

valor de Ex(2) un par de pasos temporales y luego asignarlos a Ex(1). Las condiciones de

frontera han sido implementadas en el programa ‘simple_ABC.m’ que se encuentra en el

apéndice 2. En la figura 1.7 mostramos la propagación de un pulso gaussiano. En esta

ocasión, el pulso pasa la frontera sin reflexión, que es lo esperado que ocurra.

1.4 Aplicaciones

Una de las ventajas más significativas del MDFDT es que es flexible para realizar

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

15

simulaciones tanto de ondas como de pulsos que varían en el tiempo. También es posible

simular a varias frecuencias, por ejemplo una onda sinusoidal o un pulso que contenga

varias frecuencias. A continuación veremos algunos ejemplos que ilustran la flexibilidad

del algoritmo.

1.4.1 Superposición de dos pulsos gaussianos

El primer ejemplo de modificación del programa es obtenido al analizar la

superposición de dos pulsos gaussianos. Uno de los pulsos esta ubicado en k=30 y el

segundo esta ubicado en k=70. En la Fig. 1.8 ilustramos la superposición de ambos

pulsos después de 90 pasos temporales, en el momento que se detiene la simulación los

máximos de cada pulso coinciden, aquí la amplitud total es la suma de las amplitudes

individuales, posterior a su encuentro cada pulso continua su trayecto.

Para lograr esta simulación, es necesario escribir una nueva función llamada

‘pulso_2gauss.m’ la cual mostramos en el apéndice 3. Asimismo, hay que modificar una

línea del programa ‘simple_ABC.m’. La línea 28 donde viene la instrucción

ex=pulso(ex,n,t,tao,t0,spread);

Debe ser cambiada por

ex=pulso_2gauss(ex,n,t,tao,t0,spread);

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

16

Figura 1.8 Superposición de dos pulsos gaussianos

Con esta variación estamos llamando a la nueva función ‘pulso_2gauss.m’ en vez

de a la función ‘pulso.m’. Cabe destacar que tanto el programa, como la función que

estemos llamando siempre deben de estar en el mismo subdirectorio.

1.4.2 Superposición de dos ondas sinusoidales

Este ejemplo requiere mas cambios sobre el programa ‘simple_ABC’, Ahora

vamos a requerir de una malla temporal mas grande, N_z=400 y las iteraciones

temporales van a ser N_t=350. La onda sinusoidal que vamos a simular es de la forma:

⎟⎠⎞

⎜⎝⎛= ttE

τπ21.0sin)(

Esta fuente va a ser ubicada en k=20 y k=380. Este cambio esta contenido en la

función ‘pulso_2sin.m’, la cual se llama dentro del programa modificando la línea 28

donde viene la instrucción

ex=pulso(ex,n,t,tao,t0,spread);

debe ser cambiada por

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

17

ex=pulso_2sin(ex,n,t,tao,t0,spread);

En la Fig. 1.9 mostramos la superposición de dos ondas después de 405 pasos

temporales, lo único que difiere de la anterior simulación es el tipo de fuente que se está

utilizando, la simulación se detiene cuando los máximos de las amplitudes en las ondas

coinciden.

Figura 1.9 Superposición de dos ondas sinusoidales.

1.4.3 Evolución de un grupo de ondas

Cuando dos ondas se encuentran presentes en el mismo medio estas pueden

combinarse constructivamente o destructivamente dependiendo de la coincidencia de los

máximos.

Para ilustrar este caso supongamos dos ondas sinusoidales de la forma:

⎟⎟⎠

⎞⎜⎜⎝

⎛= ttE

τπ2

3.0sin)(1

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

18

⎟⎟⎠

⎞⎜⎜⎝

⎛= ttE

τπ2

28.0sin)(2

El grupo de ondas formado por los campos anteriores obedece a la función:

⎟⎠⎞

⎜⎝⎛+⎟

⎠⎞

⎜⎝⎛= tttE

τπ

τπ 228.0sin23.0sin)(

Similarmente, la línea 28 donde viene la instrucción

ex=pulso(ex,n,t,tao,t0,spread);

debe ser cambiada por

ex=pulso_2group (ex,n,t,tao,t0,spread);

En este ejemplo consideramos una malla de tamaño N_z = 500 y en la fig. 2.10

mostramos la situación después de 1000 pasos temporales.

Figura 1.10 Evolución de la suma de dos ondas sinuosidades ligeramente desfasadas.

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

19

0 1 2 3

x 10-4

-1

-0.5

0

0.5

1

z (metros)

E(z

,t)

0 1 2 3

x 10-4

-1

-0.5

0

0.5

1

z (metros)

H(z

,t)

1.4.4 Dirección de propagación.

La dirección de propagación de la luz esta dictada por la regla de la mano derecha

según lo indica el vector de Poynting:

HES ×=π4c .

Dependiendo del signo que acompañe a los campos esto determinará la dirección del

pulso. Sea entonces la fuente del campo eléctrico para t=0 de la forma 22/

21

sin)(⎟⎟⎠

⎞⎜⎜⎝

⎛

Δ

−−

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛= z

zz

fx

f

ezz

zE π

22/21

sin)(⎟⎟⎠

⎞⎜⎜⎝

⎛

Δ

−−

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛= z

zz

fy

f

ezz

zH π

Cuando ambos campos son tomados positivos, estos se propagan en la dirección positiva,

como lo muestra la fig. 1.11.

Figura 1.11 Ilustración de una onda electromagnética propagándose a la derecha.

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

20

Cuando alguno de los campos es negativo, tomemos como ejemplo al campo eléctrico

negativo entonces la onda correspondiente se propaga en la dirección negativa como lo

ilustra la fig. 1.12.

Figura 1.12 Ilustración de una onda electromagnética propagándose a la izquierda.

1.4.5 Simulación de una cavidad

Como un último ejemplo, vamos a simular la oscilación dentro de una cavidad.

En este caso, el campo es cero en la frontera de la malla de simulación, lo cual

físicamente representa una cavidad de aire en medio de paredes perfectamente

conductoras. La función que describe el campo es

⎟⎟⎠

⎞⎜⎜⎝

⎛= z

zzE

fx

πsin)(

donde zf define el mayor valor asignado para la coordenada espacial.

0 1 2 3

x 10-4

-1

-0.5

0

0.5

1

z (metros)

E(z

,t)

0 1 2 3

x 10-4

-1

-0.5

0

0.5

1

z (metros)

H(z

,t)

Aplicación del Método de Diferencias Finitas en el Dominio del Tiempo

a la simulación del Campo Electromagnético.

21

Este campo esta definido solamente en el tiempo inicial t=0, a continuación, el algoritmo

se encarga de realizar la simulación de la onda que oscila en el tiempo. Esta cavidad se

ilustra en la fig. 1.13.

Figura 1.13 Ilustración de una oscilación dentro de una cavidad.