TF/SF源TM波在UPML条件下的2D-FDTD中的圆柱散射（三）

目录

1.题目要求

2.总体思路

3.无PML的2D-FDTD（TM波）

4.UPML的实现

5.TF/SF源

（1）在j=ja，j=jb的外部散射场设置Hx

（2）在j=ja，j=jb处总场设置Dz

（3）在i=ia，i=ib的外部散射场设置Hy

（4）设立Ezinc与Hxinc

6.圆柱散射

7.总结与演示

 参考文献

---

接前文，这部分阐述TF/SF在2D-FDTD区域的设置和圆柱的设置。

5.TF/SF源

总场散射场的原理图如下：

 具体概念就不阐述了。可以看到，需要修正的点就是在边界上的点（或边界外半个网格）。对于二维的FDTD，其总场与散射场的分界为x=ia，x=ib，y=ja，y=jb，如下图所示：

初始化代码：

% =======================================================================
% total/scattered field and Sinusoidal source
% =======================================================================

% total/scattered field boundaries
ia=12;
ib=IE-ia+1;
ja=12;
jb=JE-ja+1;

ez_inc_low_m1=0;
ez_inc_low_m2=0;
ez_inc_high_m1=0;
ez_inc_high_m2=0;

% input source
t0=20.0;                          % Center of the incident pulse
spread=6.0;                       % Width of the incident pulse
lamda=20e-6;
freq=c0/lamda*10^(-6);            % Frequency in MHz

在使用更新方程时，由于边界上的点（或差半个网格）中的更新方程涉及到了另一个场的变量，所以必须对其四个边界的点进行修正。

前文所述2D-FDTD中的得到TM波的方程用Yee网格表达为：

（1）在j=ja，j=jb的外部散射场设置Hx

代码：

% Incident Hx values
     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

（2）在j=ja，j=jb处总场设置Dz

由式 ，在j=ja，j=jb处的总场值受到ja-1/2或jb+1/2处的散射场影响，以j=ja为例，其更新系数左式化为：

代码：

% Incident DZ values
     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

（3）在i=ia，i=ib的外部散射场设置Hy

由式可知i=ia，i=ib的外部散射场边界的点受到总场的影响。同（1）理，可以将i=ia处的更新系数化作：

代码：

% Incident Hy values
     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

（4）设立Ezinc与Hxinc

    要完成上面公式的修正还需要两个源，它们相差半个时间步长和半个空间步长。由更新方程，易表达其关系。其代码为：

for j=2:JE
    ez_inc(j)= ez_inc(j)+0.5*( hx_inc(j-1)- hx_inc(j));
end

for j=1:JE-1
    hx_inc(j)=hx_inc(j)+0.5*( ez_inc(j)-ez_inc(j+1));
end

由于前文所提的固定条件c0*dt/ddx=0.5可知，波在一个时间步长里传播的距离为ddx/2。这表明，场穿过一个空间像元需要两个时间步长，由此得到波的吸收边界条件为：

为了实现上述边界条件，我们引入代码如下：

     % ABC for the incident buffer
     ez_inc(1)    =ez_inc_low_m2;
     ez_inc_low_m2=ez_inc_low_m1;
     ez_inc_low_m1=ez_inc(2);
        
     ez_inc(JE-1)  =ez_inc_high_m2;
     ez_inc_high_m2=ez_inc_high_m1;
     ez_inc_high_m1=ez_inc(JE-2);

将正弦波源加载至ez_inc(3)处：

      % insert source
      pulse=sin(2*pi*freq*10^6*dt*T);     % Sinusoidal source
      ez_inc(3)=pulse;         

6.圆柱散射

取(round(IE/2)，round(JE/2))处为圆柱中心，rad=5.0为半径，在xdist^2+ydist^2<=rad^2时，确定圆柱范围内的参数和辅助变量：

% =======================================================================
% cylinder parameters
% =======================================================================
rad=5.0;                          % radius of cylinder
icenter=round(IE/2);              % i-coordinate of cylinder's center
jcenter=round(JE/2);              % j-coordinate of cylinder's center
sigma=2.0;
epsilon=600.0;                    
mur=ones(IE,JE);
gb=zeros(IE,JE);
iz=zeros(IE,JE);
for i=1:IE
   for j=1:JE
      xdist=i-icenter;
      ydist=j-jcenter;
      dist=sqrt(xdist^2 + ydist^2);
      if dist <= rad 
         ga(i,j)=1./(epsilon+sigma*dt/epsz);
         gb(i,j)=sigma*dt/epsz;
      end
   end
end

在加入了圆柱后，其区别在于由D推导E的方程发生了变化，需要引入ga与gb变量辅助推导（之前的ga在真空中归一化后默认为1）。现推导新的更新方程与ga，gb的形式。

由前文可知，归一化后的麦克斯韦方程组在频域内表示为：

其中

则

1/jω向时域变换为对t积分，故上式变为：

Ez的更新方程变为：

 故而引入辅助变量:

方程改写为：

引入ga，gb：

代码：

 % Calculate the Ez field
     for j=2:JE-1                        
         for i=2:IE-1 
            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  

7.总结与演示

总结一下，总的算法框图为：

 演示：

UPML条件下的2D-FDTD的圆柱散射（TF/SF正弦源）

 参考文献：

[1] Sacks Z S , Kingsland D M . A perfectly matched anisotropic absorber for use as an absorbingboundary condition[J]. IEEE Transactions on Antennas & Propagation, 1995, 43(12):P.1460-1463.

[2] Electromagnetic Simulation Using the FDTD Method (IEEE2nd)-Sullivan

[3] 赵维. PML吸收边界条件下TM波的时域有限差分法分析[J]. 科技信息, 2009, 000(028):724-725.

[4] 李伟, 赵春晖, 徐娜. 基于MATLAB的时域有限差分法(FDTD)的研究[J]. 牡丹江大学学报, 2007, 016(007):91-92.

[5] J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” L Computational Phys., vol. 114, no. 2, pp. 185-200, Oct. 1994.