0   引　言

全电磁粒子模拟方法通过求解带电粒子与电磁场之间的耦合方程^[1]，能准确描述电磁场与带电粒子之间的非线性相互作用，在求解电磁场随时间的演化过程以及带电粒子的运动信息^[2]、真空电子器件的设计与优化^[3-7]、微波介质击穿模拟、高能射线产生系统电磁脉冲^[8-9]及核爆炸人工辐射带模拟^[10]等方面得到了广泛应用。传统的全电磁粒子模拟方法采用时域有限差分(finite-difference time-domain, FDTD)方法求解电磁场随时间的变化，具有算法实现相对简单、容易并行等优点。基于共形网格的FDTD方法能进一步提高全电磁粒子模拟方法的计算精度^[11]。目前，国内外均研制了多款基于FDTD方法的全电磁粒子模拟软件，主要用于模拟强场相关的等离子体现象，如电子科技大学研制的CHIPIC软件^[12]、北京应用物理与计算研究所研制的NEPTUNE3D软件^[13]、西北核技术研究所和西安交通大学联合研制的UNIPIC软件和UNIPIC-3D^[14-15]软件等，国外研制了MAGIC、KARAT及CST等粒子模拟软件。

电磁场时域间断伽辽金 (discontinuous Galerkin time-domain, DGTD) 方法是一种新发展的电磁场时域计算方法^[16-19]。该方法采用非结构网格剖分计算模型，采用基函数的线性组合描述网格单位内电磁场分布，通过引入数值通量，每个单元内电磁场的求解仅与本网格及相邻网格单元有关，避免了大型矩阵求逆。该方法具有计算精度高、易于并行等优点，国内外学者开展了大量的时域间断伽辽金-粒子模拟(discontinuous Galerkin time-domain particle-in-cell, DGTD-PIC)方法研究^[20]，有望进一步提高全电磁粒子模拟方法的计算精度，打通传统全电磁粒子模拟技术与采用非结构网格有限元方法之间的技术壁垒。近年来，DGTD-PIC方法在核心关键技术方面取得了许多突破性研究进展^[21]，为DGTD-PIC软件研制及解决实际工程问题提供了可行的技术方案。由于DGTD-PIC方法涉及的计算模型复杂，程序和软件实现难度大，同时部分适用于传统全电磁粒子模拟方法的数值算法无法有效应用于DGTD-PIC方法。为促进DGTD-PIC方法的发展及应用，结合作者所在团队在DGTD-PIC方法方面已开展的研究工作，详细介绍了国内外在DGTD-PIC方法方面的研究成果，主要包括基于修正势的DGTD电磁场计算方法^[22-25]、复频移完全匹配层(complex frequency shifted perfectly matched layer, CFS-PML)边界^[26-28]、电荷及电流分配^[29-31]等关键技术方面取得的研究进展。采用该方法研制的软件已成功应用于回旋振荡器的数值模拟研究，解决了工作在高阶模式下回旋管的波束互作用及模式转换等计算难题^[32-33]。

1   基于修正势的DGTD电磁场计算方法

在全电磁粒子模拟算法中，电磁场以及电荷分配在离散的网格单元上，由于电荷分配以及电流分配采用不同的计算公式造成数值误差，导致电场计算不准确，甚至出现计算结果发散的情况。相比基于FDTD的全电磁粒子模拟方法，DGTD-PIC方法采用非结构网格，该方法在计算电荷分配以及电流分配时更易导致离散的电场值不满足高斯定律。为此，发展了基于修正势的Maxwell方程：

式中：${\boldsymbol{E}}$和${\boldsymbol{H}}$分别为电场强度和磁场强度；${\varepsilon _0}$为真空介电常数；$\rho$和${\boldsymbol{J}}$分别为电荷密度和电流密度；${c_0}$为电磁波在自由空间中的传播速度；$\chi ,\kappa$为修正系数；${\varepsilon }_{{\mathrm{d}}}$为衰减系数，且${\varepsilon }_{{\mathrm{d}}}\ll 1$。辅助变量$\varPhi$和$\varPsi$及相关辅助方程的引入是为了减小数值计算中不满足高斯定律所导致的误差，当$\varPhi \to 0、 {\mkern 1mu} {\mkern 1mu} \varPsi \to 0$时，方程(1)退化为标准的Maxwell方程。为减小数值计算误差，方程(1)对磁场强度作了如下变换：${\boldsymbol{H}} \to \sqrt {{\mu _0}/{\varepsilon _0}} {\boldsymbol{H}}$，${\mu _0}$为真空磁导率。

从方程(1)可以推导得到如下表达式：

式(2)表明修正势$\varPhi$和$\varPsi$用于矫正电荷守恒和磁荷守恒的误差，且可以看出误差以波的形式传播，修正势的传播速度分别为$\chi {c_0}、\kappa {c_0}$。当参数$\chi > 1、\kappa > 1$时，则误差的传播速度将快于电磁波的传播速度，表明增加的辅助方程能避免电荷不守恒带来的误差对电磁场计算结果的影响。

可将方程(1)写成如下表达式：

式中，$\bar {\boldsymbol{I}}$为单位张量。将计算域$\varOmega$划分为$K$个非结构网格单元，即$\varOmega = \bigcup\limits_{k = 1}^K {{\varOmega _k}}$，采用两次分部积分，同时引入数值通量，将式(3)变成强解积分形式：

式中，${\ell _i}$为第$i$个插值点对应的基函数，其可为标量基函数或矢量基函数，通过选择高阶基函数能进一步提高计算精度；$\partial\mathrm{\varOmega}_k$为电磁场在$\mathrm{\varOmega}_k$的边界区域；$\hat {\boldsymbol{n}}$为网格单元边界的外法向矢量；${{\boldsymbol{E}}^ * }$、${{\boldsymbol{H}}^ * }$、${\varPhi ^ * }$、${\varPsi ^ * }$分别为网格单元表面处电场强度、磁场强度以及修正势的数值通量。依据上述方程的黎曼问题求解结果，可得到各个数值通量的表达式：

上述公式中$\langle {\text{ }}\rangle , [\kern-0.15em[ {\text{ }} ]\kern-0.15em]$符号的意义为：

式中，${\boldsymbol{u}}$、${{\boldsymbol{u}}^ + }$分别为本网格单元以及邻接网格单元的物理量，其可以为标量或矢量。将方程(5)代入到方程(4)可以得到如下强变分形式：

对理想导体边界，数值通量满足的关系式如下：

将式(5)代入方程(8)可以得到相邻的两个网格在边界面处电场强度、磁场强度及引入的辅助变量所满足的表达式，其他的理想磁边界以及吸收边界等边界条件可以通过相同的方式推导得到。当采用一阶基函数，$\chi、\kappa$取不同值时电场和磁场的散度误差如图1所示，可以看出，采用修正势的DGTD电磁场计算方法能够大幅降低散度误差，减少散度误差的累积，使得长时间模拟时计算结果保持稳定。

图  1  电场强度和磁感应强度的散度误差随时间的变化^[22]

Fig.  1  Numerical errors of Gauss’s law for the electric field and the magnetic flux intensity vs. time^[22]

下载: 全尺寸图片 幻灯片

2   CFS-PML吸收边界

DGTD-PIC方法只能对有限区域进行模拟，为了能模拟开域电磁场与带电粒子之间的非线性相互作用，在模拟区域的截断边界处设置吸收边界条件。各向异性单轴完全匹配层(uniaxial perfectly matched layer, UPML)由于具有实现简单、吸收效果好等优点，被广泛应用于电磁场DGTD数值模拟算法中。卷积完全匹配层(convolutional perfectly matched layer, CPML)能够有效吸收凋落波，在长时间的模拟计算中保持稳定，被广泛应用于传统的全电磁粒子模拟算法^[34]。DGTD-PIC方法中，CPML辅助变量的推进公式需要求解电磁场分量的空间导数项，导致在计算辅助变量时需要求解电磁场分量的数值通量，会大幅度增加程序实现的复杂度，同时如果数值通量计算不准确，又会导致程序的计算结果不稳定。CFS-PML与CPML具有相同的计算公式，在方程离散的过程中，通过引入特殊表达形式的辅助变量，避免在求解辅助变量时计算电磁场的空间偏导数，数值模拟结果表明该吸收边界对电磁波的吸收效果好，且经过较长步数的模拟计算依然能够保持稳定。下面对DGTD-PIC方法中采用的CFS-PML进行详细的推导，CFS-PML中电磁场方程如下式所示：

式中：$\omega$为角频率；${\varepsilon _{\text{r}}}$为相对介电常数；${\mu _{\text{r}}}$为相对磁导率；${s_\nu } = {\kappa _\nu } + \dfrac{{{\sigma _\nu }}}{{{\text{i}}\omega {\varepsilon _0}}}\;(\nu = x,y,{\textit{z}})$，${\kappa _\nu }$为坐标伸缩系数，${\sigma _\nu }$为PML的电导率。由于方程形式与CPML的表达式相同，因此可按照相同的方式设置参数大小。为便于计算，需要将CFS-PML变换成更为合理的形式。以方程(9)中第四个式子${H_x}$的计算公式为例，将${s_\nu}$的表达式代入${H_x}$的公式，得到如下表达式：

引入辅助变量$P_v^{\text{h}}$和$P_v^{\text{e}}$：

将式(11)代入式(10)，则进一步得到${H_x}$的推进公式，

进一步引入辅助变量${\tilde H_\nu }$和${\tilde E_\nu }$：

将式(13)代入式(12)，则得到辅助变量${\tilde H_x}$的推进公式如下：

通过式(11)和(13)，可进一步得到$P_x^{\text{h}}$的表达式，如式(15)所示：

将式(15)代入式(14)，并进一步化简，则得到${\tilde H_x}$、${\tilde E_{\textit{z}}}$、${\tilde E_y}$、$P_x^{\text{h}}$四个变量之间的关系如下：

通过将式(15)与(16)变换为时域的表达式，即可得到辅助变量${\tilde H_x}$、$P_x^{\mathrm{h}}$在时域的推进公式：

同样可推导磁场相关辅助变量$\tilde {\boldsymbol{H}}$、${{\boldsymbol{P}}^{\text{h}}}$在其他坐标方向的推进公式，并写成如下统一的表达形式：

式中与辅助变量${\tilde H_x}$、$P_x^{\text{h}}$相关的系数表达式如式(19)所示，其他方向的系数可按照坐标轮换的方式$\left( {x \to y,y \to {\textit{z}},{\textit{z}} \to x} \right)$计算得到。

按照电场与磁场之间的对偶原理，可推导得到与电场相关的辅助变量$\tilde {\boldsymbol{E}}$、${{\boldsymbol{P}}^{\text{e}}}$的计算公式：

式(18)和(20)是CFS-PML中电场和磁场以及相关辅助变量的推进公式，在真空区域时，这两个方程直接退化成经典的电磁场推进方程。

3   带电粒子定位及电荷、电流分配

带电粒子与电磁场之间产生非线性相互作用，电磁场作用在带电粒子上，带电粒子的运动满足牛顿-洛伦兹力方程，为保证带电粒子运动时能量守恒，全电磁粒子模拟方法通常采用Boris提出的“半加速-旋转-半加速”方法求解带电粒子运动方程。首先需要对带电粒子定位，确定带电粒子所在位置处的网格编号，通过将网格单元内定义的离散电磁场值插值计算得到带电粒子位置处电磁场值。在进行三维全电磁粒子模拟时，网格单元的数目可达到数百万，为减少数值噪声，往往需要跟踪数百万甚至上亿个带电宏粒子，因此快速定位宏粒子位置能够大幅提高计算速度。在基于FDTD方法的全电磁粒子模拟方法中，通过比较带电粒子坐标位置与结构网格单元的坐标点，能快速得到带电粒子所在的网格编号。采用均匀结构网格时，如果网格的起始位置位于坐标原点，将带电粒子的位置坐标除以网格步长，就能得到带电粒子所在网格的坐标编号以及电荷、电流分配的权重系数。

由于DGTD-PIC方法采用非结构网格，如果对每个网格单元遍历确定带电粒子所在非结构网格单元的编号，会非常耗费计算资源。为实现带电粒子在非结构网格的快速定位，国外学者提出了基于辅助背景网格耦合技术的带电粒子定位方法。在非结构网格所在区域定义一层网格尺寸远大于非结构网格尺寸的结构网格，通过比较结构网格单元的坐标与非结构网格单元的坐标构建结构网格与非结构网格之间的拓扑关系，如果非结构网格单元与某个结构网格单元所在区域重合，则表明该非结构网格从属于该结构网格单元。DGTD-PIC方法在确定带电粒子所在非结构网格单位编号时，首先快速确定带电粒子所在结构网格单元的编号，获取从属于该结构网格单元的所有非结构网格单元编号，通过判断带电粒子所在位置是否位于某个从属的非结构网格单元内，确定带电粒子所在非结构网格编号，该方法无须遍历所有非结构网格单元；此外，通过建立非结构网格单元之间的拓扑结构关系能进一步提高带电粒子定位的速度，在构建非结构网格的拓扑关系时，保留相邻网格的拓扑信息，实现网格与相邻网格之间的快速检索，当带电粒子运动导致所在非结构网格单元编号改变时，从带电粒子运动前所在网格的相邻网格开始搜索，这样可大幅提高DGTD-PIC方法定位带电粒子的速度。

带电粒子运动产生的电流作为Maxwell方程的电流源项，实现了带电粒子对电磁场的耦合作用。电磁场方程的源项$\rho$、${\boldsymbol{J}}$须从连续的带电粒子位置插值到离散的非结构网格位置得到。在采用结构网格的全电磁粒子模拟方法中，通常采用线性插值的电荷和电流分配方式，双线性插值的电流分配方法会导致电荷不守恒，该电流作用下计算得到的下一个时间步长的电场会偏离高斯定律，需要对电场分布作一定的调整，即散度校正，使之满足高斯定律。散度校正的方法按照实现方式主要有两种：Langdon-Marder校正和Boris校正方式。国外学者提出的Zigzags等算法可大幅提高计算速度，同时能够保证电流分配过程中电荷守恒^[35-36]。在DGTD-PIC方法中，电流分配时同样需要保证电场满足高斯定律，由于该方法采用了非结构网格，无法借鉴传统全电磁粒子模拟方法中线性插值以及Zigzags等相关算法。

下面以一个非结构网格常用的三角网格单元为例说明DGTD-PIC方法所采用的电荷守恒的电荷和电流分配方式，电荷分配采用面积权重的方式进行插值分配，如图2所示。其中图2(a)为带电粒子的电荷分配示意图，图中三角网格单元三个顶点的局部编号分别标示为0、1和2，带电粒子电荷量为$q$、所在位置为${P_{\text{s}}}$，该点与三角网格单元的连线将三角网格单元划分为三个小三角形，这三个小三角形的面积占该三角单元面积的比例分别为$\lambda _0^{\text{s}}$、$\lambda _1^{\text{s}}$和$\lambda _2^{\text{s}}$，且$\lambda _0^{\text{s}} + \lambda _1^{\text{s}} + \lambda _2^{\text{s}} = 1$，按照面积权重分配给三角网格单元三个顶点的电荷量分别为$q_0^{\text{s}} = q \cdot \lambda _0^{\text{s}}$、$q_1^{\text{s}} = q \cdot \lambda _1^{\text{s}}$和$q_2^{\text{s}} = q \cdot \lambda _2^{\text{s}}$。以第0号顶点为例，$q_0^{\text{s}}$为0号顶点所分配的电荷量，${P_{\text{s}}}$点与三角网格单元1号顶点和2号顶点所构成的三角形占三角网格单元的面积比例为$\lambda _0^{\text{s}}$。从图2(a)可以看出，当带电粒子与某个顶点的距离越近，则该顶点所能分配到的电荷量就越大。当经过时间步长$\Delta t$，带电粒子从${P_{\text{s}}}$运动到${P_{\text{e}}}$，如图2(b)所示，带电粒子在三个网格顶点所对应的面积分配权重更新为$\lambda _0^{\text{e}}$、$\lambda _1^{\text{e}}$和$\lambda _2^{\text{e}}$，则三角网格单元三个顶点分配的电荷量分别为$q_0^{\text{e}} = q \cdot \lambda _0^{\text{e}}$、$q_1^{\text{e}} = q \cdot \lambda _1^{\text{e}}$和$q_2^{\text{e}} = q \cdot \lambda _2^{\text{e}}$。沿着网格边线0-1、0-2、1-2方向分别产生电流${{\boldsymbol{I}}_0}$、${{\boldsymbol{I}}_1}$和${{\boldsymbol{I}}_2}$。该带电粒子运动产生的电流幅值计算公式为：

式中：$\Delta t$为时间步长；${i_0}$、${i_1}$、${i_2}$为三条边上电流${{\boldsymbol{I}}_0}$、${{\boldsymbol{I}}_1}$、${{\boldsymbol{I}}_2}$的幅值，如图2(b)所示。该电流分配方式能保证计算过程中电荷守恒，将所有宏粒子产生的电流在所在网格进行叠加，就能得到所有宏粒子运动在网格边上产生的电流。

图  2  带电粒子电荷以及电流密度分配

Fig.  2  Charge allocation and current density allocation

下载: 全尺寸图片 幻灯片

4   DGTD-PIC模拟软件计算结果

国内外学者研制了DGTD-PIC软件，并利用此软件对激光与带电粒子相互作用进行了数值模拟研究^[37]，验证了DGTD-PIC方法能够应用于模拟电磁场与带电粒子之间的非线性相互作用。一个电量为−$1.6 \times {10^{ - 15}}$ C的宏粒子初始时刻为静止状态，在时间上呈高斯分布的平面电磁波的激励下，其中平面电磁波的电场幅值为10 GV/m，脉冲持续时间为15 fs，脉冲的中心波长为1 μm，计算结果如图3所示。其中图3(a)为采用不同阶数基函数情况下计算结果与理论值的对比，可以看出软件可以很好模拟带电粒子产生的电场，同时随着基函数的阶速增加，计算精度大幅提高；图3(b)为计算误差随基函数的阶数和时间步长的变化，可以看出采用高阶基函数和缩减时间步长能够进一步提高计算的精度；图3(c)为模拟得到的辐射电场分布图。

图  3  激光加速单粒子辐射场模拟结果^[37]

Fig.  3  Results for the radiated field of one particle accelerated by a laser beam^[37]

下载: 全尺寸图片 幻灯片

DGTD-PIC软件也被用于模拟真空电子器件，回旋振荡器作为一种能够产生高频电磁波的真空电子器件，德国斯图加特大学和卡尔斯鲁厄理工大学的学者采用研制的DGTD-PIC软件对回旋管振荡器进行了数值模拟研究^[33]。回旋管的电子束电压为79 kV、电流为15 A，回旋电子束的横纵速度比值为1.5，外部约束磁场为1.16 T。采用四面体单元对三维模型进行网格剖分，生成的面体网格剖视图如图4所示，模拟得到y=0处x-z平面${E_y}$剖面图如图5所示，在z=50 mm处x-y平面${E_\varphi }$剖面图如图6所示，z=0处x-y平面电子${v_ \bot }$剖面图如图7所示。可以看出，电子束横向速度实现了调制，器件在电磁场与回旋电子束的高频互作用区的模式分布稳定，且与理论设计符合较好。

图  4  四面体网格的剖视图^[33]

Fig.  4  Cut-open view of tetrahedral grid^[33]

下载: 全尺寸图片 幻灯片

图  5  y=0处x-z平面${{\boldsymbol{E}}_{\boldsymbol{y}}}$剖面图^[33]

Fig.  5  ${{{{\boldsymbol{E}}_{\boldsymbol{y}}}}}$-field slice in the x-z plane at y = 0^[33]

下载: 全尺寸图片 幻灯片

图  6  z=50 mm处x-y平面${E_\varphi }$剖面图^[33]

Fig.  6  ${{\boldsymbol{E}}_{\boldsymbol{\varphi}} }$-field slice in the x-y plane at z = 50 mm^[33]

下载: 全尺寸图片 幻灯片

图  7  z=0处x-y平面${v_ \bot }$剖面图^[33]

Fig.  7  The tangential velocity${{\boldsymbol{v}}_ {\boldsymbol{\bot }}}$slice of particles in the x-y plane at z = 0^[33]

下载: 全尺寸图片 幻灯片

5   展　望

DGTD-PIC方法采用有限元网格生成技术，实现了复杂结构内电磁场与带电粒子之间非线性相互作用的准确计算。目前DG方法已广泛应用于弹性力学问题、热传导问题以及可压缩流体等问题的研究，通过将这些方程与全电磁粒子模拟方程耦合求解，能够实现强电磁场相关多物理场问题的耦合计算。相比传统的基于FDTD方法的全电磁粒子模拟方法，DGTD-PIC方法在提高计算精度的同时，计算量会大幅增加，主要包括电磁场的高阶算法以及带电粒子定位与推进，因此需要进一步开展深入的研究，如研制高效的数学计算模型，利用大规模并行技术，提高DGTD-PIC方法的计算性能；研制基于超算的大规模并行计算软件，针对实际工程需求开展相关的应用研究^[33,37-38]。