【科学代码】FDTD 求解三维带 UPML的电磁波传播问题

一、数学模型

所有的电磁现象都可以由 Maxwell 方程来进行描述，考虑如下频域 Maxwell 旋度方程

其中   是介电常数张量，  是磁导率张量，  是虚数单位。

这里我们考虑电磁波在真空   中的传播问题。假设真空中某一点处有一个强度随时间正弦变化的电场（称为激励源），其周围就会形成电磁波，并传播到无穷远处。但要从数值上模拟上述电磁波的传播问题，我们需要对其数学模型进行适当的修改，将其变为一个可以计算的模型。

首先，由于数值模拟只能模拟有界区域，我们需要先把求解区域限定为激励源周围的有界区域。不妨假设计算区域为一六面体区域，记为  ，而激励源就位于   的中心。

其次需要在   的边界   上设置合适的“边界条件”，保证电磁波可以畅通无阻地穿过  。但我们事先并不知道   上电场或磁场变化情况，因此我们无法设置通常意义下的边界条件。一个解决问题的思路是在有界区域边界外面人为增加一层可以把电磁波完全吸收掉的“材料”。当电磁波穿过   进入这层“材料”时，由于完全被吸收掉了，站在    里面看就好像传播出去一样的效果。

完美匹配层（Perfect Matched Layer，PML）就是这样一种人工的吸收边界层。本文主要讨论的是一种比较常用的 PML，称为单轴完美匹配层 （Uniaxial Perfect Matched Layer ，UPML）。下面具体介绍 UPML 层中介电常数和磁导率参数的具体选取办法。 

在计算区域的四周是 UPML ，UPML 厚度设为  ，且有， 、 、 、 、 、 。计算区域中的辐射源产生的外行波穿过 UPML 层的分界面，在 UPML 层中被吸收。 所有计算区域内和 UPML 层中的介质参数可以用以下的形式统一的表达：

其中， 、  和   是辅助参数，其具体形式如下，

其中， 、  和  分别是  、  和   方向上的电导率，  是光速，  和   分别是真空中的介电常数和磁导率， 、  和   定义如下：

其中,  是吸收常数，取值依赖于波的特点和入射角度， 、 、  在 UPML 层中电导率变化采用幂函数的形式，文献[1]表明，  左右时，可以达到较好的吸收效果。

当 UPML 采用以上参数时，通过 UPML 层的平面波在棱边区和角顶区的 UPML 分界面均没有反射。在 UPML 内部，电磁波会逐渐损耗掉。

将式（1.2）代入式（1.1）第一式得

即

引入中间变量

于是式（1.6）变为

应用频域到时域算子关系

将式（1.3）和式（1.9）代入式（1.8）得到

将式（1.3）代入式（1.7）得

即

应用式（1.9），上式的时域形式为

将式（1.2）代入式（1.1）第二式得

即

引入中间变量

于是式（1.15）变为

将式（1.3）和式（1.9）代入式（1.17）得到

将式（1.3）代入式（1.16）得

即

应用式（1.9），上式的时域形式为

二、FDTD 离散与程序实现

下面我们基于 FEALPy 讨论如何实现 FDTD （Finite Difference Time Domain，时域有限差分算法）方法模拟入射电磁波在区域   中以及在周围 PML 层的无反射传播过程。

计算域  ，UPML 层厚度  。假设把区域  、 和 三个方向都剖分   段，得到一个笛卡尔网格  ，相应的空间离散步长记为  。时间离散   段，时间步长记为  ，则相应的网比（或者时间稳定因子）为

进一步，一个波长划分的网格数量记为  ，则

激励源强加在   上，使用正弦波激励，第   个时刻   的入射电磁波取值为

下面我们讨论采用数组化编程实现 UPML 的算法。首先，我们把  、  、 、  和   作为程序的输入参数。

其次用 12 个三维数组来存储电场和磁场，其中电场定义在笛卡尔网格的边上，磁场定义在笛卡尔网格的面上。

• Bx、Hx ：定义在法向与  轴平行的面上，形状为 (NS+1, NS, NS)
• By、Hy ：定义在法向与  轴平行的面上，形状为 (NS, NS+1, NS)
• Bz、Hz ：定义在法向与  轴平行的面上，形状为 (NS, NS, NS+1)
• Dx、Ex ：定义在切向与  轴平行的边上，形状为 (NS, NS+1, NS+1)
• Dy、Ey ：定义在切向与  轴平行的边上，形状为 (NS+1, NS, NS+1)
• Dz、Ez ：定义在切向与  轴平行的边上，形状为 (NS+1, NS+1, NS)

注意 12 个三维数组第   轴对应   轴、第   轴对应   轴、第   轴对应   轴。

在 FEALPy 中，StructureHexMesh 中定义了一个 function 的方法，可以用来定义上面的 12 个三维数组，示例代码如下：

Bx0, By0, Bz0 = mesh.function(etype='face', dtype=np.float64)
Bx1, By1, Bz1 = mesh.function(etype='face', dtype=np.float64)
Hx0, Hy0, Hz0 = mesh.function(etype='face', dtype=np.float64)
Hx1, Hy1, Hz1 = mesh.function(etype='face', dtype=np.float64)
Dx0, Dy0, Dz0 = mesh.function(etype='edge', dtype=np.float64)
Dx1, Dy1, Dz1 = mesh.function(etype='edge', dtype=np.float64)
Ex0, Ey0, Ez0 = mesh.function(etype='edge', dtype=np.float64)
Ex1, Ey1, Ez1 = mesh.function(etype='edge', dtype=np.float64)

接下来我们讨论不同区域所定义的  、  和   如何计算，在 FEALPy 中，StructureQuadMesh 中定义了一个 interpolation 的方法，可以把一个已知函数插值到网格单元上或者边上，示例代码如下：

domain = [0, 1, 0, 1, 0, 1]  # 原来的区域
def sigma_x(p):
    x = p[..., 0]
    shape = p.shape[:-1]
    val = np.zeros(shape, dtype=np.float64)
    flag = x < 0
    val[flag] = sigma * ((0 - x[flag]) / delta) ** m
    flag = x > 1
    val[flag] = sigma * ((x[flag] - 1) / delta) ** m
    return val

def sigma_y(p):
    y = p[..., 1]
    shape = p.shape[:-1]
    val = np.zeros(shape, dtype=np.float64)
    flag = y < 0
    val[flag] = sigma * ((0 - y[flag]) / delta) ** m
    flag = y > 1
    val[flag] = sigma * ((y[flag] - 1) / delta) ** m
    return val

def sigma_z(p):
    z = p[..., 2]
    shape = p.shape[:-1]
    val = np.zeros(shape, dtype=np.float64)
    flag = z < 0
    val[flag] = sigma * ((0 - z[flag]) / delta) ** m
    flag = z > 1
    val[flag] = sigma * ((z[flag] - 1) / delta) ** m
    return val

domain = [0 - delta, 1 + delta, 0 - delta, 1 + delta, 0 - delta, 1 + delta]  # 增加了 PML 层的区域
mesh = StructureHexMesh(domain, nx=NS + 2 * NP, ny=NS + 2 * NP, nz=NS + 2 * NP)  # 建立结构网格对象

sx0 = mesh.interpolation(sigma_x, intertype='facex')
sy0 = mesh.interpolation(sigma_y, intertype='facex')
sz0 = mesh.interpolation(sigma_z, intertype='facex')

sx1 = mesh.interpolation(sigma_x, intertype='facey')
sy1 = mesh.interpolation(sigma_y, intertype='facey')
sz1 = mesh.interpolation(sigma_z, intertype='facey')

sx2 = mesh.interpolation(sigma_x, intertype='facez')
sy2 = mesh.interpolation(sigma_y, intertype='facez')
sz2 = mesh.interpolation(sigma_z, intertype='facez')

sx3 = mesh.interpolation(sigma_x, intertype='edgex')
sy3 = mesh.interpolation(sigma_y, intertype='edgex')
sz3 = mesh.interpolation(sigma_z, intertype='edgex')

sx4 = mesh.interpolation(sigma_x, intertype='edgey')
sy4 = mesh.interpolation(sigma_y, intertype='edgey')
sz4 = mesh.interpolation(sigma_z, intertype='edgey')

sx5 = mesh.interpolation(sigma_x, intertype='edgez')
sy5 = mesh.interpolation(sigma_y, intertype='edgez')
sz5 = mesh.interpolation(sigma_z, intertype='edgez')   

下面讨论每一个时间步，更新磁场和电场的 FDTD 方法。每个时间步，变量的更新顺序为

磁场更新方程

即

对应的离散代码为：

Bx1 = c1 * Bx0 - c2 * (Ez0[:, 1:, :] - Ez0[:, 0:-1, :] - Ey0[:, :, 1:] + Ey0[:, :, 0:-1])
By1 = c3 * By0 - c4 * (Ex0[:, :, 1:] - Ex0[:, :, 0:-1] - Ez0[1:, :, :] + Ez0[0:-1, :, :])
Bz1 = c5 * Bz0 - c6 * (Ey0[1:, :, :] - Ey0[0:-1, :, :] - Ex0[:, 1:, :] + Ex0[:, 0:-1, :])

Hx1 = c7 * Hx0 + c8 * Bx1 - c9 * Bx0
Hy1 = c10 * Hy0 + c11 * By1 - c12 * By0
Hz1 = c13 * Hz0 + c14 * Bz1 - c15 * Bz0

电场更新方程

即

对应的离散代码为：

Dx1[:, 1:-1, 1:-1] = c16 * Dx0[:, 1:-1, 1:-1] + c17 * (Hz1[:, 1:, 1:-1] - Hz1[:, 0:-1, 1:-1] - Hy1[:, 1:-1, 1:] + Hy1[:, 1:-1, 0:-1])
Dy1[1:-1, :, 1:-1] = c18 * Dy0[1:-1, :, 1:-1] + c19 * (Hx1[1:-1, :, 1:] - Hx1[1:-1, :, 0:-1] - Hz1[1:, :, 1:-1] + Hz1[0:-1, :, 1:-1])
Dz1[1:-1, 1:-1, :] = c20 * Dz0[1:-1, 1:-1, :] + c21 * (Hy1[1:, 1:-1, :] - Hy1[0:-1, 1:-1, :] - Hx1[1:-1, 1:, :] + Hx1[1:-1, 0:-1, :])

Ex1[:, 1:-1, 1:-1] = c22 * Ex0[:, 1:-1, 1:-1] + c23 * Dx1[:, 1:-1, 1:-1] - c24 * Dx0[:, 1:-1, 1:-1]
Ey1[1:-1, :, 1:-1] = c25 * Ey0[1:-1, :, 1:-1] + c26 * Dy1[1:-1, :, 1:-1] - c27 * Dy0[1:-1, :, 1:-1]
Ez1[1:-1, 1:-1, :] = c28 * Ez0[1:-1, 1:-1, :] + c29 * Dz1[1:-1, 1:-1, :] - c30 * Dz0[1:-1, 1:-1, :]

上述离散格式中  ，其中  ，当   取不同值时   的插值点位置不同。

• 当 时， 与 插值点位置相同；

• 当 时， 与 插值点位置相同；

• 当 时， 与 插值点位置相同；

• 当 时， 与 插值点位置相同；

• 当 时， 与 插值点位置相同；

• 当 时， 与 插值点位置相同。

实验结果如下：

全部代码见 https://github.com/weihuayi/fealpy/blob/master/example/EMWaveFDTDWithPML3d_example.py

参考文献

[1] 崔凡, 陈毅, 薛晗鹏, 等. 三角网格剖分时域有限元法的探地雷达正演模拟[J]. 地球物理学进展, 2022, 37(2): 797-809.

推荐阅读

编辑：黄琦