@Finite Difference Frequency Domain(FDFD) 算法和C++实现

Finite Difference Frequency Domain(FDFD) 算法和C++实现

中文网站上目前有不少关于Finite Different Time Domain(FDTD) 的介绍和实现，但是有关FDFD的介绍非常少。因两者的核心方法都用到了差分方法，区别仅在于差分对象分别在时域和频域上。最近正好做到一个使用FDFD模拟电磁场分布的项目，顺便记录下来。本文中对object boundary使用PEC边界条件，对space boundary使用U-PML边界条件，使用Sparse Matrix储存稀疏矩阵信息，使用BICGSTAB进行矩阵计算。

FDFD基本算法

我们从最基本的麦克斯韦方程组开始，一步一步推导出FDFD计算公式。首先，我们有频域中无源媒质里的麦克斯韦旋度方程：

$$\begin{gathered} \nabla \times \vec{E}=-j\omega \mu \vec{H}(1) \\ \nabla \times \vec{H}=j\omega ϵ\vec{E}(2) \end{gathered}$$
我们对$\vec{H}$标准化后，公式变成：

$$\begin{gathered} \nabla \times \vec{E}=k_0{\mu }_r\overrightarrow{\stackrel{~}{H}}(3) \\ \nabla \times \overrightarrow{\stackrel{~}{H}}=k_0{ϵ}_r\vec{E}(4) \end{gathered}$$
其中
$k_0=\omega /c_0$，

$\overrightarrow{\stackrel{~}{H}}=-j\sqrt{\frac{{ϵ}_0}{{\mu }_0}}\vec{H}$。

Perfectly Matched Layer

为了防止遇到外边界而产生的反射波对结果的影响，我们使用PML边界条件模拟无限大空间的波传播方式，即波逐渐被媒质吸收不产生反射波。我们使用张量$\overline{\stackrel{ˉ}{s}}$来代表PML介质，那么公式(3)(4)就变成：

$$\begin{gathered} \nabla \times \vec{E}=k_0{\mu }_r\overline{\stackrel{ˉ}{s}}\overrightarrow{\stackrel{~}{H}}(5) \\ \nabla \times \overrightarrow{\stackrel{~}{H}}=k_0{ϵ}_r\overline{\stackrel{ˉ}{s}}\vec{E}(6) \end{gathered}$$
其中：
$$\overline{\stackrel{ˉ}{s}}=\left[\begin{array}{ccc}\frac{s_y{s}_z}{s_x}& 0& 0\\ 0& \frac{s_x{s}_z}{s_y}& 0\\ 0& 0& \frac{s_x{s}_y}{s_z}\end{array}\right](7)$$

$$\begin{gathered} s_x(x)=s_{x0}[1+\frac{{\eta }_0}{j}{\sigma }_x^{\ast }(x)](8) \\ {s}_y(y)=s_{y0}[1+\frac{{\eta }_0}{j}{\sigma }_y^{\ast }(y)](9) \\ {s}_z(z)=s_{z0}[1+\frac{{\eta }_0}{j}{\sigma }_z^{\ast }(z)](10) \end{gathered}$$
其中
$s_{m0}(m)=1+s_{max}(\frac{m}{L}{)}^p$
${\sigma }_m^{\ast }(m)={\sigma }_{max}^{\ast }si{n}^2(\frac{\pi m}{2L})$

L是PML的层数，*是当前位置距离PML中心的层数。通常
$0⩽s_{max}⩽5$
${\sigma }_{max}^{\ast }\approx 1$
$0⩽p⩽5$

差分法

现在我们要通过差分法把连续的旋度公式离散化。我们知道，对于三个间距为$△$的连贯点$x_0,x_1,x_2$，它们函数值为$f_0,f_1,f_2$。那么对于$f_1$的导数我们可以近似地离散为：
$$\frac{\partial f_1}{\partial x}\approx \frac{f_2-f_0}{2△}(11)$$
同样的，对两个连贯点来说：
$$\frac{\partial f_{1.5}}{\partial x}\approx \frac{f_2-f_1}{△}(12)$$
在直角坐标系下，结合公式(7)(8)(9)(10)，我们可以把公式(5)(6)展开如下：
$$\begin{gathered} \frac{\partial E_z}{\partial y^{\prime }}-\frac{\partial E_y}{\partial z^{\prime }}={\mu }_r\frac{s_y{s}_z}{s_x}\widetilde{H_x}(13) \\ \frac{\partial E_x}{\partial z^{\prime }}-\frac{\partial E_z}{\partial x^{\prime }}={\mu }_r\frac{s_x{s}_z}{s_y}\widetilde{H_y}(14) \\ \frac{\partial E_y}{\partial x^{\prime }}-\frac{\partial E_x}{\partial y^{\prime }}={\mu }_r\frac{s_x{s}_y}{s_z}\widetilde{H_z}(15) \\ \frac{\partial \widetilde{H_z}}{\partial y^{\prime }}-\frac{\partial \widetilde{H_y}}{\partial z^{\prime }}={ϵ}_r\frac{s_y{s}_z}{s_x}{E}_x(16) \\ \frac{\partial \widetilde{H_x}}{\partial z^{\prime }}-\frac{\partial \widetilde{H_z}}{\partial x^{\prime }}={ϵ}_r\frac{s_x{s}_z}{s_y}{E}_y(17) \\ \frac{\partial \widetilde{H_y}}{\partial x^{\prime }}-\frac{\partial \widetilde{H_x}}{\partial y^{\prime }}={ϵ}_r\frac{s_x{s}_y}{s_z}{E}_z(18) \end{gathered}$$
其中($x^{\prime }=k_{0}x,y^{\prime }=k_{0}y,z^{\prime }=k_{0}z$)。
之后，把求解区域使用Yee-Grid离散化，如下图所示：

再利用差分方程离散化公式(13) - (18)得到：

其中：

使用E代替H，并加入激励电流源得到：

这便是我们期望得到的线性方程组Ax = b。此处引用WIKIPEDIA对无预处理BICGSTAB算法的介绍：

线性方程组求解

线性方程组的求解一般分为直接求解和迭代求解。对于维度较小，奇异性较弱的算子A，我们可以采用高斯消去进行直接求解；如果算子A维度较大，例如三维的情况，需要大量内存空间储存A，我们则使用数值迭代的方法求解x。此处使用Biconjugate gradient stabilized method(BICGSTAB) 算法求解。此处引用WIKIPEDIA对无预处理BICGSTAB算法的介绍：

结果

设置一个二维区域，激励源为处于正中的偶极子源，周边灰色区域为PML层。模拟结果如下：

[1]: R. C. Rumpf, “Simple Implementation Of Arbitrarily Shaped Total-Field/Scattered-field Regions In Finite Difference Frequency-domain,” Progress In Electromagnetics Research B, Vol. 36, 221-248, 2012
[2]: R. C. Rumpf, “Design And Optimization Of Nano-optical Elements By Coupling Fabrication To Optical Behavior”(2006). Electronic Theses And Dissertations, 2004-2019.1081
[3]: https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method