时域有限差分方法实现电磁场仿真-chap3

3.2完美匹配层
由Sacks等人总结得到的PML完全匹配层条件为
1、由背景媒质到PML层的阻抗满足关系
其值等于1满足我们事先规定的归一化条件；
2、垂直于边界方向上的相对介电常数和磁导率，是另一个方向 上的倒数。
电磁仿真与微波电路CAD课程
3.2完美匹配层
电磁仿真与微波电路CAD课程
3.2完美匹配层
3.2完美匹配层
对3.13c式采用相同的处理方法，可得
其中
这里需要注意的是，由于 H y 在网格中的位置位于i+ 2 ,因而围 绕该位置进行展开。
1
对3.13b式需要采用不同的处理方法，将其写成
电磁仿真与微波电路CAD课程
3.2完美匹配层
这里注意 1/ j 可以看成是对时间的积分算符，而 j 是对时间 的微分算符，电场对空间的差分写成
电磁仿真与微波电路CAD课程
Anyway, different suggestions, good ideas and any critique are welcome.
电磁仿真与微波电路CAD课程
电磁仿真与微波电路CAD课程
3.2完美匹配层
x方向上有
各参量的变化范围满足
电磁仿真与微波电路CAD课程
3.2完美匹配层
程序要点
clear ez(1:60,1:60)=0.0;hx(1:60,1:60)=0.0;hy(1:60,1:60)=0.0; ga(1:60,1:60)=1.0;dz(1:60,1:60)=0.0;ihx(1:60,1:60)=0.0;ihy(1:60,1:60)=0.0; gi2(1:60)=1.0;gi3(1:60)=1.0;gj2(1:60)=1.0;gj3(1:60)=1.0; fi1(1:60)=0.0;fi2(1:60)=1.0;fi3(1:60)=1.0; fj1(1:60)=0.0;fj2(1:60)=1.0;fj3(1:60)=1.0; ic=30; 初始值设置为主计算空间的取值 jc=30; ddx=0.01; epsz=8.85419*(1E-12); dt=ddx/(2*3*(1E8)); T=0; t0=40.0; spread=15.0; nsteps=500; npml=8;
计算z方向电场
计算x方向磁场
计算y方向磁场
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
if mod(n,3)==0 meshz(1:200,1:200,ez) zlim([0 1]) pause(0.5); end end
动画表现3维场
电磁仿真与微波电路CAD课程
3.2完美匹配层
时域有限差分方法（FDTD） 实现电磁仿真
孔祥鲲 副研究员
xkkong@nuaa.edu.cn 雷达成像与微波光子技术教育部重点实验室 电子信息工程学院
Chap.3二维电磁仿真问题求解 本章内容
3.1 二维问题中的时域有限差分法 3.2 完美匹配层
本章主要介绍二维FDTD方法仿真电磁问题。 第一节中，将介绍FDTD方如何在二维场中展开应用，并用点源 作为范例。
电磁仿真与微波电路CAD课程
3.2完美匹配层
for n=1:nsteps T=T+1 %calculate the dz field for j=2:60 for i=2:60 dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); end end %put a gaussian pulse in the middle pulse=sin(2*pi*1500*(1E6)*dt*T); dz(ic,jc)=pulse; %calculate the ez field for j=1:60 for i=1:60 ez(i,j)=ga(i,j)*dz(i,j); end end
3.1二维问题中的时域有限差分方法
for n=1:nsteps T=T+1; %calculate the Dx field 计算电位移矢量 for j=2:200 for i=2:200 dz(i,j)=dz(i,j)+0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)); end end %put a gaussian pulse in the middle 在场中心加源 pulse=exp(-0.5*((t0-T)/spread)^2); dz(ic,jc)=pulse;
TE模式
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
以TM模式为例，3.1式推导为
二维问题，将电场和磁场进行交错计算处理； 对时间和空间离散，作中心差分近似，得到
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
电场由周围的四个磁场共同计算 得到；
X方向磁场由相邻的两个电场计算 得到；
在主计算区域中， fi1 0, gi 2 gi3 1 因此，该计算式可以无缝 连接主计算空间和PML空间。
电磁仿真与微波电路CAD课程
3.2完美匹配层
到目前为止，我们仅仅考虑了x方向上的PML，同样的道理必须 在y方向上实施。代替3.13式，有
相同的方法，代替3.14式，可以获得
Y方向上 H y 需要类似于 H x 相同的处理，得到
电磁仿真与微波电路CAD课程
3.2完美匹配层
%edges to 0 %calculate the hy field for j=1:60 for j=1:60 ez(1,j)=0.0; for i=1:59 ez(60,j)=0.0; curl_e=ez(i+1,j)-ez(i,j); end ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; for i=1:60 hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j)); ez(i,1)=0.0; end ez(i,60)=0.0; end end %calculate the hx for j=1:59 for i=1:60 curl_e=ez(i,j)-ez(i,j+1); ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j)); end end
利用FDTD方法展开，可以写成
最终可以获得
电磁仿真与微波电路CAD课程
3.2完美匹配层
3.13b式的实现，通过以下一组方程
其中
在计算f和g参数的时候，没有必要真实的改变电导率，这里使 用辅助参量
它会随着PML深度的增加而增加
电磁仿真与微波电路CAD课程
3.2完美匹配层
则，该参数的计算可以写成
经验性稳定值 各参量的变化范围满足
电磁仿真与微波电路CAD课程
3.2完美匹配层
for i=1:npml xnum=npml-i; xd=npml; xxn=xnum/xd; xn=0.33*(xxn^3); gi2(i+1)=1.0/(1.0+xn); gi2(60-i)=1.0/(1.0+xn); gi3(i+1)=(1.0-xn)/(1.0+xn); gi3(60-i)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*(xxn^3); fi1(i+1)=xn; fi1(60-i)=xn; fi2(i+1)=1.0/(1.0+xn); fi2(60-i)=1.0/(1.0+xn); fi3(i+1)=(1.0-xn)/(1.0+xn); fi3(60-i)=(1.0-xn)/(1.0+xn); end for j=1:npml xnum=npml-j; xd=npml; xxn=xnum/xd; xn=0.33*(xxn^3); gj2(j+1)=1.0/(1.0+xn); gj2(60-j)=1.0/(1.0+xn); gj3(j+1)=(1.0-xn)/(1.0+xn); gj3(60-j)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*(xxn^3); fj1(j+1)=xn; fj1(60-j)=xn; fj2(j+1)=1.0/(1.0+xn); fj2(60-j)=1.0/(1.0+xn); fj3(j+1)=(1.0-xn)/(1.0+xn); fj3(60-j)=(1.0-xn)/(1.0+xn); end
电磁仿真与微波电路CAD课程
3.2完美匹配层
电磁仿真与微波电路CAD课程
3.2完美匹配层
采取的办法就是将介电常数和磁导率均设为复数，虚数部分代 表衰减程度。 对3.2式采用傅里叶变换，得到
将假设的介电常数和磁导率带入到方程组中，得到
电磁仿真与微波电路CAD课程
3.2完美匹配层
电磁仿真与微波电路CAD课程
第二节主要讨论如何利用吸收边界配合二维场展开计算。
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
从第二章中得到的归一化麦克斯韦方程出发
当我们处理三维仿真问题的时候，按照直角坐标系的分量处理 Hx , H y , Hz ； 办法，有6个参量 Ex , Ey , Ez ， 处理二维问题的时候，选取两类代表性的电磁传输模式进行研 究。 TM模式
Y方向磁场由相邻 的两个电场计算得 到；
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
采用第一章相同的操作
获得二维迭代公式
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
注意
Ez 和 Dz 的关系，采用了一维有耗媒质中的关系。 jc=100; 程序要点 ddx=0.01 ; clear epsz=8.85419*(1E-12); ez(1:200,1:200)=0.0;hx(1:200,1:200) dt=ddx/(2*3*(1E8)); =0.0;dx(1:200,1:200)=0.0;hy(1:200,1 T=0; :200)=0.0; t0=20.0; ga(1:200,1:200)=1.0; spread=6.0; dz(1:200,1:200)=0.0; nsteps=200; ic=100; 电磁仿真与微波电路CAD课程
电磁仿真与微波电路CAD课程
3.1二维问题中的时域有限差分方法
%calculate the ez field for j=1:200 for i=1:200 ez(i,j)=ga(i,j)*dz(i,j); end end %calculate the hx for j=1:199 for i=1:199 hx(i,j)=hx(i,j)+0.5*(ez(i,j)-ez(i,j+1)); end end %calculate the hy field for j=1:199 for i=1:199 hy(i,j)=hy(i,j)+0.5*(ez(i+1,j)-ez(i,j)); end end
随着电磁波向周围空间的传播，必然会遇到计算边界，如果
边界不处理，反射波将会在计算空间内产生。
程序显示二维计算边界产生的反射波
因而不能分析清楚什么是波源产生的波，什么是边界反射的波， 因此吸收边界伴随着FDTD发展的始终，因而人们发明很多方法 解决吸收边界问题。 Berenger发现的PML（perfectly matched layer ）是解决该 问题的最有效方法之一。
首先实现X方向上的PML吸收层，保留3.8式X方向上的分量
利用3.12给予的值，可得
满足PML条件2
电磁仿真与微波电路CAD课程Baidu Nhomakorabea
3.2完美匹配层
将3.13式用时域有限差分方法展开，3.13a左边可以写成
转换为时域的形式后差分展开
将上式带入到3.13a，对磁场进行空间差分，可以获得
其中
两个新参量
电磁仿真与微波电路CAD课程