0   引　言

全电磁粒子模拟方法通过模拟带电粒子与电磁场之间的非线性相互作用过程，能够给出电场、磁场随时间演化过程及带电粒子的运动轨迹，广泛应用于复杂结构真空电子器件的数值模拟研究及参数的优化设计^[1-3]。传统的全电磁粒子模拟方法采用时域有限差分(finite-difference time-domain, FDTD)法求解Maxwell方程，为准确模拟小尺寸结构的影响，通常采用空间步长较小的阶梯网格逼近计算模型的复杂边界，导致网格数量剧增，影响计算效率^[4]。采用交替隐式FDTD方法^[5]或显隐式FDTD方法^[6-7]可加速电磁场的求解。基于计算机辅助设计(computer-aided design，CAD)的结构网格剖分方法能生成共形结构网格，实现复杂模型在离散网格中重构^[8-10]，采用共形结构网格的FDTD方法能更准确描述计算模型中复杂边界的影响，大幅提高全电磁粒子模拟方法的计算效率和计算精度。国内外学者将CAD引入粒子模拟方法中^[11-16]，采用共形结构网格对计算模型剖分，并在共形结构网格的基础上构建了更高精度的电磁场计算模型及带电粒子发射模型。此外，研制了基于共形结构网格的全电磁粒子模拟软件，并将其成功应用于真空电子器件的数值模拟研究，为研制工作频率更高、高频结构更为复杂的真空电子器件提供了重要的设计工具^[17-18]。

电磁场时域有限元(fintte-element time-domain，FETD)方法对几何模型采用非结构网格剖分，能够更好地描述复杂边界^[19-21]。传统的电磁场FETD方法需要对矩阵求逆，当计算模型非常大时，所需要的计算资源非常庞大，导致数值计算效率低下，近年来随着稀疏矩阵近似求逆(sparse approximate inverse，SPAI)数值技术的发展，大幅提高了FETD求解矩阵的规模^[22-24]。随着电磁场FETD方法的深入研究，FETD方法实现了与FDTD方法相同的显式时间推进格式，同时成熟的非结构网格剖分技术也为电磁场FETD方法的应用提供了重要的技术支撑，将电磁场FETD方法应用于全电磁粒子模技术成为国内外研究的热点^[25-28]。

本文详细介绍了国内外在复杂结构共形全电磁粒子模拟研究方面的最新研究成果^[29-31]，主要包括基于共形FDTD方法的全电磁粒子模拟技术及基于FETD方法的全电磁粒子模拟技术。在基于共形FDTD方法的全电磁粒子模拟技术方面主要介绍了该技术涉及的关键技术，包括共形结构网格生成技术、共形FDTD技术及共形电子束发射技术。在基于FETD方法的全电磁粒子模拟技术方面主要介绍了该方法所采用的基函数、电磁场方程的求解以及电磁场与带电粒子耦合求解过程中涉及的电荷守恒和高斯守恒定律等关键技术。

1   基于共形FDTD方法的全电磁粒子模拟技术

1.1   复杂模型共形结构网格边界重构技术

正交网格线定义的射线、网格面定义的平面通过与复杂三维模型进行求交运算，确定模型的边界，实现模型边界的准确重构。实现的原理性框架如图1所示，基于几何建模库构建复杂几何模型，设置模型参数，根据模型参数以及网格尺寸生成全局结构网格，全局结构网格的正交网格线、网格面与几何模型求交运算确定复杂模型的边界。

图  1  生成边界数据框架图

Fig.  1  Frame for boundary data generation

下载: 全尺寸图片 幻灯片

采用射线追踪法确定模型在正交网格线上的边界。将正交网格线构建为射线，通过求解射线与几何体各个面的交点确定模型在正交网格线上的边界，交点可提供三个方面的信息：1)交点的位置信息；2)交点处射线的方向信息，表征交点位置处射线方向是指向几何体内部或几何体外部；3)交点所在模型的几何体及几何面的编号，通过上述信息准确定义边界点的位置及属性(见图2)。

图  2  网格线上边界点属性

Fig.  2  Properties of boundary point on the grid line

下载: 全尺寸图片 幻灯片

1.2   共形结构网格生成技术

基于射线追踪法求得的模型边界点信息，并利用几何模型的边界与几何模型的拓扑关系，通过在结构网格单元中引入局部非结构网格的思想，采用局部填充(partically filled cell，PFC)技术实现几何模型在正交结构网格内的共形描述。共形结构网格体系中网格元(网格顶点、网格边、网格面和网格体)的定义如图3所示。几何模型与正交结构网格相交会产生边界信息，几何模型在结构网格的边界采用网格边上的边界点、网格面上的边界线以及网格元胞内的截面描述，如图4所示。可以看出网格边上的边界点将网格边分割为多个线元；模型在网格面上的边界线将网格面分割为多个小面元。如图5所示，模型在网格体内的边界面将网格体构成的网格元胞分割为多个体元，这些网格元胞上的局部边界将网格元胞局部地分割成非规则网格体元。通过构建网格边、网格面及网格体之间的拓扑关系，并依据模型边界的连接特性和电磁场计算公式所规定的定向拓扑关系构建模型在离散网格元上的不同类型网格元之间的层次结构，可得到计算模型的共形结构网格。

图  3  网格元的定义^[25]

Fig.  3  Definition of grid element^[25]

下载: 全尺寸图片 幻灯片

图  4  结构网格面上的非结构网格面^[25]

Fig.  4  Unstructured mesh surface on structured grid surface^[25]

下载: 全尺寸图片 幻灯片

图  5  网格单元上的斜截面^[25]

Fig.  5  Oblique cross-section in grid element^[25]

下载: 全尺寸图片 幻灯片

1.3   电磁场共形FDTD技术

基于共形结构网格系统中拓扑数据结构，将扩展网格技术(enlarged cell technique, ECT)与FDTD方法结合^[32-33]，避免了传统共形FDTD方法中时间步长缩短问题。基于ECT的FDTD方法(ECT conformal FDTD，ECT-CFDTD)将所有离散方程转换为局部有限积分方程，建立局部与全局结构自洽一致的ECT-CFDTD有限积分方程，避免了大型矩阵求解，同时保持了不同区域推进方程采用计算公式的一致性。

按照电场环路积分的定义，任一面元(图6)上有向边界上的电场环路积分为

图  6  ECT-CFDTD方法中小面元

Fig.  6  Small face in ECT-CFDTD method

下载: 全尺寸图片 幻灯片

式中：V_p为标号为p的初始网格面元的电场环路积分；线元L_m为面元S_p的边界；$\partial {S _p}$为S_p的有向边界；$r_p^{(m)}$为两边界之间的相对定向关系；E_m为定义在L_m上的电场值。$r_p^{(m)}$的取值为$\pm 1$，当线元L_m的自然定向(坐标轴定向)与S_p边界的定向相反时，$r_p^{(m)} = - 1$，反之$r_p^{(m)} = + 1$。ECT-CFTD方法实现的电压计算公式为：

式中，${R_{p,m}}$为面元${S _m}$的电场环路积分借给面元${S _p}$的压降积分的比例。而对于“被借用”(intruded)的面元${S _p}$上相应的值则通过求平均的方法给出：

从式(4)～(5)可看出“被借用”面元的电场环路积分由该面元的原电场环路积分与“借用”该面元电场环路积分的小面元被放大后的电场环路积分按比例进行平均处理。

对应的小面元和“被借用”的大面元上磁场的推进分别为：

按照相同的方式，ECT-CFDTD方法中电场的推进为：

式中：$\partial {\tilde S _p}$为${\tilde S _p}$的有向边界；$\tilde V_p^{\left(n + \frac{1}{2}\right)}$为标号为p的对偶网格面元的磁场环路积分。

1.4   共形电子束发射技术

全电磁粒子模拟方法须准确模拟复杂边界产生带电粒子的物理过程，对计算模型的电子发射阴极表面进行三角网格剖分，每个小三角形单元构成一个独立的电子发射阴极，如图7所示。采用ECT-CFDTD方法准确给出边界处电磁场值后，基于阴极发射面的三角网格建立阴极电子发射模型，解决复杂阴极表面的电子发射问题。

图  7  采用SCL定律计算阴极表面发射电流模型

Fig.  7  Surface emission current of cathode by SCL law

下载: 全尺寸图片 幻灯片

真空电子器件一般采用高电压产生强流相对论电子束，当在器件的阴阳极间加载高电压时，阴极表面局部尖端会由于强电场而产生电子，阴极表面在电流欧姆加热作用下发生气化，形成覆盖在电极表面的等离子体层，并由该等离子体层继续发射电子。由于形成的等离子体密度非常大，具有很强的带电粒子发射能力，最终在阴极表面强电场作用下形成空间电荷限制流，可采用Child-Langmuir公式描述，如式(10)所示：

式中：$\varepsilon$为介电常数；$\eta$为电荷的荷质比；$\phi$为阴阳极之间的电势；$d$为阴阳极之间的间距。

式(10)给出的空间电荷限制流是在无限大平板的理想状态下推导的计算公式。实际应用中，由于阴极形状多种多样，并且二极管可能存在多个电极，用于束流的聚焦，使得电子发射的物理过程不再满足无限大平板的理想状态，导致阴阳极之间的电荷和电势分布无法由理论公式直接给出。全电磁粒子模拟方法通过将复杂的三维阴极电子发射面转变成一系列一维平板二极管模型，如图7所示。通过在阴极面附近设置一个面，假定该面到阴极面的电势在没有空间电荷作用下呈线性变化，并假定该面与阴极面所构成的区域构成一维平板二极管模型，阴极面产生的电流密度遵循Child-Langmuir定律。

令阴极小三角面元的面积为$\Delta S$，定义单个宏粒子的电量为q，根据式(10)计算得到三角面元需要产生的电流密度| j |，发射前通过发射面元法向的电流密度为$|{{\boldsymbol{j}}_{{\text{old}}}}|$，可计算得到当前时间所产生的电流密度为

根据电荷守恒定理，可计算得到当前时刻该三角面元需发射的宏粒子数为

粒子运动的牛顿-洛伦兹方程可采用Boris推进方法求解^[2-3]。开发边界的截断可采用卷积完全匹配层实现^[34]。

2   基于FETD方法的全电磁粒子模拟技术

在FETD方法中，非结构网格单元上物理量的分布采用基函数的线性组合描述。由于在全电磁粒子模拟方法中，电磁场方程与带电粒子推进公式耦合计算，电荷以及电流的分配方式会导致电场和电荷不满足高斯定律。因此，在基函数的选取及电流分配方面需慎重，确保计算过程中算法能满足电荷守恒定律、电场满足高斯定律^[35-38]。

2.1   电磁场、电流密度基函数

在基于FETD方法的全电磁粒子模拟中，电场、磁场、电流以及电荷均可采用Whitney形式的基函数描述。电场和磁场采用时间交错的方式推进，在整数时间步长时推进电场的计算公式，在半整数时间步长时推进磁场的计算公式。$n$时刻电场强度${\boldsymbol{E}}^{({{n}})}$及$n+\dfrac{1}{2}$时刻磁感应强度${\boldsymbol{B}}^{\left(n+\frac{1}{2}\right)}$的插值函数为：

式中：$N_{1}$为非结构网格单元棱边的数目；$N_{2}$为非结构网格单元小面元的数目；$e_{i}$、$b_{i}$分别为电场强度和磁感应强度离散后与基函数相对应的系数；${\boldsymbol{W}}_{i}^{1}(r)$为1阶Whitney基函数；${\boldsymbol{W}}_{i}^{2}(r)$为2阶Whitney基函数。电流密度和电荷密度的定义如下：

式中：$N_{0}$为非结构网格单元上网格顶点的数目；${{W}}_{i}^{0}({r})$为0阶Whitney基函数；$j_i$、$q_i$分别为电流密度和电荷密度采用基函数离散后与基函数相对应的系数。

2.2   电场和磁感应强度的推进公式

根据式(13)～(15)定义的电场强度、磁感应强度及电流的基函数，同时利用对偶网格原理，离散后的电场强度和磁感应强度的推进公式为：

式中：$\boldsymbol{b}=\left[b_{1}, b_{2}, \cdots, b_{N_{2}}\right]^{\mathrm{T}}$、$\boldsymbol{e}=\left[e_{1}, e_{2}, \cdots, e_{N_1}\right]^{\mathrm{T}}$分别为离散后的磁感应强度、电场强度所对应基函数的系数向量；$\boldsymbol{j}=\left[j_{1}, j_{2}, \cdots, j_{N_{1}}\right]^{\rm{T}}$为棱边上的电流密度向量；$\left[\star_{\varepsilon }\right]$、$\left[\star_{\mu^{-1}}\right]$矩阵中参数的表达式为：

${\boldsymbol{C}}$为电场强度和磁感应强度在非结构网格中旋度关联矩阵；$\tilde{{\boldsymbol{C}}}$为电场强度和磁感应强度在非结构网格所对应的对偶网格上的旋度关联矩阵，且$\tilde{{\boldsymbol{C}}}={\boldsymbol{C}}^{\mathrm{T}}$，即两个矩阵互为转置关系。${\boldsymbol{C}}$矩阵中每个元素的值可以为–1、0和1，该值的取值与定义的电场强度及磁感应强度基函数的方向、电场强度及磁感应强度之间的拓扑结构有关。以图8所示的非结构网格为例，电场强度和电流密度定义在非结构网格的棱边上，共24个数值；磁感应强度定义在三角形网格单元面上，共13个数值；电荷定义在网格顶点上，共12个数值。以图8中编号为6的三角形网格单元面为例，编号为8的棱边定义的电场强度方向与磁感应强度的旋度方向相反，${C}_{6,8}=-1$；编号为9的棱边定义的电场强度方向与电场强度的旋度方向相同，${C}_{6,9}=1$；编号为20的棱边定义的电场强度方向与磁感应强度的旋度方向相反，${C}_{6,20}=-1$，其他的${C}_{6, k}=0$。

图  8  非结构网格示例^[30]

Fig.  8  Example (primal) unstructured mesh^[30]

下载: 全尺寸图片 幻灯片

2.3   算法中电荷守恒和高斯守恒定律

FETD全电磁粒子模拟方法在计算电流时需要保证电荷守恒，即满足：

式中，${\boldsymbol{S}}$为电场强度或电流与电荷在非结构网格上的散度关联矩阵，而$\tilde{{\boldsymbol{S}}}$为电场强度或电流与电荷在非结构网格所对应的对偶网格上的散度关联矩阵，且$\tilde{\boldsymbol{S}}={\boldsymbol{S}}^{\mathrm{T}}$，即两个矩阵互为转置关系。${\boldsymbol{S}}$矩阵中每个元素的值可以为–1、0或1，该值的取值与定义的电场强度及电流基函数的方向、电场强度及电荷密度之间拓扑结构有关。在图8所示的非结构网格中，以编号为5的三角形网格顶点为例，定义电流以从顶点沿着棱边流出的方向为正，可得到$\tilde{S}_{5, 7}=1$、$\tilde{S}_{5,10}=-1$、$\tilde{S}_{5,12}=1$、$\tilde{S}_{5,13}=1$、$\tilde{S}_{5,14}=1$、$\tilde{S}_{5,15}=1$，其余的$\tilde{S}_{5, k}=0$。

将式(21)在时间上离散，可得到：

下面以一个带电粒子的运动为例证明采用式(15)、(16)所示的电流以及电荷基函数时，式(22)成立，即满足电荷守恒定律。假设带电粒子所带电荷为$q$，从$n$时刻的位置$r_{p}^{(n)}$，经过了一个时间步长$\Delta t$，在$n+1$时刻运动到相同三角形的其他位置$r_{p}^{(n+1)}$。三角形三个网格顶点的编号分别为${i}$，$l$，$k$。带电粒子运动产生沿$i 、k$顶点所构成棱边的电流如式(23)所示，沿着$i、 l$顶点所构成棱边的电流如式(24)所示。

利用$\lambda_i\left(r_p^{(n)}\right) + \lambda_l\left(r_p^{(n)}\right) + \lambda_k\left(r_p^{(n)}\right)=1$，$\lambda_i\left(r_p^{(n+1)}\right) + \lambda_l\left(r_p^{(n+1)}\right) + \lambda_k\left(r_p^{(n+1)}\right)=1$，可得到电荷$q$从$r_{p}^{(n)}$运动到$r_{p}^{(n+1)}$时，该电荷在${i}$号顶点产生的流出电流为

式(22)中${i}$号顶点电荷的变化如式(26)所示：

从式(25)和(26)可以看出式(22)成立，表明采用式(15)、(16)给出的电流以及电荷基函数能够保证计算过程中的电荷守恒。

将式(18)两边同时乘以${{\boldsymbol{S}}}$并化简，得到：

由于对偶网格中旋度和散度矩阵满足$\tilde{{\boldsymbol{S}}} \cdot {\boldsymbol{C}}^{\mathrm{T}}=0$，将式(22)代入式(27)，可进一步简化为

式(28)为如下方程的离散形式：

式(29)表明如果初始时刻电磁场方程满足高斯定律守恒$\tilde{\boldsymbol{S}}\cdot \left[\star_{t}\right]\cdot {\boldsymbol{e}}^{(0)}=q^{(0)}$，则在后续推进过程中方程仍可保证电场和电荷满足高斯守恒定律。

3   复杂真空电子器件的全电磁粒子模拟

3.1   基于共形FDTD方法的全电磁粒子模拟软件计算结果

西北核技术研究所研制了基于ECT-CFDTD-PIC技术的三维并行粒子模拟软件^[13]，软件集成了CAD功能，能够进行三维复杂几何模型建模，同时采用射线追踪算法得到正交直角网格线与真空电子器件模型的交点和交线等模型边界信息，通过边界信息及网格点构建能够准确描述三维复杂模型的共形结构网格。采用该软件在超算平台上对磁绝缘线振荡器进行数值模拟研究，所建三维模型如图9所示；对几何模型进行了三维共形结构网格剖分，磁绝缘线振荡器模型的三个正交剖面的共形结构网格如图10所示。软件可对模型的部分波-束互作用结构进行网格细节的显示，如图11所示，其中红色部分为位于真空区域的正交网格线，可以看出共形结构网格剖分模块能准确给出三维复杂模型与正交网格线的交点。在国产超算平台上模拟得到的磁绝缘线振荡器中带电粒子相空间图如图12所示，可以看出，电子束得到了较好的调制。模拟得到的电场随时间变化如图13所示。模拟结果表明采用共形正交网格的全电磁粒子模拟软件的计算结果与2.5维粒子模拟软件UNIPIC计算结果吻合较好^[39-42]。

图  9  磁绝缘线振荡器的三维模型^[13]

Fig.  9  3D model of mangnetically insulated transmission line oscillator (MILO)^[13]

下载: 全尺寸图片 幻灯片

图  10  磁绝缘线振荡器剖面共形结构网格示意图^[13]

Fig.  10  Conformal grids of three different sections of MILO^[13]

下载: 全尺寸图片 幻灯片

图  11  共形结构网格的局部示意图^[13]

Fig.  11  Details of the conformal structured grid^[13]

下载: 全尺寸图片 幻灯片

图  12  磁绝缘线振荡器内带电粒子的速度空间分布图^[13]

Fig.  12  Velocity space of particles inside MILO^[13]

下载: 全尺寸图片 幻灯片

图  13  电场随时间的变化^[13]

Fig.  13  Electric fields vs. time^[13]

下载: 全尺寸图片 幻灯片

3.2   基于FETD方法的全电磁粒子模拟软件计算结果

俄亥俄州立大学基于FETD方法研制了2.5维全电磁粒子模拟软件^[35]，并采用该软件模拟了相对论返波管。返波管三维模型如图14所示，它由阴极区、电磁波与电子束产生非线性相互作用的慢波结构区以及同轴输出结构等部分组成。由于三维模型为轴对称结构，且器件工作在轴对称模式，因此只须模拟其中一个剖面，模拟的二维剖面结构如图15所示，电子束在外部引导磁场的约束下轰击到收集极。将模型进一步放大可观察到慢波结构区域内非结构网格分布及电子束的调制情况，如图16所示。可以看出，非结构网格可以很好地描述慢波结构区域内存在的细小结构。模拟得到的剖面电场分布如图17所示，模拟得到的电子束相空间如图18所示，模拟得到的电场随时间的变化如图19所示，电场的频谱分布如图20所示。可以看出，电子束沿着轴向的调制稳定，同时最后输出电场的频谱较为单一。

图  14  返波管模型图^[35]

Fig.  14  Schematic of a backward wave oscillator (BWO)^[35]

下载: 全尺寸图片 幻灯片

图  15  相对论返波管模型以及粒子分布图^[35]

Fig.  15  Schematic of relativistic BWO (RBWO) with instant particle distribution snapshots^[35]

下载: 全尺寸图片 幻灯片

图  16  相对论返波管慢波结构最右段的放大图像^[35]

Fig.  16  A zoomed-in region of rightmost corrugations of RBWO^[35]

下载: 全尺寸图片 幻灯片

图  17  模拟得到的稳定自洽场分布^[35]

Fig.  17  Snapshot of steady-state self-fields^[35]

下载: 全尺寸图片 幻灯片

图  18  电子束的3维速度空间分布图^[35]

Fig.  18  3D velocity plots for an electron beam^[35]

下载: 全尺寸图片 幻灯片

图  19  产生的电场随时间的变化^[35]

Fig.  19  Output electric field vs. time^[35]

下载: 全尺寸图片 幻灯片

图  20  电场的频谱分布图^[35]

Fig.  20  Distribution of electric field frequency spectrum^[35]

下载: 全尺寸图片 幻灯片

4   展　望

复杂结构内强电磁现象的数值模拟研究需要准确模拟电磁场与带电粒子之间的非线性相互作用。借鉴传统全电磁粒子模拟技术，采用共形结构网格及在此基础上发展的共形FDTD方法、共形电子束发射等相关数值模拟技术，形成了基于共形结构网格的全电磁粒子模拟方法。共形和非共形全电磁粒子模拟方法构成一个整体，可用于模拟复杂结构真空电子高功率微波器件和太赫兹器件设计及参数优化^[40-44]、复杂结构在高能射线作用下产生系统电磁脉冲^[45-47]、强场击穿产生等离子体及气体放电等物理^[48]、高空核爆炸等离子体环境^[49]等。电磁场FETD方法采用非结构网格描述复杂模型，采用Whitney基函数描述网格单元内电磁场分布，同时采用SPAI方法求解方程，可大幅度降低FETD方法的计算量，同时有效解决了FETD方法中电流、电荷分配引起的电荷不守恒问题，确保了计算过程中电场和电荷满足高斯守恒定律，实现了基于FETD方法的全电磁粒子模拟的长时间稳定计算。在该算法基础上发展新的物理计算模块及高效算法，利用大规模并行技术，可进一步提高基于FETD方法的全电磁粒子模拟方法的计算性能。近年来，电磁场时域间断有限元数值模拟技术不断发展，该方法不需要对大型矩阵求逆，基于时域间断有限元的全电磁粒子模拟技术有望能更好地解决电大尺寸复杂结构中等离子体现象模拟难题^[50-51]。