FDTD-二维圆柱散射

基于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E1、二维FDTD基本原理对时间和空间差分后，迭代公式1/21/200001,,1/2,1/2,1,1/2,1/2nnnnzzyynnxxtijijijijxtijijyDDHHHH11/21/2001,,,1,nnnnxxzztijijijijyHHEE11/21/2001,,1,,nnnnyyzztijijijijxHHEE1/21/2,,0.51/2,1/2,0.5,1/2,1/2nnnnzzyynnxxijijijijijijDDHHHH11/21/2,,0.5,1,nnnnxxzzijijijijHHEE11/21/2,,0.51,,nnnnyyzzijijijijHHEE完美匹配层（PerfectlymatchedLayer,PML）是由Berenger提出的，使用最为灵活、广泛的一种ABC（AbsorbingBoundaryConditions）。其基本原理是：电磁波的反射量由两个介质的波阻抗决定：ABAB其中2、完全匹配层基本原理2.1在X方向上实现PML（仅仅保留与X方向有关的εF*、μF*）*001yxzFzHHjxxyD*0zxFxjHxcyE*0zyFyjHxcxE将εF*、μF*带入到左式00()1yxDzHHxjcjxyD100()1DzxxjHcjyE00()1DzyxjHcjxE注意到：Hx方向上的磁导率与Hy方向上的磁导率互成倒数。因此，满足了PML的第二个条件。(1)(2)(3)00()1yxzHHcxjjyxD00()()1DzzzxxjjjDDD时域0()zDzxtDD1/21/21/21/20,,,,()2nnnnzzzzDijijijijitDDDD1/21/200()()11,1,122nnDDzzititijijttDD对(1)式左边进行差分00()1yxzHHcxjjyxD1/21/2,3,20.51/2,1/2,,1/2,1/2nnnnzzyynnxxijgiiijgiiHijHijHijHijDD其中012tcx0121/2Dgiiit001/231/2DDitgiiit00()1DzyxjHcjxE111221111,3,20.51,,2222nnnnyyzzijfiiijfiiEijEijHH其中012tcx0112121/22Dfiiit0011/2123121/22DDitfiiit同理(2)式100()1DzxxjHcjyE11/2111,,_22211,2nnxxnHxijijcurlefiiIijHH1/21/2_,,1nnzxcurleijEijE1/21/211,,_22nnHxHxijijcurleII00.5ctx012Dtifii其中(3)式引进辅助参数02Dtxn随着电磁波进入到PML，该参数是增大的。30.333_ixnilengthpml121giixni131xnigiixnii=1,2,…,length_pml注意到参变量i/length_pml的变化范围是从0到1，而权值0.333是保持稳定状态的最大值1,0.751,0.50,0.3331fii2.2在y方向上实现PML000()()11yxDDzHHxyjcjjxyD1000()()11DDzxxyjHcjjyE1000()()11DDzyxyjHcjjxE(1)(2)(3)1/21/2,33,220.51/2,1/2,,1/2,1/2nnzznnyynnxxijgiigjjijgiigjjHijHijHijHijDD000()()11yxDDzHHxyjcjjxyD(1)1000()()11DDzxxyjHcjjyE11/211111,3,2_2222211,2nnxxnHxijfjjijfjjcurlefiiIijHH00.5ctx012Dtifii其中1/21/211,,_22nnHxHxijijcurleII1/21/2_,,1nnzxcurleijEijE(2)1000()()11DDzyxyjHcjjxE1/12111,3,222120.5_21110.5,22nnyynHyijfiiijfiicurlefiiIijHH其中1/21/2_,,1nnxxcurleijEijE1/21/211,,_22nnHyHyijijcurleII(3)源程序clearall；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;%初始化系统变量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%初始化变量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;%越大衰减越大forj=ja:jbfori=ia:ibxdist=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;endendend%内圆radius1=10;epsilon1=1000;sigma1=10;forj=ja:jbfori=ia:ibxdist=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;endendend%输入PMLCell个数，即PML有多少个单元网格，在此，x方向和y方向上的PML网格相同npml=input('PleaseinputthenumberofPMLCell:');%x方向用*i*表示，一方向用*j*表示%x方向上的PML参数设置fori=1:npmlxnum=npml-i+1;%从npml到0xd=npml;xxn=xnum/xd;%辅助变量xxn，从1到0xn=0.33*xxn^3;%成立方衰减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30.333_ixnilengthpmlxn=0.33*xxn^3;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参数设置forj=1:npml+1xnum=npml-j+1;xd=npml;xxn=xnum/xd;xn=0.33*xxn^3;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*xxn^3;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('Pleaseinputthenumberofnsteps');