一种基于测地线的用于粒子CFD应用的复杂几何重建方法

《Computer Physics Communications》：A geodesic based complex geometry reconstruction method for particle CFD applications

　　

基于粒子的等值面重构技术对于复杂几何体在科学计算、工程应用和计算机图形学等多个领域具有至关重要的作用。本文建立了一套系统的理论框架，用于在任意复杂几何体上准确重构粒子。该方法有效地解决了在重构具有复杂结构的表面时常见的重要挑战，如粒子聚集或稀疏分布问题，尤其是在包含锐边或其他局部特征的区域。所提出的方法将粒子重分配从有约束的欧几里得优化重构为离散流形上无约束的黎曼优化。基于测地线理论，该方法将基于欧几里得距离的传统邻近粒子搜索和能量最小化策略，从传统的三维欧几里得空间扩展到基于测地线距离的表示，使梯度下降可以直接在复杂表面上进行。这使得粒子在整个运动过程中能够在明确定义的几何表面上得到全面的表征。从计算几何的角度来看，该方法确保粒子的移动严格限制在表面上，从而在整个重构过程中保持原有拓扑结构和所有几何细节。

本文选取了一系列多边形表面数字模型——包括CAD设计模型和3D扫描模型——作为测试案例，以验证所提出方法的有效性。结果表明，该方法能够实现复杂几何体的高精度重构。数值模拟进一步证明，该方法为应用于具有复杂几何结构的真实三维工程问题的基于粒子的CFD方法提供了高精度且适定的初始条件，例如均匀且物理一致的粒子分布。这一进展为在几何条件具有挑战性的情况下开发高效且稳定的基于粒子的求解器奠定了坚实的基础。

引言：基于粒子的方法[1]、[2]、[3]、[4]、[5]已成为模拟复杂工程问题（如大变形、断裂、动态界面和结构失效）的重要数值工具。与传统基于网格的方法不同，这些方法使用相互作用的粒子离散化物理域，无需计算网格，从而具有更好的适应性。然而，尽管具有这些优势，这些方法在收敛性和一致性方面仍面临重大挑战[6]、[7]，因为这些性质高度依赖于粒子分布。具体而言，不规则的粒子分布无法满足够约束条件[6]，从而破坏了零阶一致性[8]、[9]、[10]、[11]、[12]。此外，这种不规则性违反了单元划分性质，导致收敛率降低[12]。因此，生成均匀的粒子分布已成为粒子方法中的关键步骤。该任务涉及两项关键要求：(1) 域均匀性，即粒子必须在计算域的欧几里得空间内均匀分布；(2) 边界均匀性，即粒子必须沿着嵌入在欧几里得空间中的边界流形上均匀排列，以准确表示几何边界。对于具有复杂几何结构的计算域，同时满足这两项要求尤为具有挑战性。

在当前的基于粒子的方法中，现有的粒子初始化技术可以基于其生成策略分为两大类：直接生成方法和迭代生成方法。直接生成方法直接在网格上确定粒子位置。最早的直接方法在规则立方晶格结构上生成粒子。例如，Mayrhofer等人[13]开发了一个名为SPHysicsGen的预处理工具，专门用于开源SPH代码SPHysics[14]、[15]、[16]。该工具利用规则晶格结构构建初始粒子配置；然而，它仅限于创建基本的二维几何结构，并且仅支持有限范围的几何形状和物体类型。在此基础上，Domínguez等人[17]开发了一个名为GenCase的预处理工具，集成到支持GPU加速的开源代码DualSPHysics[18]、[19]、[20]、[21]中。GenCase通过在规则笛卡尔网格上遍历整个计算域来填充粒子，仅将粒子放置在目标几何所需的节点位置上，并根据几何布尔运算删除域外的粒子。该方法能够以非常简单、近乎实时的方式快速初始化数百万个粒子。此外，生成的粒子分布自然地满足域均匀性要求，因此在域内表现出各向同性的分布特征。然而，这种基于晶格的方法在表面附近生成的边界粒子呈阶梯状分布。这导致边界曲率失真和边界法向量不连续，无法满足边界均匀性要求，从而对数值结果产生不利影响。

除了基于晶格的生成方法外，粒子也可以直接在体网格上生成。这些方法利用有限元方法(FEM)或有限体积方法(FVM)中成熟的预处理工具[22]、[23]，将三维体单元的中心、顶点或高斯积分点映射到粒子位置。该方法广泛应用于鸟撞模拟中[24]、[25]、[26]、[27]、[28]。典型的工作流程包括：从鸟的CAD模型出发，将其离散化为四面体或六面体网格，直接在单元几何中心生成粒子，并为其分配相应的单元体积。最近，Wang等人[29]提出了一种基于移动粒子半隐式(MPS)方法的数值框架，用于模拟离心泵中的流动。在他们的方法中，墙边界粒子通过将原始3D CAD模型导入网格化软件，将其离散化为表面壳网格，并利用壳网格节点作为边界粒子位置来生成。尽管目前主流商业软件支持这种方便的转换，并有效地避免了晶格方法相关的阶梯效应，但生成的粒子分布存在严重的非均匀性。因此，它们无法满足域内和边界上的均匀性要求。

除了直接生成方法之外，迭代生成方法近年来也得到了快速发展。第一类涉及粒子打包算法，Colagrossi等人[30]引入了由分布非均匀性驱动的斥力，以驱动粒子迭代移动来填充域内的空隙。同时，施加阻尼以确保平滑收敛，使内部粒子演化至静态平衡位置。该过程在简单的二维区域内生成规则、非重叠的初始分布，满足内部均匀性。随后，Jiang等人[31]设计了一个具有蓝噪声特征的三维粒子采样工具，利用粒子间斥力迭代地将随机分布转变为均匀状态。上述方法主要解决计算域内粒子的均匀性问题，而不考虑边界处理或边界均匀性。因此，边界条件依赖于虚拟粒子；然而，在复杂的三维几何结构中构建虚拟粒子层是困难的，这严重限制了它们在复杂场景中的应用。在此基础上，Negi等人[32]基于开源PySPH框架[33]、[34]进一步开发了粒子打包算法。通过围绕计算域外围布置虚拟粒子而不是直接将其放在边界上，他们降低了构建虚拟粒子层的难度。在演化过程中，近壁粒子被投影到几何表面上作为边界粒子，而其余粒子则被分类为内部或外部。尽管生成的边界粒子紧贴表面，但它们并非均匀排列，因此无法满足边界均匀性要求。最近，Liu等人[35]提出了一种基于边界积分的粒子初始化算法。该方法首先在边界附近排列并固定一层粒子以替代虚拟粒子来施加边界条件，然后调整域内的粒子分布。虽然该方法生成的初始粒子配置满足内部均匀性，但边界均匀性问题仍未解决。

除了粒子打包算法之外，加权沃罗诺伊镶嵌(WVT)方法[36]、[37]通过引入斥力在域内重新分配粒子。它根据粒子间距和局部空间分辨率动态调整粒子位置，从而使粒子更紧密地贴合复杂几何边界。Diehl和Siemann[38]、[39]等先后将WVT方法应用于天体物理中的SPH模拟和鸟撞问题。随后，Vela等人[40]将WVT扩展到支持任意分辨率的低噪声粒子生成。在此基础上，Fu等人[41]将质心沃罗诺伊镶嵌(CVT)与沃罗诺伊粒子动力学相结合，开发了一种在简单的二维几何域内生成均匀粒子分布的方法。然而，与粒子打包算法类似，上述方法依赖于虚拟粒子来施加边界条件，因此不适合用于复杂的三维几何结构。

除了上述方法之外，以粒子松弛方法为代表的迭代预处理技术仍在继续发展。该方法已被集成到开源代码SPHinXsys[42]、[43]中。该方法最初由Nikolaus等人[44]提出，将粒子视为受短程斥力作用的质点。在边界约束下，这些粒子进行动态演化，直到达到静态平衡状态。松弛后，内部粒子之间的间距变得均匀，靠近边界的粒子紧密贴合复杂几何表面。然而，与粒子打包算法和WVT方法类似，粒子松弛方法仅解决域内粒子的均匀分布问题，而未解决边界均匀性。因此，必须手动布置虚拟粒子来施加边界条件，以确保边界附近的内部粒子获得完整的核函数支撑。

为解决这些挑战，提出了一系列改进方法。Ji等人[45]在松弛框架中引入了特征边界修正项，初步缓解了近壁区域核支撑不足的影响；然而，这种方法仍然需要精确地在几何表面上生成一层边界粒子。随后，Zhu等人[46]将水平集场纳入松弛技术中，以实时检测粒子到几何边界的距离，从而消除对虚拟粒子的依赖。他们还对试图穿越边界的内部粒子施加了约束力。2023年，Zhang等人[47]直接在CAD表面上生成边界粒子以替代虚拟粒子，利用粒子松弛技术优化内部粒子分布，并采用水平集为边界粒子提供法向量。同年，针对由非解析几何形状引起的粒子松弛失效问题，Yu等人[48]使用水平集场来识别和光滑锐利几何特征。2024年，Kang等人[49]通过动态地在边界附近添加粒子来减少总体积误差，并结合松弛方法，实现了体积守恒的均匀内部粒子分布，同时确保粒子不穿透边界。尽管该方法生成的边界粒子紧贴几何表面，但它们并非均匀排列。同样在2024年，Zheng等人[50]提出了一种无参数粒子松弛技术。该方法摒弃了先前方法中对水平集背景网格的依赖，转而采用零阶一致性准则和梯度下降方法来修正粒子位置。这在减少积分误差的同时实现了更均匀的内部粒子分布。

此外，在大型三维问题中，粒子松弛方法通常需要数千次迭代才能达到稳定状态。针对高计算成本的问题，Yang等人[51]（2025年）在共享内存架构下实现了松弛过程的完全并行化，并集成了布尔运算，以高效处理包含多组件装配的复杂工业模型。同年，Zhao和Zhang等人[52]、[53]将粒子松弛技术扩展到多体系统，并明确指出均匀初始粒子对于提高Eulerian SPH的精度也至关重要。最近，Fan等人[54]将现有的粒子松弛方法重新表述为一个优化问题，从理论上揭示了均匀粒子分布与零阶一致性优化之间的等价性。总体而言，尽管现有的松弛方法可以实现域内粒子均匀、贴合物体的排列，但松弛过程仅发生在几何体内部。因此，这些方法无法保证边界粒子分布的质量，无法满足边界均匀性要求。

最近，为满足粒子初始化中自适应分辨率的需求，Liu等人[55]（2025年）提出了一种无网格自适应粒子生成方法。该方法首先通过粒子分裂操作生成初始粒子，然后采用粒子移动算法使粒子自动填充并贴合几何形状。它在复杂域内实现了具有自适应分辨率的均匀粒子分布，同时保持了几何体的体积。同年，Villodi等人[56]也提出了一种无需背景网格的粒子初始化方法。通过将恢复力驱动的动力学与移动技术相结合，他们的方法能够迅速将粒子松弛为均匀各向同性的分布，同时结合粒子分裂和合并机制以支持局部自适应性。

然而，尽管这些方法能够实现域内粒子的自适应分辨率，却未能解决边界粒子的自适应分辨率问题。综上所述，从基于规则格点和体网格的直接生成方法，到粒子堆积算法、加权Voronoi划分、粒子松弛技术、水平集以及自适应分辨率方法等迭代生成技术，近期预处理方面的研究显著提高了复杂几何体粒子初始化的质量。然而，大多数成果集中在实现域内均匀性上。相比之下，关于边界均匀性的研究仍相对匮乏。这种差异源于这两种均匀性在实现空间上的根本不同。实现内部均匀性仅需在欧氏空间中均匀分布粒子，因此可直接使用欧氏距离计算排斥力。相反，边界均匀性需要在嵌入于欧氏空间中的曲面流形上操作。因此，应采用能更好反映真实曲面距离的测地距离来计算排斥力。此外，施加于粒子的运动学约束在两种情形下也存在本质差异。内部均匀性要求施加边界条件以确保粒子留在域内，涉及简单的直线空间位移。然而，对于边界均匀性，粒子被约束在表面边界上，必须严格沿着流形移动。因此，它们的运动必须遵循更为严格的"最直"测地路径，以确保不偏离曲面。总体而言，满足边界均匀性的要求提出了显著更大的挑战。然而，现有的边界粒子生成方法通常将近壁粒子直接投射到表面上作为边界粒子。虽然这确保了粒子紧贴曲面，但无法实现均匀分布。或者，一些方法，如Zhang等人[57]的工作，将边界均匀性表述为欧氏空间中的约束优化问题。这些方法在欧氏空间内建立粒子与曲面之间的约束关系，通过欧氏距离和直线位移迭代实现边界粒子的均匀分布。然而，在实践中，由于缺乏严格的几何保证，仅由欧氏空间中的约束控制的粒子在迭代后不可避免地会偏离或穿透曲面。为应对在复杂几何曲面上实现均匀边界粒子分布的挑战，本文在分析现有方法，特别是Zhang等人[57]工作的基础上，提出了一种创新方法。我们将粒子重分布从欧氏空间中的约束优化扩展为离散流形上的无约束黎曼优化。该方法将边界均匀性问题重构为直接在曲面流形上执行的梯度下降，确保粒子均匀化操作在流形本身而非欧氏空间中执行。这一转变是必要的，因为现有方法所用的欧氏距离无法准确反映表面粒子之间的约束距离。此外，欧氏空间中的简单直线移动不可避免地导致粒子偏离或穿透曲面，导致迭代过程中累积误差。因此，本文将邻域搜索和能量梯度计算均用测地距离进行重构。通过使用精确单源测地距离(SSGD)和最直测地线理论计算粒子的黎曼梯度和测地轨迹，我们确保粒子运动严格贴合曲面。随后，应用迭代下降算法最小化相互作用能，直至达到平衡状态。该方法在任意复杂曲面上实现了均匀且几何保真的边界粒子分布。本文的主要贡献总结如下：•建立了用于任意复杂几何体粒子重分布的黎曼优化框架，将问题从约束欧氏优化重构为无约束流形优化。这有效解决了重建过程中因几何曲面复杂性导致的粒子聚集和稀疏分布等关键问题。•通过将邻域搜索和梯度下降从三维欧氏空间扩展到复杂曲面上基于测地的表示，提出了一个完整的离散流形上能量最小化的计算框架。这确保了优化过程中表面约束的精确满足。•所提出的算法在多边形曲面上生成高质量表面粒子分布，同时完全保留原始几何边界的拓扑和形状细节，为复杂几何体中的粒子CFD模拟方法提供了鲁棒且适定初始粒子构型。本文其余部分组织如下。第2.1节和2.2节总结了粒子附着于曲面的理论及相关几何约束。第2.3节提出了一种在任意复杂曲面上生成边界粒子的一般化方案，为后续重分布过程提供鲁棒初始化。第3节建立了等值面粒子重分布的物理模型。基于精确测地理论，第4.1节和4.2节分别详细分析了邻域粒子搜索和曲面约束粒子运动。第4.3节提出了在保持几何保真度的同时实现均匀表面粒子分布的算法。第5节和第6节分别通过一系列几何测试案例和典型数值算例进行了验证。第7节概述了结论和展望。各章节摘要粒子重建曲面基本理论回顾了基于欧氏距离和优化理论的典型粒子重分布方法[57], [58]。该方法将粒子和曲面的运动表述为约束动力系统，解决了建模受可变形曲面约束的n个运动粒子的一般性问题。曲面上移动粒子的轨迹表示为g(x, t)，而曲面的时变形状参数表示为q(t)。则曲面可表示为等值面粒子重分布模型本节将等值面粒子重分布问题[57], [58]重构为离散流形上的无约束能量最小化问题。与在R3中执行约束优化的欧氏方法不同，我们的方法直接在多边形曲面Γ上执行梯度下降，曲面约束在内禀上得到满足。粒子动力学由相互作用核的能量梯度控制，所有运动严格限制在曲面流形上。多边形曲面上的邻域粒子搜索不同于基于欧氏的粒子搜索方法，在曲面Γ上识别相互作用粒子需要确定它们之间的测地路径和距离。由于任何平面多边形都可以三角剖分，采用三角网格作为背景几何来说明粒子搜索过程。对于两个分别位于宿主三角形f?和f?上的表面粒子i和j，如图2所示，存在三种不同的情形。因此，通过将局部邻域映射验证案例与分析本节使用一系列三维几何模型对所提出的算法进行验证。此外，我们定量分析了基于测地重分布前后的粒子分布，并与最新的基于欧氏距离的方法[57]的结果进行了比较。SPH数值算例本节通过分析静水压和水入水场景，考察初始粒子分布质量对数值结果的影响，其中不同的分布导致不同的构型和初始条件。所有模拟均使用开源代码DualSPHysics进行。结论本研究首次系统地将粒子重分布技术从平面三维欧氏空间通过黎曼优化适应到复杂曲面。关键创新在于将粒子重分布从约束欧氏优化重构为直接在曲面上执行的无约束梯度下降，通过测地路径确保精确的曲面约束，如算法4中所述。基于离散微分几何，该方法保证拓扑CRediT作者贡献声明Jiatao Zhang: 撰写—审阅与编辑、可视化、验证、软件、方法学、形式分析、数据整理。Xiaohu Guo: 撰写—审阅与编辑、撰写—初稿、可视化、验证、监督、方法学、调查、形式分析、数据整理、概念化。Xiufang Feng: 撰写—审阅与编辑、撰写—初稿、验证、监督、资源、项目管理、方法学、调查、资金获取、形式声明利益冲突声明作者声明以下财务利益/个人关系，可能被视为潜在的利益冲突：通讯作者目前为编辑委员会委员。鉴于其在HPC、网格与新型计算领域的专业编辑身份，他未参与本文的同行评审，也未获取关于其同行评审的信息。本文编辑过程的全部责任已委托给另一位期刊编辑。如果Jiatao Zhang|Xiaohu Guo|Xiufang Feng|Li Zhu宁夏大学数学与统计学院, 750021, 中国