FDTD 二维圆柱散射

1、基于fdtd的PML的TM波散射和mur边界的TE波散射组员：樊家伟、黄登祥、袁一粟、 江亚男、陈永炜我们的任务1、没有吸收边界条件下的TM波同轴圆柱散射；2、PML吸收边界条件下的TM波同轴圆柱散射；3、Mur吸收边界条件下的TE波圆柱散射。重写Maxswell方程组00D1t H001t HE *rD=E二维Maxswell方程组标量方程001yxzHHtxy D001xzHty E *zrzD=E001yzHtx E 1、二维FDTD基本原理对时间和空间差分后，迭代公式对时间和空间差分后，迭代公式1/21/200001,1/ 2,1/ 2,1,1/ 2,1/ 2nnnnzzyynnxxt

2、i ji jijijxti ji jy DDHHHH11/21/2001,1,nnnnxxzzti ji ji ji jy HHEE11/21/2001,1,nnnnyyzzti ji jiji jx HHEE1/21/2,0.51/ 2,1/ 2,0.5,1/ 2,1/ 2nnnnzzyynnxxi ji jijiji ji jDDHHHH11/21/2,0.5,1,nnnnxxzzi ji ji ji jHHEE11/21/2,0.51,nnnnyyzzi ji jiji jHHEE完美匹配层（Perfectly matched Layer, PML）是由Berenger 提出的，使用最为

3、灵活、广泛的一种ABC（Absorbing Boundary Conditions）。其基本原理是：电磁波的反射量由两个介质的波阻抗决定：ABAB其中2、完全匹配层基本原理2.1 在X 方向上实现PML（仅仅保留与X 方向有关的 F* 、 F* ） *001yxzFzHHjxxy D *0zxFxj Hxcy E *0zyFyj HxcxE将 F* 、 F*带入到左式00( )1yxDzHHxjcjxyD100( )1Dzxxj Hcjy E00( )1Dzyxj HcjxE注意到： H x方向上的磁导率与 H y方向上的磁导率互成倒数。因此，满足了PML 的第二个条件。(1)(2)(3)00

4、( )1yxzHHcxjjyxD00( )( )1DzzzxxjjjDDD时域0( )zDzxtDD1/21/21/21/20,( )2nnnnzzzzDi ji ji ji jitDDDD1/21/200( )( )11,1,122nnDDzzititi ji jttDD对(1)式左边进行差分00( )1yxzHHcxjjyxD 1/21/2,3,20.51/ 2,1/ 2,1/ 2,1/ 2nnnnzzyynnxxi jgiii jgiiHijHijHi jHi jDD其中012tcx 0121/ 2Dgiiit 001/ 231/ 2DDitgiiit00( )1Dzyxj HcjxE1

5、11221111,3,20.51,2222nnnnyyzzijfiiijfiiEijEi jHH其中012tcx0112121/ 22Dfiiit0011/ 2123121/ 22DDitfiiit同理(2)式100( )1Dzxxj Hcjy E 11/2111,_22211,2nnxxnHxi ji jcurlefii Ii jHH1/21/2_,1nnzxcurlei jEi jE1/21/211,_22nnHxHxi ji jcurleII00.5ctx 012Dtifii其中(3)式引进辅助参数02Dtxn随着电磁波进入到 PML，该参数是增大的。 30.333_ixn ilengt

6、hpml 121giixn i 131xn igiixn ii=1,2, length_pml注意到参变量 i/length_pml的变化范围是从0 到1，而权值0.333 是保持稳定状态的最大值1,0.751,0.50,0.333 1fii2.2 在y 方向上实现PML000( )( )11yxDDzHHxyjcjjxyD1000( )( )11DDzxxyj Hcjjy E1000( )( )11DDzyxyjHcjjxE(1)(2)(3) 1/21/2,33,220.51/ 2,1/ 2,1/ 2,1/ 2nnzznnyynnxxi jgii gjji jgii gjjHijHijHi

7、jHi jDD000( )( )11yxDDzHHxyjcjjxyD(1)1000( )( )11DDzxxyj Hcjjy E 11/211111,3,2_2222211,2nnxxnHxi jfjji jfjjcurlefii Ii jHH00.5ctx 012Dtifii其中1/21/211,_22nnHxHxi ji jcurleII1/21/2_,1nnzxcurlei jEi jE(2)1000( )( )11DDzyxyjHcjjxE1/12111,3,222120.5_21110.5,22nnyynHyijfiiijfiicurlefiiIijHH其中1/21/2_,1nnxx

8、curlei jEi jE1/21/211,_22nnHyHyi ji jcurleII(3)源程序clear all；clc;%设置网格数IE = 101; % x 方向网格100,实际有101个点JE = 101; % y 方向网格多少%设置圆柱中心ic = round(IE / 2);jc = round(JE / 2);%总场区区域ia=15;ib=IE-ia+1;ja =15;jb=JE-ja+1;%初始化设置ddx = .01; %空间网格大小，x 方向和y 方向网格大小相同dt = ddx / 6e8; %时域网格大小epsz = 8.8e-12; %介电常数epsilon=30

9、; sigma=0.3;pi = 3.13159; %常数PIc=3e8;%初始化系统变量dz = zeros(IE,JE); %电场通量Dhx = dz; %x 方向磁场强度hy = dz; %y 方向磁场强度ihy = dz; %用于中间变量ihy，是用来计算磁场强度的，是curl_e 的积分ihx = dz; %用于中间变量ihy，是用来计算磁场强度的，是curl_e 的积分ga = ones(IE,JE); %电场强度与D 之间的关系矩阵，gb=zeros(IE,JE);iz=gb;real_pt=gb;imag_pt=gb;real_in=0;imag_in=0;amp=zeros(

10、JE,1);%由于做了归一化，同时在真空中，因此ga 为1%初始化变量ez_inc=zeros(JE,1);hx_inc=ez_inc;hy_inc=ez_inc;ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;radius = 15;epsilon=10;%越大衰减越大sigma=0.3;%越大衰减越大for j=ja:jb for i=ia:ib xdist=ic-i; ydist=jc-j; dist=sqrt(xdist*xdist+ydist*ydist); if(dist=radius) ga(i,

11、j)=1./(epsilon+(sigma*dt/epsz); gb(i,j)=sigma*dt/epsz; end endend%内圆radius1 = 10;epsilon1=1000;sigma1=10;for j=ja:jb for i=ia:ib xdist=ic-i; ydist=jc-j; dist=sqrt(xdist*xdist+ydist*ydist); if(dist=radius1) ga(i,j)=1./(epsilon1+(sigma1*dt/epsz); gb(i,j)=sigma1*dt/epsz; end endend%输入PML Cell 个数，即PML 有

12、多少个单元网格，在此，x 方向和y 方向上的PML 网格相同npml = input(Please input the number of PML Cell: );%x 方向用*i*表示，一方向用*j*表示%x 方向上的PML 参数设置for i = 1:npmlxnum = npml -i + 1; %从npml 到0 xd = npml;xxn = xnum / xd; %辅助变量xxn，从1 到0 xn = 0.33 * xxn3; %成立方衰减gi2(i) = 1 / ( 1 + xn );gi2(IE-i+1) = 1 / ( 1 + xn );gi3(i) = (1 - xn) /

13、 (1 + xn);gi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn );xxn = ( xnum - 0.5 ) / xd;if(xxn0) break;end 30.333_ixn ilengthpmlxn = 0.33 * xxn3;fil(i) = xn;fil(IE-i+1) = xn;fi2(i) = 1 / ( 1 + xn );fi2(IE-i+1) = 1 / ( 1 + xn );fi3(i) = ( 1-xn ) / ( 1 + xn );fi3(IE-i+1) = ( 1 - xn ) / ( 1 + xn );end;%y 方向上的PML 参数设

14、置for j = 1:npml+1xnum = npml - j + 1;xd = npml;xxn = xnum / xd;xn = 0.33 * xxn3;gj2(j) = 1 / ( 1 + xn );gj2(JE+1-j) = 1 / ( 1 + xn );gj3(j) = ( 1 - xn ) / ( 1 + xn );gj3(JE-j+1) = ( 1 - xn ) / ( 1 + xn );xxn = ( xnum - 0.5 ) / xd;if(xxn0) breakendxn = 0.33 * xxn3;fj1(j) = xn;fj1(JE-j+1) = xn;fj2(j)

15、= 1 / ( 1 + xn );fj2(JE-j+1) = 1 / ( 1 + xn );fj3(j) = ( 1 - xn ) / ( 1 + xn );fj3(JE-j+1) = ( 1 - xn ) / ( 1 + xn );end;%高斯脉冲变量设置t0 = 40;spread = 10;T = 0;%输入nsteps 必须为正整数nsteps = input(Please input the number of nsteps );%设置循环次数，从常数可以得到空间和时间上的传输情况for n = 1:nsteps T = T+1; for j=2:JE ez_inc(j)=ez_i

16、nc(j)+0.5*(hx_inc(j-1)-hx_inc(j); end real_in=real_in+cos(2*pi*frequency*T*dt)*ez_inc(ja); imag_in=imag_in-sin(2*pi*frequency*T*dt)*ez_inc(ja); %边界条件 ez_inc(1)=ez_inc_low_m2; ez_inc_low_m2=ez_inc_low_m1; ez_inc_low_m1=ez_inc(1); ez_inc(JE)=ez_inc_high_m2; ez_inc_high_m2=ez_inc_high_m1; ez_inc_high_m

17、1=ez_inc(JE-1); %计算D for j = 2:JE for i = 2:IE 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; 1/21/2,33,220.51/2,1/2,1/2,1/2nnzznnyynnxxi jgii gjji jgii gjjHijHijHi jHi jDD% 激励元 在此采用正弦波，频率为1.5G，可选为高斯脉冲% hard source%pulse_choice

18、= input(Enter 1 for Gaussian, 2 for sine wave pulse: );%if pulse_choice = 1 pulse = exp( -0.5 * ( t0-T ) / spread ) 2 ) ; %clear j; %pulse=1;%exp(-j*omiga*ddx*T/c);%dz(IE/2,JE/2) = pulse + dz(IE/2,JE/2);%else %pulse = sin( 2*pi*1500*1e6*dt*T); ez_inc(10) = pulse;%end;for i=ia:ib dz(i,ja)=dz(i,ja)+0.

19、5*hx_inc(ja-1); dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);end%计算E for j = 1:JE for i = 1:IE ez(i,j) = ga(i,j) * (dz(i,j)-iz(i,j); iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j); end; end;%Set the Ez edges to 0, as part of the PML for j=1:JE ez(1,j) = 0; ez(IE,j) = 0; end; for i = 1:IE ez(i,1) = 0; ez(i,JE) = 0; end; for j=

20、1:JE-1 hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1); end for j=1:JE for i=1:IE real_pt(i,j)=real_pt(i,j)+cos(2*pi*frequency*T*dt)*ez(i,j); imag_pt(i,j)=imag_pt(i,j)-sin(2*pi*frequency*T*dt)*ez(i,j); end end %计算Hx for j = 1:JE-1 for i = 1:IE curl_e = ez(i,j) - ez(i,j+1); ihx(i,j) = ihx(i,j) + fi1(i

21、) * curl_e; hx(i,j) = fj3(j) * hx(i,j) + fj2(j) * 0.5 * ( curl_e + ihx(i,j) ); end; end;%修正场区域 for i=ia:ib hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb); end %计算Hy for j = 1:JE for i = 1:IE -1 curl_e = ez(i+1,j) - ez(i,j); ihy(i,j) = ihy(i,j) + fj1(j) * curl_e; hy(i,j) = fi3

22、(i) * hy(i,j) + fi2(i) * 0.5 * ( curl_e + ihy(i,j) ); end; end; %修正场区域 for j=ja:jb hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)=hy(ib,j)+0.5*ez_inc(j); end %画图 figure(1); surfc(1:IE,1:JE,ez) %surfc(ia:ib,ja:jb,ez(ia:ib,ja:jb) axis(0 IE 0 JE -0.2 1.) pause(.01);end;amp_in=sqrt(real_in2+imag_in2);for

23、 j=ja:jb if(ga(ic,j)=0右侧区域同时存在入射波和反射波叠加，此区域中有：式中设 左行波右行波Engquist-Majda吸收边界条件n带入得定义微分算子 注：保留对x的导数n算子L做因式分解 左行波算子右行波算子Engquist-Majda吸收边界条件n将左行波算子作用在平面波上得若在 截断边界处设置条件 ，相当于使截断截面处右行波（反射波）成分等于零。n将算子具体表示代入上式得n作频域到时域的转换Engquist-Majda吸收边界条件nEngquist-Majda吸收边界条件 （截断边界位于区域左侧） （截断边界位于区域右侧）n通过波动方程的因子分解获得直角坐标系下FD

24、TD吸收边界的单向波动方程一阶和二阶近似吸收边界n算子中含有根式运算，限制了数值实现。做近似处理：n利用Taylor级数展开n重写左行波算子n保留级数第一项做近似这就是x=0平面作为区域左侧界面不产生反射波的一阶近似吸收边界条件。一阶和二阶近似吸收边界n保留级数至第二项这就是x=0平面作为区域左侧界面不产生反射波的二阶近似吸收边界条件。n二阶近似比一阶近似吸收边界条件有所改善，残留反射波小。nMur对表中的吸收边界条件引入了一种简单有效的差分数值算法，即对时间和空间的偏微分取二阶中心差分近似，将单向波方程离散化，便形成了Mur的一阶二阶吸收边界条件，其总体虚假反射在1%5%，能满足许多工程需求

25、。Mur 吸收边界条件n对于二维电磁场问题，Mur指出二阶近似吸收边界可降低为只含E,H分量的一阶导数，从而使数值计算简化。nTM波 二维直角坐标TM波将上式对t积分，并设初始场为0Mur 吸收边界条件n对于TE波，同理 Mur二阶近似吸收边界条件比一阶近似多出一项一阶导数。nMur吸收边界条件的FDTD离散式n先看TM波的一阶近似 Mur 吸收边界条件n将上式在 作离散 n利用下式线性插值关系消去tntji)2/1()2/1(及，ZEMur 吸收边界条件n将上式带入，再带入得n整理后得到TM波边界条件的一阶近似n对照Yee元胞图可知位于左截断边界上的 节点值是用区域内部节点值及前一时刻边界上

26、的节点值来表示，不涉及截断边界以为的场量。ZEMur 吸收边界条件n对于TM波的Mur吸收边界条件二阶形式，比一阶多一项 离散式 线性插值n整理后得吸收边界的离散式：Mur 吸收边界条件nTE波一阶二阶吸收边界条件可类似推导得到Mur 吸收边界条件n对于二维TM波情况，电磁场除了 还有 、 。由图可知，在用FDTD计算边界处的TM波元胞 时并不会涉及截断边界以外E或H的节点。只有 涉及截断边界外侧的H节点。因此只需给出边界处切向场分量 的吸收边界条件。n对TE波同理只需 的 边界。ZEHxHyZEzHMur 吸收边界条件n以上只讨论了x=0的左侧界面的吸收边界条件。对于其余三边有相似结果。Mur程序 网格划分% Define the