用于不可压缩两相流问题的半拉格朗日无单元伽辽金方法

《Engineering Analysis with Boundary Elements》：A semi-Lagrangian element-free Galerkin method for incompressible two-phase flow problems

　　

• 开发了一种 Galerkin 型无网格方法用于不可压缩两相流问题。• 该方法融合了半拉格朗日方法和无单元 Galerkin 方法。• 该方法通过半拉格朗日时间离散获得了良好的稳定性。• 该方法在两相流问题求解中展现了良好的性能。• 无网格特性简化了半拉格朗日方法的实施。引言流体流动是自然界和工程设施中最普遍存在的物理现象之一。对流机理的深入研究可以为工程技术的升级提供关键的理论指导。实际工程中的流动问题通常表现出强非线性，并包含多组分流体输运、对流和扩散等复杂的物理机制。因此，开发复杂流动问题的数值模拟方法具有重要的实际意义 [1], [2]。两相流的数值模拟是流体力学中的一个热门且具有挑战性的课题，界面捕获/追踪技术是两相流问题研究中的关键组成部分 [3]。界面捕获方法是目前最广泛采用的技术，其中成熟的方法包括体积分数法（VOF）[4]、水平集法（LS）[5]和相场法 [6], [7]。VOF 方法由 Hirt 和 Nichols [8], [9]提出。其核心思想是定义一个体积分数函数，表示每个计算单元内目标流体的比例，并根据流场速度直接对函数进行平流。该方法具有良好的质量守恒特性，但在界面重构方面需要相当高的计算成本。LS 方法由 Osher 等人 [10]开发。它构建了一个具有符号距离函数特性的 LS 函数，该函数在速度场的控制下随时间演化。每一步中函数的零等值线对应于界面的位置。LS 方法避免了显式界面追踪，能够在不进行界面重构的情况下捕获具有复杂变形的移动界面。它还具备公式简单的特点，可以方便地扩展到高维问题。由于这些优点，LS 方法已被广泛应用于多个领域，如图形和图像处理、结构拓扑优化以及计算流体力学等 [11]。LS 函数在初始时刻表现为符号距离函数。然而，经过若干时间步迭代后，它将失去这一特性，导致界面识别不准确。1994 年，Sussman 等人 [5]提出了重初始化技术以维持 LS 函数的符号距离特征，这在一定程度上改善了质量守恒。在过去几十年中，已开发出大量方法来弥补标准 LS 方法的缺点。典型的方法包括粒子水平集法（PLS）[12]、耦合 LS 和 VOF 法（CLSVOF）[13]以及守恒 LS 法（CLS）[14]。本文旨在提出一种新的数值方法来模拟两相流。所采用的两相流模型由不可压缩 Navier–Stokes 方程、LS 方程和 LS 重初始化方程组成。由于所有这些方程都由纯对流或对流–扩散效应控制，数值求解的主要挑战在于对流主导问题。为了获得此类问题的稳定数值结果，通常需要对对流项进行特殊处理。在过去几十年中，已开发了各种有效的数值技术，包括迎风格式 [15]、人工粘性方法 [16]、稳定化 Galerkin 方法 [17]、变分多尺度方法 [18], [19] 和特征 Galerkin 方法 [20]，以及半拉格朗日（SL）方法 [21]。值得注意的是，SL 方法连接了欧拉和拉格朗日框架。它以拉格朗日方式处理对流–扩散方程中的对流项，通过追踪流体粒子实现，并在每个时间步内在固定的欧拉网格上研究粒子运动。一方面，SL 方法保留了欧拉方法固定计算网格的优点，消除了受 Courant–Friedrichs–Lewy（CFL）条件约束的缺点，这使得它们在采用大时间步时难以保持数值精度和稳定性。另一方面，SL 方法继承了拉格朗日方法的优点，包括没有显式对流项和无条件数值稳定性，同时避免了由对流引起的严重网格畸变。因此，SL 方法可以在固定欧拉网格上克服 CFL 限制，即使在大时间步下也能保持良好的数值稳定性和精度 [22]。SL 方法可以追溯到 20 世纪 50 年代。它最初被用于求解对流方程，并避免数值积分过程中的网格变形 [23], [24]。在 20 世纪 60–70 年代，该方法被应用于各种数值预测模型中，SL 格式的特性也得到了广泛研究。20 世纪 80 年代，Robert [25]证明了基于原始方程的 SL 格式适用于气象、海洋和水文预报。与欧拉方法相比，SL 方法允许采用大得多的时间步同时保持数值稳定性，这为其在数值天气预报中的应用奠定了坚实的基础。此后，研究人员将 SL 方法与各种数值方法相结合，开发了谱元 SL 方法 [26]、SL 有限元方法 [27]、SL 有限体积方法 [28] 等技术。这些方法已被广泛用于各种流体流动问题的数值模拟，并能够产生稳定的数值结果。在每个时间步中，SL 方法从固定网格节点反向追踪流体粒子以确定其起始位置。由于这些起始点通常不与网格节点重合，因此需要进行函数插值以获取这些位置的物理量。相应地，反向轨迹追踪和函数插值构成了 SL 方法的两个核心组成部分 [22]。作为 SL 方法的关键技术，轨迹追踪长期以来一直是主要研究焦点。现有的反向追踪方法包括后向欧拉法 [22]、中点法 [29]、Runge–Kutta 格式 [30]、Adams–Moulton 格式 [31] 和 Adams–Bashforth 格式 [31] 等。所有这些方法都依赖于迭代计算，不可避免地降低了计算效率。为解决这一问题，人们提出了各种非迭代反向追踪策略。一个典型例子是基于 Taylor 展开的方法，该方法无需迭代，只需遍历计算域内的所有网格节点 [32]。作为 SL 方法的另一个核心组成部分，函数插值也经历了持续的发展。广泛采用的插值方案包括 Lagrange 插值、Hermite 插值和 WENO 插值。对于这些基于网格的插值方法，必须设计特殊的单元搜索算法来首先定位包含起始点的单元，然后在该位置执行插值以评估解。这些操作导致较高的计算成本，尤其对于非结构化网格 [29], [33]。由于网格限制是 SL 方法中插值问题的主要原因，消除网格依赖提供了一种直接有效的方式来提高插值质量。无网格方法是一种仅依赖离散节点而不需要进行单元划分的数值方法 [34]。它们不受网格和单元的限制，在处理涉及网格畸变的问题时避免了网格重构。在各种无网格方法中，无单元 Galerkin（EFG）方法因其高精度和稳定性而脱颖而出，成为目前研究与应用最为广泛的无网格方法之一 [35], [36]。尽管 EFG 方法能够产生可靠的数值结果，但在处理对流主导问题时稳定性较差 [37]。将无网格 Galerkin 方法的网格灵活性和高计算精度与 SL 方法的良好稳定性相结合，为解决对流主导流动问题提供了一条有前景的途径。同时，其无网格特性还可以简化 SL 方法中的函数插值过程。在我们之前的工作中，已成功将 EFG 方法与 SL 方法相结合，开发了用于求解对流主导对流–扩散方程 [22] 和不可压缩 Navier–Stokes 方程 [38] 的 SL-EFG 方案。数值结果表明，所提出的方法继承了两种方法的优点。它保持了 SL 方法的强稳定性和大时间步能力，同时保留了 EFG 方法的高精度。此外，它还简化了传统 SL 方法中固有的反向追踪和函数插值过程。本文的目标是将 SL-EFG 方法扩展到不可压缩两相流。采用耦合不可压缩 Navier–Stokes 方程和 LS 方程的模型来描述两相流 [39]。数值解由三个部分组成：不可压缩两相 Navier–Stokes 方程、LS 方程和 LS 重初始化方程。第一部分属于对流–扩散方程，而后两部分是纯对流方程。为了保证良好的数值稳定性，采用 SL 方法进行时间离散，将对流项的效果转化为沿特征线的反向追踪。同时，对扩散项应用 Galerkin 离散以保持计算精度。特别是，借助无网格特性简化了 SL 方法的实施。对控制方程的三个部分采用了时间解耦策略。首先，使用时间分裂算法解耦不可压缩两相 Navier–Stokes 方程中的速度和压力。这使得在每个时间步内可以分别求解流动变量，绕开了速度–压力插值的 Ladyzhenskaya–Babu?ka–Brezzi（LBB）约束，并产生稳定的压力场 [38]。其次，在应用 SL 方法后，LS 方程及其重初始化方程都被转化为常微分方程。相应地，无需进行空间离散，只需沿特征线执行反向追踪。所提出的方法充分融合了 SL 方法和 EFG 方法的优点。SL 方法为对流项提供无条件稳定性，Galerkin 方法为扩散项提供最优逼近。分裂技术实现了速度–压力解耦，SL 方法和 LS 方法的结合确保了界面捕获的良好质量守恒 [40]。此外，无网格特性简化了 SL 方法的实施。因此，所开发的 SL-EFG 方法有望成为求解不可压缩两相流问题的有效算法。本文其余部分组织如下。第 2 节简要介绍不可压缩两相流的控制方程，包括两相流方程、界面演化的 LS 方程和重初始化方程。第 3 节详细阐述基于 SL 方法的控制方程时间离散。首先，在 SL 框架内为不可压缩两相 Navier–Stokes 方程建立分数步方法。其次，推导 LS 方程及其重初始化方程的 SL 离散方案。最后，描述 SL 离散中涉及的反向追踪过程。第 4 节聚焦于将 EFG 方法应用于解耦后 Navier–Stokes 方程的空间离散。第 5 节总结所提出的用于不可压缩两相流的 SL-EFG 方法并给出详细的算法实施。第 6 节采用五个数值算例全面验证 SL-EFG 方法在界面捕获中的质量守恒及其在两相流模拟中的整体性能。最后，第 7 节总结本文并总结主要成果和研究意义。章节摘要两相流控制方程采用不可压缩两相流模型研究两种不相容流体在接触或混合状态下共同流动时的流动特性。它通常由描述流动状态的流动控制方程和表征界面演化的输运方程组成。基于半拉格朗日方法的时间离散由上述控制方程可以看出，不可压缩两相流问题的数值模拟主要涉及求解两相 Navier–Stokes 方程 (8)、LS 界面演化方程 (10) 和界面重初始化方程 (13)。这些方程的求解将遇到对流主导问题。为了获得良好的数值稳定性，采用 SL 方法对所有这些方程进行时间离散。基于无单元 Galerkin 方法的空间离散在 SL 框架内，界面演化方程和重初始化方程均呈现常微分方程的形式，新时间步的未知变量更新可以通过逐节点赋值来完成。然而，对于两相 Navier–Stokes 方程，即使在该框架内，它们仍然保持偏微分方程的形式。

为了计算新时间步下的速度和压力，还需要对空间导数进行离散化。第3节基于SL-EFG方法的计算流程，SL-EFG方法针对前述不可压缩两相流问题包含三个组成部分：不可压缩两相Navier-Stokes方程的求解、界面演化方程的求解以及界面重新初始化方程的求解。这三组方程均包含对流项；特别是后两者为纯对流方程，应用SL方法可以获得良好的数值稳定性。完全离散化数值算例在本节中，给出了五个数值算例，以全面验证SL-EFG方法求解不可压缩两相流问题的有效性。其中，前两个算例的速度场是预先给定的，主要用于检验该方法在捕捉移动界面时的质量守恒性能。其余三个算例涉及界面与流场之间的耦合相互作用，旨在进一步验证结论本研究开发并扩展了一种SL-EFG方法，用于求解不可压缩两相流问题，旨在克服传统基于网格的方法和基于Galerkin的无网格方法在对流主导流动和界面捕捉模拟中的固有局限。所提出的数值框架将用于流场演化的不可压缩Navier-Stokes方程与用于两相界面演化的包含重新初始化方程的LS方程进行耦合。作者贡献声明Ying Liu：撰写——初稿，验证，软件，获取经费，概念构思。Xiaodong Wang：撰写——审稿与编辑，监督，获取经费，概念构思。利益冲突声明作者声明，没有已知的 competing financial interests 或可能影响本文报告工作的个人关系。致谢本工作受陕西省自然科学基金基本研究计划（基金编号2026JC-YBMS-0091）、陕西省高校青年创新团队（基金编号24JP019）以及国家自然科学基金（基金编号12471389）资助。Ying Liu | Xiaodong Wang