计算电磁学：FDTD算法总结

JiNan.YouQuan.Soft

已于 2024-06-15 00:35:45 修改

阅读量3.5k

点赞数 10

分类专栏： CAx HPC 文章标签：算法

于 2024-02-24 01:09:21 首次发布

JiNan YouQuan Software Co. Ltd.

本文链接：https://blog.csdn.net/qq_26221775/article/details/136264673

版权

CAx 同时被 2 个专栏收录

65 篇文章 20 订阅

订阅专栏

HPC

7 篇文章 6 订阅

订阅专栏

本文介绍了计算电磁学中的FDTD（FiniteDifferenceTimeDomain）算法，涵盖了物理模型（如库仑定律和洛伦兹力）、数学模型（Maxwell方程组及其离散化）、时间与空间离散方法、边界条件处理、源项处理以及收敛性和稳定性条件。同时提到了相关理论书籍和开源资源。

摘要由CSDN通过智能技术生成

计算电磁学(Computational Electromagnetics, CEM)是通过数值计算来研究电磁场的交叉学科。

数值求解电磁学问题的方法可以分成频域(Frequency Doamin, FD)、时域(Time Domain, TD)等两类。

频域法基于时谐微分，通过对多个采样值的傅里叶逆变换得到所需的脉冲响应，使用这种方法，每次计算只能求得一个频率点上的响应。这类方法又可进一步分成低频算法、高频算法等两类。低频算法包括矩量法(Method of Moment, MoM)、频域有限差分(Finite Difference Frequency Doamin, FDFD)等；高频算法包括几何光学法、物理光学法等。

时域法按时间步进求得有关场量，一次求解可以获得很宽频带范围内的解。这类方法包括时域有限差分(Finite Difference Time Domain，FDTD)、时域有限单元(Finite Element Time Domain，FETD)等。

在时域法中，最为著名的就是FDTD。因此，本文拟对FDTD算法涉及的数理模型、数值模型等内容进行简要介绍。

注1：限于研究水平，分析难免不当，欢迎批评指正。

注2：文章内容会不定期更新。

零、预修：物理模型

0.1 Coulomb's law

1785年由法国科学家C,-A.de库仑根据实验得出，真空中两个静止的点电荷之间的相互作用力与它们的电荷量的乘积成正比，与它们的距离的二次方成反比，作用力的方向在它们的连线上，同名电荷相斥，异名电荷相吸。

$\boldsymbol{f}=k\frac{q_{1}q_{2}}{r^{2}}\boldsymbol{e}_{r}$

其中， $\boldsymbol{e}_{r}$ 是从 $q_{1}$ 到 $q_{2}$ 的单位矢量， $k$ 为库伦常数，大小为 $9\times 10^{9}Nm^{2}/C^{2}$ 。

0.2 Lorentz Force

1895年，荷兰物理学家H·A·洛伦兹建立经典电子论时,作为基本假定首先提出了运动电荷产生磁场和磁场对运动电荷有作用力的观点。

洛伦兹力是运动电荷在磁场中所受到的力，即磁场对运动电荷的作用力。

$\boldsymbol{f}=q\left ( \boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B} \right )$

其中，q是带电粒子的电荷量； $\boldsymbol{E}$ 是电场强度； $\boldsymbol{v}$ 是带电粒子的速度； $B$ 是磁感应强度。

洛伦兹力的方向可以通过左手定则来判断，即伸开左手，使拇指与其余四个手指垂直，并且都与手掌处于同一水平面，让磁感线从掌心进入，四指指向正电荷运动的方向，拇指指的方向即洛伦兹力的方向。

0.3 Ampere circuital theorem

0.4 Faraday law of electromagnetic induction

一、数学模型：Maxell方程组

1864年，Maxwell在前人工作的基础之上，建立了统一的电磁场理论，并用一组数学方程揭示了宏观电磁场的基本规律，这就是著名的Maxwell方程组。

Maxwell方程组有四个方程组成：描述电荷如何产生电场的高斯定律；论述磁单极子不存在的高斯磁定律；描述电流和时变电场怎样产生磁场的安培环路定律；描述时变磁场如何产生电场的法拉第感应定律。

$\left\{\begin{matrix} \nabla \cdot \boldsymbol{D}=\rho\\ \nabla\cdot \boldsymbol{B}=0\\ \nabla\times \boldsymbol{H}=\boldsymbol{J}+\frac{\partial \boldsymbol{D}}{\partial t}\\ \nabla\times \boldsymbol{E}=-\frac{\partial \boldsymbol{B}}{\partial t}-\boldsymbol{M} \end{matrix}\right.$

可以看出，电流密度 $\boldsymbol{J}$ 、磁流密度 $\boldsymbol{M}$ 是以源项的形式出现在磁场旋度方程、电场旋度方程。

本构关系的一般表达式为，

$\left\{\begin{matrix} \boldsymbol{D}=\boldsymbol{\epsilon}\boldsymbol{E}+\boldsymbol{\xi }\boldsymbol{H}\\ \boldsymbol{B}=\boldsymbol{\mu} \boldsymbol{E}+\boldsymbol{\zeta} \boldsymbol{H} \end{matrix}\right.$

若 $\boldsymbol{\epsilon},\boldsymbol{\xi },\boldsymbol{\mu},\boldsymbol{\zeta}$ 是角频率 $\omega$ 的函数，则称介质是色散介质。

对于均匀各向同性介质，本构方程简化为，

$\left\{\begin{matrix} \boldsymbol{D}=\varepsilon\boldsymbol{E}\\ \boldsymbol{B}=\mu \boldsymbol{H}\end{matrix}\right.$

在介质交界面处，场变量可能会发生突变，基于积分形式的Maxwell方程组，可得以下关系式，

$\left\{\begin{matrix} \boldsymbol{n}\times \left ( \boldsymbol{E}_{1}-\boldsymbol{E}_{2} \right )=\boldsymbol{0}\\ \boldsymbol{n}\times \left ( \boldsymbol{H}_{1}-\boldsymbol{H}_{2} \right )=\boldsymbol{J}_{s} \\ \boldsymbol{n}\cdot \left ( \boldsymbol{D}_{1}-\boldsymbol{D}_{2} \right )=\rho _{s} \\ \boldsymbol{n}\cdot \left ( \boldsymbol{B}_{1}-\boldsymbol{B}_{2} \right )=0 \end{matrix}\right.$

在上述公式中，各符号含义如下：

$\boldsymbol{E}$	电场强度	$V/m$
$\boldsymbol{D}$	电通量密度	$C/m^{2}$
$\epsilon$	介电常数	在真空中， $\epsilon=\epsilon_{0}=8.8542\times 10^{-12} F/m$
$\rho$	电荷密度	$C/m^{3}$
$\mathbf{J}$	电流密度	$A/m^{2}$
$\sigma$	电导率	$S/m$
$\boldsymbol{H}$	磁场强度	$A/m$
$\boldsymbol{B}$	磁通量密度	$Wb/m^{2}$ ， $T$
$\mu$	磁导率	在真空中， $\mu =\mu_{0}=4\pi \times 10^{-7} H/m$
$\boldsymbol{M}$	磁流密度	$V/m^{2}$
$\boldsymbol{J}_{s}$	交界面电流i密度	$A/m$
$\rho _{s}$	交界面电荷密度	$C/m^{2}$

不考虑电流密度 $\boldsymbol{J}$ 与磁流密度 $\boldsymbol{M}$ ，在三维笛卡尔坐标系下，有

$\frac{\partial E_{x}}{\partial x}+\frac{\partial E_{y}}{\partial y}+\frac{\partial E_{z}}{\partial z}=\frac{\rho }{\varepsilon }$

$\frac{\partial H_{x}}{\partial x}+\frac{\partial H_{y}}{\partial y}+\frac{\partial H_{z}}{\partial z}=0$

$\begin{matrix} \frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}=\epsilon \frac{\partial E_{x}}{\partial t}\\ \frac{\partial H_{x}}{\partial z}-\frac{\partial H_{z}}{\partial x}=\epsilon \frac{\partial E_{y}}{\partial t}\\ \frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}=\epsilon \frac{\partial E_{z}}{\partial t} \end{matrix}$

$\begin{matrix} \frac{\partial E_{y}}{\partial z}-\frac{\partial E_{z}}{\partial y}=\mu \frac{\partial H_{x}}{\partial t}\\ \frac{\partial E_{z}}{\partial x}-\frac{\partial E_{x}}{\partial z}=\mu \frac{\partial H_{y}}{\partial t}\\ \frac{\partial E_{x}}{\partial y}-\frac{\partial E_{y}}{\partial x}=\mu \frac{\partial H_{z}}{\partial t} \end{matrix}$

若电磁场与 $z$ 无关时，两个旋度方程可简化为，

$TM\left\{\begin{matrix} \epsilon \frac{\partial E_{x}}{\partial t}=\frac{\partial H_{z}}{\partial y}\\ \epsilon \frac{\partial E_{y}}{\partial t}=-\frac{\partial H_{z}}{\partial x}\\ \mu \frac{\partial H_{z}}{\partial t}=\frac{\partial E_{x}}{\partial y}-\frac{\partial E_{y}}{\partial x}\end{matrix}\right.$

$TE\left\{\begin{matrix} \mu \frac{\partial H_{x}}{\partial t}=-\frac{\partial E_{z}}{\partial y}\\ \mu \frac{\partial H_{y}}{\partial t}=\frac{\partial E_{z}}{\partial x}\\\epsilon \frac{\partial E_{z}}{\partial t}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y} \end{matrix}\right.$

二、数学模型：Maxwell波动方程

在无源情况下，对于均匀各向同性介质，Maxwell方程组简化为，

$\left\{\begin{matrix} \nabla \cdot \boldsymbol{E}=0 \\ \nabla\cdot \boldsymbol{H}=0\\ \nabla\times \boldsymbol{H}=\epsilon\frac{\partial \boldsymbol{E}}{\partial t}\\ \nabla\times \boldsymbol{E}=-\mu\frac{\partial \boldsymbol{H}}{\partial t} \end{matrix}\right.$

对于电场旋度方程，有，

$\mu^{-1}\nabla\times \boldsymbol{E}=-\frac{\partial \boldsymbol{H}}{\partial t} \Rightarrow \mu^{-1}\nabla\times \nabla\times \boldsymbol{E}=-\epsilon\frac{\partial^{2} \boldsymbol{E}}{\partial t^{2}}$

考虑到恒等式 $\nabla\times \nabla\times \boldsymbol{E}=\nabla\nabla\cdot \boldsymbol{E}-\nabla^{2}\boldsymbol{E}$ ，即有

$\frac{\partial^{2} \boldsymbol{E}}{\partial t^{2}}=\frac{1}{\epsilon\mu}\nabla^{2}\boldsymbol{E}$

对于磁场旋度方程，亦可得出，

$\frac{\partial^{2} \boldsymbol{H}}{\partial t^{2}}=\frac{1}{\epsilon\mu}\nabla^{2}\boldsymbol{H}$

三、数值模型

FDTD是美籍华人K.S. Yee 于1966 年提出的求解微分形式Maxwell方程组的数值方法。FDTD在时域上直接求解微分形式的Maxwell方程组。

3.1 空间离散

在FDTD中，Yee使用Yee Cell在空间上交叉布置电场、磁场分量。电场分量放置到Yee单元各棱的中间，方向平行于各棱；磁场分量放置到Yee单元各面中心，方向平行于各面法线。

这样，每个电场分量由四个磁场分量环绕；每个磁场分量也由四个电场分量环绕。

对于安培环路定律与法拉第电磁感应定律，采用中心差分来离散磁场旋度与电场旋度。

为对于编号为 $\left ( i,j,k \right )$ 的Yee Cell，其中 $i\in \left [ 0, n_{x}\right )$ ， $j\in \left [ 0, n_{y}\right )$ ， $k\in \left [ 0, n_{z}\right )$ ， $n_{x}$ 、 $n_{y}$ 、 $n_{z}$ 为计算域在三个坐标轴方向的网格数，

可得法拉第感应定律、安培环路定律在空间域上的半离散格式，

$\begin{matrix}\mu \frac{\partial H_{x}\left ( i,j,k \right )}{\partial t}= \frac{E_{y}\left ( i,j,k \right )-E_{y}\left ( i,j,k-1 \right )}{\Delta z}-\frac{E_{z}\left ( i,j,k \right )-E_{z}\left ( i,j-1,k \right )}{\Delta y}\\ \mu \frac{\partial H_{y}\left ( i,j,k \right )}{\partial t}=\frac{E_{z}\left ( i,j,k \right )-E_{z}\left ( i-1,j,k \right )}{\Delta x}-\frac{E_{x}\left ( i,j,k \right )-E_{x}\left ( i,j,k-1 \right )}{\Delta z}\\ \mu \frac{\partial H_{z}\left ( i,j,k \right )}{\partial t}=\frac{E_{x}\left ( i,j,k \right )-E_{x}\left ( i,j-1,k \right )}{\Delta y} -\frac{E_{y}\left ( i,j,k \right )-E_{y}\left ( i-1,j,k \right )}{\Delta x}\end{matrix}$

$\begin{matrix} \epsilon \frac{\partial E_{x}\left ( i,j,k \right )}{\partial t}=\frac{H_{z}\left ( i,j,k \right )-H_{z}\left ( i,j-1,k \right )}{\Delta y}-\frac{H_{y}\left ( i,j,k \right )-H_{y}\left ( i,j,k-1 \right )}{\Delta z}\\ \epsilon \frac{\partial E_{y}\left ( i,j,k \right )}{\partial t}=\frac{H_{x}\left ( i,j,k \right )-H_{x}\left ( i-1,j,k \right )}{\Delta z}-\frac{H_{z}\left ( i,j,k \right )-H_{z}\left ( i-1,j,k \right )}{\Delta x}\\ \epsilon \frac{\partial E_{z}\left ( i,j,k \right )}{\partial t}= \frac{H_{y}\left ( i,j,k \right )-H_{y}\left ( i-1,j,k \right )}{\Delta x}-\frac{H_{x}\left ( i,j,k \right )-H_{x}\left ( i,j-1,k \right )}{\Delta y}\end{matrix}$

3.2 时间离散

在时间域上，FDTD采用蛙跳格式，由此可得法拉第感应定律、安培环路定律在时间域上的半离散格式，

$\begin{matrix} \mu \frac{H_{x}^{n+\frac{1}{2}}-H_{x}^{n-\frac{1}{2}}}{\Delta t}=\frac{\partial E_{y}^{n}}{\partial z}-\frac{\partial E_{z}^{n}}{\partial y}\\ \mu \frac{H_{y}^{n+\frac{1}{2}}-H_{y}^{n-\frac{1}{2}}}{\Delta t}=\frac{\partial E_{z}^{n}}{\partial x}-\frac{\partial E_{x}^{n}}{\partial z}\\ \mu \frac{H_{z}^{n+\frac{1}{2}}-H_{z}^{n-\frac{1}{2}}}{\Delta t}=\frac{\partial E_{x}^{n}}{\partial y}-\frac{\partial E_{y}^{n}}{\partial x} \end{matrix}$

$\begin{matrix}\epsilon \frac{E_{x}^{n+1}-E_{x}^{n}}{\Delta t}=\frac{\partial H_{z}}{\partial y}-\frac{\partial H_{y}}{\partial z}\\ \epsilon \frac{E_{y}^{n+1}-E_{y}^{n}}{\Delta t}=\frac{\partial H_{x}}{\partial z}-\frac{\partial H_{z}}{\partial x}\\ \epsilon \frac{E_{z}^{n+1}-E_{z}^{n}}{\Delta t}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y} \end{matrix}$

3.3 FDTD差分方程组

结合法拉第感应定律、安培环路定律在空间域、时间域上的离散格式，可得最终差分方程组，

$\begin{matrix}\mu \frac{H_{x}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{x}^{n-\frac{1}{2}}\left ( i,j,k \right )}{\Delta t}= \frac{E_{y}^{n}\left ( i,j,k \right )-E_{y}^{n}\left ( i,j,k-1 \right )}{\Delta z}-\frac{E_{z}^{n}\left ( i,j,k \right )-E_{z}^{n}\left ( i,j-1,k \right )}{\Delta y}\\ \mu \frac{H_{y}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{y}^{n-\frac{1}{2}}\left ( i,j,k \right )}{\Delta t}=\frac{E_{z}^{n}\left ( i,j,k \right )-E_{z}^{n}\left ( i-1,j,k \right )}{\Delta x}-\frac{E_{x}^{n}\left ( i,j,k \right )-E_{x}^{n}\left ( i,j,k-1 \right )}{\Delta z}\\ \mu \frac{H_{z}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{z}^{n-\frac{1}{2}}\left ( i,j,k \right )}{\Delta t}=\frac{E_{x}^{n}\left ( i,j,k \right )-E_{x}^{n}\left ( i,j-1,k \right )}{\Delta y}-\frac{E_{y}^{n}\left ( i,j,k \right )-E_{y}^{n}\left ( i-1,j,k \right )}{\Delta x} \end{matrix}$

$\begin{matrix} \epsilon \frac{E_{x}^{n+1}\left ( i,j,k \right )-E_{x}^{n+1}\left ( i,j,k \right )}{\Delta t}=\frac{H_{z}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{z}^{n+\frac{1}{2}}\left ( i,j-1,k \right )}{\Delta y}-\frac{H_{y}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{y}^{n+\frac{1}{2}}\left ( i,j,k-1 \right )}{\Delta z}\\ \epsilon \frac{E_{y}^{n+1}\left ( i,j,k \right )-E_{y}^{n}\left ( i,j,k \right )}{ \Delta t}=\frac{H_{x}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{x}^{n+\frac{1}{2}}\left ( i-1,j,k \right )}{\Delta z}-\frac{H_{z}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{z}^{n+\frac{1}{2}}\left ( i-1,j,k \right )}{\Delta x}\\ \epsilon \frac{E_{z}^{n+1}\left ( i,j,k \right )-E_{z}^{n}\left ( i,j,k \right )}{\Delta t}= \frac{H_{y}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{y}^{n+\frac{1}{2}}\left ( i-1,j,k \right )}{\Delta x}-\frac{H_{x}^{n+\frac{1}{2}}\left ( i,j,k \right )-H_{x}^{n+\frac{1}{2}}\left ( i,j-1,k \right )}{\Delta y}\end{matrix}$

3.4 边界条件

由于数值模拟只能选取有限区域作为计算域进行仿真计算，因此，在计算域边界处，需要给出吸收边界条件。

目前效果最好的吸收边界层是Berenger于1994年提出的完美匹配层(Perfect Matched Layer, PML)。完全匹配层是数学上在FDTD区域截断边界处虚拟设置一种特殊介质层,通过合理地选择PML的本构参数,能够使FDTD计算域的外行电磁波无反射地穿过分界面,并在吸收层内被迅速吸收分解。

3.5 源项处理

根据源项随时间的变化，源项可分为周期性源项和非周期性源项。

根据源项几何形状，可分为点源、线源、面源等。

3.6 收敛性、稳定性条件

设 $c$ 为计算空间内电磁波最大传播速度， $\Delta t$ 为时间步长， $\delta$ 为网格尺寸，则有

对于三维均匀网格， $c\Delta t\leq \frac{\delta }{\sqrt{3}}$ 。

对于二维均匀网格， $c\Delta t\leq \frac{\delta }{\sqrt{2}}$ 。

对于一维均匀网格， $c\Delta t\leq \delta$ 。

三、实现与测试

参考文献

王长清. 现代计算电磁学基础. 2005.
K S Yee .Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[J].IEEE Transactions on Antennas & Propagation, 1966, 14(5):302-307.DOI:10.1109/TAP.1966.1138693.
赖生建. 计算电磁学讲义.
梁铨廷. 物理光学(第五版). 2018
王元明. 数学物理方程与特殊函数.
陆庆乐. 复变函数.
张元林. 积分变换.
老大中. 变分法基础.