FDTD二维圆柱散射.ppt

《FDTD二维圆柱散射.ppt》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、基于fdtd的PML的TM波散射和mur边界的TE波散射组员：樊家伟、黄登祥、袁一粟、江亚男、陈永炜我们的任务1、没有吸收边界条件下的TM波同轴圆柱散射；2、PML吸收边界条件下的TM波同轴圆柱散射；3、Mur吸收边界条件下的TE波圆柱散射。重写Maxswell方程组二维Maxswell方程组标量方程1、二维FDTD基本原理对时间和空间差分后，迭代公式完美匹配层（PerfectlymatchedLayer,PML）是由Berenger提出的，使用最为灵活、广泛的一种ABC（AbsorbingBoundaryCondit

2、ions）。其基本原理是：电磁波的反射量由两个介质的波阻抗决定：其中2、完全匹配层基本原理2.1在X方向上实现PML（仅仅保留与X方向有关的εF*、μF*）将εF*、μF*带入到左式注意到：Hx方向上的磁导率与Hy方向上的磁导率互成倒数。因此，满足了PML的第二个条件。(1)(2)(3)时域对(1)式左边进行差分其中其中同理(2)式其中(3)式引进辅助参数随着电磁波进入到PML，该参数是增大的。i=1,2,…,length_pml注意到参变量i/length_pml的变化范围是从0到1，而权值0.333是保持稳定状态的

3、最大值2.2在y方向上实现PML(1)(2)(3)(1)其中(2)其中(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;sig

4、ma=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

5、=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

6、+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/

7、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

8、);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(xxn<0)break;endxn=0.33*xxn^3;fil(i)=xn;fil(IE-i+1)=xn;fi2(i)=1/(1+xn);fi2(I