技术领域
[0001] 本
发明属于粒子模拟(Particle-in-cell，简写为PIC)的数值模拟领域，具体涉及一种应用于PIC的推动粒子运动有限元算法。
背景技术
[0002] PIC方法是一种被广泛应用于带电粒子与电
磁场相互作用物理问题中的数值模拟方法，它通过
跟踪大量带电粒子在外加及自洽
电磁场中的运动并统计平均而得到宏观特性及运动规律。经过几十年的发展，PIC模拟方法已经成为研究带电粒子与电磁场相互作用物理问题的一种强有
力的数值手段，广泛应用于带电粒子与电磁场相互作用所涉及的许多领域，如
磁约束聚变
等离子体、
惯性约束聚变等离子体、核爆、空间等离子体、人造等离子体(包括
电子枪、离子源等)、电推进、自由电子激光以及电
真空器件等。
[0003] PIC求解的核心步骤如下：
[0004] 1、电磁场求解，即通过求解电磁场所满足的麦克斯韦方程组(静电模型下退化为泊松方程)，得出所有网格点上的电磁场；
[0005] 2、粒子受力求解，即通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在
位置处
电场，然后求解受力；
[0006] 3、推动粒子运动，即通过求解离散粒子运动方程，更新粒子的速度及位置等运动信息，进而更新粒子所属网格；
[0007] 4、源的分配，即根据粒子所在的位置求得其对周
围网格点电荷与
电流的贡献，然后将所有粒子对网格点上的电荷与电流贡献累加得到网格点上的电荷
密度与电流密度；
[0008] 不断循环如上过程，直到计算结果收敛或人为设置的时间为止。
[0009] 其中步骤3推动粒子运动是PIC的核心步骤之一，该步骤的精确且高效的求解对于PIC的整体求解
精度和效率的控制十分重要。截止目前为止，PIC中推动粒子运动主要有两种方法，分别是有限差分(FD)方法和嵌入式有限元(IFE)方法。
[0010] FD方法：在PIC推动粒子运动的应用中，FD方法是通过采用结构化网格离散求解区域的，因此求解离散粒子运动方程和判断粒子所属网格等算法形式简单、易于理解，但是在PIC推动粒子运动的应用中存在如下缺点：
[0011] 1、FD方法采用的是由
正交线划分而成的结构化网格，对于复杂曲形边界的拟合较差，从而使得复杂曲形边界附近的推动粒子运动的求解精度较低；
[0012] 2、由于FD方法对网格尺寸均匀性的要求比较高，因此受限于模拟系统中细小物理结构的限制，必须划分足够小的网格才能满足计算精度要求，从而使得总网格数巨大，而模拟粒子数目正比于总网格数，这就导致推动粒子运动求解的计算量十分庞大；
[0013] 3、FD方法严格受限于数值
稳定性条件的限制，即在对PIC推动粒子运动的数值模拟中，如果空间网格尺寸很小，则时间步长也会随之取的很小，这会进一步增加对推动粒子运动时间循环求解的FD数值模拟负担。
[0014] 针对FD方法在PIC推动粒子运动应用中出现的不足，Kafafy和Wang在2003年，提出了可应用于PIC推动粒子运动中的IFE方法。
[0015] IFE方法：在PIC推动粒子运动的应用中，IFE方法通过采用侵入式非结构化网格离散求解区域。侵入式非结构化网格划分情况如图1所示，可以看出此网格划分相当于有2重网格，其中1重网格是结构化网格，2重网格是将1重网格中的每一个结构化网格进一步划分成五个四面体的侵入式非结构化网格。IFE方法在PIC的应用中，场求解采用第2重侵入式非结构网格，而推动粒子运动是在第1重结构化网格内进行的，因此IFE方法并没有解决FD方法在PIC推动粒子运动应用中的缺点。
发明内容
[0016] 针对上述存在的问题和不足，为解决FD和IFE方法在推动粒子运动求解中对边界匹配度不高、求解精度不高以及数值模拟负担大的问题，本发明提供了一种应用于PIC的推动粒子运动有限元(FEM)算法。具体技术方案如下：
[0017] 步骤1、电磁场求解。
[0018] 通过求解电磁场所满足的离散麦克斯韦方程组(静电模型下退化为静电离散泊松方程)，得出所有网格点上的电磁场；
[0019] 步骤2、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力。
[0020] 步骤3、推动粒子运动
[0021] 采用全局非结构化网格，如图2所示。通过求解离散粒子运动方程，更新粒子的速度及位置等运动信息，进而更新粒子所属网格。详细的求解过程分为如下两个步骤：
[0022] 步骤3.1、更新粒子的速度和位置
[0023] 在直
角坐标系下，考虑相对论效应得出的粒子的运动方程如(1)式所示：
[0024]
[0025] 其中：上标一个点表示该变量对时间的一次导数，两个点表示该变量对时间的二次导数；η为粒子的荷质比；γ为相对论因子；Ex、Ey、Ez、Bx、By、Bz为粒子所在位置处的电磁场分量，由步骤2求得；
[0026] 可利用Boris或Runge-Kutta等方法求解式(1)，更新粒子的速度和位置。
[0027] 步骤3.2、更新粒子所属网格
[0028] 每次更新粒子的速度和位置过后更新粒子所属网格是PIC中不可缺少的环节，只有更新粒子所属网格，才能进行源的计算(即步骤4)、网格内发生的特征物理过程的模
块添加以及相关数值诊断等等。
[0029] 由于PIC中的模拟粒子数目巨大且此时采用的是非结构化网格(并非像结构化网格情形下可直接得出粒子所属网格)，因此对所有粒子所属网格的计算是影响整体PIC计算精度及效率的关键步骤。为此，本发明提出一种基于非结构化网格的粒子快速
定位算法，用于准确高效地计算粒子所属网格，具体实现方式如下：
[0030] 1)在第一步时间步长时，根据粒子出现位置或发射面设定，初始化粒子所处网格编号。
[0031] 2)从第二步时间步长开始，首先判断粒子是否仍然在当前网格内。如果在，则结束搜索；
[0032] 如果不在，则进入3)。
[0033] 3)连线前一时刻粒子所在网格PreId的
重心P1和当前时刻粒子所在位置P2，得到线段P1P2，并求出线段P1P2与PreId的相交面FaceId，同时得到与PreId共用该相交面FaceId的新网格，将新网格设为PreId，回到2)。
[0034] 对于总是选择网格的重心而不是其他网格点作为计算点的主要原因有三个：1、重心一定位于网格内部；2、重心计算方便，只需要将重心所在网格的各个
顶点的x,y,z坐标分量分别求和除以网格顶点总数即可得到该重心坐标的x,y,z方向分量；3、避免粒子在特殊位置时程序判断时出现死循环，如某一网格顶点正好处于线段P1P2上时。
[0035] 上述基于非结构化网格的粒子快速定位算法，结合粒子上一时间步长所处网格，按深度优先算法进行搜索，按照明确的方向从起始位置开始，逐步逼近粒子所在网格，最终实现准确快速定位。
[0036] 算法具体示意图如图3所示。为了更好地描述该粒子快速定位算法，示意图采用二维表现形式，在图中，二维的三角形网格相当于三维的四面体网格，二维的线相当于三维的面。
[0038] 步骤4、源的分配。
[0039] 根据粒子所在的位置求得其对周围网格点电荷与电流的贡献，然后将所有粒子对网格点上的电荷与电流贡献累加得到网格点上的电荷密度与电流密度；
[0040] 步骤1、3至4的求解采用结构化、浸入式非结构化或完全非结构化网格。
[0041] 循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。
[0042] 本发明适用于二维及三维结构，适用于二维时，网格划分从四面体网格变为三角形网格。
[0043] 相对于PIC推动粒子运动的FD方法和IFE方法，本发明的有益效果体现在：
[0044] 1、使用非结构化网格，该网格能够更好的拟合模型边界的形状，使得在复杂边界情况下PIC推动粒子运动算法具有更高的计算精度；
[0045] 2、将FEM推动粒子运动的方法结合到典型的PIC方法中，在保持典型PIC方法的计算简单、快速的优良特性的同时，利用FEM得到更高的有限元计算精度；
[0046] 3、由于FEM方法既可以很好地匹配复杂边界，又可以根据模拟需要使用非均匀网格，且不受数值稳定性条件的限制，因此可以在保持计算精度的条件下，优化空间网格和时间步长，从而大幅地提高模拟效率。
附图说明
[0047] 图1为PIC求解的IFE网格示意图；
[0048] 图2为PIC求解的FEM网格示意图；
[0049] 图3为粒子快速定位算法示意图；
[0050] 图4为粒子快速定位算法流程图；
[0051] 图5为七孔双栅离子光学系统的PIC静电模型计算实例示意图；
[0052] 图6为七孔双栅离子光学系统的PIC静电模型计算实例网格划分示意图。
具体实施方式
[0053] 下面通过实施实例对本发明作进一步详细说明。
[0054] 以离子
推进器七孔双栅离子光学系统为例，其示意图如图5所示。采用本发明中算法对这一实例进行PIC静电模拟的具体实施步骤如下：
[0055] 步骤1、电场求解。
[0056] 通过求解电场所满足的静电离散泊松方程，得出所有网格点上的电场；
[0057] 步骤2、离子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
[0058] 步骤3、推动离子运动。
[0059] 采用全局非结构化网格，通过求解离散离子运动方程，更新离子的速度及位置等运动信息，进而更新离子所属网格。详细的求解过程分为如下两个步骤：
[0060] 步骤3.1、更新离子的速度和位置
[0061] 此处不考虑磁场，直角坐标系下的离子运动方程如(2)式所示：
[0062]
[0063] 其中：上标一个点表示该变量对时间的一次导数，两个点表示该变量对时间的二次导数；η为离子的荷质比；γ为相对论因子；Ex、Ey、Ez为离子所在位置处的电场分量，由步骤2求得；
[0064] 利用Runge-Kutta方法求解式(2)，得到一个时间步长后离子新的速度和位置。
[0065] 步骤3.2、更新离子所属网格
[0066] 每次更新离子的速度和位置过后更新离子所属网格，为后续的步骤4做准备。
[0067] 这里采用基于非结构化网格的离子快速定位算法，用于准确高效地计算离子所属网格，具体实现方式如下：
[0068] 1)在第一步时间步长时，根据离子出现位置或发射面设定，初始化离子所处网格编号。
[0069] 2)从第二步时间步长开始，首先判断离子是否仍然在当前网格内。如果在，则结束搜索；
[0070] 如果不在，则进入3)。
[0071] 3)连线前一时刻离子所在网格PreId的重心P1和当前时刻离子所在位置P2，得到线段P1P2,并求出线段P1P2与PreId的相交面FaceId，同时得到与PreId共用该相交面FaceId的新网格，将新网格设为PreId，回到2)。
[0072] 步骤4、电荷分配。
[0073] 根据离子所在的位置求得其对周围网格点电荷的贡献，然后将所有离子对网格点上的电荷贡献累加得到网格点上的电荷密度；
[0074] 步骤1、3至4的求解可采用结构化、浸入式非结构化以及完全非结构化等网格。
[0075] 循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。采用本发明中算法对这一实例进行PIC静电模拟的结果如图5所示。
