CN111695281A - 一种四面体网格划分有限元粒子模拟的粒子快速定位方法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

粒子模拟(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_p,y_p,z_p)，当前搜索四面体网格的四个顶点按右手定则依次为A(x_A,y_A,z_A)，B(x_B,y_B,z_B)，C(x_C,y_C,z_C)和D(x_D,y_D,z_D)。

首先，依次取顶点A、B、C、D相对的三角形面，各顶点相对的三角形面的重心点依次为M_A,M_B,M_C,M_D，分别连接粒子点P与各顶点相对三角形面的重心点，得到向量

再分别求得三角形面指向对应顶点的法向量，依次为

规定θ_A为向量

与

之间的夹角，同理θ_B,θ_C,θ_D分别为

与

与

和

与

的夹角，如图2所示。

然后，依次计算θ_A,θ_B,θ_C,θ_D的角度，此时若四个角度都小于或等于90°则粒子位于当前网格内，返回粒子所属网格单元索引值；否则选择大于90°且角度最大的夹角θ，取其下标对应顶点相对的三角形面Ω，作为下一步搜索的方向，进行步骤4。

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

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

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

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

附图说明

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

图2为四面体有限元网格公共面选取策略示意图；

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

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

具体实施方式

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

以基于有限元四面体网格划分的三维计算区域内的PIC为例，离散化的计算区域的有限元四面体网格如图3所示。采用本发明中的粒子快速定位方法来确定这一实例中运动更新后粒子所属的网格索引值。操作流程如图4所示，具体实施步骤如下：

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

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

步骤3、确定进行下一步搜索的方向。本实例为三维计算区域，采用特定策略选择四面体网格单元的一个三角形面，详细方法如下：已知运动更新后的粒子位置为点P₁(x_p,y_p,z_p)，当前搜索的四面体网格的四个顶点按右手定则依次为A(x_A,y_A,z_A)，B(x_B,y_B,z_B)，C(x_C,y_C,z_C)和D(x_D,y_D,z_D)；首先，依次取顶点A、B、C、D相对的三角形面，各顶点相对的三角形面的重心点依次对应为M_A,M_B,M_C,M_D，分别连接粒子点与各顶点相对三角形面的重心点，得到向量

再分别求得三角形面指向对应顶点的法向量，依次为

规定θ_A为向量

与

之间的夹角，同理θ_B,θ_C,θ_D分别为

与

与

和

与

的夹角。

然后，依次计算θ_A,θ_B,θ_C,θ_D的角度，此时若四个角度都小于或等于90°则粒子位于当前网格内，返回粒子所属网格单元索引值；否则选择超过90°且角度最大的夹角θ，取其下标对应顶点相对的三角形面，作为下一步搜索的方向，在本实例中，选取的面为节点A、C、D构成的三角形面Ω₁，然后进行步骤4。

步骤4、进行搜索，即由步骤3得到当前四面体网格的一个三角形面Ω₁，根据步骤1存储的数据结构，通过网格拓扑关系得到与当前体四面体网格单元共用三角形面Ω₁的网格单元，本实例中与a₁号网格单元共用面Ω₁的网格单元记为a₂网格，将粒子运动更新前所属单元索引值更新为a₂，在新得到的a₂号网格单元内继续进行步骤2。

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

在本实例中，通过不断递归迭代，依次经过a₁,a₂,a₃,…,a_n号网格单元，在n次搜索操作后最终得到运动更新后的粒子所属的网格单元索引值，在本实例中为a_n号网格单元。

Claims (1)

1.一种四面体网格划分有限元粒子模拟的粒子快速定位方法，具体步骤如下：

步骤1、对通过非结构化四面体有限元网格划分产生的离散网格建立数据结构，用于存储网格的基本信息和拓扑关系；

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

步骤3、确定进行下一步搜索的方向；

对于四面体网格划分的三维计算区域，已知运动更新后的粒子位置为点P(x_p,y_p,z_p)，当前搜索四面体网格的四个顶点按右手定则依次为A(x_A,y_A,z_A)，B(x_B,y_B,z_B)，C(x_C,y_C,z_C)和D(x_D,y_D,z_D)；

首先，依次取顶点A、B、C、D相对的三角形面，各顶点相对的三角形面的重心点依次对应为M_A,M_B,M_C,M_D，分别连接粒子点P与各顶点相对三角形面的重心点，得到向量

再分别求得三角形面指向对应顶点的法向量，依次为

规定θ_A为向量

与

之间的夹角，同理θ_B,θ_C,θ_D分别为

与

与

和

与

的夹角；

然后，依次计算θ_A,θ_B,θ_C,θ_D的角度，此时若四个角度都小于或等于90°则粒子位于当前网格内，返回粒子所属网格单元索引值；否则选择大于90°且角度最大的夹角θ，取其下标对应顶点相对的三角形面Ω，作为下一步搜索的方向，进行步骤4；

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

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

Images