时域有限差分方法发展至今日臻成熟，在瞬态电磁场模拟中有广泛应用，特别是高频瞬态电磁场。不仅如此，其在声波方程和薛定谔方程求解中也有广泛应用，形成较为规范的处理流程：1)基本控制方程的建立，FDTD方法可用于多种偏微分方程的求解，在电磁领域，主要是直接在时间域求解麦克斯韦旋度方程组。电磁波与不同的介质相互作用，麦克斯韦方程略有不同；2)空间剖分，即模拟空间的三维建模，根据所研究问题的环境和边界对模拟空间在直角坐标、柱坐标、极坐标等坐标系中进行剖分，网格尺寸(空间网格步长)和形状根据所模拟的环境和使用的电磁波波长来定，必须满足数值色散条件，对于瞬态脉冲，空间剖分后网格就是一个低通数字滤波器，一般最大网格小于最小波长的1/12，即最大的网格可以让要求的最高频率通过，这样引起的数值误差较小，引起的数值色散在可接受的范围。如，海洋电磁中使用的频率为0.1~10 Hz，真空中的波长以1 Hz为例为$3 \times {10^8}$，海水中波速远比光速小，波长也比真空中小的多，模拟中网格步长可取10 m，甚至更大；地面瞬变电磁中，电磁脉冲的频带很宽，以最大1 MHz为例，模拟中网格步长需小于25 m，可取1 m；超宽带探地雷达探测中，雷达发射ns级电磁脉冲，以1 GHz为例，模拟中网格步长需小于2.5 cm，可取1 cm；而在回旋管瞬态场模拟中，如微波频率为94 GHz，则空间步长必须小于0.026 mm，实际可取0.02 mm。根据实际，网格步长可以是均匀的，也可以是非均匀的(连续渐变或者按一定规则突变)；3)时间剖分，时间步长，可以是均匀的，也可以是非均匀的，在确定了空间步长后，时间步长要满足Courant稳定性条件，否则求解过程就不收敛，不能得到正确的结果；4)在剖分后的模拟空间网格单元填充介质参数，主要指介质的介电常数、电导率和磁导率等；5)满足稳定性要求，这是数值求解微分方程的基本要求，即求解过程必须收敛，要逐步稳定的向解析结果靠近，否则求解失败，没有意义。FDTD中最基本的稳定性要求是波动方程离散后在时间步长和空间步长之间应当满足Courant稳定性条件；6)注意数值色散，对于单色波，即使介质本身无损耗，无色散，对时空做差分离散后，会导致波数和频率的非线性关系，因而使波的相速度跟频率有关，造成波形畸变，是为数值色散。可以通过减小空间网格和使用高阶差分的方法降低数值色散；7)源的处理，麦克斯韦方程中，不论是电流源、磁流源还是运动的粒子源，最后都表现为电流密度。在回旋管瞬态场模拟中与其他不同的是，源的处理相当复杂，电流密度由电子在互作用腔中螺旋运动引起，要结合电子的相对论动力学方程和弗拉索夫电子密度分布方程来求解回旋电子引起的电流密度，还要用PIC方法将电子放到网格点上，并对非网格点进行插值；8)边界条件的处理，PML是FDTD模拟中公认较好的边界条件处理方法，要求不高的情况下可以用简单的Mur吸收边界^[4-6]，回旋管模拟中用PML边界条件。

麦克斯韦方程组的旋度(非铁磁介质)形式为：

式中，$\mathit{\boldsymbol{H}}$，$\mathit{\boldsymbol{E}}$分别为磁场强度，电场强度；$\mu $和$\varepsilon $为介质中磁导率和介电常数；$\mathit{\boldsymbol{J}}$为电流密度，通常有$\mathit{\boldsymbol{J}} = {\rm{ \mathsf{ σ} }}\mathit{\boldsymbol{E}}$，其中${\rm{ \mathsf{ σ} }}$为介质电导率。

式(1)、(2)整理变换可得到：

在经典瞬态场模拟中如果介电常数起主要作用，即电磁波频率较高，位移电流起主要作用时，传导电流可忽略，即满足${\rm{ \mathsf{ σ} }} \ll {\rm{ \mathsf{ ε}}} {\rm{ \mathsf{ ω}}} $，则式(3)、式(4)变成波动方程；如果电磁波的频率较低，考察与导电介质的相互作用，电导率起主要作用时，由此引起的传导电流起主要作用，位移电流可忽略，即满足${\rm{ \mathsf{ σ} }} \gg {\rm{ \mathsf{ ε}}} {\rm{ \mathsf{ ω}}} $，则式(3)、式(4)变为扩散方程；如果二者都不能忽略，则式(3)、式(4)的解皆具有波动和扩散的双重特性。回旋管瞬态场模拟将都情况较为特殊，其求解的是真空中的波动方程，但是有由回旋电子注引起的电流密度项。

海洋可控源电磁法是用船拖拽一条距海底50~100 m的电缆作为电偶极子，向偶极子发送几百到几千安培的电流信号，在水中激发0.1~10 Hz的低频电磁波，同时在海底放置或拖拽一定数量的电磁探测传感器，通过其数据处理以了解海底矿藏、油气藏分布及其他异常的电磁探测方法。海洋可控源电磁法作为海洋地震法的辅助手段，弥补了海洋地震法的不足，降低了海洋油气开发中钻井的风险，2000年以来受到世界各大石油公司的青睐。用FDTD法对海洋电磁场进行模拟可以了解低频电磁场在水中的传播及其与介质的作用规律。海洋电磁场满足准静态场的条件，即${\rm{ \mathsf{ σ} }} \gg {\rm{ \mathsf{ ε}}} {\rm{ \mathsf{ ω}}} $(空气中介电常数${\varepsilon _0} = 8.85 \times {10^{ - 12}}$，相对介电常数和角频率的乘积不超过2个数量级，按10^-10算，海水电导率0~10，通常取平均电导率3.3 s/m，电导率起主要作用，位移电流可忽略，方程退化为扩散方程。式(2)变为：

式中没有电场对时间的导数项，电场和磁场无法随时间交错递推。对式(5)差分后做一些数学变换，

构造一个虚拟介电常数${\rm{ \mathsf{ γ} }} \geqslant \frac{{3\Delta {t^2}}}{{\mu {\Delta ^2}}}$，其中$\Delta t$为时间步长，有经验公式$\Delta {t_{{\rm{max}}}} = {\rm{ \mathsf{ α} }}{\Delta _{{\rm{min}}}}\sqrt {\frac{{\mu {\rm{ \mathsf{ σ} }}t}}{6}} $，${\Delta _{{\rm{min}}}}$为最小空间步长，${\rm{ \mathsf{ α} }} \in (0.1,0.2)$，可见γ是模拟时间的函数，变换后式(5)变为：

式(6)与式(1)结合，整理后得到的微分方程与原来的扩散方程是相容的。在偶极子源处，式(6)需再加一个源电流密度项，构造了虚拟介电常数后，FDTD方法得以实现。图 1所示为导电全空间的电导率3.3 s/m，计算网格100×100×100，在中心有一个指向x正方向的电偶极子源，网格步长取均匀网格50 m，最大时间步长12.9 ms，则得出1.16 s时各坐标面切片上Ex的等值线分布。与用解析公式算出的结果相似，证明FDTD计算的正确性。

图  1  1.16 s时Ex空间分布的等值线图

该模拟中用的边界条件是自然边界，即无穷远处场值衰减为0，这比用Mur吸收边界略差一点，从图 1中可以看到边界的回波，但是计算速度却加快了近一倍。减小边界反射的另一个方法就是采用渐变网格，即在源附近采用细网格，而远离源靠近边界的区域采用逐渐增大的网格。

以上公式在处理好海水和空气边界后就可以模拟海洋电磁中更复杂的情形^[7]。

图 2为三维FDTD模拟的海底低阻体地质模型在不同时刻的电磁响应快照。

图  2  海底低阻体的电磁响应FDTD模拟

地面瞬变电磁法具有对低电阻率目标反映灵敏的特性。它作为一种有效的电磁勘探方法在近地表地球物理勘查中发挥了重要的作用，如考古、地热、地下水勘探，城市大型建筑物施工前的地质勘查、桥梁隧道施工前的地质勘查和地质灾害预报等。用FDTD方法进行正演模拟可以对地质体三维建模，快速求解地质模型的电磁响应，为探测数据的反演解释提供理论依据。所用的理论模型与2.1节相同，只是目前地面瞬变电磁法多数情况用的是回线源，所激发的电磁脉冲所含频率比海洋电磁要高，如1 K~1 MHz。文献[8-9]对隧道中充水断层的灾害勘测做了较为全面的三维FDTD模拟。一方面验证了FDTD方法在地面瞬变电磁探测模拟中的有效性；另一方面为地面时频电磁勘探提供了强有力的工具。如图 3a为隧道XZ断面图，图 3b为隧道掌子面图，右手坐标系，X垂直向上，Y水平向右，Z垂直纸面向里为正，6 m的方形隧道掌子面中央放置3 m宽的方形回线源，掌子面后10 m处有一充水断层，图 3c为不同时刻充水断层的电磁响应在XZ断面上Ex电场空间分布。

图  3  隧道模型及其电磁响应的FDTD模拟

UWB探地雷达用于城市灾害救助，如生命搜救、反恐探测，工程验收时材料、结构及厚度的检测，市政金属或非金属管道、下水道、缆线、缆线管道、孔洞、基础层、混凝土中的钢筋及其他地下埋件的位置，坝体渗漏探测中，渗漏部位及范围的反演成像，工程施工前作业面的地质调查等。用FDTD方法对UWB探地雷达正演模拟，其所用理论模型为式(1)、式(2)，既要考虑位移电流的作用，因为雷达电磁脉冲所含频谱宽，频率高，又要考虑传导电流的作用，因为电磁波在地下不可避免的会碰到导电介质。

图 4为UWB雷达穿墙探测人的侧面模型图，图 5为人体的正面模型，建立坐标系，矩形面天线馈入半带宽为1.22 GHz的截断三正弦电流信号，观察不同时刻目标的电磁响应。用FDTD方法进行模拟，图 6a~6d分别为8 ns时刻的天线面、墙面、人体正面以及垂直于墙面中央的YZ面上的Ex分布。

图  4  UWB穿墙探测模型侧面图

图  5  UWB穿墙探测模型人体截面图

图  6  8 ns时刻各个面上的Ex分布

从图 6可以看出，超宽带雷达就像一个面光源照亮了整个探测空间。图 6d显示，人体反射波即将到达墙面，此时对雷达扫描数据进行处理就可以对墙后的物体进行成像。

图 7为文献[10]做的UWB雷达探测地下钢筋混凝土FDTD模拟，单频1.5 GHz探地雷达和中心频率为4 GHz的UWB雷达探测正演模拟对比图，可见该雷达可以清晰分辨混凝土中间距为5 cm的钢筋。此模型说明，混凝土中有两组间距不同的钢筋，钢筋直径为2.5 cm，中心埋深为5 cm，上层介质为混凝土，相对介电常数为6，电导率为0.005 S/m；下层为干沙，相对介电常数为3，电导率为0.000 1 S/m。SBS表示钢筋间距。实验体现了UWB雷达的高分辨率成像特性。

图  7  UWB探地雷达对混凝土中钢筋分辨的FDTD模拟

基于相对论电子回旋脉塞理论的回旋管内电磁场是电子的回旋受激辐射，而高能电子在电磁场中运动时基于相对论效应而产生角向群聚，并在波频率和电子的回旋频率满足回旋谐振产生换能^{[3, 11]}。提高工作频率、功率、互作用效率和稳定性是回旋管的研究主题，而这些都被如起振电流、电压、电子数、回旋半径、磁场、互作用腔结构等多种因素影响。能动态了解回旋管内电子的相位和电磁场的相位，对回旋管的理论研究和工程设计有益。三维FDTD-PIC模拟是在固定腔体空间中研究运动的电子和变化的电磁场互作用的规律，从而检验所设计腔体和输出电磁场频率、模式、功率等参数的合理性，目的是降低实际制管的风险，节约人力、物力和财力。

目前，国外有MAGIC和MAGY等回旋管瞬态场模拟软件，国内回旋管瞬态场的模拟软件还存在局限性，特别是三维模拟还不够，FDTD方法较为成熟，但涉及到高能电子的运动及与场的互作用，就使得理论模型变的复杂。虽然基本方程仍然用式(1)、式(2)，采用导电边界，但电流密度的表达式确定却比较困难，电流密度由电子在腔体中的分布和运动的速度等因素决定，因此模拟的关键在于确定电流密度。如图 8所示，用现有的MAGIC模拟软件做94 GHz回旋管TE6, 2模的输出过程，由于篇幅所限，腔体参数暂未给出。工作参数如下：电子注电压为42 kV，电流为2 A，横纵速度比为1.3，引导中心半径为3 mm，工作磁场取3.556 8 Tesla.

图  8  94 GHz回旋管TE6, 2工作过程FTDT-PIC模拟

所使用的模拟软件为单线程MAGIC，过程耗时，40 ns的工作时间需要计算13天。经过3个月的优化，最后得到如图 8所示的模拟结果，从8a、8b可知30 ns后有稳定输出，但功率并不高，从8c可知互作用效率低，电子交出能量不充分，从8d可知，有两个输出，87.4 GHz和95.4 GHz，存在模式竞争，TE6, 2模还没有出现，模拟结果仍然不理想。有待进一步优化参数，提高互作用效率，抑制模式竞争。

柱坐标系下，r、φ、z方向满足Courant稳定性条件的最大步长分别为0.02、${\rm{ \mathsf{ π} }}/12$、0.02 mm。按照给定尺寸，剖分网格数最大为400×24×4 000，内存使用小于500 MB，40 ns的模拟时间耗时300 h，对现有计算机资源造成极大的浪费。如果能将上述程序做一个8线程并行，不到两天就可以完成一轮测试，将会提高工作效率。因此，回旋管瞬态场三维并行FDTD-PIC模拟软件的开发很有必要。

时域有限差分方法非常适合并行计算，虽然时间维无法并行，但可以在空间维进行分块并行，也可以针对各个场量进行任务并行，每一个时间步协调刷新全模拟空间各个任务的电磁场值。可以从多方面对程序加速：1)基于多节点的网络并行，即多处理器集群。从空间剖分和任务划分的角度，由主计算节点把整个空间分成尽量一致的若干块，按计算任务块分配到不同的从计算节点，每一个时间步由主节点与从节点进行数据交换与衔接，以保持整个模拟空间数据的连续性和时间的一致性。该种并行技术可以组建超大型的并行集群或计算网格，适合大尺度、大规模的计算模拟，由于节点间的通信开销较大，小规模的计算任务并不适合在集群上计算；2)基于共享存储的多线程并行技术OpenMP，针对多核计算机，可以充分利用目前多核计算机的资源，既可以基于任务，把块分给不同线程计算，也可以对某一个循环进行分割。三维FDTD计算中，用HPZ820工作站进行计算，网格规模200×200×160，开20个线程，对循环进行并行，并行加速比可达14；开10个线程，并行加速比可达7.5。相同条件下，并行加速比与所开线程数、计算规模以及任务的分配都有关系，并不一定开的线程越多，加速比越高；3)基于SSE指令集的核内细粒度并行技术。在多核处理器中，每个核内都有自己的缓存、浮点运算单元(FPU)和向量算术逻辑单元(VALU)。FPU一次处理1对数据，得到1个计算结果，而VALU一次可以处理4对数据，得到4个计算结果，是典型的单指令多数据操作。一般的高级语言调用的是FPU，无法调用VALU，而SSE指令集却可以调用VALU进行运算。理论上其加速比为4，文献[12]的实验中，网格规模是200×200×200，得到加速比为2.44；4)基于GPU技术的协同加速计算。GPU特别适合高度并行以及数据密集型任务的计算，FDTD计算方法空间高度并行，模型数据密集，因此可以实现CPU+GPU的协同计算，但需要注意是部分任务是无法并行的，需要解决好CPU与GPU的通信开销与效率提高的矛盾。当然，也可以将全部运算放到GPU上进行，这不失为一个很好的解决方案，CUDA架构中，被抽象化的线程并行操作一个或若干个Yee格点的迭代计算，这利用了GPU硬件中含有上万个算术逻辑单元的优势，但是这须有足够的显存^[13]。虽然GPU的加速计算在许多领域已经有成熟的应用，但还有较多技术问题需要在不同的实践中加以解决。

本文首先总结了时域有限差分方法在瞬态电磁场模拟中的基本流程，然后给出了其在海洋电磁探测模拟、地面瞬态电磁探测模拟、探地雷达探测模拟以及回旋管设计等方面的应用案例，并分别作了分析，最后研究了FDTD程序并行策略。目前，国内FDTD方法已经广泛用于各种瞬态电磁场模拟场景中，但三维FDTD-PIC方法在回旋管瞬态电磁场的数值模拟中还有研究空间。相对论电子学是回旋管瞬态电磁场模拟的理论基础，理论难度大是制约数值模拟发展的主要因素。下一步，可以从3个方面开展研究工作：1)研究相对论电子学理论，推导相关公式，建立完善的数理模型；2)引进国外先进模拟软件，提高设计水平，开发和完善面向对象的数值模拟软件模块。参照MAGIC，MAGY，VSIM等软件，结合国内外在粒子模拟方面已取得的成就和积累的经验，开发具有自主知识产权的FDTD-PIC并行模拟程序；3)加强数值模拟指导下的实验测试。有效的数值模拟可使实验次数减少，但无法取代实际测试。一方面，回旋管在交付使用前必须经过严格的实测检验，确保其稳定可靠；另一方面，强流电子束在谐振腔中的相对论效应，有待在实验中进一步观察，证实已有理论和模拟结果的正确性，进一步推动理论的发展和模拟质量的提高。

将理论研究、数值模拟和实验测试三者结合起来是开发出优秀软件的必由之路。

本文的研究工作得到了警察大学自然基金培育项目(ZKJJPY201717)的资助，在此表示感谢！