一种应用于PIC静电模型的电位有限元求解算法的制作方法

本发明属于粒子模拟(particle-in-cell，简写为pic)的数值模拟领域，具体涉及一种应用于pic静电模型的电位有限元求解算法。

背景技术：

pic方法是一种被广泛应用于带电粒子与电磁场相互作用物理问题中的数值模拟方法，它通过跟踪大量带电粒子在外加及自洽电磁场中的运动并统计平均而得到宏观特性及运动规律。经过几十年的发展，pic模拟方法已经成为研究带电粒子与电磁场相互作用物理问题的一种强有力的数值手段，广泛应用于带电粒子与电磁场相互作用所涉及的许多领域，如磁约束聚变等离子体、惯性约束聚变等离子体、核爆、空间等离子体、人造等离子体(包括电子枪、离子源等)、电推进、自由电子激光以及电真空器件等。

pic方法按照求解电磁场方程形式的不同而分为静电模型、电磁模型和静磁模型，其中静电模型主要适用于静电分离为主要物理矛盾的带电粒子与时变静电场相互作用问题，如电推进系统中的离子引出过程、朗缪尔振荡、电子枪和收集极中电子的运动轨迹演化过程等。

静电模型求解的核心步骤如下：

1、电位求解，即通过求解静电场满足的离散泊松方程，得出所有网格点上的电位；

2、粒子受力求解，即通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；

3、推动粒子运动，即通过求解离散粒子运动方程，更新粒子的动量及位置等运动信息；

4、电荷分配，即根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；

不断循环如上过程，直到计算结果收敛或人为设置的时间为止。

其中步骤1电位求解是pic静电模型必不可少的核心步骤之一，该步骤的精确且高效的求解对于pic静电模型的整体求解精度和效率的控制十分重要。截止目前为止，pic静电模型中电位求解主要有两种方法，分别是有限差分(fd)方法和嵌入式有限元(ife)方法。

fd方法：在pic静电模型电位求解的应用中，fd方法是通过采用结构化网格离散求解区域和静电泊松方程，然后采用有限差分的方法求解泊松方程来得到网格点上的电位。

fd方法完全基于结构化网格，因此它的优点是算法形式简单、易于理解，但是在pic静电模型的应用中存在如下缺点：

1、fd方法采用的是由正交线划分而成的结构化网格，对于复杂曲形边界的拟合较差，从而使得复杂曲形边界附近电位求解精度较低，进而降低全部区域电位求解精度；

2、由于fd方法对网格尺寸均匀性的要求比较高，因此受限于模拟系统中细小物理结构的限制，必须划分足够小的网格才能满足计算精度要求，这就导致计算量十分庞大。

3、fd方法严格受限于数值稳定性条件的限制，即在对pic静电模型的数值模拟中，如果空间网格尺寸很小，则时间步长也会随之取的很小，这会进一步增加对电位时间循环求解的fd数值模拟负担。

针对fd方法在pic静电模型电位求解应用中出现的不足，kafafy和wang在2003年，提出了可应用于pic静电模型电位求解的ife方法。

ife方法：在pic静电模型电位求解的应用中，ife方法是通过采用侵入式非结构化网格离散求解区域和静电泊松方程，然后采用有限元的方法求解泊松方程来得到网格点上的电位。侵入式非结构化网格划分情况如图1所示，可以看出此网格划分是将一个结构化的六面体又进一步划分成五个四面体，整体还是结构化网格。

ife方法相比于fd方法能在一定程度上解决pic静电模型电位求解中的复杂边界问题，但是，依然存在如下缺陷：

1、ife方法虽然使用侵入式非结构化网格，但是整体还是基于结构化网格的划分，由于结构化网格的边界规整，对于不规则模型的边界匹配度不高，产生的误差较大；

2、从图1中可以看出，ife方法使用的四面体网格的角度为直角，而根据有限元基本理论，四面体越接近正四面体，网格质量越好，计算精度越高，而嵌入式有限元网格与正四面体相差甚远，因此采用嵌入式有限元网格在同等网格规模的情况下，计算精度远远低于正四面体网格。另外，网格质量越好，最终形成的系数矩阵条件数越好，求解速度越快。

技术实现要素：

针对上述存在的问题和不足，为解决fd和ife方法在电位求解中边界匹配度不高、求解精度不高的问题，本发明提供了一种应用于pic静电模型的电位有限元(fem)求解方法。具体技术方案如下：

步骤1、电位求解。

采用全局非结构化网格，其三维网格划分实例如图2所示。

pic静电模型中电位的求解是二阶偏微分方程定义的边值问题。

首先将该边值问题等价变换为如下变分问题：

其中：

其中α，β，γ是于与区域物理性质有关的已知参量，f是源或激励函数。然后将计算区域体v离散为m个四面体，在每个网格单元内，未知函数φ表示为：

φ^e(x,y,z)＝a^e+b^ex+c^ey+d^ez(3)

其中上标e代表某一网格单元。将φ在四个顶点处的值带入(3)式，可由克莱姆法则解得系数a^e,b^e,c^e,d^e，并将其带回(3)式，整理可得：

其中下标j表示e网格单元中的j号顶点，网格单元的插值函数为：

然后将(4)式带入(2)式，得：

分别取f^e对每个顶点的偏导，并写成矩阵形式：

将(7)式给出的单元方程对所有单元求和，然后强加驻点条件，可得到方程组，并结合边界条件，解之即得每个网格顶点上的电位。

步骤2、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力。

步骤3、推动粒子运动

通过求解离散运动方程，更新粒子的动量及位置等运动信息；

步骤4、电荷分配

根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；

步骤2至4的求解可采用结构化、浸入式非结构化或完全非结构化网格。

循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。

本发明适用于二维及三维结构，适用于二维时，网格划分从四面体网格变为三角形网格。

相对于pic静电模型电位求解的fd方法和ife方法，本发明的有益效果体现在：

1、使用完全非结构化网格，该网格能够更好的拟合模型边界的形状，使得在复杂边界情况下pic静电模型的电位求解具有更高的计算精度、更快的求解速度；

2、将用于求解无源电磁场分布、热分析、力学分析等无粒子源问题的fem求解电位的方法结合到典型的pic方法中，在保持典型pic方法的计算简单、快速的优良特性的同时，利用fem得到更高的有限元计算精度；

3、由于fem方法既可以很好地匹配复杂边界，又可以根据模拟需要使用非均匀网格，且不受数值稳定性条件的限制，因此可以在保持计算精度的条件下，优化空间网格和时间步长，从而大幅地提高模拟效率。

附图说明

图1为pic静电模型电位求解的ife网格示意图；

图2为实施例中pic静电模型电位求解的fem网格示意图；

图3为七孔双栅离子光学系统的pic静电模型计算实例示意图；

图4为七孔双栅离子光学系统的pic静电模型计算实例网格划分示意图。

具体实施方式

下面通过实施例和附图对本发明作进一步详细说明。

以离子推进器七孔双栅离子光学系统为例，其示意图如图3所示。采用本发明中算法对这一实例进行pic静电模拟的具体实施步骤如下：

步骤1、电位求解：采用基于非结构化网格的fem算法。

求解电位满足的possion方程：

式中，ρ为电荷密度，ε0为真空介电常数。

变换为等价变分问题为：

因为求解的是possion方程，所以αx＝1，αy＝1，αz＝1，β＝0，则可得：

然后将三维计算区域体v离散为m个四面体，网格划分情况如图4所示。在每个网格单元内，未知函数φ为：

φ^e(x,y,z)＝a^e+b^ex+c^ey+d^ez(11)

将φ在四个顶点处的值带入(11)式，并由克莱姆法则可以解得

其中v^e为四面体体积，进而得：

同样地也可以解出b^e,c^e,d^e。并将解出的a^e,b^e,c^e,d^e带入(11)式，合并同类项，可得

其中下标j表示e网格单元中的j号顶点，网格单元的插值函数为：

然后将(14)式带入(10)式，得：

取f^e对每个顶点的偏导：

分别取f^e对每个顶点的偏导，并对f应用驻点条件，则针对网格单元e可构成线性方程组：

其中，

最后，将m个如式(18)这样的线性方程组组合成一个线性方程组，求解方程组即得每个网格顶点上的电位。

步骤2、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；

步骤3、推动离子运动

通过求解离散运动方程，更新离子的动量及位置等运动信息；

步骤4、电荷分配

根据离子所在的位置求得其对周围网格点电荷的贡献，然后将所有离子对网格点上的电荷贡献累加得到网格点上的电荷密度；

步骤2至4的求解可采用结构化、浸入式非结构化以及完全非结构化等网格。

循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。采用本发明中算法对这一实例进行pic静电模拟的结果如图3所示。