FDTD-二维圆柱散射市公开课一等奖百校联赛获奖课件.pptx

基于fdtdPMLTM波散射和mur边界TE波散射组员：樊家伟、黄登祥、袁一粟、江亚男、陈永炜第1页我们任务1、没有吸收边界条件下TM波同轴圆柱散射；2、PML吸收边界条件下TM波同轴圆柱散射；3、Mur吸收边界条件下TE波圆柱散射。第2页重写Maxswell方程组二维Maxswell方程组标量方程 1、二维FDTD基本原理第3页对时间和空间差分后，迭代公式对时间和空间差分后，迭代公式第4页第5页完美匹配层（Perfectly matched Layer,PML）是由Berenger 提出，使用最为灵活、广泛一个ABC（Absorbing Boundary Conditions）。其基本原理是：电磁波反射量由两个介质波阻抗决定：其中2、完全匹配层基本原理第6页2.1 在X 方向上实现PML（仅仅保留与X 方向相关 F*、F*）将 F*、F*带入到左式注意到：H x方向上磁导率与 H y方向上磁导率互成倒数。所以，满足了PML 第二个条件。(1)(2)(3)第7页时域对(1)式左边进行差分第8页其中第9页其中同理(2)式第10页其中(3)式第11页引进辅助参数伴随电磁波进入到 PML，该参数是增大。i=1,2,length_pml注意到参变量 i/length_pml改变范围是从0 到1，而权值0.333 是保持稳定状态最大值第12页2.2 在y 方向上实现PML(1)(2)(3)第13页(1)第14页其中(2)第15页其中(3)第16页源程序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;sigma=0.3;pi=3.13159;%常数PIc=3e8;第17页%初始化系统变量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(JE,1);%因为做了归一化，同时在真空中，所以ga 为1第18页%初始化变量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,j)=1./(epsilon+(sigma*dt/epsz);gb(i,j)=sigma*dt/epsz;end endend第19页%内圆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第20页%输入PML Cell 个数，即PML 有多少个单元网格，在此，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)/(1+xn);gi3(IE-i+1)=(1-xn)/(1+xn);xxn=(xnum-0.5)/xd;if(xxn0)break;end第21页xn=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 参数设置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;第22页if(xxn0)breakendxn=0.33*xxn3;fj1(j)=xn;fj1(JE-j+1)=xn;fj2(j)=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);%设置循环次数，从常数能够得到空间和时间上传输情况第23页for n=1:nsteps T=T+1;for j=2:JE ez_inc(j)=ez_inc(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_m1=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;第24页%激励元 在此采取正弦波，频率为1.5G，可选为高斯脉冲%hard source%pulse_choice=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.5*hx_inc(ja-1);dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);end第25页%计算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=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第26页%计算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)*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第27页%计算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(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第28页%画图 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 j=ja:jb if(ga(ic,j)=0右侧区域同时存在入射波和反射波叠加，此区域中有：式中设 左行波右行波第31页Engquist-Majda吸收边界条件n带入得定义微分算子 注：保留对x导数n算子L做因式分解 左行波算子右行波算子第32页Engquist-Majda吸收边界条件n将左行波算子作用在平面波上得若在 截断边界处设置条件 ，相当于使截断截面处右行波（反射波）成份等于零。n将算子详细表示代入上式得n作频域到时域转换第33页Engquist-Majda吸收边界条件nEngquist-Majda吸收边界条件 （截断边界位于区域左侧）（截断边界位于区域右侧）n经过波动方程因子分解取得直角坐标系下FDTD吸收边界单向波动方程第34页一阶和二阶近似吸收边界n算子中含有根式运算，限制了数值实现。做近似处理：n利用Taylor级数展开n重写左行波算子n保留级数第一项做近似这就是x=0平面作为区域左侧界面不产生反射波一阶近似吸收边界条件。第35页一阶和二阶近似吸收边界n保留级数至第二项这就是x=0平面作为区域左侧界面不产生反射波二阶近似吸收边界条件。n二阶近似比一阶近似吸收边界条件有所改进，残留反射波小。nMur对表中吸收边界条件引入了一个简单有效差分数值算法，即对时间和空间偏微分取二阶中心差分近似，将单向波方程离散化，便形成了Mur一阶二阶吸收边界条件，其总体虚假反射在1%5%，能满足许多工程需求。第36页Mur 吸收边界条件n对于二维电磁场问题，Mur指出二阶近似吸收边界可降低为只含E,H分量一阶导数，从而使数值计算简化。nTM波 二维直角坐标TM波将上式对t积分，并设初始场为0第37页Mur 吸收边界条件n对于TE波，同理 Mur二阶近似吸收边界条件比一阶近似多出一项一阶导数。nMur吸收边界条件FDTD离散式n先看TM波一阶近似 第38页Mur 吸收边界条件n将上式在 作离散 n利用下式线性插值关系消去第39页Mur 吸收边界条件n将上式带入，再带入得n整理后得到TM波边界条件一阶近似n对照Yee元胞图可知位于左截断边界上 节点值是用区域内部节点值及前一时刻边界上节点值来表示，不包括截断边界认为场量。第40页Mur 吸收边界条件n对于TM波Mur吸收边界条件二阶形式，比一阶多一项 离散式 线性插值n整理后得吸收边界离散式：第41页Mur 吸收边界条件nTE波一阶二阶吸收边界条件可类似推导得到第42页Mur 吸收边界条件n对于二维TM波情况，电磁场除了 还有 、。由图可知，在用FDTD计算边界处TM波元胞 时并不会包括截断边界以外E或H节点。只有 包括截断边界外侧H节点。所以只需给出边界处切向场分量 吸收边界条件。n对TE波同理只需 边界。第43页Mur 吸收边界条件n以上只讨论了x=0左侧界面吸收边界条件。对于其余三边有相同结果。第44页Mur程序 网格划分%Define the FDTD grids_smaller=1;t_smaller=sqrt(2);del_s=lambda/20/s_smallerdel_t=del_s/c/sqrt(2)/t_smallerdim_s=5;s_range=dim_s*lambda;s_range=ceil(s_range/del_s)*del_s;s_cells=s_range/del_s;cells=s_cells;nodes=cells+1;第45页主循环while(done=0),%Store previous values of E Ex_o=Ex;Ey_o=Ey;Hz_o=Hz;第46页Mur边界%Update boundary conditions(Ex,Ey)Ex(1,1:nodes)=Ex_o(2,1:nodes)+c_MUR*(Ex(2,1:nodes)-Ex(1,1:nodes);Ey(1:nodes,nodes)=Ey_o(1:nodes,nodes-1)+c_MUR*(Ey(1:nodes,nodes-1)-Ey(1:nodes,nodes);Ex(nodes,1:nodes)=Ex_o(nodes-1,1:nodes)+c_MUR*(Ex(nodes-1,1:nodes)-Ex(nodes,1:nodes);Ey(1:nodes,1)=Ey_o(1:nodes,2)+c_MUR*(Ey(1:nodes,2)-Ey(1:nodes,1);第47页作图%Plot responsefigure(5);mesh(abs(Ex+i*Ey);axis(0 nodes 0 nodes 0 1);end;第48页TE波MUR效果第49页5、仿真结果及分析TM波在圆柱体内幅度伴随相对介电常数不一样而展现不一样程度衰减第50页TM波在圆柱体内幅度伴随电导率不一样而展现不一样程度衰减第51页没有吸收边界时散射情况：第52页PML吸收边界条件下散射情况：第53页