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

最新推荐文章于 2023-09-20 17:01:40 发布

miro生

最新推荐文章于 2023-09-20 17:01:40 发布

阅读量1.8k

点赞数 7

分类专栏：光电数值仿真文章标签： matlab 差分法

本文链接：https://blog.csdn.net/mirosheng/article/details/108142218

版权

光电数值仿真专栏收录该内容

2 篇文章 3 订阅

订阅专栏

4.UMPL的实现

UPML基于各向异性的材料属性(ε，σ)来描述吸收层，通过选择合适的介质的材料特性，使得各向异性介质和自由空间之间的平面界面完全无反射，即反射系数R为0。这种方法不需要修改麦克斯韦方程组，实现比较简单，特别是在频域中。

（1）UPML条件下的麦克斯韦旋度方程

对于UPML的引入部分，详见刚刚提到的文献中，这里简单概述一下。

由阻抗匹配条件

并考虑均为复对角张量的情况，可以改写为：

对于在xoy平面上的二维FDTD，Sz=1，UPML边界如下图所示：

将[S]代入麦克斯韦频域卷积方程组：

经电场归一化条件

得：

与前文一样，假定只有对角张量，利用各方向分量展开卷积方程，有：

此时UPML边界如下图所示：

（2）x方向的UPML条件下的更新方程

Dz：

引入辅助变量

Hy:

注意fi2、fi3和gi2、gi3的唯一差别是在于它们有半个网格差。

Hx:

其中Lx为x方向PML的厚度，x为其距PML内边界的距离。系数为0.33是因为这是在稳定条件下能取得的最大值。如果按照PPT上的公式推导的话，则此处系数应该为0.25。对于这里的x，在代码取值时要注意取nmpl-i+1，另外还有右侧的PML区域，其与左侧区域i相对应的点应在IE-i+1。如图所示：

（3）二维的UPML条件下的更新方程

在y方向的UPML条件下的更新方程过程类似于x方向，这里略去。其中，X方向的Hx类似于Y方向的Hy，Y方向的Hx类似于X方向的Hy。最后将二者合在一起，得到二维的更新方程为：

Dz:

Hy:

Hx:

图示：

（4）代码

设立UPML范围及辅助变量（以x向为例）：

% Define the PML for the x direction
for i=1:npml
    xnum=npml-i+1;
    xd=npml;                       
    xxn=xnum/xd;
    xn=0.33*xxn^3;                 
    gi2(i)=1.0/(1.0+xn);           % left side
    gi3(i)=(1.0-xn)/(1.0+xn);
    gi2(IE-i+1)=1.0/(1.0+xn);          % Right side
    gi3(IE-i+1)=(1.0-xn)/(1.0+xn);
    
    xxn=(xnum-0.5)/xd;
    xn=0.33*xxn^3;
    fil(i)=xn;
    fi2(i)=1.0/(1.0+xn);
    fi3(i)=(1.0-xn)/(1.0+xn);
    fi1(IE+1-i)=xn;
    fi2(IE+1-i)=1.0/(1.0+xn);
    fi3(IE+1-i)=(1.0-xn)/(1.0+xn);
end

主循环中的更新方程：

     % Calculate the Dz field
     for j=2:JE-1             
         for i=2:IE-1 
            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            
     
     % Calculate the Ez field
     for j=2:JE-1                        
         for i=2:IE-1 
            Ez(i,j)=ga(i,j)*Dz(i,j);
         end
     end  
     
     % Calculate the Hx field
     for j=1:JE-2                        
         for i=1:IE-1 
            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/mur(i,j)*(curl_e+ihx(i,j));
         end
     end 

      % Calculate the Hy field
     for j=1:JE-1                       
         for i=1:IE-2 
            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/mur(i,j)*(curl_e+ihy(i,j));
         end
     end

注：Ez的更新方程在后文加入圆柱后会有所变化。

至此已实现了UPML边界。TF/SF的设置、圆柱散射、总结与演示请见下一篇文章。