基于时域有限差分法的FDTD的计算电磁学算法（含Matlab代码）-YEE网格下的更新公式推导

参考书籍：The finite-difference time-domain method for electromagnetics with MATLAB simulations（国内翻译版本：MATLAB模拟的电磁学时域有限差分法）
代码推荐：The finite-difference time-domain method for electromagnetics with MATLAB simulations的附件代码
我最初也是基于这个代码学习的

FDTD算法：采用差分直接离散时域Maxwell方程，电磁场的求解基于时间步的迭代，无需存储全空间的电磁场信息，内存消耗较小，同时采用立方体网格和差分算法，网格形式和算法均十分简单，计算速度快，基于时域算法，特别适合“宽带问题”的求解。但是，简单的立方体方体网格带来的弊端就是模型拟合精度较低，对于含有精细结构的模型，计算精度较低，同时基于“微分方程”，计算区域需要设置截断。
详细对比参考：常用计算电磁学算法特性与电磁软件分析

1、从麦克斯韦开始的FDTD时域有限差分法

1.1 麦克斯韦方程

FDTD叫时域有限差分法，显然，其依赖的麦克斯韦方程也是时域的。麦克斯韦时域微分方程为：
$$\begin{array}{c}\nabla \times \mathbf{H}=\frac{\partial \mathbf{D}}{\partial t}+\mathbf{J}\\ \nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t}-\mathbf{M}\\ \nabla \cdot \mathbf{D}={\rho }_{\mathrm{e}}\\ \nabla \cdot \mathbf{B}={\rho }_m\end{array}$$
式中，E为电场强度(V/m);D为电位移(C/m);H为磁场强度(A/m);B为磁通量密度(Wb/m°);J为电流密度(A/m);M为磁流密度(V/m);${\rho }_e$为电荷密度(C/m);${\rho }_m$为磁荷密度(Wb/m)。

依稀记得当时老师说，麦克斯韦方程有其直观理解，分别是：
1. 变化的电场和电流会产生磁场
2. 变化的磁场和磁荷会产生电场（自然界无磁荷，一般是等效出来）
3. 电流源产生电场
4. 磁流源产生磁场

1.2 本构关系

本构关系对补充麦克斯韦方程和描述媒质的特性是必要的,本构关系对线性、各向同性和非色散媒质可以写成：
$$\begin{array}{rl}D& =\epsilon E\\ B& =\mu H\end{array}.$$
其中，$\epsilon$为媒质的介电常数;$\mu$为媒质的磁导率。在自由空间，有:
$$\begin{array}{rl}\epsilon =& {\epsilon }_0=8.854\times 10^{-12}F/m\\ \mu =& {\mu }_0=4\pi \times 10^{-7}H/m\end{array}$$
在常规的电磁学表述中，我们更多的使用相对介电常数。比如说耳熟能详的FR4板材，其相对介电常数大概是${\epsilon }_r=4.2$。 这就代表其实际的介电常数为${\epsilon }_{FR4}={\epsilon }_r{\epsilon }_0$。但是，还有一个重要参数和本构关系相关，那就是损耗角正切$tan\delta$。

对于FR4板材，一般认为其损耗角正切为$tan\delta =0.02$，根据微波工程1.3小节的公式：
$$ϵ={ϵ}^{\prime }-j{ϵ}^{\prime \prime }={ϵ}^{\prime }(1-j\tan \delta )={ϵ}_0{ϵ}_r(1-j\tan \delta )$$，其对应的介电常数应该是：
$${\epsilon }_{FR4}={\epsilon }_r(1-j\tan \delta ){\epsilon }_0=(4.2-j0.02){\epsilon }_0$$
其对应的相对介电常数为：4.2-j0.02

在进行FDTD的推导时，因为在 FDTD 的更新方程的过程中满足散度方程，所以只需要考虑两个旋度方程即可。麦克斯韦中的电流密度$\mathbf{J}$等于导体电流密度${\mathbf{J}}_{\mathbf{c}}$与施加电流密度${\mathbf{J}}_{\mathbf{i}}$之和，即：
$$\mathbf{J}={\mathbf{J}}_{\mathrm{c}}+{\mathbf{J}}_{\mathrm{i}}$$
对于磁流密度，也类似：
$$\mathbf{M}={\mathbf{M}}_{\mathrm{c}}+{\mathbf{M}}_{\mathrm{i}}$$
因此，对原来的麦克斯韦方程拆分一下，就是：
$$\nabla \times \mathbf{H}=\epsilon \frac{\partial \mathbf{E}}{\partial t}+{\sigma }^e\mathbf{E}+{\mathbf{J}}_{\mathbf{i}}$$
和：
$$\nabla \times \mathbf{E}=-\mu \frac{\partial \mathbf{H}}{\partial t}-{\sigma }^m\mathbf{H}-{\mathbf{M}}_{\mathbf{i}}$$

旋度的计算公式大家还记得不：
$$\begin{array}{rl}& \nabla \times \mathbf{F}(x,y,z)=\left|\begin{array}{ccc}\hat{\mathbf{i}}& \hat{\mathbf{j}}& \hat{\mathbf{k}}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ F_x& F_y& F_z\end{array}\right|\\ & =\left(\frac{\partial F_z}{\partial y}-\frac{\partial F_y}{\partial z}\right)\hat{\mathbf{i}}+\left(\frac{\partial F_x}{\partial z}-\frac{\partial F_z}{\partial x}\right)\hat{\mathbf{j}}+\left(\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}\right)\hat{\mathbf{k}}\end{array}$$
把麦克斯韦旋度方程按照三个方向x，y，z全部展开，就可以得到6个方程：
$$\begin{array}{c}\frac{\partial {\mathbf{E}}_x}{\partial t}=\frac{1}{{\epsilon }_x}(\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-{\sigma }_x^e{E}_x-J_{ix})\\ \frac{\partial E_y}{\partial t}=\frac{1}{{\epsilon }_y}(\frac{\partial H_x}{\partial z}-\frac{\partial H_z}{\partial x}-{\sigma }_y^e{E}_y-J_{iy})\\ \frac{\partial E_z}{\partial t}=\frac{1}{{\epsilon }_z}(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}-{\sigma }_z^e{E}_z-J_{iz})\\ \frac{\partial H_x}{\partial t}=\frac{1}{{\mu }_x}(\frac{\partial E_y}{\partial z}-\frac{\partial E_z}{\partial y}-{\sigma }_x^m{H}_x-M_{ix})\\ \frac{\partial H_y}{\partial t}=\frac{1}{{\mu }_y}(\frac{\partial {\mathbf{E}}_x}{\partial x}-\frac{\partial {\mathbf{E}}_x}{\partial \mathbf{z}}-{\mathbf{\sigma }}_y^{\mathfrak{m}}{H}_y-{\mathbf{M}}_{iy})\\ \frac{\partial H_z}{\partial t}=\frac{1}{{\mu }_z}(\frac{\partial {\mathbf{E}}_x}{\partial y}-\frac{\partial {\mathbf{E}}_y}{\partial x}-{\sigma }_z^{\mathfrak{m}}{H}_z-{\mathbf{M}}_{iz})\end{array}$$

2、空间差分与时间差分

2.1、非常简单的差分方程

FDTD是在离散网格中进行迭代的，上面的麦克斯韦公式有大量的求导计算，这该如何解决呢？答案是差分近似。大家学高数都学过导数的近似吧：
$$f^{{}^{\prime }}(x)=\underset{\Delta x\to 0}{\lim }\frac{f(x+\Delta x)-f(x)}{\Delta x}$$
如果$\Delta x$非常小，那么：
$$f^{{}^{\prime }}(x)\approx \frac{f(x+\Delta x)-f(x)}{\Delta x}$$
但是为了实现更高的精度，所以采用FDTD都会采用双向差分公式：
$$f^{{}^{\prime }}(x)\approx \frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$$
实际上，此处使用的是近似，也存在高阶的FDTD的算法，对于此近似考虑了更多项，精度会更高(参考“基于高阶时域有限差分法平面波及完全匹配层的研究”等)：
$$f^{\prime }(x)=\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}-\frac{(\Delta x^2)}{6}+...=\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}+O((\Delta x{)}^2)$$

2.2、差分方程的运用

在FDTD算法中，网格被剖分为YEE网格的形式，电场和磁场元胞差半个身位，其更新的时间步也是差$0.5\Delta t$：

具体来讲，实际的电场网格和磁场网格的位置是：
$$\begin{array}{r}{E}_x(i,j,k)\Rightarrow \left((i-0,5)\Delta x,(j-1)\Delta y,(k-1)\Delta z\right)\\ E_y(i,j,k)\Rightarrow \left((i-1)\Delta x,(j-0.5)\Delta y,(k-1)\Delta z\right)\\ E_z(i,j,k)\Rightarrow \left((i-1)\Delta x,(j-1)\Delta y,(k-0.5)\Delta z\right)\\ H_x(i,j,k)\Rightarrow \left((i-1)\Delta x,(j-0.5)\Delta y,(k-0.5)\Delta z\right)\\ H_y(i,j,k)\Rightarrow ((i-0.5)\Delta x,(j-1)\Delta y,(k-0.5)\Delta z)\\ H_z(i,j,k)\Rightarrow ((i-0.5)\Delta x,(j-0.5)\Delta y,(k-1)\Delta z)\end{array}$$

更新的时间步也是差$0.5\Delta t$：FDTD算法在离散的时间瞬间取样和计算场值,但是电场和磁场取样计算并不是在相同的时刻。对时间步$\Delta t$,电场E的取样时刻为:0,$\Delta t$,2$\Delta t$,3$\Delta t$，…,n$\Delta t$；而磁场H取样时刻为:0.5$\Delta t$,1.5$\Delta t$,2.5$\Delta t$,…(n+0.5)$\Delta t$。即电场取样在时间的整数步长时刻,而磁场取样时刻为半整数时间步时刻。它们之间的时间差为半个时间步。

因此，考虑一个上面得到的麦克斯韦的方程（以Ex方向为例）：
$$\frac{\partial E_x}{\partial t}=\frac{1}{{\epsilon }_x}\left(\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-{\sigma }_x^e{E}_x-J_{ir}\right)$$
观察其导数项，分别有时间的差分项$\frac{\partial E_x}{\partial t}$和空间的差分项$\frac{\partial H_z}{\partial y}$和$\frac{\partial H_y}{\partial z}$。

方程中的导数可以用中心差分来近似，此时$E_x^n(i,j,k)$的位置为中心差分公式的中心点，而时间上应以$(n+0.5)\Delta t$作为中心点（因为电场E的取样时刻为:0,$\Delta t$,2$\Delta t$,3$\Delta t$，…,n$\Delta t$，而$(n+0.5)\Delta t$差分后可以得到n和n+1，符合取样时刻）。因此，第一项$\frac{\partial E_x}{\partial t}$可以写成如下的差分形式：
$$E_x^{n+0.5}(i,j,k)=\frac{E_x^{n+1}(i,j,k)-E_x^n(i,j,k)}{\Delta t}$$
而空间的差分项$\frac{\partial H_z}{\partial y}$可以写成：
$$\frac{\partial H_z}{\partial y}=\frac{H_z^{n+\frac{1}{2}}(i,j,k)-H_z^{n+\frac{1}{2}}(i,j-1,k)}{\Delta y}$$

2.3、得到差分方程

把所有项都写成差分形式，就可以得到3D的FDTD更新方程：
$$\begin{array}{rl}{E}_x^{n+1}\left(i,j,k\right)& =C_{exe}(i,j,k)\times E_x^n(i,j,k)\\ & +C_{exhz}(i,j,k)\times (H_z^{n+\frac{1}{2}}(i,j,k)-H_z^{n+\frac{1}{2}}(i,j-1,k))\\ & +C_{\mathrm{exhy}}(i,j,k)\times (H_y^{n+\frac{1}{2}}(i,j,k)-H_y^{n+\frac{1}{2}}(i,j,k-1))\\ & +C_{exj}\left(i,j,k\right)\times J_{ix}^{n+\frac{1}{2}}(i,j,k)\end{array}$$
C开头的都是系数，为了书写方便，其实际的值为：
$$\begin{array}{c}{C}_{exe}(i,j,k)=\frac{2{\epsilon }_z(i,j,k)-\Delta t{\sigma }_z^e(i,j,k)}{2{\epsilon }_z(i,j,k)+\Delta t{\sigma }_z^e(i,j,k)}\\ C_{exhy}(i,j,k)=\frac{2\Delta t}{(2{\epsilon }_z(i,j,k)+\Delta t{\sigma }_z^e(i,j,k))\Delta x}\\ C_{exhy}(i,j,k)=-\frac{2\Delta t}{(2{\epsilon }_z(i,j,k)+\Delta t{\sigma }_z^e(i,j,k))\Delta y}\\ C_{exj}\left(i,j,k\right)=-\frac{2\Delta t}{2{\epsilon }_z(i,j,k)+\Delta t{\sigma }_z^e(i,j,k)}\end{array}$$
当然，这只是6个方程中的一个，更加详细的方程参考：
MATLAB模拟的电磁学时域有限差分法的1.3。看看对应的matlab代码是怎么写的（没有电流就可以省略Cexj）：

current_time  = current_time + dt/2;

Ex(1:nx,2:ny,2:nz) = Cexe(1:nx,2:ny,2:nz).*Ex(1:nx,2:ny,2:nz) ...
                     + Cexhz(1:nx,2:ny,2:nz).*...
                     (Hz(1:nx,2:ny,2:nz)-Hz(1:nx,1:ny-1,2:nz)) ...
                     + Cexhy(1:nx,2:ny,2:nz).*...
                     (Hy(1:nx,2:ny,2:nz)-Hy(1:nx,2:ny,1:nz-1));   

% General electric field updating coefficients
% Coeffiecients updating Ex
Cexe  =  (2*eps_r_x*eps_0 - dt*sigma_e_x) ...
    ./(2*eps_r_x*eps_0 + dt*sigma_e_x);
Cexhz =  (2*dt/dy)./(2*eps_r_x*eps_0 + dt*sigma_e_x);
Cexhy = -(2*dt/dz)./(2*eps_r_x*eps_0 + dt*sigma_e_x);