一种三角面网格划分有限元粒子模拟的粒子快速定位方法

摘要

本发明属于有限元粒子模拟领域，具体涉及一种三角面网格划分有限元粒子模拟的粒子快速定位方法。本发明通过评估粒子运动后所在点与三角形网格单元各个顶点形成的向量的夹角来选取网格中的一条边作为网格进行下一步搜索的方向，使用网格现有的数据，避免了保存其他额外信息带来的存储占用开销，同时也避免了计算并判断线段与三角形哪条边相交之类带来的繁重计算负担，方法简单快速，有效降低了实现成本。使得非结构化三角面有限元PIC程序具有更高的运行效率。这样，在利用非结构化三角面有限元网格更好地拟合非规则非正交物理模型边界、提高计算精度的同时，仍然不失PIC方法计算简单、快速的优良特性。

说明书

技术领域

本发明属于有限元粒子模拟领域，具体涉及一种三角面网格划分有限元粒子模拟的粒子快速定位方法。

背景技术

粒子模拟(Particle-in-Cell，简记为PIC)方法是用计算机跟踪大量微观粒子的运动，再对组成物理系统的微观粒子进行统计平均，由此得到物理系统的宏观特性和运动规律的数值模拟方法。对于物理规律还不十分清楚的问题，PIC方法可以帮助建立明确的物理图像，促进新的理论模型的萌芽。现在，PIC方法已经成为研究有源复杂物理问题非线性效应的强有力手段，在可控热核聚变、空间物理、自由电子激光和一般等离子体问题的研究中都有广泛的应用。PIC方法的基本思路是给定大量带电粒子的初始位置和速度，对它们统计平均求出空间的电荷和电流密度分布，再通过麦克斯韦方程组求出空间中的电场和磁场，进而求得每个粒子所受到的洛伦兹力，下一时刻每个粒子的位置和速度则可以通过运动方程更新，如图1所示。如此循环进行，通过跟踪计算大量带电粒子的运动，再根据所感兴趣的问题对这些大量带电粒子的某些物理量作统计平均，即可得到宏观的物理特性和运动规律。

在上述场与粒子相互作用过程中，有时结构和物理特性是满足轴对称性的，或者某一个方向的物理特性可以通过理论解析的方法得到，则对该类问题的数值模拟可以采用二维PIC方法来实现。

传统的PIC方法采用结构化的正交网格对计算区域和计算方程进行离散，在规则正交求解区域其求解精度较高。但是，当求解区域整体或部分为非规则非正交的区域时，结构化网格大多采用阶梯近似来离散非规则非正交的求解区域，这种阶梯近似网格对模型的拟合效果较差，从而使得非规则非正交区域的计算求解精度较低。如果要提高非规则非正交区域的求解精度，则必须在整个计算区域划分足够小的网格，而这样会大大增加网格数量，进而导致繁重的计算负担。

针对结构化正交网格PIC方法的上述缺陷，鉴于有限元网格对非规则非正交区域强大的拟合能力，可以采用基于有限元网格的PIC方法，即使用非结构化的有限元网格(二维情况下为三角形面网格)对计算区域进行离散，并在非结构化的有限元网格上完成PIC方法的计算流程。

有限元PIC方法需要基于完全非结构化的网格进行场量求解、粒子处场插值、粒子运动更新以及电荷、电流源分配的一系列操作，其中在进行粒子处场插值和电荷、电流源分配操作时，需要确定运动更新后粒子所属的网格单元。在结构化网格下，由于网格坐标和网格索引之间存在简单有序的映射关系，因此可以通过粒子坐标除以网格步长直接求得粒子所属网格索引。而在非结构化网格下，网格坐标和网格索引之间不存在上述简单有序的映射关系，因此无法实现如结构化网格下的粒子所属网格的直接快速求解，必须通过遍历对计算区域内的网格进行搜索。在具有大规模网格和粒子数目的情况下，粒子所属网格定位操作十分频繁且相当耗时，这极大地限制了PIC方法的运行效率。

确定了粒子上一步所在的网格单元，根据粒子运动的局部性，按照某种策略搜索粒子运动后所在的位置，将会比大海捞针的方式遍历确定所属网格要快得多。当前主要有两类搜索方法，即深度优先搜索方法和广度优先搜索方法。

按广度优先的方式进行搜索，按层次搜索围绕在粒子上一步所在三角形网格PreID周围的网格单元，称为“扫描法”。如果在第一层范围内没能找到则接着按上述方法在PreID网格周围的第二层网格单元继续搜索。由于限制了粒子运动的步长，因此保证一定能在有限次数内定位粒子或者判断粒子出界。按深度优先的方式进行搜索，按照明确的方向从起始位置开始，逐步逼近粒子所在的网格，最终实现定位，称为“追踪法”。具体确定搜索方向的策略，可根据粒子运动前后位置形成的射线与三角形网格PreID相交的某一边，找到共用该边的网格并判断粒子是否在该网格单元内，若不在则在该网格内确定一个点，构成新的射线，重复上述步骤，直到最终找到粒子所在的网格。也可以根据计算三角形有限元网格各个顶点的插值基函数，比较数值大小来确定搜索的方向。

此外，还可以通过坐标变换将非结构化网格转化为正交的结构化网格，进而确定粒子所属网格单元。

上述方法都能够在一定程度上有效解决粒子定位的问题，但是，其中“扫描法”在搜索第二层的时候会出现重复判断的情况，需要额外设置一个状态数组来记录网格被搜索的情况。而“追踪法”中，基于判断射线和三角形面的边相交的方法将带来一定的计算负担，基于有限元插值基函数计算的方法则需要额外存储所有有限元网格各个顶点的插值系数，将带来可观的存储开销。对于坐标变换则涉及繁杂的数学公式推导，在计算机实现上成本过高。

发明内容

针对上述存在的问题和不足，为解决有限元PIC方法中粒子所属网格定位过程计算效率低或存储占用高的问题，本发明提供了一种三角面网格划分有限元粒子模拟的粒子快速定位方法。

具体技术方案如下：

步骤1、首先对通过非结构化三角面有限元网格划分产生的离散网格建立数据结构，用于存储网格的基本信息和拓扑关系，包括网格节点的坐标、局部索引(即某一网格单元内各节点的索引)和全局索引以及网格单元的全局索引和具有公共边的相邻网格单元信息。

步骤2、利用粒子属性获得运动更新前所属网格索引值，将其设置为当前搜索网格，并判断粒子是否仍处于当前网格内。如果是，则搜索结束，返回当前网格索引值；如果否，则进行步骤3。

步骤3、确定进行下一步搜索的方向。详情如下：

对于三角面网格划分的二维计算区域，已知运动更新后的粒子位置为点P(x

依次判断θ

步骤4、进行搜索，即由步骤3得到当前三角形网格的其中一条边α，根据步骤1存储的数据结构，通过网格拓扑关系得到与当前三角面网格单元共用边α的网格单元，将粒子运动更新前所属单元索引值设置为所得网格的索引值，在得到的新的网格单元内继续进行步骤2。

步骤5、循环步骤2至4，直至搜索到粒子所属网格单元索引值返回；或者，当前网格单元在步骤3得到的边的方向上不存在共边网格时，则说明粒子此时已经运动到计算区域外，返回-1表示粒子出界，并存储当前网格索引值以及对应的边，以便后续进行粒子出界的处理。

上述步骤适用于基于非结构化三角面有限元网格PIC的粒子快速定位求解。

本发明通过评估粒子运动后所在点与三角形网格单元各个顶点形成的向量的夹角来选取网格中的一条边作为网格进行下一步搜索的方向，使用网格现有的数据，避免了保存其他额外信息带来的存储占用开销，同时也避免了计算并判断线段与三角形哪条边相交之类带来的繁重计算负担，方法简单快速，有效降低了实现成本。使得非结构化三角面有限元PIC程序具有更高的运行效率。这样，在利用非结构化三角面有限元网格更好地拟合非规则非正交物理模型边界、提高计算精度的同时，仍然不失PIC方法计算简单、快速的优良特性。

附图说明

图1为粒子模拟方法的流程图；

图2为三角形有限元网格公共边选取策略示意图；

图3为粒子快速定位操作示意图；

图4为粒子快速定位方法的流程图。

具体实施方式

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

以基于有限元三角面网格划分的二维计算区域内的PIC为例，离散化的计算区域的有限元三角形网格如图3所示。采用本发明中的粒子快速定位方法来确定这一实例中运动更新后粒子所属的网格索引值。

操作流程如图4所示，具体实施步骤如下：

步骤1、首先对通过非结构化三角面有限元网格划分产生的离散网格建立数据结构，用于存储网格的基本信息和拓扑关系，包括网格节点的坐标、局部索引和全局索引以及网格单元的全局索引和具有公共边的相邻网格单元信息等。

步骤2、利用粒子属性获得运动更新前所属网格索引值，本实例中记为a

步骤3、确定进行下一步搜索的方向。本实例为二维计算区域，采用特定策略选择三角形网格单元的一个边，详细方法如下：

已知运动更新后的粒子位置为P

依次判断θ

步骤4、进行搜索，即由步骤3得到当前三角形网格的一条边l

步骤5、循环步骤2至4，直至搜索到粒子所属网格单元索引值返回，或者，当前网格单元在步骤3得到的边不存在共边面网格单元时，则说明粒子此时已经运动到计算区域外，返回-1表示粒子出界，并存储当前网格索引值以及对应的边，以便后续进行粒子出界的处理。

在本实例中，通过不断递归迭代，依次经过a