FDTD笔记

FDTD

FDTD简介

时域有限差分方法基本方程

麦克斯韦方程组

本构关系

$\begin{aligned}D&=_\varepsilon E\\B&=_\mu H\end{aligned}$

把电流密度$J$写成导体电流密度$J_c$和施加电流密度$J_i$之和，同样对磁流密度也就行分解，重写麦克斯韦方程：

$\begin{aligned}\nabla\times H&=\varepsilon\frac{\partial E}{\partial t}+\sigma^{e}E+J_{i}\\\nabla\times E&=-\mu\frac{\partial H}{\partial t}-\sigma^{m}H-M_{i}\end{aligned}$

上面的方程仅设计电磁场$E$和$H$，而为涉及通量密度矢量$D$和$B$。公式中出现了四个本构关系参量，所以只能描述任意线性、各向同性的媒质。虽然FDTD公式中未出现散度方程，但是可以用来检测仿真结果。

在三维空间中把方程分解为3个标量方程。得到：

$\begin{gathered} \frac{\partial\boldsymbol{E}_{x}}{\partial t}= \frac{1}{\varepsilon_{x}}\Big(\frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}-\sigma_{x}^{e}E_{x}-J_{ix}\Big) \\ \frac{\partial E_{y}}{\partial t}= \frac{1}{\varepsilon_{y}}\Big(\frac{\partial H_{x}}{\partial z}-\frac{\partial H_{z}}{\partial x}-\sigma_{y}^{e}E_{y}-J_{iy}\Big) \\ \frac{\partial E_{z}}{\partial t}= \frac{1}{\varepsilon_{z}}\Big(\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}-\sigma_{z}^{e}E_{z}-J_{iz}\Big) \\ \frac{\partial H_{x}}{\partial t}= \frac{1}{\mu_{x}}\Big(\frac{\partial E_{y}}{\partial z}-\frac{\partial E_{z}}{\partial y}-\sigma_{x}^{m}H_{x}-M_{ix}\Big) \\ \frac{\partial H_{y}}{\partial t}= \frac{1}{\mu_{y}}\Big(\frac{\partial\boldsymbol{E}_{x}}{\partial x}-\frac{\partial\boldsymbol{E}_{x}}{\partial z}-\boldsymbol{\sigma}_{y}^{\mathfrak{m}}H_{y}-\boldsymbol{M}_{iy}\Big) \\ \frac{\partial H_{z}}{\partial t}= \frac{1}{\mu_{z}}\Big(\frac{\partial E_{x}}{\partial y}-\frac{\partial E_{y}}{\partial x}-\sigma_{z}^{m}H_{z}-M_{iz}\Big) \end{gathered}$

导数的近似差分

前向差分公式$f'(x)\approx\frac{f(x+\Delta x)-f(x)}{\Delta x}$

后向差分公式$f^{'}(x)\approx\frac{f(x)-f(x-\Delta x)}{\Delta x}$

中心差分公式$f^{\prime}(x){\approx}\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$

其中前向差分和后向差分公式具有一阶精度，但是中心差分公式具有二阶精度。

二阶导数的近似差分

中心差分公式$f^{\prime\prime}(x)\approx\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\left(\Delta x\right)^{2}}$

三维问题的FDTD更新方程

FDTD 技术将三维问题的几何结构分解为单元,以构成相应的网格。下图示出由（$N_x\times N_y\times N_z$）个Yee 单元构成的网格。使用矩形 Yee 单元，以单元的大小作为分辨率,用阶跃或阶梯形式来近似表面和内部几何结构。

下图示出了标志为$(i,j,k)$的 Yee 单元方式中离散空间各场分量位置。电场分量放置在 Yee 单元各棱的中间,平行于各;磁场分量放置在 Yee 单元各面的中心,平行于各面的法线。

令Yee单元的尺寸分别为$\Delta x，\Delta y，\Delta z$，可以得到坐标：

$\begin{aligned} & E_x(i, j, k) \Rightarrow((i-0.5) \Delta x,(j-1) \Delta y,(k-1) \Delta z), \\ & E_y(i, j, k) \Rightarrow((i-1) \Delta x,(j-0.5) \Delta y,(k-1) \Delta z), \\ & E_z(i, j, k) \Rightarrow((i-1) \Delta x,(j-1) \Delta y,(k-0.5) \Delta z), \\ & H_x(i, j, k) \Rightarrow((i-1) \Delta x,(j-0.5) \Delta y,(k-0.5) \Delta z), \\ & 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{aligned}$

FDTD算法在离散时间瞬间取样和计算场值，但是电场取样时刻为整数步长时刻，磁场取样为半整数时间步时刻。由于场分量和时间和空间都有关，故用上标表示时间标志。例如：$((i-1)\Delta x,(j-1）\Delta y,(k-0.5)\Delta z)$记为$E_{z}^{n}\left(i,j,k\right)$。其他以此类推。

物质参量(介电常数、磁导率、电导率和导磁性)分布在整个 FDTD 网格上，并且与场量相关,因此它们的标记与专场量相同。下图给出了介电常数和磁导率的标记。电导率分布在整个网格，其标记同介电常数。同样,导磁性的标记也与磁导率相同。对离散取样的场分量，在空间和时间上都具备适当的标记方式，这样麦克斯韦方程的旋度方程就可以标量方程的差分形式给出。

考虑方程：$\frac{\partial E_x}{\partial t}=\frac1{\varepsilon_x}\left(\frac{\partial H_x}{\partial y}-\frac{\partial H_y}{\partial z}-\sigma_x^eE_x-J_{i\pi}\right)$

可以写出差分形式

$\begin{aligned} \frac{E_{x}^{n+1}\left(i,j,k\right)-E_{x}^{n}\left(i,j,k\right)}{\Delta t}=& \frac{1}{\varepsilon_{x}(i,j,k)}\frac{H_{z}^{n+\frac{1}{2}}(i,j,k)-H_{z}^{n+\frac{1}{2}}(i,j-1,k)}{\Delta y}- \\ &\frac1{\varepsilon_{x}(i,j,k)}\frac{H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i,j,k-1)}{\Delta y}- \\ &\frac{\sigma_{x}^{e}(i,j,k)}{\varepsilon_{x}(i,j,k)}E_{x}^{n+\frac12}(i,j,k)-\frac{\int_{ix}^{n+\frac12}(i,j,k)}{\varepsilon_{x}(i,j,k)} \end{aligned}$

上式右边包含的电场是半整数时刻的，，这些电场可写为两个整数时刻的平均。$E_x^{n+\frac12}(i,j,k)=\frac{E_x^{n+1}(i,j,k)+E_x^n(i,j,k)}2$

由上面两式可以推出FDTD更新方程：

$\begin{aligned} E_{x}^{n+1}\left(i,j,k\right)& =\frac{2\varepsilon_{x}\left(i,j,k\right)-\Delta t\sigma_{x}^{\mathrm{e}}\left(i,j,k\right)}{2\varepsilon_{x}\left(i,j,k\right)+\Delta t\sigma_{x}^{\mathrm{e}}\left(i,j,k\right)}E_{x}^{n}\left(i,j,k\right)+ \\ &\frac{2\Delta t}{(2\varepsilon_x(i,j,k)+\Delta t\sigma_x^e(i,j,k))\Delta y}(H_z^{n+\frac12}(i,j,k)-H_z^{n+\frac12}(i,j-1,k))- \\ &\frac{2\Delta t}{(2\varepsilon_r(i,j,k)+\Delta t\sigma_r^e(i,j,k))\Delta z}(H_y^{n+\frac12}(i,j,k)-H_y^{n+\frac12}(i,j,k-1))- \\ &\frac{2\Delta t}{2\varepsilon_{x}\left(i,j,k\right)+\Delta t\sigma_{x}^{\epsilon}\left(i,j,k\right)}J_{\dot{x}}^{n+\frac12}\left(i,j,k\right) \end{aligned}$

上式表明了如何用用过去时刻的电场/磁场分量和激励分量来计算下一时间步的电场分量。同样可以推出其他分量的更新方程。

对于磁场考虑方程$\frac{\partial H_x}{\partial t}=\frac1{\mu_x}\left(\frac{\partial E_y}{\partial z}-\frac{\partial H_z}{\partial y}-\sigma_x^mH_x-M_{ix}\right)$

写出差分形式：

$\begin{aligned} \frac{H_x^{n+\frac12}(i,j,k)-H_x^{n-\frac12}(i,j,k)}{\Delta t}=& \frac1{\mu_x(i,j,k)}\frac{\boldsymbol{E}_y^n(i,j,k+1)-\boldsymbol{E}_y^n(i,j,k)}{\Delta z}- \\ &\frac1{\mu_x(i,j,k)}\frac{E_z^n(i,j+1,k)-E_z^n(i,j,k)}{\Delta y}- \\ &\frac{\sigma_{x}^{\mathfrak{m}}(i,j,k)}{\mu_{x}(i,j,k)}H_{x}^{n}(i,j,k)-\frac1{\mu_{x}(i,j,k)}M_{i_{x}}^{n}(i,j,k) \end{aligned}$

整理之后同样可以得到磁场的更新公式

$\begin{aligned} H_{x}^{n+\frac{1}{2}}\left(i,j,k\right)& =\frac{2\mu_x(i,j,k)-\Delta t\sigma_x^{\mathfrak{m}}(i,j,k)}{2\mu_x(i,j,k)-\Delta t\sigma_x^{\mathfrak{m}}(i,j,k)}H_x^{n-\frac12}(i,j,k) \\ &+\frac{2\Delta t}{(2\mu_x(i,j,k)+\Delta t\sigma_x^{\mathfrak{m}}(i,j,k))\Delta z}(E_y^n(i,j,k+1)-E_y^n(i,j,k)) \\ &-\frac{2\Delta t}{(2\mu_x(i,j,k)+\Delta t\sigma_x^m(i,j,k))\Delta y}(E_z^n(i,j+1,k)-E_z^n(i,j,k)) \\ &-\frac{2\Delta t}{2\mu_x(i,j,k)+\Delta t\sigma_x^n(i,j,k)}M_{ix}^n\left(i,j,k\right) \end{aligned}$

同样可以推出其他分量的更新方程。

最后，引入系数项，可以得到

$\begin{aligned} E_x^{n+1}(i,j,k)& =C_{exe}(i,j,k)\times E_x^n(i,j,k) \\ &+C_{exhz}(i,j,k)\times\left(H_z^{n+\frac12}(i,j,k)-H_z^{n+\frac12}(i,j-1,k)\right) \\ &+C_{exhy}(i,j,k)\times\left(H_y^{n+\frac12}(i,j,k)-H_y^{n+\frac12}(i,j,k-1)\right) \\ &+C_{exj}(i,j,k)\times J_{ix}^{n+\frac12}(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_{exe}(i,j,k)&& =\frac{2\varepsilon_x(i,j,k)-\Delta t\sigma_x^e(i,j,k)}{2\varepsilon_x(i,j,k)+\Delta t\sigma_x^e(i,j,k)}, \\ &C_{exhz}(i,j,k)&& =\frac{2\Delta t}{\left(2\varepsilon_x(i,j,k)+\Delta t\sigma_x^e(i,j,k)\right)\Delta y}, \\ &C_{exhy}(i,j,k)&& =-\frac{2\Delta t}{\left(2\varepsilon_x(i,j,k)+\Delta t\sigma_x^e(i,j,k)\right)\Delta z}, \\ &C_{exj}(i,j,k)&& =-\frac{2\Delta t}{2\varepsilon_x(i,j,k)+\Delta t\sigma_x^e(i,j,k)}. \end{aligned}$

$\begin{aligned} E_{y}^{n+1}(i,j,k)& =C_{eye}(i,j,k)\times E_{y}^{n}(i,j,k) \\ &+C_{eyhx}(i,j,k)\times\left(H_x^{n+\frac12}(i,j,k)-H_x^{n+\frac12}(i,j,k-1)\right) \\ &+C_{eyhz}(i,j,k)\times\left(H_z^{n+\frac12}(i,j,k)-H_z^{n+\frac12}(i-1,j,k)\right) \\ &+C_{eyj}(i,j,k)\times J_{iy}^{n+\frac12}(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_\text{еуе}{ ( i , j , k )}&& =\frac{2\varepsilon_y(i,j,k)-\Delta t\sigma_y^e(i,j,k)}{2\varepsilon_y(i,j,k)+\Delta t\sigma_y^e(i,j,k)}, \\ &C_{\textit{еуh}x}(i,j,k)&& =\frac{2\Delta t}{\left(2\varepsilon_y(i,j,k)+\Delta t\sigma_y^e(i,j,k)\right)\Delta z}, \\ &C_{\textit{еу}hz}(i,j,k)&& =-\frac{2\Delta t}{\left(2\varepsilon_y(i,j,k)+\Delta t\sigma_y^e(i,j,k)\right)\Delta x}, \\ &C_{eyj}(i,j,k)&& =-\frac{2\Delta t}{2\varepsilon_y(i,j,k)+\Delta t\sigma_y^e(i,j,k)}. \end{aligned}$

$\begin{aligned} E_z^{n+1}\left(i,j,k\right)& =C_{eze}(i,j,k)\times E_z^n(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_x^{n+\frac12}(i,j,k)-H_x^{n+\frac12}(i,j-1,k)\right) \\ &+C_{ezj}(i,j,k)\times J_{iz}^{n+\frac12}(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_{eze}(i,j,k)&& =\frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)}, \\ &C_{ezhy}(i,j,k)&& =\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)\right)\Delta x}, \\ &C_{ezhx}(i,j,k)&& =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)\right)\Delta y}, \\ &C_{ezj}(i,j,k)&& =-\frac{2\Delta t}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)}. \end{aligned}$

$\begin{aligned} H_x^{n+\frac12}(i,j,k) =&C_{hxh}(i,j,k)\times H_x^{n-\frac12}(i,j,k) \\ &+C_{hxey}(i,j,k)\times\left(E_{y}^n(i,j,k+1)-E_{y}^n(i,j,k)\right) \\ &+C_{hxez}(i,j,k)\times\left(E_z^n(i,j+1,k)-E_z^n(i,j,k)\right) \\ &+C_{hxm}(i,j,k)\times M_{ix}^n(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_{hxh}(i,j,k) =\frac{2\mu_x(i,j,k)-\Delta t\sigma_x^m(i,j,k)}{2\mu_x(i,j,k)+\Delta t\sigma_x^m(i,j,k)}, \\ &C_{hхеу} ( i , j , k ) =\frac{2\Delta t}{\left(2\mu_x(i,j,k)+\Delta t\sigma_x^m(i,j,k)\right)\Delta z}, \\ &C_{hxez}(i,j,k) =-\frac{2\Delta t}{\left(2\mu_x(i,j,k)+\Delta t\sigma_x^m(i,j,k)\right)\Delta y}, \\ &C_{hxm}(i,j,k) =-\frac{2\Delta t}{2\mu_x(i,j,k)+\Delta t\sigma_x^m(i,j,k)}. \end{aligned}$

$\begin{aligned} H_{y}^{n+\frac12}(i,j,k)& =C_{hyh}(i,j,k)\times H_y^{n-\frac12}(i,j,k)+C_{hyez}(i,j,k) \\ &\times\left(E_z^n(i+1,j,k)-E_z^n(i,j,k)\right)+C_{hyex}(i,j,k) \\ &\times\left(E_x^n(i,j,k+1)-E_x^n(i,j,k)\right)+C_{hym}(i,j,k)\times M_{iy}^n(i,j,k) \end{aligned}$

式中

$\begin{aligned} &C_{hyh}(i,j,k) =\frac{2\mu_y(i,j,k)-\Delta t\sigma_y^m(i,j,k)}{2\mu_y(i,j,k)+\Delta t\sigma_y^m(i,j,k)}, \\ &C_{\textit{hye}z}(i,j,k) =\frac{2\Delta t}{\left(2\mu_y(i,j,k)+\Delta t\sigma_y^m(i,j,k)\right)\Delta x}, \\ &C_{\textit{bуех}}(i,j,k) =-\frac{2\Delta t}{\left(2\mu_y(i,j,k)+\Delta t\sigma_y^m(i,j,k)\right)\Delta z}, \\ &C_{hym} ( i , j , k ) =-\frac{2\Delta t}{2\mu_y(i,j,k)+\Delta t\sigma_y^{m}(i,j,k)}. \end{aligned}$

$\begin{aligned} H_z^{n+\frac12}(i,j,k)& =C_{hzh}(i,j,k)\times H_z^{n-\frac12}(i,j,k) \\ &+C_{hzex}(i,j,k)\times\left(E_x^n(i,j+1,k)-E_x^n(i,j,k)\right) \\ &+C_{hzey}(i,j,k)\times\left(E_y^n(i+1,j,k)-E_y^n(i,j,k)\right) \\ &+C_{hzm}(i,j,k)\times M_{iz}^n(i,j,k), \end{aligned}$

式中

$\begin{gathered} C_{hzh}(i,j,k) =\frac{2\mu_y(i,j,k)-\Delta t\sigma_z^m(i,j,k)}{2\mu_z(i,j,k)+\Delta t\sigma_z^m(i,j,k)}, \\ C_{hzex}(i,j,k) =\frac{2\Delta t}{\left(2\mu_z(i,j,k)+\Delta t\sigma_z^m(i,j,k)\right)\Delta y}, \\ C_{hzey}(i,j,k) =-\frac{2\Delta t}{\left(2\mu_z(i,j,k)+\Delta t\sigma_z^m(i,j,k)\right)\Delta x}, \\ C_{hzm}(i,j,k) =-\frac{2\Delta t}{2\mu_z(i,j,k)+\Delta t\sigma_z^m(i,j,k)}. \end{gathered}$

应该注意的是,在所有系数项中,第一、二个系数下标与更新的场量相关,而对三个下标系数,第三个下标与所乘的场或源相关,对四下标系数,第三个和第四个下标与所乘的场量相关。

推导出FDTD的更新公式之后，就可以构造出时间步进算法。流程图如下：

二维问题的FDTD迭代方程

在二维情况下，场仅在一维方向分布时，可以得到简化的二维方程：

$\begin{aligned} &\frac{\partial E_x}{\partial t}= \frac{1}{\varepsilon_{x}}\big(\frac{\partial H_{z}}{\partial y}-\sigma_{x}^{\prime}E_{x}-J_{ix}\big) \\ &\frac{\partial E_y}{\partial t}=\frac1{\varepsilon_y}\left(-\frac{\partial H_z}{\partial x}-\sigma_y^eE_y-J_{iy}\right) \\ & \frac{\partial E_{z}}{\partial t}=\frac{1}{\varepsilon_{z}}\Big(\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}-\sigma_{z}^{e}E_{z}-J_{iz}\Big) \\ &\frac{\partial H_x}{\partial t}= \frac{1}{\mu_{x}}\Big(-\frac{\partial\boldsymbol{E}_{z}}{\partial y}-\sigma_{x}^{m}H_{x}-\boldsymbol{M}_{ix}\Big) \\ & \frac{\partial H_{y}}{\partial t}=\frac{1}{\mu_{y}}\Big(\frac{\partial E_{z}}{\partial x}-\sigma_{y}^{m}H_{y}-M_{iy}\Big) \\ &\frac{\partial H_{z}}{\partial t}=\frac{1}{\mu_{z}}\left(\frac{\partial E_{z}}{\partial y}-\frac{\partial E_{y}}{\partial x}-\sigma_{z}^{\mathfrak{m}}H_{z}-M_{iz}\right) \end{aligned}$

可将上面的方程根据场量对$z$方向的关系分为两组，即TMz模式和TEz模式。可以对这两类问题分别求解，问题的全解可以由这两类解的和来得到。

TEz模式

FDTD更新方程：

$\begin{aligned} E_{x}^{n+1}(i,j)& =C_{exe}(i,j)\times E_{x}^{n}(i,j)+C_{exhz}(i,j)\times\left(H_{z}^{n+\frac{1}{2}}(i,j)-H_{z}^{n+\frac{1}{2}}(i,j-1)\right) \\ &+C_{exj}(i,j)\times J_{ix}^{n+\frac{1}{2}}(i,j), \end{aligned}$

式中

$\begin{aligned} &C_{exe}(i,j) =\frac{2\varepsilon_{x}(i,j)-\Delta t\sigma_{x}^{e}(i,j)}{2\varepsilon_{x}(i,j)+\Delta t\sigma_{x}^{e}(i,j)}, \\ &C_{exhz}(i,j) =\frac{2\Delta t}{\left(2\varepsilon_x(i,j)+\Delta t\sigma_x^e(i,j)\right)\Delta y}, \\ &C_{exj}(i,j) =-\frac{2\Delta t}{2\varepsilon_{x}(i,j)+\Delta t\sigma_{x}^{e}(i,j)}. \end{aligned}$

$\begin{aligned} E_{y}^{n+1}(i,j)& =C_{eye}(i,j)\times E_{y}^{n}(i,j)+C_{eyhz}(i,j)\times\left(H_{z}^{n+\frac12}(i,j)-H_{z}^{n+\frac12}(i-1,j)\right) \\ &+C_{eyj}(i,j)\times J_{iy}^{n+\frac{1}{2}}(i,j), \end{aligned}$

式中

$\begin{aligned} &C_{eye}(i,j) =\frac{2\varepsilon_y(i,j)-\Delta t\sigma_y^e(i,j)}{2\varepsilon_y(i,j)+\Delta t\sigma_y^e(i,j)}, \\ &C_{eyhz}(i,j) =-\frac{2\Delta t}{\left(2\varepsilon_y(i,j)+\Delta t\sigma_y^e(i,j)\right)\Delta x}, \\ &C_{eyj}(i,j) =-\frac{2\Delta t}{2\varepsilon_{y}(i,j)+\Delta t\sigma_{\nu}^{e}(i,j)}. \end{aligned}$

$\begin{aligned} H_{z}^{n+\frac12}(i,j)& =C_{hzh}(i,j)\times H_{z}^{n-\frac12}(i,j)+C_{bzex}(i,j)\times\left(E_{x}^{n}(i,j+1)-E_{x}^{n}(i,j)\right) \\ &+C_{hzey}(i,j)\times\left(E_{y}^{n}(i+1,j)-E_{y}^{n}(i,j)\right)+C_{hzm}(i,j)\times M_{iz}^{m}(i,j), \end{aligned}$

式中

$\begin{aligned} &C_{hzh}(i,j) =\frac{2\mu_z(i,j)-\Delta t\sigma_z^m(i,j)}{2\mu_z(i,j)+\Delta t\sigma_z^m(i,j)}, \\ &C_{hzex}(i,j) =\frac{2\Delta t}{\left(2\mu_z(i,j)+\Delta t\sigma_z^m(i,j)\right)\Delta y}, \\ &C_{hzey}(i,j) =-\frac{2\Delta t}{\left(2\mu_z(i,j)+\Delta t\sigma_z^m(i,j)\right)\Delta x}, \\ &C_{hzm}(i,j) =-\frac{2\Delta t}{2\mu_z(i,j)+\Delta t\sigma_z^m(i,j)}. \end{aligned}$

TMz模式

$\begin{aligned} E_{z}^{n+1}(i,j)& =C_{eze}(i,j)\times E_{z}^{n}(i,j)+C_{ezhy}(i,j)\times\left(H_{y}^{n+\frac{1}{2}}(i,j)-H_{y}^{n+\frac{1}{2}}(i-1,j)\right) \\ &+C_{ezhx}(i,j)\times\left(H_{x}^{n+\frac{1}{2}}(i,j)-H_{x}^{n+\frac{1}{2}}(i,j-1)\right)+C_{ezj}(i,j)\times J_{iz}^{n+\frac{1}{2}}(i,j), \end{aligned}$

式中

$\begin{aligned} &C_{eze}(i,j)&& =\frac{2\varepsilon_{z}(i,j)-\Delta t\sigma_{z}^{e}(i,j)}{2\varepsilon_{z}(i,j)+\Delta t\sigma_{z}^{e}(i,j)}, \\ &C_{ezhy}(i,j)&& =\frac{2\Delta t}{\left(2\varepsilon_z(i,j)+\Delta t\sigma_z^e(i,j)\right)\Delta x}, \\ &C_{ezhx}(i,j)&& =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j)+\Delta t\sigma_z^e(i,j)\right)\Delta y}, \\ &C_{ezj}(i,j)&& =-\frac{2\Delta t}{2\varepsilon_z(i,j)+\Delta t\sigma_z^e(i,j)}. \end{aligned}$

$\begin{aligned} H_{x}^{n+\frac{1}{2}}(i,j)& =C_{h\mathrm{x}h}(i,j)\times H_x^{n-\frac12}(i,j)+C_{h\mathrm{x}ez}(i,j)\times\left(E_z^n(i,j+1)-E_z^n(i,j)\right) \\ &+C_{hxm}(i,j)\times M_{ix}^{n}(i,j), \end{aligned}$

式中

$\begin{aligned} &C_{hxh}(i,j) =\frac{2\mu_x(i,j)-\Delta t\sigma_x^m(i,j)}{2\mu_x(i,j)+\Delta t\sigma_x^m(i,j)}, \\ &C_{hxez}(i,j) =-\frac{2\Delta t}{\left(2\mu_x(i,j)+\Delta t\sigma_x^m(i,j)\right)\Delta y}, \\ &C_{hxm}(i,j) =-\frac{2\Delta t}{2\mu_{x}(i,j)+\Delta t\sigma_{x}^{m}(i,j)}. \end{aligned}$

一维FDTD问题的更新方程

一维情况下，场分布与两个坐标变量无关。此时麦克斯韦方程可以重写为：

$\begin{gathered} \frac{\partial E_{x}}{\partial t}= \frac1{\varepsilon_{x}}\big(-\sigma_{x}^{e}E_{x}-J_{ix}\big), (1.39\mathbf{a}) \\ \frac{\partial E_{y}}{\partial t}= \frac{1}{\varepsilon_{y}}\biggl(-\frac{\partial H_{z}}{\partial x}-\sigma_{y}^{e}E_{y}-J_{iy}\biggr), (1.39\mathbf{b}) \\ \frac{\partial E_z}{\partial t}= \frac{1}{\varepsilon_{z}}\biggl(\frac{\partial H_{y}}{\partial x}-\sigma_{z}^{e}E_{z}-J_{iz}\biggr), (1.39\mathfrak{c}) \\ \frac{\partial H_{x}}{\partial t}= \frac1{\mu_{x}}\left(-\sigma_{x}^{m}H_{x}-M_{ix}\right), (1.39\mathrm{d}) \\ \frac{\partial H_{y}}{\partial t}= \frac{1}{\mu_{y}}\left(\frac{\partial E_{z}}{\partial x}-\sigma_{y}^{m}H_{y}-M_{iy}\right), (1.39\mathbf{e}) \\ \frac{\partial H_{z}}{\partial t}= \frac{1}{\mu_{z}}\left(-\frac{\partial E_{y}}{\partial x}-\sigma_{z}^{m}H_{z}-M_{iz}\right). (1.39\mathrm{f}) \end{gathered}$

上式第一个和第四个方程不包含对空间的导数的项，因此这两式子不代表传播场。而其他方程对于$x$轴属于横向场，因此，在一维的情况下存在TEM波，其传播如同平面波的传播。

分析一维场的时候采用和二维场同样的方法把波分解为两类分离的类型。因为方程(1.39b)和方程(1.39f)仅含$E_z$和$H_x$两项,它们与方程(1.39c)方程(1.39e)是去耦合的。后者只含有E,和H两项。对方程(1.39b)和方程(1.39f),场量的空间位置如下图所示。

根据中心差分方法可得到FDTD的更新方程：

$E_y^{n+1}(i)=C_{eye}(i)\times E_y^n(i)+C_{eyhz}(i)\times\left(H_z^{n+\frac{1}{2}}(i)-H_z^{n+\frac{1}{2}}(i-1)\right)+C_{eyj}(i)\times J_{iy}^{n+\frac{1}{2}}(i),$

式中

$\begin{aligned} &C_{eye}(i) =\frac{2\varepsilon_{y}(i)-\Delta t\sigma_{y}^{e}(i)}{2\varepsilon_{y}(i)+\Delta t\sigma_{y}^{e}(i)}, \\ &C_{eyhz}(i) =-\frac{2\Delta t}{\left(2\varepsilon_y(i)+\Delta t\sigma_y^e(i)\right)\Delta x}, \\ &C_{eyj}(i) =-\frac{2\Delta t}{2\varepsilon_y(i)+\Delta t\sigma_y^e(i)}, \end{aligned}$

$H_z^{n+\frac{1}{2}}(i)=C_{hzh}(i)\times H_z^{n-\frac{1}{2}}(i)+C_{hzey}(i)\times\left(E_y^n(i+1)-E_y^n(i)\right)+C_{hzm}(i)\times M_{iz}^n(i),\quad(1.41)$

式中

$\begin{aligned} &C_{hzh}(i) =\frac{2\mu_{z}(i)-\Delta t\sigma_{z}^{m}(i)}{2\mu_{z}(i)+\Delta t\sigma_{z}^{m}(i)}, \\ &C_{hzey}(i) =-\frac{2\Delta t}{\left(2\mu_z(i)+\Delta t\sigma_z^m(i)\right)\Delta x}, \\ &C_{hzm}(i) =-\frac{2\Delta t}{2\mu_z(i)+\Delta t\sigma_z^m(i)}. \end{aligned}$

对于方程(1.39c)和(1.39e)，场量的空间位置如图：

FDTD更新方程：

$E_z^{n+1}(i)=C_{eze}(i)\times E_z^n(i)+C_{ezhy}(i)\times\left(H_y^{n+\frac{1}{2}}(i)-H_y^{n+\frac{1}{2}}(i-1)\right)+C_{ezj}(i)\times J_{iz}^{n+\frac{1}{2}}(i),$

其中

$\begin{aligned} &C_{eze}(i) =\frac{2\varepsilon_{z}(i)-\Delta t\sigma_{z}^{e}(i)}{2\varepsilon_{z}(i)+\Delta t\sigma_{z}^{e}(i)}, \\ &C_{ezhy}(i) =\frac{2\Delta t}{\left(2\varepsilon_z(i)+\Delta t\sigma_z^e(i)\right)\Delta x}, \\ &C_{ezj}(i) =-\frac{2\Delta t}{2\varepsilon_z(i)+\Delta t\sigma_z^e(i)}, \end{aligned}$

$H_y^{n+\frac{1}{2}}(i)=C_{hyh}(i)\times H_y^{n-\frac{1}{2}}(i)+C_{hyez}(i)\times\left(E_z^n(i+1)-E_z^n(i)\right)+C_{hym}(i)\times M_{iy}^n(i),$

式中

$\begin{aligned} &C_{hyh}(i) =\frac{2\mu_{y}(i)-\Delta t\sigma_{y}^{m}(i)}{2\mu_{y}(i)+\Delta t\sigma_{y}^{m}(i)}, \\ &C_{hyez}(i) =\frac{2\Delta t}{\left(2\mu_y(i)+\Delta t\sigma_y^m(i)\right)\Delta x}, \\ &C_{hym}(i) =-\frac{2\Delta t}{2\mu_{y}(i)+\Delta t\sigma_{y}^{m}(i)}. \end{aligned}$

数值稳定性和色散

在FDTD中，取样的周期应遵循一定的规则，确保解的稳定性。

时域算法中的稳定性

对于一个简单的波方程$\frac{\partial u(x,t)}{\partial t}+\frac{\partial u(x,t)}{\partial x}=0,\quad u(x,t=0)=u_0(x)$，其解析解为$u(x,t)=u_0(x-t)$。离散化后得到时域数值算法。

$\frac{u_i^{n+1}-u_i^{n-1}}{2\Delta t}+\frac{u_{i-1}^n-u_{i+1}^n}{2\Delta x}=0$整理后得到$u_i^{n+1}=u_i^{n-1}+\lambda(u_{i+1}^n-u_{i-1}^n)=0,\quad\lambda=\frac{\Delta t}{\Delta x}$

下图是一系列表示时间和空间的网格，假定对$u$,在时间标记为 $n$ 和$(n-1)$的时刻是已知的。为简单起见,假定所有这些值都为零。还假定,在时间标注为$n$ 时,空间下标为时$i$,有一小的误差量,用$\varepsilon$来表示。当时间向前推进时，$u$在$(n+1)$的值,用上面的递推式可以由低两行的相应值算出。

因此，上述递推式的稳定是有条件的。

FDTD方法的CFL稳定条件

FDTD方法的数值稳定由CFL条件确定。要求时间增量$\Delta t$相对于空间网格小于一特定值

$\Delta t\leqslant\frac1{c\sqrt{\frac1{\left(\Delta x\right)^2}+\frac1{\left(\Delta y\right)^2}+\frac1{\left(\Delta z\right)^2}}}$其中，$c$为自由空间光速。

对立方体网格有，$\Delta x=\Delta y=\Delta z$，则CFL条件变为$\Delta t\leqslant\frac{\Delta x}{c\sqrt{3}}$。

一维情况下，CFL条件简化为$\Delta t\leqslant\Delta x/c$。表明，在一个时间步长内，不容许波行进距离超过一个网格尺寸。

稳定性问题的存在表明，当FDTD迭代进行时，在问题空间可能激励出虚假的收敛场。

数值色散

用FDTD数值方法得到的相速与实际的相速之间的差别称为数值色散。

解析解的色散关系$k_x^2=\omega^2\mu_0\varepsilon_0=\left(\frac\omega c\right)^2$

在一维情况下，离散后可以推的色散关系$\left[\frac1{c\Delta t}\sin\left(\frac{\omega\Delta t}2\right)\right]^2=\left[\frac1{\Delta x}\sin\left(\frac{k_x\Delta x}2\right)\right]^2$

两者是不一样的，这意味着存在与实际解的偏差，这就是由数值解带来的误差。

二维色散关系是$\left[\frac1{c\Delta t}\sin\left(\frac{\omega\Delta t}2\right)\right]^2=\left[\frac1{\Delta x}\sin\left(\frac{k_x\Delta x}2\right)\right]^2+\left[\frac1{\Delta y}\sin\left(\frac{k_y\Delta y}2\right)\right]^2$

三维色散关系是

在Yee网格中创建目标

目标的定义

在FDTD算法中,开始 FDTD 时间步进迭代之前,第一步应该创建问题。问题的创建需要运行两个子程序:一是定义问题的子程序;二是初始化问题空间和 FDTD参量的子程序。问题的定义过程包括创建数据结构和存储与要解的电磁问题相关的构成信息。其中信息部分包括:在问题空间中的目标几何形状;构成材料的类型;材料的电磁特性;激励源的类型和波形;由 FDTD 所得仿真结果的类型;其他定义 FDTD算法的参量。如时间步数、问题边界的类型。初始化处理将结构数据嵌入到 FDTD媒质网格中,并构造和初始化数据。这包括 FDTD更新参量数组、场的系数数组和在 FDTD 迭代中及选代后用到的数据。简而言之,它为 FDTD计算和后处理准备数据结构。

定义问题空间参量

在问题空间中定义目标

媒质近似

近似策略：假定，每一种媒质位于一个网格单元的中心。这些网格与FDTD网格之间有一定的位置偏移，称为媒质网格。

点留驻方法：粗略近似加权体积平均。Yee单元中的每一媒质单元划分为8个子网格。

可以得到$\varepsilon_z(i,j,k)=\frac{\varepsilon_1+\varepsilon_2+\varepsilon_3+\varepsilon_4+\varepsilon_5+\varepsilon_6+\varepsilon_7+\varepsilon_8}8$

切向和法向分量的子网格平均方案

计算媒质分量平行于两种媒质的界面的等效介电常数$\begin{aligned}\varepsilon_1A_1+\varepsilon_2A_2&=\varepsilon_{\mathrm{eff}}\times(A_1+A_2)\Rightarrow\varepsilon_{\mathrm{eff}}\\&=\varepsilon_{\mathrm{eff}}(i,j,k)=\frac{\varepsilon_1A_1+\varepsilon_2A_2}{A_1+A_2}\end{aligned}$

计算媒质分量垂直于两种媒质的界面的等效介电常数$\varepsilon_{eff}=\varepsilon_{z}(i,j,k)=\frac{\varepsilon_{1}\varepsilon_{2}(l_{1}+l_{2})}{l_{1}\varepsilon_{2}+l_{2}\varepsilon_{1}}$

可以证明，其他物质参量仍满足这种关系，例如等效磁导率$\mu$小电导率$\sigma^e$和小电导率$\sigma^m$

定义目标

定义零厚度PEC目标

任意环绕着并与PEC平板共线的电导分量都可以指定为PEC电导率，以便模拟似零厚度的PEC平板。

创建媒质网格

改善8个网络平均

有源和无源集总参数电路

FDTD中集总参数元件的更新公式

电压源

节点$(i,j,k)$和$(i,j,k+1)$之间极化方向为$z$正向的电压源的FDTD更新方程

$\begin{aligned} E_{z}^{n+1}\left(i,j,k\right)& =C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac{1}{2}}(i,j,k)-H_{y}^{n+\frac{1}{2}}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac{1}{2}}(i,j,k)-H_{x}^{n+\frac{1}{2}}(i,j-1,k)\right) \\ &+C_{ezs}(i,j,k)\times V_{s}^{n+\frac{1}{2}}(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_{eze}(i,j,k)&&= \begin{aligned}\frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)-\frac{R_s\Delta x\Delta y}{R_s\Delta x\Delta y}}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}},\end{aligned} \\ &C_{ezhy}(i,j,k)&&=\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)+\frac{\Delta t\Delta z}{R_{s}\Delta x\Delta y}\right)\Delta x}, \\ &C_{ezhx}(i,j,k)&& =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}\right)\Delta y}, \\ &C_{ezs}(i,j,k)&& =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}\right)(R_s\Delta x\Delta y)}, \end{aligned}$

硬激励电压源

无内阻电压源

$C_{eze}(i,j,k)=-1,\quad C_{ezhy}(i,j,k)=0,\quad C_{ezhx}(i,j,k)=0,\quad C_{ezs}(i,j,k)=-\frac{2}{\Delta z}.$

电流源

$\begin{aligned} E_{z}^{n+1}\left(i,j,k\right)& =C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac{1}{2}}(i,j,k)-H_{y}^{n+\frac{1}{2}}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac{1}{2}}(i,j,k)-H_{x}^{n+\frac{1}{2}}(i,j-1,k)\right) \\ &+C_{ezs}(i,j,k)\times I_{s}^{n+\frac{1}{2}}(i,j,k), \end{aligned}$

式中

$\begin{aligned} &C_{eze}(i,j,k) =\frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)-\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}}, \\ &C_{ezhy}(i,j,k) =\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}\right)\Delta x}, \\ &C_{ezhx}(i,j,k) =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}\right)\Delta y}, \\ &C_{ezs}(i,j,k) =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{\Delta t\Delta z}{R_s\Delta x\Delta y}\right)\Delta x\Delta y}, \end{aligned}$

电阻的FDTD建模

$\begin{aligned} E_{z}^{+1}(i,j,k)& =C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left[H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i-1,j,k)\right] \\ &+C_{ezhz}\left(i,j,k\right)\times\left[H_{z}^{n+\frac12}\left(i,j,k\right)-H_{z}^{n+\frac12}\left(i,j-1,k\right)\right] \end{aligned}$

式中

$\begin{gathered} C_{\epsilon\epsilon\epsilon}(i,j,k)=\frac{2\varepsilon_{\epsilon}(i,j,k)-\Delta t\sigma_{\epsilon}^{\epsilon}(i,j,k)-\frac{\Delta t\Delta z}{R\Delta x\Delta y}}{2\varepsilon_{\epsilon}(i,j,k)+\Delta t\sigma_{\epsilon}^{*}(i,j,k)+\frac{\Delta t\Delta z}{R\Delta x\Delta y}} \\ C_{\epsilon thy}\left(i,j,k\right)=\frac{2\Delta t}{\left[2\varepsilon_{\epsilon}\left(i,j,k\right)+\Delta t\sigma_{\epsilon}^{\epsilon}\left(i,j,k\right)+\frac{\Delta t\Delta z}{R\Delta x\Delta y}\right]\Delta\dot{x}} \\ C_{\epsilon*h*}(i,j,k)=\frac{2\Delta t}{\left[2\varepsilon_{\epsilon}(i,j,k)+\Delta t\sigma_{\epsilon}^{\epsilon}(i,j,k)+\frac{\Delta t\Delta x}{R\Delta x\Delta y}\right]\Delta y} \end{gathered}$

电容的FDTD建模

$\begin{aligned} E_{z}^{n+1}(i,j,k)& =C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac{1}{2}}(i,j,k)-H_{y}^{n+\frac{1}{2}}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac{1}{2}}(i,j,k)-H_{x}^{n+\frac{1}{2}}(i,j-1,k)\right), \end{aligned}$

式中

$\begin{aligned} &C_{eze}(i,j,k) =\frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)+\frac{2C\Delta z}{\Delta x\Delta y}}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{2C\Delta z}{\Delta x\Delta y}}, \\ &C_{ezhy}(i,j,k) =\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{2C\Delta z}{\Delta x\Delta y}\right)\Delta x}, \\ &C_{ezhx}(i,j,k) =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)+\frac{2C\Delta z}{\Delta x\Delta y}\right)\Delta y}. \end{aligned}$

电感的FDTD建模

$\begin{aligned} E_{z}^{n+1}(i,j,k)& =C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k) \\ &+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac{1}{2}}(i,j,k)-H_{x}^{n+\frac{1}{2}}(i,j-1,k)\right) \\ &+C_{ezj}(i,j,k)\times J_{iz}^{n+\frac{1}{2}}(i,j,k), \end{aligned}$

式中

$J_{iz}^{n+\frac{1}{2}}(i,j,k)=J_{iz}^{n-\frac{1}{2}}(i,j,k)+\frac{\Delta t\Delta z}{L\Delta x\Delta y}E_{z}^{n}(i,j,k).$

$\begin{gathered} C_{eze}(i,j,k)= \frac{2\varepsilon_z(i,j,k)-\Delta t\sigma_z^e(i,j,k)}{2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)}, \\ C_{ezhy}(i,j,k)= =\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)\right)\Delta x}, \\ C_{ezhx}(i,j,k) =-\frac{2\Delta t}{\left(2\varepsilon_z(i,j,k)+\Delta t\sigma_z^e(i,j,k)\right)\Delta y}, \\ C_{ezj}(i,j,k) =-\frac{2\Delta t}{2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)}. \end{gathered}$

位于表明或体积内的集总参数元件

电阻$R^{'}=R\times\frac{(\text{rie}-\text{ris}+1)\times(\text{rje}-\text{rjs}+1)}{\text{rke}-\text{rks}}$

电感$L^{'}=L\times\frac{(\mathrm{ie}-\mathrm{is}+1)\times(\mathrm{je}-\mathrm{js}+1)}{\mathrm{ke}-\mathrm{ks}}$

电容$C^{\prime}=C\times\frac{\mathrm{ke}-\mathrm{ks}}{(\mathrm{ie}-\mathrm{is}+1)\times(\mathrm{je}-\mathrm{js}+1)}$

二极管的FDTD模拟

正向放置的二极管

$A\mathrm{e}^{\mathrm{Br}}+x+C=0$

式中

$\begin{aligned} &C=C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k)+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i-1,j,k)\right. \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac12}(i,j,k)-H_{x}^{n+\frac12}(i,j-1,k)\right)+C_{ezd}(i,j,k), \\ &C_{eze}(i,j,k)=-\frac{2\varepsilon_{z}(i,j,k)-\Delta t\sigma_{z}^{e}(i,j,k)}{2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)}, \\ &C_{ezhy}(i,j,k)=-\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)\right)\Delta x}, \\ &C_{ezhx}(i,j,k)=\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)\right)\Delta y}, \\ &C_{ezd}(i,j,k)=-\frac{2\Delta tI_{0}}{2\varepsilon_{z}(i,j,k)\Delta\mathrm{x}\Delta y+\Delta t\sigma_{z}^{e}(i,j,k)\Delta\mathrm{x}\Delta y}. \end{aligned}$

解此方程可用Newton-Raphson法，由于在连续时间步长内变化很小，所以只需要几次迭代就可以接近目标值。

反向放置的二极管

$A\mathrm{e}^{Bx}+x+C=0$

式中

$\begin{aligned} &C=C_{eze}(i,j,k)\times E_{z}^{n}(i,j,k)+C_{ezhy}(i,j,k)\times\left(H_{y}^{n+\frac12}(i,j,k)-H_{y}^{n+\frac12}(i-1,j,k)\right) \\ &+C_{ezhx}(i,j,k)\times\left(H_{x}^{n+\frac12}(i,j,k)-H_{x}^{n+\frac12}(i,j-1,k)\right)+C_{ezd}(i,j,k), \\ &C_{eze}(i,j,k)=-\frac{2\varepsilon_{z}(i,j,k)-\Delta t\sigma_{z}^{e}(i,j,k)}{2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)}, \\ &C_{ezhy}(i,j,k)=-\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)\right)\Delta x}, \\ &C_{ezhx}(i,j,k)=\frac{2\Delta t}{\left(2\varepsilon_{z}(i,j,k)+\Delta t\sigma_{z}^{e}(i,j,k)\right)\Delta y}, \\ &C_{ezd}(i,j,k)=\frac{2\Delta tI_{0}}{2\varepsilon_{z}(i,j,k)\Delta\mathrm{x}\Delta\mathrm{y}+\Delta t\sigma_{z}^{e}(i,j,k)\Delta\mathrm{x}\Delta\mathrm{y}}. \end{aligned}$

集总参数元件的定义、初始化和模拟

集总参数元件的定义

FDTD参量和数组的初始化

集总参数元件的初始化

更新系数的初始化

电压源更新系数的初始化

电流源更新系数的初始化

电感更新系数的初始化

电阻更新系数的初始化

电容更新系数的初始化

二极管更新系数的初始化

电场和磁场以及电压和电流的取样

计算取样电压

计算取样电流

输出参数的定义与初始化

初始化输出参量

运行时间中的图形化显示的初始化

运行FDTD模拟：时进循环

磁场更新

电场更新

与电压源相关的电场更新

与电流源相关的电场更新

与电感相关的电场更新：注意需要使用之前时间的参数来更新。

与二极管相关的电场更新

磁场分量取样的获取

取样电流的获取

取样电场分量的获取

取样电压的获取

显示取样参量

显示FDTD模拟结果

激励源的波形与从时域到频域的变换

常用FDTD仿真波形

正弦波形

如果观察电磁问题对正弦波的激励响应,则需要运行足够长的时间,以使瞬态响应中，由于源的开始而激励出的高频信号完全逝去，只留下主正弦信号。

高斯波形

对于大多数应用,在合理的FDTD模拟中,设最高频率波长大于20倍网格尺寸是充分的。

高斯波形$g(t)=\mathrm{e}^{-(\frac{t}{r})^{2}}$，$\tau$为一参数决定高斯脉冲在时域和频域的宽度。其傅里叶变换是$G(\omega)=\tau\sqrt{\pi}\mathrm{e}^{-\frac{(\pi\omega)}{4}2}$

在FDTD计算中，最高频率由精度参量number of cells wavelength决定：$f_{\mathrm{max}}=\frac{c}{\lambda_{\mathrm{min}}}=\frac{c}{n_{\mathrm{~c}}\Delta s_{\mathrm{max}}}$式中，$c$自由空间光速，$\Delta s_{max}$为$\Delta x,\Delta y,\Delta z$中的最大网格尺寸，$\lambda_{min}$为自由空间中最高工作频率波长。

将高斯波形的有效最高频率考虑为其幅度是最高高斯频谱幅度的10%的频率，可以计算出$\tau=\frac{\sqrt{2.3}n_{\mathrm{c}}\Delta s_{\mathrm{max}}}{\pi c}\approx\frac{n_{\mathrm{c}}\Delta s_{\mathrm{max}}}{2c}$

时域移动高斯波形$g(t)=\mathrm{e}^{-(\frac{r-t_{0}}{r})^{2}}$，$t_0$为时间移动。

高斯波形的导数归一化

该封面图片由Roland Mey在Pixabay上发布