Problemas De Vector De Poynting 653d4n

142

4.5.

Problemas

4-1. Calcule la energ´ıa almacenada en un condensador plano sometido a una diferencia de potencial entre placas V . El area de las placas es A y la distancia entre las mismas d. Exprese el resultado en funci´ on de (C , Q), (C , V ) y E. Desprecie los efectos de bordes. 4-2. Calcule la energ´ıa almacenada en un solenoide recto de secci´ on circular que est´ a recorrido por una intensidad I. Su longitud es l, su secci´ on A y el n´ umero de espiras por unidad de longitud n. Exprese el resultado en funci´ on de (L , I), (L , Φ) y B. Desprecie los efectos de bordes. Soluci´ on: Tenemos dos alternativas, calcularla en funci´ on de la densidad de energ´ıa Wm =

1 2µ0

o en funci´ on del flujo

Z

B 2 dv

V→∞

1 Wm = IΦ 2 El campo magn´ etico en el interior del solenoide que la primera alternativa nos da Wm =

14

es constante, por lo

1 V B2 2µ0

donde V = A l es el volumen del solenoide.

Por otra parte, teniendo en cuenta que la autoinducci´ on se ha definido 15 Φ como L ≡ I Wm =

1 2 1 L I2 = Φ 2 2L

4-3. Sea una l´ınea coaxial, figura 4.11, compuesta por dos conductores cil´ındricos. El interno de radio a, es hueco y el externo de radio interior b, ambos de espesor despreciable. La longitud de ambos es L ≫ a, b. Considere una peque˜ na secci´ on de la misma de longitud ∆z ≪ L y alejada de los extremos. Por los conductores circulan intensidades I, uniformemente distribuidas, iguales y en sentido contrario y la diferencia de potencial entre ellos, en dicha secci´ on, es V . Halle la energ´ıa electromagn´etica almacenada, la capacidad y la autoinducci´ on por unidad de longitud. 14 15

V´eanse los problemas del cap´ıtulo segundo. V´eanse los problemas del cap´ıtulo segundo.

143

I + V −

2a I

2b

L>>a,b Figura 4.11: 4-4. Una onda plana se propaga en el vac´ıo en la direcci´ on del eje z. El campo electrico en z = 0 es

~ =→ E

 b [e−t − e−2t ] , , t > 0  E0 x 

0,, t < 0

donde E0 = 100 µV · m−1 . Halle, para z > 0, a) Campo magn´etico. b) Vector de Poynting. c) Energ´ıa total que atraviesa a un casquete hemisf´erico, de radio a = 1 m, cuyo eje es paralelo al z. SOLUCION: (a) - Para encontrar la expresi´ on del campo electromagn´ etico en cualquier posici´ on z, basta con tener en cuenta que la onda viaja, a trav´ es del vac´ıo, en la direcci´ on positiva del eje z, por lo que t = (t− zc )z=0 y i h z z ~ b E(z, t) = E0 e −(t− c ) − e −2(t− c ) x

Haciendo uso de la relaci´ on de estructura con ~n = zb 1 1 ~ ~ B(z, t) = zb ∧ E(z, t) = E(z, t) yb c c

(b) - el vector de Poynting es

1 ~ ~ E ∧ B = c ω0 zb = c ε0 E 2 zb S~0 = µ0

(c) - Puesto que el vac´ıo no es disipativo, la energ´ıa que atraviesa el casquete es la misma que atraviesa un disco, del mismo radio y con el

Scan109, abril 02, 2008.max

143

I + V −

2a I

2b

L>>a,b Figura 4.11: 4-4. Una onda plana se propaga en el vac´ıo en la direcci´ on del eje z. El campo electrico en z = 0 es

~ =→ E

 b [e−t − e−2t ] , , t > 0  E0 x 

0,, t < 0

donde E0 = 100 µV · m−1 . Halle, para z > 0, a) Campo magn´etico. b) Vector de Poynting. c) Energ´ıa total que atraviesa a un casquete hemisf´erico, de radio a = 1 m, cuyo eje es paralelo al z. SOLUCION: (a) - Para encontrar la expresi´ on del campo electromagn´ etico en cualquier posici´ on z, basta con tener en cuenta que la onda viaja, a trav´ es del vac´ıo, en la direcci´ on positiva del eje z, por lo que t = (t− zc )z=0 y i h z z ~ b E(z, t) = E0 e −(t− c ) − e −2(t− c ) x

Haciendo uso de la relaci´ on de estructura con ~n = zb 1 1 ~ ~ B(z, t) = zb ∧ E(z, t) = E(z, t) yb c c

(b) - el vector de Poynting es

1 ~ ~ E ∧ B = c ω0 zb = c ε0 E 2 zb S~0 = µ0

(c) - Puesto que el vac´ıo no es disipativo, la energ´ıa que atraviesa el casquete es la misma que atraviesa un disco, del mismo radio y con el

144 t=t 0

E 0.25

E

0.2

tangentes iniciales

1

0.15 0.1

^z

0.05

pulso

0.5

0.5

1

1.5

2

2.5

3

z=c t0 (b)

t 0.5

1

1.5

2

2.5

3

x^

-0.5

∆S -1

^z y^

(c)

(a)

Figura 4.12: mismo eje, situado, como se muestra en la figura 4.12-c, en z = 0. As´ı, dicha energ´ıa puede calcularse de la forma Z ∞ Z ∞ Z −t 2 2 2 ~ · d~s dt = ε0 c E πa W = S0 e − e −2t dt t=0

z=0

∆S

4-5. Verifique que

t=0

2 2 2 2 ~ =x E b e(−a t −b z +2 a b z t)

puede corresponder a una onda que viaja a trav´es del vac´ıo. En caso afirmativo: a) Determine las condiciones que deben cumplir las constantes a y b. b) Calcule el campo magn´etico. c) Suponga que dicha onda incide, desde el vac´ıo, sobre un conductor ideal que ocupa el semiespacio z > 0 y halle el campo electromagn´etico total resultante. ´mico, se dice que un conductor es ideal si su conducEn el caso dina tividad es infinita. Como se comprobar´ a m´ as adelante, dicho conductor anula las componentes tangenciales del campo el´ectrico. SOLUCION : (a) - Obviamente, si de la forma

a b

=c=

√1 , µ0 ε0

el campo el´ ectrico puede escribirse

2 2 ~ E(z, t) = x b e −b (z−ct)

145 que pone de manifiesto que E = f (z − ct) es una onda que viaja en el sentido positivo del eje z. (b) - El campo magn´ etico viene dado por la relaci´ on de estructura de las onda electromagn´ eticas en el vac´ıo 1 1 ~ ~ t) = E(z, t) yb B(z, t) = zb ∧ E(z, c c

(c) - Si la onda incide sobre el conductor, debe existir, adem´ as de la incidente, una onda reflejada, de forma que, para z < 0 h i i r −b2 (z−ct)2 −b2 (z+ct)2 ~ ~ ~ Et (z, t) = E (z, t) + E (z, t) = x b e + Ae

El argumento de la onda reflejada lleva un signo + puesto que ´ esta viaja en sentido contrario al de la incidente. La constante A se determinar´ a mediante las condiciones de frontera en z = 0. Et (z = 0) = 0 ⇒ A = −1 Aplicando por separado la relaci´ on de estructura a las ondas incidente y reflejada, a la primera con ~n = zb y a la segunda con ~n = −b z, i h ~ t (z, t) = B ~ i (z, t) + B ~ r (z, t) = 1 yb e −b2 (z−ct)2 + e −b2 (z+ct)2 B c

4-6. Para una onda plana monocrom´ atica que se propaga en un medio con: ~ t) = (b E(z, y + 2b x) cos(109 t + 30z) V · m−1 Calcule: a) La direcci´ on de propagaci´ on. b) El vector unitario en la direcci´ on de la polarizaci´ on. c) La frecuencia f , la angular ω y el periodo T de la onda. d) El n´ umero de onda k y la longitud de onda λ. e) la velocidad de fase ¿Se propaga esta onda en el vac´ıo? f) El campo magn´etico.

4-7. Considere una onda electromagn´etica que se propaga en la direcci´ on x b en el vac´ıo, la cual se descompone en la suma de otras dos: (a) -

~ = yˆ E1 cos(ωt − kx) + yˆ E2 cos(ωt − kx + α) E

~ y hPi. ~ Los valores obtenidos ¿son la suma de los correspondientes a Obtenga P cada componente de la onda ?

146 (b) - Repita el problema si ~ = yˆ E1 cos(ωt − kx) + zˆ E2 cos(ωt − kx + α) E SOLUCION: (a) - Aplicando la relaci´ on de estructura ~ = 1x ~ = E zb = zb (E1 cos ϕ + E2 cos (ϕ + α)) , ϕ = ωt + kx B b∧E c c c

y el vector de Poynting

1 2 ~ ∧B ~ =x ~ = 1 E E1 cos2 ϕ + E22 cos2 (ϕ + α) + 2 E1 E2 cos ϕ cos (ϕ + α) b P µ0 µ0 c

Dado que cos2 a = 21 (1 + cos 2a) , cos a cos b = 21 [cos (a + b) + cos (a − b)] 



(c)

}| { z b  2  2 2 2 ~ =1 x P E1 + E2 + 2 E1 E2 cos α + E1 cos 2ϕ + E2 cos (2ϕ + 2α) + 2 E1 E2 cos (2ϕ + α) {z } | {z } 2 µ0 c | (a)

(b)

Como puede verse, los t´ erminos (b) oscilan en el tiempo sobre un valor medio nulo, por lo que el valor medio viene dado por los (a). Las dos componentes de la onda no constituyen modos independientes, por lo que aparece el sumando (c) que representa a la interferencia entre las dos componentes   ~ = hPi

1 x b  2  2 E1 + E2 + |2 E1 E{z 2 cos α } 2 µ0 c (c)

(b) - En este segundo caso las dos componentes est´ an polarizadas en direcciones distintas, por lo que son independientes. Aplicando la relaci´ on de estructura se tiene que el valor medio del vector de Poynting es ~ = hPi

b 2 1 x E1 + E22 2 µ0 c

Comprobamos que en este caso no existe t´ ermino de interferencia y cada una aporta independientente su contribuci´ on al flujo de potencia. ~ = E0x cos (ω t − k z) x 4-8. Sea el campo E b + E0y cos (ω t − k z + ϕ) yb = x x b + y yb. Dibuje ~ la trayectoria del extremo del vector E en el plano z = 0, a lo largo del tiempo. ´ polarizada linealmente si la trayectoria Se dice que una onda esta es una recta, el´ıpticamente , si es una elipse, y circularmente si es una circunferencia.¿Qu´e condiciones deben cumplir los par´ ametros E0x , E0y y ϕ para que la polarizaci´ on sea de cada uno de los tipos arriba mencionados?

147 4-9. Demuestre que los siguientes campos el´ectricos pueden corresponder a ondas electromagn´eticas planas en el vac´ıo, que se propagan seg´ un el eje z, y calcule los campos magn´eticos asociados. a)

~ = E0 e ±j ω z/c e j ω t x E b

b)

~ = E0 cos (ω z ) cos (ω t) yb E c

Calcule el vector de Poynting y analize en que forma transporta la energ´ıa cada una de estas ondas. ~ = E0 cos k(z − ct) x 4-10. Demuestre que la onda plana E b puede ser descrita con s´ olo un potencial vector que, adem´ as, cumple la condici´ on de Culomb. Soluci´ on:

El vector unitario de propagaci´ on es ~n = zb. Escribiendo α = k(z − ct) y haciendo uso de la relaci´ on de estructura ~ = E0 cos α x ~ = E0 cos α yb E b , B c

Debemos comprobar que se llega a este mismo resultado siguiendo el enunciado del problema: ~ ~ ~ ~ = −∂ A = −d A ∂ α = k c d A E ∂t dα ∂ t dα ~ = Ax Luego A by

Z E0 E0 sen α cos α dα + cte = kc kc La constante de integraci´ on se ha tomado como nula. A=

Dado que el campo no depende de x ni de y ~ ~ ~ =∇∧A ~ = yb ∂ A = yb d A ∂ α = E0 cos α yb B ∂z dα ∂ z c

Por u ´ ltimo, dado que hemos supuesto que, en este caso, V puede tomarse como nulo, el contraste de Coulomb y el de Faraday se expresan de la misma forma. 4-11. Demuestre que las ondas f (u) (u = x − ct), que se propagan en la direcci´ on ~n = +b x, cumplen la ecuaci´ on

∂ 1 ∂ + ∂x c ∂t



f (u) = 0

148 y las g(w) (u = x + ct), que lo hacen en la ~n = −b x, cumplen la ecuaci´ on

.

∂ 1 ∂ − ∂x c ∂t



g(w) = 0

4-12. El ritmo medio al que la energ´ıa solar incide sobre la Tierra es aproximadamente 1400 W · m−2 . a) Calcule el valor del campo el´ectrico en la superficie terrestre considerando a la luz solar como monocrom´ atica y linealmente polarizada. b) Si el Sol radia isotr´ opicamente, ¿con que potencia lo hace ? La distancia de la Tierra al Sol es 1,49 × 108 km.

c) Calcule la potencia total recibida por la Tierra sabiendo que su radio es 6,37× 103 km.

Se dice que un cuerpo radia de forma is´ otropa cuando la intensidad ´ n es independiente de la direccio ´ n. de radiacio 4-13. Considere un transmisor que colocado en la Luna radia isotr´ opicamente a la frecuencia de 5 GHz y con una potencia de 1 W . Sabiendo que la distancia TierraLuna es 3,8 × 105 km, calcule: a) El valor de los campos E y B en la superficie de la Tierra. b) El valor medio del vector de Poynting en la superficie terrestre. c) La densidad media de energ´ıa. d) El tiempo que tarda una se˜ nal en alcanzar la Tierra. 4-14. Sean dos cargas puntuales ± q(t), en el vac´ıo, situadas respectivamente en los puntos ± d2 zb. a) Halle el potencial retardado producido en un punto arbitrario ~r.

b) Haga lo mismo, en la zona lejana (r ≫ λ0 y r ≫ d), para q(t) = q0 cos ω0 t.

c) Aproxime el resultado anterior para un dipolo el´ectricamente peque˜ no (d << λ0 ).

SOLUCION: (a) - El potencial retardado viene dado por la integral 1 V (~r, t) = 4πε0

Z

V′

[ρ] ′ dv R

Si las dos cargas se situan sobre el eje z, figura 4.13 la densidad puede describirse como

149 ^z

∆V

P R1

+ d ^z 2

r

11 00

^ n

θ

R2

^y -

d ^ z 2

11 00 x^

Figura 4.13: 1 1 ′ ′ ρ(z ) = q(t) δ(~r − d zb) − δ(~r + d zb) 2 2 ′

La densidad retardada es, por lo tanto, R [ρ] = q(t − ) c Realizando la integral

1 1 ′ ′ δ(~r − d zb) − δ(~r + d zb) 2 2

1 V (~r, t) = 4πε0 donde,

(

q(t − Rc1 ) q(t − Rc2 ) − R1 R2

)

(4.46)

~ 2 = ~r + 1 d zb ~ 1 = ~r − 1 d zb , R R 2 2

(b) - En el caso de que q(t) sea una funci´ on arm´ onica la expresi´ on 4.46 toma la forma q0 V (~r, t) = 4πε0 Los t´ ermino

(

cos ω0 (t − R1

R1 c )

cos ω0 (t − − R2

R2 c )

)

(4.47)

1 podemos aproximarlos en la zona lejana como R

1 1 ≃ R r pero este tipo de aproximaci´ on no debemos hacerlo en el argumento del coseno porque el potencial resultante ser´ıa nulo. Haremos una aproximaci´ on que nos de un resultado significativo.

150 Definamos δ =r−R con lo que

cos ω0 (t −

R ) = cos (ω0 t − k0 r + k0 δ) c = cos (ω0 t − k0 r) cos (k0 δ) − sen (ω0 t − k0 r) sen (k0 δ)

2π ω0 = c λ0 Si R ≫ d, seg´ un la figura 4.13, podemos hacer la aproximaci´ on

donde k0 =

δ = ~n · ~r ′ =

  δ1 = 21 d cos θ 

para ~r ′ = 12 d

δ2 = − 12 d cos θ para ~r ′ = − 21 d

de acuerdo con ´ esto la ecuaci´ on 4.47 se aproxima por V (~r, t) ≃ −

d q0 sen ω0 (t − rc ) sen (k0 cos θ) 2πε0 r 2

(c) - En este caso d << λ0 , luego y x=

1 d k0 d cos θ = π cos θ << 1 ⇒ sen x ≃ x ⇒ 2 λ0

El potencial retardado producido por un dipolo corto es, en consecuencia V (~r, t) ≃ −

q0 d sen ω0 (t − rc ) cos θ 2ε0 λ0 r

4-15. La figura 4.13 muestra dos cargas puntuales, de igual magnitud y signo opuesto, unidas mediante un conductor de longitud d por el que circula una corriente tal que q = q0 cos ωt. Bajo la condici´ on de campo lejano, r >> d, y suponiendo que el dipolo es el´ectricamente corto, d << λ, halle La potencia radiada. 16 Soluci´ on: Para este apartado seguiremos los pasos dados en la secci´ on 4.4 Previamente calcularemos la intensidad aplicando la ley de conservaci´ on de la carga: La intensidad que entra en el volumen V es igual a la raz´ on de incremento temporal de la carga almacenada en el mismo. 16

Esta u ´ltima condici´ on, como puede verse en la secci´ on E.2, implica que la intensidad es cuasiestacionaria, es decir, que I = I(t) es aproximadamente independiente de z.

151

I=

dq = −q0 ω sen ω t dt

Puesto que se trata de un conductor filiforme, el potencial vector viene dado por

Aproximando R ≃ r

donde C =

~ r, t) = zb µ0 A(~ 4π

Z

1 d 2

− 21 d

I(t −

R ′ c ) dz

R

~ r, t) = −C zb 1 sen u A(~ r

µ0 ω p0 r , p0 = q0 d , u = ω (t − ). 4π c

1 1 Si despreciamos el t´ ermino proporcional a ∇ ∼ 2 , y tenemos en r r rb cuenta que ∇ u = −ω c µ0 ω 2 p0 1 r ~ =∇∧A ~ = −C 1 ∇ ∧ (b z sen u) = − cos ω (t − ) sen θ ϕ b B r 4π c r c 2

~ = c B rb sobre la superficie de una Integrando el vector de Poynting P µ0 esfera de radio r obtenemos la potencia radiada P (r, t) =

r µ0 ω 4 p20 cos2 ω (t − ) 6π c c

La potencia radiante de un dipolo el´ ectrico crece con la cuarta potencia de la frecuencia, por lo que las oscilaciones de baja frecuencia radian muy poca energ´ıa. 4-16. Halle la potencia que radiar´ıa un electr´ on cl´ asico girando alrededor de un n´ ucleo de hidr´ ogeno, en una o ´rbita de radio a y con velocidad angular uniforme ω. 4-17. Halle la intensidad radiada en la colisi´ on frontal de dos part´ıculas, con carga de distinto signo, que se mueven en las proximidades del origen. Suponga que una de las part´ıculas tiene una masa M muy superior a la de la otra, m, y que el movimiento de ambas es lento.

152

4.6.

Resoluci´ on de las ecuaciones de Maxwell unidimensionales mediante el m´ etodo FD–TD: F DT D 1D − vacio.nb

17

Lo que expondremos a continuaci´ on ilustra la soluci´ on num´erica de la ecuaci´ on de ondas mediante el uso del m´etodo de las Diferencias Finitas en el Dominio del Tiempo (FD-TD) 18 .

4.6.1.

La ecuaci´ on de onda unidimensional

Las ecuaciones rotacionales de Maxwell, para ondas que est´an polarizadas en la direcci´ on del eje y y se propagan, en el vac´ıo, a lo largo del eje x, pueden expresarse como

∂ ey (x, τ ) ∂ Bz (x, τ ) =− ∂τ ∂x

(4.48)

∂ ey (x, τ ) ∂ Bz (x, τ ) =− ∂τ ∂x

(4.49)

para lo cual se han normalizado el campo el´ectrico y el tiempo

ey ≡

19 .

Ey , τ ≡ ct c

Eliminando el campo magn´etico, o el el´ectrico, de las ecuaciones anteriores se obtienen las ecuaci´ ones de onda para estos campos: ∂ 2 ey (x, τ ) ∂ 2 ey (x, τ ) ∂ 2 Bz (x, τ ) ∂ 2 Bz (x, τ ) = , = ∂ x2 ∂ τ2 ∂ x2 ∂ τ2 La soluci´ on general compatible con las ecuaciones de Maxwell tiene la forma

ey (x, τ ) = f (x − τ ) + g(x + τ ) , Bz (x, τ ) = f (x − τ ) − g(x + τ )

(4.50)

en la que f y g son funciones arbitrarias, aunque de buen comportamiento, que se propagan respectivamente en el sentido positivo y negativo del eje x. 17

Basado en un programa de S. Gonz´ alez. Finite Difference–Time Domain.V´ease [Taflove]. 19 Estas definiciones son ventajosas para el c´ alculo num´erico ya que, en el vac´ıo, para una onda que se propaga en un sentido determinado, e = B. 18

153 Condiciones iniciales: Para obtener la soluci´ on temporal a partir de las condiciones iniciales ey (x, 0) = f (x) + g(x) , Bz (x, 0) = f (x) − g(x) basta con substituir los argumentos ’x’ de dichas funciones por ’x ± τ ’, seg´ un se trate de f (x − τ ) o de g(x + τ ). Por ejemplo, si se excita inicialmente el espacio con un campo el´ectrico de perfil gausiano

1

ey (x, 0) = e − a2 (x−x0 )

2

, Bz (x, 0) = 0

(4.51)

la soluci´ on a lo largo del tiempo toma la forma ey (x, τ ) =

Bz (x, τ ) =

1 2 1 − 12 (x−τ −x0 )2 e a + e − a2 (x+τ −x0 ) 2 1 2 1 − 12 (x−τ −x0 )2 − e − a2 (x+τ −x0 ) e a 2

que representa a dos ondas de igual amplitud y que viajan en sentidos opuestos. Puede comprobarse que cada una de ´estas cumple su respectiva relaci´ on de estructura y que sus vectores de Poyntyng tienen las direcciones ±b x. Por otra parte, para τ = 0, cumplen con las con las condiciones iniciales 4.51. Por el contrario, si la excitaci´ on es del tipo

1

ey (x, 0) = e − a2 (x−x0 )

2

, Bz (x, 0) = +ey (x, 0)

(4.52)

la soluci´ on temporal es una sola onda que viaja en la direcci´ on +b x.

1

ey (x, τ ) = e − a2 (x−τ −x0 )

4.6.2.

2

1

, Bz (x, τ ) = e − a2 (x−τ −x0 )

2

Soluci´ on num´ erica

Para resolver num´ericamente estas ecuaciones, se definir´a un dominio num´erico en cuyos nudos se tomar´ an muestras del valor de los campos y sobre el cual se definir´ an operadores num´ericos de diferenciaci´ on espacial y temporal basados en diferencias centradas. Con estas herramientas, se dise˜ nar´ an algoritmos de avance temporal de los campos que permiten, en sucesivas iteraciones, actualizar sus valores en cada nudo a lo largo del tiempo.

154

n=τ/δτ

campo magnético

iteración ni campo eléctrico

2 ni (i-1,n+1)

iteración 1

condiciones iniciales

A

2

c

+ A

1

1

1

0 1

celda 1

2

-A

1

+

a 0

condición absorbente

(i,n)

-A

b

3

2 nc-2

celda 2

2nc

i=x/δx

celda nc celda nc+1

condiciones de contorno Figura 4.14: Red num´erica FD–TD Dominio num´ erico: Como se muestra en la figura 4.14, el espacio y el tiempo se discretizan a intervalos δx y δτ de lo que resulta un dominio num´erico constituido por los nudos

x = iδx ,

i = 0, · · · 2nc

τ = nδτ ,

n = 0, · · · 2ni + 1

(4.53)

Los nudos se agrupan por parejas, en celdas, segun la direcci´ on espacial, y en iteraciones seg´ un la temporal. El campo el´ectrico se evalua en los nudos pares i = 0, 2, · · · 2nc y en instantes pares n = 0, 2, · · · 2ni (en la figura ’⋄’) y el magn´etico en los impares i = 1, 3, · · · 2nc − 1 y en instantes impares i = 1, 3, · · · 2ni + 1 (en la figura ’◦’). nc es el n´ umero de celdas completas; en la u ´ltima, 2nc + 1, solo se toma la muestra del campo el´ectrico. Los nudos de coordenada temporal n = 0, 1 corresponden a la iteraci´ on 0 y contienen las condiciones iniciales. El resto contiene a los campos en las sucesivas iteraciones (n = 2, 4 · · · , 2ni. En cada una de las iteraciones se calcula, mediante el algoritmo de avance temporal, primero el campo el´ectrico y posteriormente el magn´etico. Desarrollo del programa:

155 El programa que proponemos a continuaci´ on puede dividirse esquem´ aticamente en las siguientes etapas: 1. Establecimiento de las condiciones iniciales. 2. C´ alculo de los campos en iteraciones sucesivas mediante los algoritmos de avance temporal y las condiciones de contorno: a) C´ alculo del campo el´ectrico de los nudos interiores y en el instante (n). b) C´ alculo del campo magn´etico en el instante (n + 1). c) C´ alculo del campo el´ectrico en los nudos extremos y en el instante (n). 3. Generaci´ on de los fotogramas correspondientes a las distintas iteraciones. Diferencias finitas centradas: Es posible resolver num´ericamente las ecuaciones rotacionales de Maxwell mediante la aproximaci´ on de los operadores diferenciales por otros num´ericos en diferencias centradas. Con este fin, se divide el espacio α en intervalos de longitud h y se define el operador derivada centrada de una funci´ on f (α) de la forma Dα [f (α)] =

1 {f (α + h) − f (α − h)} 2h

Desarrollando en serie de Taylor f (α + h) y f (α − h), restando, despejando y despreciando t´erminos de orden O(h2 ), se tiene que d f (α) ≃ Dα [f (α)] dα Este operador aproxima al operador derivada anal´ıtica hasta el segundo orden en h. Avance temporal: El algoritmo que permite avanzar temporalmente al campo en los distintos puntos de la red se obtiene substituyendo a los operadores diferenciales por los correspondientes en diferencias. Con este f´ın, en cada iteraci´ on n = 2, 4 · · · 2ni seguiremos los siguientes pasos: C´ alculo del campo electrico. El circulo (b) de la figura 4.14 est´ a centrado en el punto (i, n − 1), marcado por una cruz, y contiene una estrella de cuatro puntas cuyo nudo superior es el (i, n), donde i y n son pares. Para calcular el campo el´ectrico en este u ´ltimo nudo, evaluamos la ecuaci´ on 4.48 en el centro de la estrella y despejamos ey (i, n). De esta forma se obtiene una expresi´ on que relaciona a este campo con ´el mismo y con el magn´etico en instantes anteriores.

156 Simplificaremos la notaci´ on escribiendo: ey ≡ e , Bz ≡ B , f (i, n) ≡ fin – Ecuaci´ on 4.48 evaluada en (i,n-1). n−1 n−1 , i = 2, 4 · · · 2nc − 2 + A Bi−1 − Bi+1 eni = en−2 i

(4.54)

donde se excluyen los nudos frontera, en los que hay que aplicar las conciciones de contorno. C´ alculo del campo magn´etico. El circulo (a) est´ a centrado en el punto (i − 1, n) y el nudo superior de la estrella correspondiente es el (i − 1, n + 1). Para calcular el campo magn´etico en este u ´ltimo nudo, evaluamos la ecuaci´ on 4.49 n+1 en el centro de la estrella y despejamos Bi−1 . – Ecuaci´ on 4.49 evaluada en (i-1,n). n+1 n−1 Bi−1 = Bi−1 + A eni−2 − eni , i = 2, 4 · · · 2nc − 2

(4.55)

En el programa siguiente, en cada iteraci´ on se aplicar´ an los algoritmos de avance anteriores, para ambos campos y para i = 2, 4 · · · 2nc − 2, con lo que, v´ease la n+1 figura 4.14, no se cubre a la componente B2nc−1 la cual habr´ a de calcularse posteriormente. La estructura de los algoritmos se muestran dentro de los c´ırculos (a) y (b) de la figura 4.14 Condici´ on de Courant: El coeficiente que aparece en el segundo miembro es

A=

δτ δx

Puede demostrarse que, para que los algoritmos num´ericos sean estables, es necesario que se cumpla la condici´ on de Courant δx ≥1 δτ En adelante se tomar´ a A = 1 ⇒ δx = δτ = δ

(4.56)

157 Condiciones iniciales: Las condiciones iniciales se imponen en los dos primeros instantes (n = 0, 1) (figura 4.14)

e0i = e0 (i) , Bi1 = B1 (i) ,

i = 0, 2, · · · , 2nc i = 1, 3, · · · , 2nc − 1

(4.57)

Condiciones de contorno: Dado que el dominio num´erico es finito, X = [0 ≤ x ≤ L], en sus extremos deben aplicarse condiciones de contorno adecuadas. En el primer ejemplo supondremos X limitado por planos conductores ideales x = 0 yx=L´ o i = 0 e i = 2 nc. Obs´ervese en la figura 4.14 que este dominio se ha configurado de forma tal que en los nudos extremos s´ olo est´ a definido el campo el´ectrido. Dado que ´este es tangencial a los conductores, y estos son perfectos, debe anularse en la superficie. Impondremos, por lo tanto, condiciones reflectantes

Condiciones reflectantes →

 n  e0 = 0 

(4.58)

en2nc = 0

En el segundo ejemplo se impondr´ an condiciones absorbentes, con las cuales se simula la existencia en los puntos extremos de superficies perfectamente absorbentes (no reflectantes), de forma que una onda incidente sobre las mismas contin´ ua propag´ andose, sin reflejarse, como si el medio correspondiente se extendiese m´ as all´ a del la frontera num´erica. Las ondas que inciden desde el interior de X sobre el punto x = 0 son las f (w = x+τ ), que viajan enel sentido negativo del eje x, y su velocidad de fase normalizada es dx v− = = −1. De forma an´ aloga, las que inciden sobre el punto ν− ≡ c d τ w=cte x = L son las f (u = x − τ ), que viajan en el sentido positivo del eje x, y cuya velocidad de fase normalizada es ν+ = 1. En cualquiera de estos dos casos, tomando δx = δτ = δ, las ondas recorren un espacio 2 δ en el tiempo 2 δ y las condiciones absorbentes toman la siguiente forma

Condiciones absorbentes →

 n n−2  e0 = e2 

en2nc

=

(4.59)

en−2 2nc−2

La curva (c) de la figura 4.14 delimita el esquema de la condici´ on absorbente en 2nc.

158 Programa Mathematica F DT D 1D − vacio.nb: En este programa se hace uso de funciones discretas f (i, n). Ejemplo 1o En este primer ejemplo se simula una onda de perfil gausiano que viaja inicialmente hacia la derecha y que se refleja sucesivamente sobre dos planos conductores. Remove[”Global‘ ∗ ”]; Off [General :: ”spell1”]; $TextStyle = {FontFamily → ”Courier”, FontSize → 12}; f (i, n) la escribiremos de la forma. 1 ∗ (i − n − nc)2 ]; na2 en la que el punto central del pulso se ha situado en i = nc en el instante inicial n = 0. Se asignan valores concretos al n´ umero de celdas nc, al de iteraciones ni y al factor que define la anchura del pulso na. f [i , n ] := Exp[−

nc = 60; ni = 120; na = 8; Los valores de los campos en cada posici´ on (i) y en cada instante (n) se almacenan en las funciones ey (i, n), con (i, n) pares, y Bz (i, n) 20 , con (i, n) impares. Para cada posici´ on (i), estas funciones se incluyen en listas con dos elementos {i, ey (i, n)} y {i, Bz (i, n)}. Se comienza por definir las condiciones iniciales en funci´ on de f (i, n). te0 = Table[{i, ey[i, 0] = N[f [i, 0]]}, {i, 0, 2 ∗ nc, 2}]; ey[0, 0] = 0; ey[0, 2 ∗ nc] = 0 Se anula el campo el´ectrico en los nudos frontera porque las condiciones en los mismos son reflectantes. tB1 = Table[{i, Bz[i, 1] = N[f [i, 1]]}, {i, 1, 2 ∗ nc − 1, 2}]; Se hace uso de la orden N[] para que la tabla se llene con los valores num´ericos y no de expresiones anal´ıticas. En caso contrario el c´ alculo se ralentiza y el ordenador puede saturarse. Los representamos gr´ aficamente, sin mostrarlos, ey(i, 0) en rojo y Bz(i, 1) en azul: grey0 = ListPlot[te0, PlotStyle → {PointSize[0.015], RGBColor[1, 0, 0]}, DisplayFunction → Identity, PlotRange → {0, 1}];

20 Esto supone guardar en memoria los campos correspondientes a todas las iteraciones que se realicen a lo largo de la ejecuci´ on del programa, lo que es poco eficiente y resulta prohibitivo para problemas de dos y tres dimensiones. El esquema FDTD permite guardar en memoria u ´nicamente los datos m´ as 1 recientes para cada posici´ on i, lo que reduce la necesidad de memoria en el factor 2nc . Haremos uso de esta posibilidad en la secci´ on B.2.1.

159

grBz1 = ListPlot[tB1, PlotStyle → {PointSize[0.015], RGBColor[0, 0, 1]}, DisplayFunction → Identity, PlotRange → {0, 1}];

y los mostramos conjuntamente. Show[grey0, grBz1, DisplayFunction → $DisplayFunction]; La figura 4.15 presenta solamente la parte significativa de las funciones. Puede verse como las muestras de cada campo se toman en nudos alternados de la red. Dado que 1 Ey = f (x) y Bz = g(x) = f (x), este pulso corresponde a una sola onda cuyo vector c de Poynting tiene la direcci´ on x b. 1

0.8 0.6 0.4 0.2

45

50

55

60

65

70

75

80

Figura 4.15: Para calcular el resto de los valores se emplea un lazo Do que aplica las reglas de avance temporal. Dentro de este lazo se hace uso de otro Do para determinar los campos en los nudos interiores, salvo el campo magn´etico de los nudos 2nc − 1 21 y el el´ectrico de los de la frontera que se calculan fuera de este u ´ltimo lazo.

Do[{Do[{ey[i, n] = ey[i, n − 2] − Bz[i + 1, n − 1] + Bz[i − 1, n − 1],

Bz[i − 1, n + 1] = Bz[i − 1, n − 1] − ey[i, n] + ey[i − 2, n]}, {i, 2, 2 ∗ nc − 2, 2}],

Bz[2 ∗ nc − 1, n + 1] = Bz[2 ∗ nc − 1, n − 1] − ey[2 ∗ nc, n] + ey[2 ∗ nc − 2, n]},

ey[0, n] = 0, ey[2 ∗ nc, n] = 0{n, 2, 2 ∗ ni, 2}];

Para visualizar la propagaci´ on del pulso electromagn´etico, con los colores empleados en la gr´ afica anterior, se representa a los campos en cada instante temporal y se activa la pel´ıcula de la forma ya descrita en un programa anterior. Para que las se˜ nales no se solapen totalmente, Bz se multiplica por el factor 0,9. 21

N´ otese que el algoritmo de avance temporal no puede aplicarse a los nudos de la frontera porque, en ese caso, necesita datos que estar´ıan fuera del dominio num´erico. Compruebe, figura 4.14, que Bz[2*nc - 1, n + 1] queda fuera del lazo.

160

Table[Show[ListPlot[Table[{i, ey[i, n]}, {i, 0, 2 ∗ nc, 2}], Ploted → True, PlotRange → {−2.1, 2.1},

DisplayFunction → Identity, PlotStyle → {RGBColor[1, 0, 0]}],

ListPlot[Table[{i, 0.9 ∗ Bz[i, n + 1]}, {i, 1, 2 ∗ nc − 1, 2}], Ploted → True, PlotRange → {−2.1, 2.1},

DisplayFunction → Identity, PlotStyle → {RGBColor[0, 0, 1]}],

DisplayFunction → $DisplayFunction], {n, 0, 2 ∗ ni, 2}];

La figura 4.16 corresponde al fotograma no 35 de la pel´ıcula. En ´el se ve como la onda se ha reflejado ya en el extremo derecho, donde se simula al conductor ideal. El campo el´ectrico, que es tangencial, se anula en ese punto, mientras que el magn´etico se refuerza. Al mismo tiempo, la onda se propaga hacia la izquierda: el Bz reflejado no cambia de signo mientras que Ey lo inviete; el vector de Poynting cambia su sentido a −b x. 0.75 0.5 0.25 90

100

110

120

-0.25 -0.5 -0.75 -1

Figura 4.16: Ejemplo 2o En este caso se simula, junto con la onda del ejemplo anterior, otra que viaja en sentido contrario y se hace uso de las condiciones de frontera absorbentes. Remove[”Global‘ ∗ ”]; Se elijen los valores de los par´ ametros. nc = 60; ni = 40; na = 8; n0 = 40; 1 ∗ (i − n − nc)2 ]; na2 La segunda onda se define con un perfil que resulta de modular el pulso gausiano del primer ejemplo con una funci´ on senoidal de periodo T = n0 δ f [i , n ] := Exp[−

ge[i , n ] := f [i, n] ∗ 4 ∗ Sin[

π ∗ (i + n − nc) ]; n0

161 Se generan las gr´ aficas sin mostrarlas. f0 = Table[{i, f [i, 0]}, {i, 0, 2 ∗ nc, 2}]; ge0 = Table[{i, ge[i, 0]}, {i, 0, 2 ∗ nc, 2}];

grf = ListPlot[f0, Ploted → True, PlotRange → {−1.2, 1.2},

PlotStyle → {RGBColor[1, 0, 0]}, DisplayFunction → Identity];

grge = ListPlot[ge0, Ploted → True, PlotRange → {−1.2, 1.2}, PlotStyle → {RGBColor[0, 0, 1]}, DisplayFunction → Identity];

La figura 4.17 muestra conjuntamente ambas funciones. Show[grf , grge, DisplayFunction → $DisplayFunction]; 1 0.5

20

40

60

80

100

120

-0.5 -1

Figura 4.17: A partir de estas funciones se obtienen los valores iniciales de ambos campos ey0b = Table[ey[i, 0] = N[f [i, 0] + ge[i, 0]], {i, 0, 2 ∗ nc, 2}]; Bz1b = Table[ Bz[i, 1] = N[f [i, 1] − ge[i, 1]], {i, 1, 2 ∗ nc − 1, 2}]; se generan las gr´ aficas

grey0b = ListPlot[Table[{i, ey[i, 0]}, {i, 0, 2 ∗ nc, 2}],

PlotStyle → {PointSize[0.015], RGBColor[1, 0, 0]}, PlotRange → {−1, 2},

DisplayFunction → Identity];

162

grBz1b = ListPlot[Table[{i, Bz[i, 1]}, {i, 1, 2 ∗ nc − 1, 2}],

PlotStyle → {PointSize[0.015], RGBColor[0, 0, 1]}, PlotRange → {−1, 2},

DisplayFunction → Identity];

En la figura 4.18 se muestra la suma de los pulsos iniciales. Dado que ey = f (x) y Bz = g(x) 6= f (x) estos campos se desglosan en dos ondas que se propagan en sentidos opuestos. Show[grey0b, grBz1b, DisplayFunction → $DisplayFunction]; 2 1.5 1 0.5 20

40

60

80

100

120

-0.5 -1

Figura 4.18: Por u ´ltimo, se forman las tablas ey(i, n) y Bz(i, n), se generan las gr´ aficas temporales, como en el primer ejemplo, y se animan. Do[{Do[{ey[i, n] = ey[i, n − 2] − Bz[i + 1, n − 1] + Bz[i − 1, n − 1],

Bz[i − 1, n + 1] = Bz[i − 1, n − 1] − ey[i, n] + ey[i − 2, n]}, {i, 2, 2 ∗ nc − 2, 2}],

Bz[2 ∗ nc − 1, n + 1] = Bz[2 ∗ nc − 1, n − 1] − ey[2 ∗ nc, n] + ey[2 ∗ nc − 2, n],

ey[0, n] = ey[2, n − 2], ey[2 ∗ nc, n] = ey[2 ∗ nc − 2, n − 2]}, {n, 2, 2 ∗ ni, 2}] Table[Show[ListPlot[Table[{i, ey[i, n]}, {i, 0, 2 ∗ nc, 2}], Ploted → True, PlotRange → {−2.2, 2.2}, DisplayFunction → Identity,

PlotStyle → {RGBColor[1, 0, 0]}],

ListPlot[Table[{i, 0.9 ∗ Bz[i, n + 1]}, {i, 1, 2 ∗ nc − 1, 2}], Ploted → True,

PlotRange → {−2.2, 2.2}, DisplayFunction → Identity,

PlotStyle → {RGBColor[0, 0, 1]}],

DisplayFunction → $DisplayFunction], {n, 0, 2 ∗ ni, 2}]; La figura 4.19a corresponde al fotograma no 20. Cada uno de los pulsos se ha propagado en sentido opuesto pero a´ un no ha llegado al extremo correspondiente. La figura 4.19b corresponde al fotograma no 30. Como puede observarse, los pulsos traspasan los l´ımites sin reflejarse.

163

1

1

0.5

0.5

20

40

60

80

100

120

20

-0.5

-0.5

-1

-1

(a)

40

60

(b)

Figura 4.19:

80

100

120