JCP：最新多相界面数值处理方法 — 多相流尖锐界面追踪方法综述

12天前浏览799

本文摘要：（由ai生成）

本文综述了模拟多相流和自由表面流动的数值方法，重点介绍了虚流体和Voronoi界面等界面处理方法，并探讨了自适应网格改进技术和并行计算方法的应用。这些方法已成功应用于多个领域，包括艺术领域的计算机图形学。此外，文章还展望了机器学习和深度学习在多相流模拟中的潜在应用。多相流模拟涉及空化、金属凝固等复杂现象，本文所提方法为解决这些问题提供了有力工具，对推动多相流研究具有重要意义。

本文综述了模拟多相流和自由表面流动的数值方法，特别关注能保持跨界面相不连续性质的数值方法。本文讨论了虚流体（Ghost-Fluid）方法和Voronoi界面方法，这两种方法在处理表面张力时，可以避免严格的时间步长限制。此外，本文还探讨了自适应网格改进技术的效率以及并行计算方法的应用。最后本文展示了使用这些方法在二维和三维空间中获得的一些模拟结果。后续还将讨论机器学习和深度学习技术在多相流中的应用，提出利用深度学习来加强多相流研究和模拟的一些未来潜在的研究方向。

引言

多相流在物理和生命科学中无处不在，其模拟也广泛应用于多个领域。其中，关键的模拟案例包括研究船舶工业中的空化现象、了解大多数金属的凝固过程、分析核电站中的沸腾流动、探索新的表面活性剂和液体进入肺部的机制、预测石油和天然气操作中的多孔介质流动等。此外，在艺术领域，模拟也被广泛应用于计算机图形学中，以产生逼真的流动效果。

模拟多相流面临诸多挑战。首先，要追踪不同相之间的界面，这一界面随时间不断变化，并可能经历拓扑结构的演变。其次，一些关键物理量，如流体的密度、速度或压力，在介观和宏观尺度上可能会发生剧烈变化，因此数值方法必须强制实施相应的表示。最后，考虑到多相流的多尺度特性，在某些流体区域可能会出现精细的特征，同时计算资源也有限，因此我们需要考虑在自适应网格和并行架构上采用数值方法。

早期关于多相流处理的研究可以追溯到Peskin的浸入边界法工作[103、102]。尽管这些方法并没有直接考虑多相流本身，但它们利用离散的δ-函数来求解浸没弹性界面中的血液流动，这为将偏微分方程在矩形区域的离散化推广到任意区域奠定了基础。Unverdi等人在此基础上将该方法与界面追踪技术相结合，提出了多相不可压缩流求解器。这种方法允许相与相之间的界面发生大变形，包括拓扑结构的变化。在[129]和[19]中，作者采用了水平集方法来捕捉界面运动，而不是传统的追踪方法，从而减轻了对拓扑变化处理的复杂性。然而，虽然前向追踪和水平集方法能够以尖锐的界面表示方式呈现，与相场模型中将相间界面表示为数值糊状区不同，使用δ-函数公式仍然会平滑掉整个界面上的物理量。在这项工作中，我们重点关注“尖锐”界面方法，即那些在数值上能够保留不连续量中不连续性的方案。

下面考虑一个计算域Ω=Ω^-∪Ω⁺，其中Ω^-和Ω⁺被一个Γ界面隔开。流体的运动由不可压缩的Navier-Stokes方程来模拟：

其中，t是时间，u=(u,v,w)是速度场，p是压力，f包含了例如重力的表面力。本例中在子区域Ω^-,Ω⁺，我们认为流体具有均匀的粘度μ和密度ρ。在计算域边界上和边界上施加适当的边界条件被施加在边界Γ以及计算域的边界

。这些将在本文的剩余部分进行详细的讨论。

求解这些方程的框架是Chorin在1967年提出的标准投影方法[20]，它由以下三个步骤组成。利用动量方程(1)并忽略压力梯度项等计算一个中间速度场，同时时间导数用向前欧拉步离散:

第二步是基于Helmoltz-Hodge分解，即分解一个二次连续可微的有界向量场分解成一个无散度分量

和无旋度分量

：

其中，指Hodge变量，对于

施加不可压缩性约束可以得到一个的方程：

边界条件以及如何施加它们将在整个手稿中描述。

在第2节中，我们回顾了水平集和粒子水平集方法；在第3节中，我们讨论了水平集方法在四叉树/八叉树笛卡尔网格上的并行以及如何求解自由表面流；在第4节中，我们讨论了求解多相流的数值方法；第5节聚焦于自由表面流动，最后给出一个简短的结论。

表面追踪与捕获的对比

利用无质量粒子可以实现相间界面的显式表示(如前向追踪)[42,43,57,58,135]。由于粒子坐标可以高精度地更新，显式方法提供了非常精确的界面演化。另一方面，拓扑结构的变化也需要显式地处理，这增加了这些方法的复杂性。内隐表征以直截了当的方式处理拓扑结构的变化。在多相流背景下，使用最多的隐式方法是流体体积法（Volume of Fluid，VOF）[92,10,9,24,51]、流体矩量法[27]和水平集方法[95,94,120,121,39]。隐式方法将相间界面表示为连续函数(在VOF情况下为一相所占单元格的分数,或在水平集方法情况下为Lipschitz连续水平集函数)的一个等式。VOF方法通过构造保证了质量守恒，但它对界面的不连续描述使得几何量(界面法向,界面曲率等)。的计算变得困难，尽管我们指出Sussman [128]的工作试图耦合水平集和VOF方法。水平集表示法可以精确地计算几何量，但不能保证质量守恒，除非使用精细的网格。为了解决这个问题，引入了两种互补的方法：使用无质量粒子来补充水平集表示和使用自适应网格。在这篇综述中，我们关注粒子水平集方法[31]和自适应八叉树笛卡尔网格上的水平集方法[75,76,85]。

01 水平集方法

给定计算域Ω和分离开内部Ω⁻和外部Ω⁺子域的界面Γ：Ω=Ω⁻∪Ω⁺。水平集借助Lipschitz连续函数φ来描述不同的区域：

根据该表示，界面的外法线n和平均曲率κ表示为:

两个主要方程是:在给定速度场u下演化水平集函数；

并且将任意函数变换成带符号的距离函数，即满足|∇φ|=1的函数。虽然最后一个方程可以通过快速方法（包括其并行对应物）直接求解，例如快速推进法[119，136，18]或快速扫描法[143，25]，但也可以通过求解以下方程将任意水平集函数φ⁰转换为带符号的距离函数:

在伪时间τ内。在最后一个等式中，Sign是指signum函数。

在数值上，由于水平集函数在附近是光滑的，所以我们一般使用中心差分来逼近法线和平均曲率。在一般情况下，方程(6)在时间上用一阶导数的HJ-WENO近似的Godunov格式[56]和三阶TVD Runge-Kutta格式[122]逼近。在u与φ无关的情况下，不会发生激波和稀疏现象，我们往往首先选择使用半拉格朗日方法，因为此时只有精度限制时间步长。重新初始化方程总是用Godunov格式求解。

02 粒子水平集方法

粒子水平集方法[31]的中心思想是引入分别位于φ＞0和φ＜0区域的正负粒子，并使用它们来校正水平集方法中固有的质量损失。每个粒子都有一个半径r_p，定义如下:

其中s_p对于正粒子为+1，对于负粒子为-1，r_min=.1min(Δx,Δy,Δz), r_max=.5min(Δx,Δy,Δz).它们的动态特性由dx/dt=u描述，使用多线性插值在时间上向前积分（通常使用三阶TVD龙格-库塔方案）以在正确位置定义u。

一旦水平集函数和粒子迭代了一个时间步长，逃逸的粒子（例如φ＜0区域中的正粒子）被用于校正水平集函数中的任何错误。具体而言，对于每个粒子p，人们将球面水平集函数与半径r_p相关联:

在当前计算单元的八个角中的任何一个角上，φ与φ_p的任何差异都表明水平集解中可能存在错误。通过将φ变为φ_p及其自身的最小绝对值来进行校正。

实现细节：通常，在两个空间维度中每个栅格单元使用16个粒子，在三个空间维度中使用32个粒子。此外，为了保持平滑的界面，从界面的适当侧删除半径大于其半径1.5倍的逃逸粒子。图1中给出了一个例子，感兴趣的读者可以参考[31，32]，以获得关于粒子水平集方法的更多细节。

自适应笛卡尔八叉树/四叉树网格

一般情况下（湍流要求网格大小与Re^9/4成比例，Re为雷诺数）特别是多相流，在出现不同长度尺度的问题时，需要自适应网格来进行数值仿真，因为相间界面可以延伸成很薄的特征尺寸。图2给出了案例，其中一种流体推动另一种粘度更高的流体，导致粘性指进现象。

基于四叉树/八叉树数据结构的自适应网格特别方便，因为它们可以实现连续的变化细化。此外，在捕捉界面方法（如水平集方法）的背景中，网格的构建是简单、快速的，因此可以在每个时间步进行适应，而不需要太多的计算量。自适应性也有效地解决了水平集方法中的质量损失问题，因为在传播前沿附近可以使用更高的分辨率，而不需要在整个计算域中使用这种精细程度。图1（右）给出了三维空间中的Enright测试示例，说明质量损失可以忽略不计，并且可以分辨薄板流动结构[85]。

在自适应笛卡尔网格而不是均匀网格上离散Navier-Stokes方程，由于存在T型网格点（即在笛卡尔方向之一上缺少一个邻居的网格点），因此会增加一些复杂性。此外，为了保证数值稳定性，有必要定义梯度和散度算子的离散化，使其在离散过程中保持∇=−(∇·)T关系。标准的MAC网格排列提供了一种保证此类约束的直接方法，尽管其他一些工作已经考虑以基于节点的方式求解Navier-Stokes方程[84]。

例如，Popinet通过引入压力方程（5）的离散化，使用八叉树网格离散化Navier-Stokes方程[104]，这导致了一个非对称的线性系统。Losasso等人引入了一个Poisson方程的对称离散化，并将其应用于无粘性Navier-Stokes方程[76]。图3给出了一个模拟的示例，其中使用Navier-Stokes方程对流标量量（烟雾密度）进行平流，可以进行渲染，并应用于自由表面流的情况。在这项工作中，仅使用数值粘性来模拟粘性效应。后来，Losasso等人[75]在八叉树上引入了一个Poisson求解器，该求解器可以产生二阶精确解，同时保持线性系统的对称正定性。最后，Guittet等人使用该求解器求解Navier-Stokes方程，并引入了一种基于Voronoi剖分的考虑粘性应力张量的方法。

01 一阶精确泊松求解器

考虑一般的单元配置，其中一个计算单元与较小的计算单元相邻。压力方程（5）的有限体积离散化导致:

其中指数f指的是当前计算单元的面，A_f指的是它们的面积。在面上定义∇Φ分量是很自然的，例如，我们需要在与u*相同的位置定义Φ_x。在[76]中，作者利用变量位置的一阶扰动导致泊松方程[40]的收敛解的概念，提出了一个简单的∇Φ公式。例如，考虑一个大计算单元（索引为L）位于几个小单元（索引为s）的右侧，则小单元的Φ的x导数由下式给出:

其中，Δ_s是大单元L的中心和小单元s的中心之间在x方向上的距离。∇Φ的所有组件都是为所有计算单元类似构建的。这种离散化产生了一个对称正定线性系统，可以使用标准的PCG求解器有效求解。将该求解器与粒子水平集方法结合起来，[76]获得了烟雾和自由表面流的真实模拟（见图3）。

02 二阶精确泊松求解器

上述离散化仅针对Hodge变量Φ产生一阶精确解。此外，这种离散化在每个单元面上定义了不同的压力梯度。在[75]中，作者采用了不同的方法，并构建了一个大单元的所有相邻单元的Φ_x离散化等于大单元的离散化。因此有:

此外，为了使线性系统对称，Δs被替换为

[75]中给出的数值示例表明，精度阶数为二阶（表1）。后来，[48]证明了精度阶数确实是二阶，表明精度的提高来自于横向方向上的误差相消，这与[86,85]中T型结点值的高阶精确定义构造相似。

03 八叉树的稳定粘性

在自适应网格的情况下，由于MAC网格布局中速度场的分量是交错的，因此使用Voronoi剖分上的有限体积方法对粘性应力张量进行隐式离散化是方便的。在文献[50]中，作者为Vlasov-Poisson和Vlasov-Maxwell方程引入了一个隐式的二阶精确Voronoi求解器。在文献[29,30]中，该求解器被用于在Chimera网格的背景下计算无散度流体的速度。在文献[48]中，Guittet等人利用三维的Voro++[117]，在Quad-/Oc-tree笛卡尔网格的背景下，对应力张量进行了隐式离散化。图4给出了Quadtree网格上的Voronoi划分的一个示例；在这种情况下，速度场x分量之一u0的离散化表达式为：

其中，Voro(u₀)是u₀的Voronoi单元中涉及的邻居集合，d_i是u₀和u_i之间的距离，s_i是与线段[u₀,u_i]正交的Voronoi单元边的长度。将这种离散化应用于所有自由度会产生一个对称正定的线性系统，并且具有二阶精确的解（参见表2，其中考虑了精确解u(x,y)=cos(x)sin(y)）。

04 平行八叉树

在[88]中，Mirzadeh等人介绍了一种用于解决Octree森林中自由边界问题的计算框架。他们没有考虑单一的Octree网格，而是考虑了任意数量的Octree网格，这些网格分布在可用的处理器之间。为了确保负载平衡，库p4est被用来在处理器集上统一分配自由度。具体来说，使用Z-ordering对连续的网格点进行排序和分组；然后直接分布生成的数组。采用第二种算法在每个处理器上构建局部八叉树网格，确保局部八叉树的叶与数组的局部元素相匹配。然后对每个处理器进行八叉树的后续离散化，并使用ghost padding在处理器之间交换数据。

该框架已在Guittet等人的[46]中被用于模拟Octree森林上的单相流体流动。图5给出了在1024个核心上计算的一个球体流过模拟的快照，描绘了Q-criterion的Q = 0.006等位线[54]。在此模拟中，雷诺数为Re=300，八叉树具有(6,9)级，这相当于600万个叶子。作者还在4000个核心上展示了并行缩放。

多相流

01 带跳跃的泊松求解器

多相流问题的解决方案允许在各相之间的界面处发生跳跃。特别是，压力跳跃与表面张力力有关。在第1节概述的投影方法的上下文中，当求解Hodge变量Φ时，会出现这种边界条件，即需要求解以下形式的方程：

其中β=1/ρ＞0在每个子域中被视为常数，但在整个界面上不连续。函数f∈L²(Ω)、g_Γ、h_Γ和k都是给定的。跨界面Γ的量q的跳跃定义为[q]=q⁺_Γ- q^-_Γ。在下文中，我们描述了用于解决上述方程组的虚流体方法及其扩展，即Voronoi界面方法。

1.1 虚流体法

首先考虑一维情况，即试图在网格点x_i处离散以下方程（图6）:

带跳跃的泊松方程[72]的虚流体方法在界面上引入了一种虚拟解，以反映跳跃条件，同时避免在近似导数时出现O(1/Δx)的误差。具体来说，其符合标准中心差分公式:

其中(q)^±R指真实流体的Ω^±区域中的量q，修改为:

其中(q)^±G指虚流体的Ω^±区域中的量q。[Φ_x]中的跳转条件可以用来构建以下方程:

这个简单的变换以一种急剧的方式在Φ_x中强加跳跃，并且只影响线性系统的右侧。同样的程序用于解释Φ的跳跃，给出x_i的最终离散化:

这种方法产生一个对称正定线性方程组，其中跳跃条件只影响系统的右侧。感兴趣的读者可以参考[72]了解更多细节。

在两个空间维度中，等式（10）的原始系统可以被改变如下:

其中n₁和n₂是法向量n的x-和y-分量。虽然该公式与原始问题不同，但它提供了一维情况的简单扩展，并且仍然保留了界面法线方向上的跳跃。事实上，Liu和Sideris证明了这个解是一阶精确的[74]。

1.2 Voronoi界面方法

泊松方程的虚流体方法忽略了切线方向上的跳跃，这导致解梯度通常不收敛。Guittet等人[47]通过修改界面局部的网格解决了该问题，以使自由度与通量正交的线对齐。然后，离散化可以使用一维的虚流体思想。具体来说，删除界面附近的自由度，并确定它们在Γ上的投影。然后在界面的每一侧沿其法线方向构建两个额外的自由度，然后用整个自由度集构建Voronoi分区（因为网格是笛卡尔坐标的，所以只有界面附近的小带发生了变化）。图7说明了这个过程。

然后用有限体积法在每个自由度i上对方程组（10）进行离散化:

其中C是i周围的控制体，Φ_ij是i和j中间的Φ值。由于通过构造Voronoi单元的边缘符合界面，因此可以认为Φ_ij=Φ_Γ.使用一维虚流体方法的思想，网格点i的每个相邻点j对线性系统的左侧的贡献量为:

在它的右手边:

其中是和之间的调和平均值。

正如虚流体方法的情况一样，每个阶段的解与其他阶段解耦，这使得当β、Φ和∇Φ·n中的跳跃任意大时能够计算解。图8说明了(10)的典型解，表3给出了该方法在L^∞-范数下的精度顺序（来自[47]）。本例考虑了以下问题的精确解:

其中θ是（x，y）与x轴之间的角度。系数β取为:

02 多相界面

考虑多种流体（两种以上）的主要困难之一是需要跟踪多种界面。如果考虑每种流体的不同水平集函数，并根据该流体的速度场演化该水平集函数，则数值误差可能会引入无效区域。这在图9中示出，其中产生了重叠和/或真空。必须引入称为预测的其他方法来解决这些情况[81，116，124]。计算机视觉社区经常使用的一种替代方法是使用n个水平集函数的不同符号组合来表示多达2n个区域。然而，这种表示容易产生偏差。

在[77]中，洛萨索等人通过引入投影方法解决了该问题，该方法消除了重叠和真空，保留了带符号距离函数的属性，并保留了准确的界面位置。策略是考虑n个水平集函数来跟踪n种不同的流体，并根据它们的速度场来演化水平集函数。然后通过一个简单的规则更新水平集函数:考虑两个水平集函数φ₁和φ₂（见图10）已经根据它们各自的速度场更新的情况，投影步骤只是逐点减去φ₁和φ₂的平均值。在多个水平集函数φ_i的情况下，投影相当于从所有φ_i中减去最小的两个水平集函数的平均值。图11给出了使用该技术作为主干的瑞利–泰勒不稳定性模拟示例。

03 表面张力

表面张力是重要物理现象的原因，如液滴破裂或马兰戈尼效应。在数值上，显式处理导致标准的O(Δx^3/2)时间步长限制[14]。已经提出了隐式方法来减少时间步长限制。例如，[52]关注曲率随时间的变化，[55]通过使用拉普拉斯–贝尔特拉米算子简化了这种方法（参见[106]中的有限体积实施）。观察到拉普拉斯–贝尔特拉米可以分解为拉普拉斯算子和非线性项的和，[123，142，144]通过隐式离散拉普拉斯算子和显式离散非线性项提出了半隐式离散化。VOF社区中的其他工作已经提出了提高表面张力方法的稳定性[127，70，126]。

在[145]中，Zheng等人介绍了一种混合粒子/网格方法，该方法允许完全隐式离散化，从而取消了时间步长限制。表面由三角形网格近似，表面张力计算为每个元素上边界段的总和。然后可以计算每个元素上的力导数，并合成一个对称矩阵D，该矩阵可以为负半定矩阵，可写成D=-C^TC，C与表面张力系数σ相关（详情见[145]）。

Navier-Stokes方程的时间离散化表示为:

其中是体积加权梯度矩阵，是质量矩阵，每个粒子的质量定义为密度加权控制体积，是Δt-weighted压力矢量，是Δt-weighted表面张力矢量，Δt-weighted重力矢量。这个非线性方程组可以写成:

Newton-Raphson方法用于求解这个非线性系统:

其中是的雅可比矩阵，为了计算近似雅可比矩阵，每次迭代时都将和视为常数。

相对于

的导数是−Δt²D，因此近似雅可比矩阵为:

最后，剩下的是求解对称正定线性方程组:

以及:

从中更新粒子位置:

此外，[145]增加了对体积的限制，因此提出了一种可以隐式处理表面张力的体积守恒方法。这是通过将等式（13）替换为:

离散地加强体积守恒。

图12中给出了隐式处理表面张力的示例，模拟了液滴振荡。特别地，时间步长取得足够大，以便在单个时间步长内获得稳态。图13示出了液滴分裂成几个卫星液滴的模拟。

04 表气泡动力学的整体质量跟踪

当考虑流体中的气泡时，困难之一是简单的数值处理可能导致气泡崩溃。对抗这种人为塌陷的方法已经考虑通过跟踪气泡的封闭体积并使用给定的状态方程来计算气泡内部的压力[125]。然而，由于拓扑结构的变化，气泡云的体积在整个模拟过程中可能会发生剧烈变化，因此跟踪体积非常微妙。在[1]中，Aanjaneya等人提出跟踪质量，这消除了上述困难。具体来说，他们定义了每个气泡内部的质量，并使用[65]的保守方案将其与空气速度一起平流，确保质量包含在[64]之后的零水平集内。此外，从Caiden等人[17]的耦合不可压缩/可压缩流体解算器开始，他们开发了一个整体式压力解算器，该解算器类似于为流固耦合[112，45，41]开发的解算器，并避免了分区方法的标准稳定性问题。该求解器完全耦合了气泡内压力（假设为常数）和不可压缩流体中压力的自由度。然后在标准投影方法中使用产生的压力来更新流体的速度场。气泡中的速度场使用具有诺依曼边界条件的第二泊松解算器更新，以确保气泡速度与流体速度一致。这种处理允许气泡对流体引起的压缩和膨胀力做出反应。

作者在[1]中还包括了表面张力的影响，并表明这种方法可以产生稳定的模拟。图14给出了[1]中的一个例子，其中单个气泡在环境流体中上升，与球形固体物体网络相互作用，迫使其分裂成几个更小的气泡。

自由表面流动

在气体动量效应可以忽略不计的应用中，要考虑自由表面流。在这种情况下，对精度影响最大的步骤是自由表面的压力边界处理。当求解方程(5)中的Φ时应用该边界条件，即求解系统:

并通过边界条件Φ_|Γ=σκ来控制表面张力的影响。在[40]中，引入了带有狄利克雷边界条件的泊松方程的二阶精确离散化，并随后应用于自由表面流的情况[33]。

[40]中介绍的方法修改了自由表面附近网格点的标准隐式中心差分格式。由于它使用了逐维方法，我们仅描述一维常系数泊松方程情况下的处理方法，Φ_xx=u_x*，注意变系数情况的扩展很简单。假设内部网格点x_i与界面相邻，位于x_Γ，x_i＜x_Γ＜x_i+1。在这种情况下，通过定义Φ^G_i+1的虚值并考虑以下修改的离散化

来填充Φ的线性方程组中的相应行，来施加边界条件Φ=Φ_Γ。Φ^G_i+1的定义是使用Φ_Γ作为数据点之一通过外推得到的。所得线性系统的性质和解决方案的准确性取决于外推的程度[37，89]:对于常数或线性外推，相应的线性系统是对称正定的，否则是非对称的（但对角占优）。对于一个常数外推，在解中获得一阶精度，对于每一个额外的外推程度，获得一个额外的阶。在[33]中，作者使用线性外推法来定义Φ^G_i+1:

x_i的离散化变为:

由于线性系统是对称正定的，[33]使用了一种带有修正的不完全Cholesky预条件子的PCG方法[118]。这些作者证明了压力处理对振荡气泡振幅的影响，并表明压力处理对自由表面求解器的精度至关重要。

由于Φ的离散化是隐式的，必须在时间tⁿ⁺¹施加边界条件Φ=σκ，这意味着在计算平均曲率κ时使用水平集函数φⁿ⁺¹。水平集函数的运动由时间tⁿ时的流体速度场给出，为了发展水平集函数，有必要在界面附近的带中定义该速度的有效值。由于流体区域外的速度场未知，因此使用以下方程的几次迭代，在界面的法线方向上不断外推速度场: