Universidad Nacional Autónoma de México

Introducción a la Física Computacional

Implementación del Perfectly Matched Layer (PML) al código de FDTD

Autor: Lic. Favio Vázquez

Supervisor: Dr. Yevgeniy Kolokoltsev

14 de noviembre de 2016

Introducción teórica al

PML

Introducción teórica al PML

¿Qué es un PML?

En algunas situaciones necesitamos modelar un problema electrodinámico en el cual, entre comillas, no existan fronteras, es decir, un problema de fronteras abiertas.

Para hacer esto comúnmente se utilizan técnicas y métodos que absorven las ondas en los bordes de una simulación FDTD.

Introducción teórica al PML

¿Qué es un PML?

El modo en que el PML resuelve este problema es introduciendo una capa de grosor finito a las fronteras del dominio de simulación, la cual absorve ondas a través de conductividades eléctricas y magnéticas, no físicas, causando un decaimiento exponencial de la amplitud del campo.

Las reflexiones adentro del espacio de simulación con PML se eliminan apareando perfectamente las impedancias de los dos medios, y calculando con cuidado las conductividades de PML.

Introducción teórica al PML

El primer modo efectivo de aplicar un PML fue introducido por J. P. Bérenger en 1994; el cual consistía en la reflexión y refracción de ondas planas en la interfaz entre un medio dieléctrico sin pérdidas y otro con pérdidas.

Para eso es necesario separar la onda plana en dos modos ortogonales, y encontrar la condición para la no reflexión, y luego lo único que queda es separar estas ecuaciones dentro de las ecuaciones de Maxwell.  

H_z
HzH_z

Esta posibilidad fue encontrada por Bérenger en este artículo, y está sujeta a la introducción de dos conductancias independientes para el mismo vector de campo magnético, lo cual no es físico, ya que         solo tiene una dirección, pero es un artificio numérico

Introducción teórica al PML

Ecuaciones relevantes (modo TM):

\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
ϵ2∂Ex∂t+σyEx=∂Hz∂y\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
ϵ2∂Ey∂t+σxEy=−∂Hz∂x\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
\mu_2 \frac{\partial H_{z}}{\partial t} + \sigma_{m} H_{z} = \frac{\partial E_y}{\partial x} - \frac{\partial E_y}{\partial x}
μ2∂Hz∂t+σmHz=∂Ey∂x−∂Ey∂x\mu_2 \frac{\partial H_{z}}{\partial t} + \sigma_{m} H_{z} = \frac{\partial E_y}{\partial x} - \frac{\partial E_y}{\partial x}

Caso isotrópico

Caso anisotrópico

H_z = H_{zx} + H_{zy}
Hz=Hzx+HzyH_z = H_{zx} + H_{zy}
\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
ϵ2∂Ex∂t+σyEx=∂Hz∂y\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
ϵ2∂Ey∂t+σxEy=−∂Hz∂x\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
\mu_2 \frac{\partial H_{zx}}{\partial t} + \sigma_{m,x} H_{zx} = - \frac{\partial E_y}{\partial x}
μ2∂Hzx∂t+σm,xHzx=−∂Ey∂x\mu_2 \frac{\partial H_{zx}}{\partial t} + \sigma_{m,x} H_{zx} = - \frac{\partial E_y}{\partial x}
\mu_2 \frac{\partial H_{zy}}{\partial t} + \sigma_{m,y} H_{zy} = \frac{\partial E_x}{\partial y}
μ2∂Hzy∂t+σm,yHzy=∂Ex∂y\mu_2 \frac{\partial H_{zy}}{\partial t} + \sigma_{m,y} H_{zy} = \frac{\partial E_x}{\partial y}

Introducción teórica al PML

Ecuaciones relevantes (modo TM):

H_z = H_{zx} + H_{zy}
Hz=Hzx+HzyH_z = H_{zx} + H_{zy}
\frac{\sigma_x}{\epsilon_1} = \frac{\sigma_{m,x}}{\mu_1}
σxϵ1=σm,xμ1\frac{\sigma_x}{\epsilon_1} = \frac{\sigma_{m,x}}{\mu_1}
H_z = H_0 (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})
Hz=H0(e−jk1xx+k1yy)(eσxxη1cosθi)H_z = H_0 (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})
E_x = - H_0 \eta_1 \sin{\theta_i} (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})
Ex=−H0η1sinθi(e−jk1xx+k1yy)(eσxxη1cosθi)E_x = - H_0 \eta_1 \sin{\theta_i} (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})
E_y = H_0 \eta_1 \cos{\theta_i} (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})
Ey=H0η1cosθi(e−jk1xx+k1yy)(eσxxη1cosθi)E_y = H_0 \eta_1 \cos{\theta_i} (e^{-jk_{1_x}x+k_{1_y}y})(e^{\sigma_x x \eta_1 \cos{\theta_i}})

Término de atenuación

\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
ϵ2∂Ex∂t+σyEx=∂Hz∂y\epsilon_2 \frac{\partial E_x}{\partial t} + \sigma_y E_x = \frac{\partial H_z}{\partial y}
\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
ϵ2∂Ey∂t+σxEy=−∂Hz∂x\epsilon_2 \frac{\partial E_y}{\partial t} + \sigma_x E_y = - \frac{\partial H_z}{\partial x}
\mu_2 \frac{\partial H_{zx}}{\partial t} + \sigma_{m,x} H_{zx} = - \frac{\partial E_y}{\partial x}
μ2∂Hzx∂t+σm,xHzx=−∂Ey∂x\mu_2 \frac{\partial H_{zx}}{\partial t} + \sigma_{m,x} H_{zx} = - \frac{\partial E_y}{\partial x}
\mu_2 \frac{\partial H_{zy}}{\partial t} + \sigma_{m,y} H_{zy} = \frac{\partial E_x}{\partial y}
μ2∂Hzy∂t+σm,yHzy=∂Ex∂y\mu_2 \frac{\partial H_{zy}}{\partial t} + \sigma_{m,y} H_{zy} = \frac{\partial E_x}{\partial y}

Implementación

Implementación

Se creó una clase que ​hereda las funciones y estructuras de EMField pero que agrega las variables necesarias para dividir el campo. Para hacer esto se modificaron los archivos:

  • emfield.h
  • emfield.cpp

Utilizando DF se modificaron las ecuaciones de campo, para así implementar las ecuaciones mostradas.

  • optics.h
  • optics.cpp

Se creó la caja de cómputo, haciendo la geometría de la frontera en la cual se implementaría el PML.

  • visual_field.cpp

Se modificó la forma en la cual se calcula la energía de los campos para que se tomara en cuenta la división del campo magnético.

Implementación

void EMFieldPML::eval_cell_E(int i, int j){

    TECell& cij = field[i][j];
    TECell& cij1 = c(i,j-1);
    TECell& ci1j = c(i-1,j);

    double Hij = cij.Hzx + cij.Hzy;
    double Hij1 = cij1.Hzx + cij1.Hzy;
    double Hi1j = ci1j.Hzx + ci1j.Hzy;

cij.Ex += (dt/dl*(1/epsilon))*( (Hij/cij.epsilon_mult - 
Hij1/cij1.epsilon_mult) ) - cij.Ex*cij.sigmax*dt/(cij.epsilon_mult*epsilon);

cij.Ey += (dt/dl*(-1/epsilon))*( (Hij/cij.epsilon_mult - 
Hi1j/ci1j.epsilon_mult)  - cij.Ex*cij.sigmay*dt/cij.epsilon_mult*epsilon);
}
void EMFieldPML::eval_cell_H(int i, int j){

    TECell& cij = field[i][j];
    TECell& cij1 = c(i,j+1);    
    TECell& ci1j = c(i+1,j);

    cij.Hzx += (dt/dl*(1/mu))*( - (ci1j.Ey/ci1j.mu_mult -
    cij.Ey/cij.mu_mult) - cij.Hzx*cij.sigmam_x*dt/cij.mu_mult*mu);

    cij.Hzy += (dt/dl*(1/mu))*( (cij1.Ex/cij1.mu_mult -
    cij.Ex/cij.mu_mult) - cij.Hzy*cij.sigmam_y*dt/cij.mu_mult*mu);
}

Demo

Comentario: Matemáticamente hablando, el PML es simplemente un dominio que tiene una permitividad y permeabilidad anisotrópica y compleja.  Aunque los PMLs son teóricamente no reflejantes, sí exhiben alguna reflexión debido a la discretización numérica: la malla. Para minimizar esta reflexión, se debe utilizar una malla en el PML que se ajuste (lo más perfecto posible) a las anisotropías en las propiedades del material.

Made with Slides.com