JCP | E Olsson经典之作：可用于两相流动的守恒型Level-set方法

5月前浏览2736

本期推送2005年由瑞典皇家理工学院数值分析与计算机科学系的E Olsson和G Kreiss发表在《Journal of Computational Physics》上的论文《A conservative level set method for two phase flow》。该文章提出了一个重要的数值方法，用于处理两相流动问题。涉及两种不同流体介质（如水和气体）的流动，在工程和科学研究中非常常见，涉及到的应用包括但不限于气液混合、油水分离、化学反应器等。

E Olsson和G Kreiss提出的方法有效地解决了界面捕捉在多相流动模拟中的一些关键问题，特别是界面厚度的控制和数值稳定性。多相流界面捕捉涉及到如何准确模拟和跟踪不同流体介质之间的界面。传统方法如VOF（Volume of Fluid）和标准的Level set方法各有优缺点，例如在处理界面合并和断裂时的挑战。Olsson和Kreiss通过他们的守恒型Level set方法提供了改进的数值技术，使得界面的模拟更为精确和稳定。

E Olsson和G Kreiss提出的方法能够在复杂流动的模拟中更精确地追踪界面，特别是在涉及表面张力和拓扑变化的情况下。他们的工作被视为在计算流体动力学领域的一个重要贡献，尤其是在数值模拟技术方面的进步，对后续的研究和应用具有长远的影响。

SIMPOP

一、引言

涉及移动边界和界面的问题在众多应用领域中存在，例如多相流、晶体生长、图像处理、界面传播、流体-结构的相互作用等。已经开发出不同的方法来模拟这些问题。其中一些较常用的方法包括界面跟踪方法、Level-set方法以及用于不可压缩流的VOF方法。

在不可压缩两相流动的模拟中，流体体积（VOF）方法已被广泛使用。在这些方法中，界面由一个color函数隐式给出，该函数定义为每个单元中一个流体的体积分数。从color函数出发，对界面进行重构，然后通过更新color函数隐式确定界面传播。VOF方法是守恒的，并且可以处理界面的拓扑变化。然而，它们通常相对不够精确，由于color函数的不连续性，很难达到高阶精度。据我们所知，没有任何VOF方法的输运数值方案的阶数高于二阶。此外，界面的属性如法线和曲率也难以精确计算。尽管如此，良好的守恒特性仍具有吸引力，并且已开发出相当复杂的方法。关于VOF的早期工作，我们参考Noh和Woodward的研究，对这类方法的综述则可参考Scardovelli和Zaleski。

自由边界问题的另一种处理方法是通过在界面上均匀分布的标记来显式追踪边界，然后输运这些标记。通过这种方式，界面可以被清晰明确地表示。这类方法通常被称为界面跟踪方法。然而，标记可能会聚集或分散，这使得重新分布标记变得必要。在处理拓扑变化时需要特别小心。此外，如果标记彼此独立移动，界面可能会出现振荡。另一个困难是界面与固定的欧拉网格的相互作用，这通常是必需的。所有这些特点使得界面跟踪方法在一般情况下难以实现。

近年来，Level-set方法已变得流行，并已在多种应用中使用，如可压缩和不可压缩的两相流、图像处理和火焰传播等。通常情况下，界面由一个带符号的距离函数的零轮廓——level-set函数来表示。界面的移动由Level-set函数的微分方程控制。一般而言，输运通过(weighted) essentially non-oscillatory (WENO, ENO)来完成。为了保持Level-set函数是一个带符号的距离函数，需要一个重新初始化的过程。这个过程也受控于一个微分方程。Level-set方法自动处理拓扑变化，并且通常很容易获得高阶精度，只需选择具有所需精度等级的ENO或WENO方法即可。

Level-set方法的一个缺点是它们不具有守恒性。对于不可压缩的两相流，可能会发生质量的增加或损失，这在物理上是不正确的。与界面跟踪方法相比，在有限元框架中Level-set方法的质量守恒性较差，这一点已经在前人的工作中被明确指出。很多前人也进行了多次尝试来改进Level-set方法的质量守恒性。譬如为了获得流体体积（VOF）方法的良好质量守恒性，可以尝试使用Level-set方法和VOF方法的组合，但使用Level-set函数来获得更好的曲率近似。需要一个color函数，并且必须像标准VOF方法中一样对其进行输运。由于这个函数在界面处是不连续的，输运此函数时必须特别小心。因此，在不引入振荡的情况下获得高于二阶的输运方案可能很困难。原始Level-set方法的简易性也因此丧失。也有很多学者提到Level-set方法的质量守恒问题，同时提出了一种混合Level-set-标记粒子方法来提高精度，特别是在分辨率不足的区域。然而，在这两种情况下，源自Level-set方法本源的简易性都或多或少地丧失了。

我们的目标是找到一种替代的Level-set函数，以及一个输运方案，从而实现界面所围区域的面积（在三维中为体积）的守恒。速度场被假设为无散度的。为了实现我们的目标，我们使用一个模糊的阶跃函数作为Level-set函数，即在一种流体中为0，在另一种流体中为1。在界面上，它从0平滑过渡到1。Level-set函数的输运通过守恒方案进行，该方案包括一个中间步骤，用以保持界面横截面的形状和宽度不变。此外，我们的Level-set函数将是光滑的，这使得我们的方法容易扩展到更高阶，与不连续的color函数相反。

二、Level-set函数的选择

在标准的Level-set方法中，Level-set函数Φ被定义为一个带符号的距离函数：

其中I是界面，

在界面的一侧大于0，在另一侧小于0。

的输运包括一个重新初始化步骤，以保持

作为距离函数，但这种操作不是保守的，即使在无散度的速度场中也是如此。这意味着由零Level-set界定的区域面积不被保守。这是Level-set方法的一个缺点。

为了表示界面上密度和粘度的不连续性，需要使用Heaviside函数：

在计算中，为了达到数值稳定性，常使用模糊的Heaviside函数。例如：

其中ε对应界面的一半厚度。

如果我们可以选择

我们就无需从Φ计算Hsm。更重要的是，假设我们有一种保守的数值方法来输运Φ(x)，以保持Φ的平滑轮廓。由于该方法是保守的，∫Φ将精确地保证守恒。这意味着我们也可以期待对由Φ=0.5界定的区域AΦ=0.5有良好的保守，因为AΦ=0.5≈∫Φ。如果我们使用sharp Heaviside函数，我们可以精确确保面积的守恒性。然而，在离散网格上，界面的位置通过平滑函数的Level-set更好地近似。因此，界面厚度通常取决于网格大小，以便平滑轮廓可以被网格解析。

我们也可以选择定义一个在 Φ(x)=0.5处的尖锐界面或在 0<Φ(x)<1 范围内的模糊界面。在接下来的讨论中，我们将把Level-set函数称为相场函数，记为Φ。

法线和曲率可以容易地从我们的相场函数中获得，计算公式如下：

三、相场函数的输运

3.1 输运步骤

我们现在转向寻找一个输运Φ的方法，该方法应守恒且不改变界面处Φ的轮廓。假设界面通过给定的速度场u被输运。这对应于接下来的常微分方程，对界面上每个点x都成立：

如同标准Level-set方法中的处理，我们可以改为解决

在整个域上。这将移动Φ的0.5等值层（依据上述方程(2)）。对于不可压缩流，u始终无散度，即∇⋅u=0。于是方程(3)相当于质量守恒定律

在选择合适的数值方法解决这个问题时，必须考虑以下因素：

1.方法应当是守恒的；

2.不应引入虚假振荡；

3.界面的厚度和Φ的轮廓应保持恒定。

使用均匀网格，我们定义网格点

和网格函数

。速度 u=(u,v)假设给定在交错网格上，即u在

网格点上，v在

网格点上。守恒方法可以表示为

其中

在交错网格上近似通量

。使用中心平均计算通量，

这通常会在界面附近引入振荡，如图1(b)所示，显示了几个时间步后的结果。

图1. 0.05、0.5 和 0.95 的Φ轮廓线在初始状态和进行一次旋转后，使用不同数值方法得到的结果。其中(c)-(f)为二阶TVD方法，(b)为二阶但非TVD方法。(a) 初始状态；(b) 中心差分；(c) 逆风差分结合Minmod限制器；(d) 逆风差分结合Van Albada限制器；(e) 逆风差分结合Van Leer限制器；(f) 逆风差分结合Superbee限制器。

大量采用总变差减小（TVD total variation diminishing) 方法，这些方法可以近似质量守恒定律而不在不连续附近引入振荡。线性 TVD 方法在非线性界面的守恒近似会更好。一个非线性 TVD 方法将不会抹去不连续。一个典型的 TVD 方法使用逆风格式和

的分段线性重构。为了时间上的守恒性，可以使用二阶龙格-库塔离散化。使用不同限制器的结果，即不同 Lim(x,y) 在 (A.1) 中显示在图1(c)-(f)中。这里，旋转速度场 u(x,y)=(−y,x)用来测试，确保一个完整的旋转后应精确恢复初始状态。我们使用了 minmod, van Albada, van Leer 和 Superbee 限制器。在 TVD 方案中，使用 Superbee 限制器的方案最为守恒，即产生最小的模糊。从图1(c)-(f)可以清楚看出，这种限制器比其他限制器更好地保持了界面的厚度。

3.2. 中间步骤

即使使用Superbee限制器（见图1(f)），界面的轮廓和厚度仍然不保持守恒。此外，轮廓的形状似乎取决于界面的法线与运动方向的相对关系。为了解决这个问题，我们在每个时间步之后增加了一个中间步骤，以确保界面保持其厚度和形状。正如Harten最初提出的，可以添加人工压缩来保持接触间断的分辨率。这可以被视为一个解决守恒律的中间步骤。

其中

对应于压缩通量。在我们的情况下，我们希望人工压缩通量在0<Φ<1并且在界面的法向方向上发挥作用。为了实现这一点，我们选择

。我们用变量τ来表示时间，以强调这是一个人工时间，不等同于实际时间t。我们注意到（6）是一个双曲微分方程。随着τ的增加，界面上将发展出静止的shocks。

为了避免界面处的不连续性，我们引入了少量的粘性，即我们通过以下方式修改守恒律：

或者以守恒形式表示为：

其中，

通过解方程(8)至稳态，界面厚度将保持恒定并与ε成比例。这在空间中可以通过以下方式近似：

其中F和G是网格单元面上的数值通量。我们选择：

这里的f和g是

的x和y分量。界面的法线

用中心差分法近似计算：

由于每次对流步骤后只计算一次

，因此在达到(8)的稳态后，这个值保持不变。在时间上我们使用相同的Runge-Kutta近似法来对Φ的对流进行近似。所获得的方法在空间和时间上都是二阶精确的。

因为我们使用了显式时间步进，因此由于粘性项，我们在Δτ上得到稳定性限制，通常为

实验中，我们发现通过选择C=1/4获得稳定性。我们选择ε，它决定了界面厚度，取决于网格尺寸如下方式：

如果d=0，ε和因此界面的宽度将与Δx成比例。在这种情况下，即使进行网格细化，平滑的界面轮廓也不会增加。尽管如此，我们还是获得了二阶精确性，即使对于d=0在旋转圆的测试计算中也是如此。对于更复杂的流场，如第5.2节中的涡旋测试，我们必须选择d>0才能获得收敛。但是，一个小的值，d=0.1，就足够了。在所有计算中，我们使用了：

作为达到稳态的标准，我们使用了：

对于某个指定的容差TOL。数值测试表明，实际上只需几个时间步就可以达到稳态。

使用不同的对流方法以及人工压缩后，进行一次旋转的结果在图2.1中展示，即使是使用中心差分的方法也是如此。每种情况下界面的厚度都保持恒定。

图2. 使用不同方法以及人工压缩步骤后，Φ的轮廓线在最初和一个周期后（t = 2π）。(a) 初始状态；(b) 中心差分法；(c) Van Albada 方法；(d) 顺风向 Minmod 方法；(e) 顺风向 Van Leer 方法；(f) 顺风向 Superbee 方法。

3.3 Φ的边界和初始条件

Φ的适当边界和初始条件也必须指定。在我们的计算中，要么在边界上设定Φ=0，要么假设边界为与接触角θ=π的墙。在墙的边界上，应满足：

其中，

是界面的法线，

是墙的法线。这可以转化为Φ的均质Neumann条件：

在t=0时，必须对Φ进行初始化。一种方法是在界面的一侧将Φ设为1，在另一侧设为0。然后数值求解方程（7）直到达到稳态。得到的函数用作Φ的初始数据。然而，在某些情况下，（7）的稳态可以解析找到。对于稳态，应满足

，并且Φ的0.5等高线应与所需界面重合。如果

满足x≠xc时的Φ对应于以

为中心、半径为r的圆。对应于y=yint的水平界面，我们可以设置

以同时满足稳态和在y=yint获得0.5等高线。

四、不可压缩两相流

下一步是将我们的输运方法与两相不可压流求解器耦合。我们采用一种扩散界面模型，在此模型中，表面张力转化为在几层单元内分布的体积力。界面上的密度和粘度不连续性也被平滑处理。界面上每一点

的表面张力为：

选择在任一点

的体积力为：

这将产生与

相同的总力，但在有限的界面宽度上分布。在此，保持过渡层厚度恒定非常重要。只有在界面厚度趋向于零时，

才与

相等。如果界面变宽，则

将不再是一个好的近似，因为

的精确数值计算将因需要到二阶导数而变得困难。

无量纲化的不可压缩纳维-斯托克斯方程与表面张力和重力为：

其中，

是雷诺数，

是弗劳德数，

是韦伯数。

是重力方向的单位向量，ρ和μ是非量纲化的密度和粘度。密度和粘度在界面上平滑变化，定义为：

其中ρ1，ρ2和μ1，μ2是两种流体的密度和粘度。

我们使用交错网格进行离散化，即Φ和p定义在网格点

处，而u定义在

处，v定义在

处。方程(10)和(11)使用Marker和Cell方法的扩展数值解，结合我们的保守型形状保持对流方案求解方程(12)。压力通过求解具有变系数的泊松方程隐式更新。该线性系统通过直接带状求解器求解。然后显式更新速度，并最终将Φ与计算出的速度一起输运。我们的方法在空间上具有二阶精度，在时间上具有一阶精度。

请注意，由于我们使用守恒型方法更新Φi,j，即

被精确保守，因此密度也可以确保守恒，即

也被我们的方法确保精确守恒。

五、结果

我们在四个不同的工况下测试了我们的方法。首先，我们对给定无散度速度场的输运方法进行了两次测试。包括对收敛性的研究，评估了顺序准确性。然后，输运被耦合到第4节描述的纳维-斯托克斯求解器中进行测试，考虑了界面因流体流动而发生的物理变化。最后，一个测试案例是考虑通过空气的水滴下落并撞击水面。

在最后两个模拟中，针对所有边界的切向分量施加了齐次诺依曼条件。所有边界对Φ也使用齐次诺依曼条件，即接触角为π。从这些边界条件和纳维-斯托克斯方程得到，压力梯度 ∂x=0在垂直边界上，∂p/∂y=ρ/Fr2在水平边界上。

5.1 旋转圆

为了测试我们的输运方案，我们在常速度场(u,v)=(y,−x)中旋转了一个圆。在一次旋转后，即t=2π时，不同网格上的解进行了比较。在（7）式中，粘度参数ε=Δx/2，以使得随着网格变细，过渡层的厚度逐渐减小。在四种不同的网格上，Δx=0.08，0.04，0.02，0.01，以及Δt=Δτ=Δx/2，解在t=2π时对应于Φ=0.05，0.5和0.95的等值线在图3(a)中显示。在每个时间步后执行了四次人工压缩步骤。在边界上我们使用了Φ=0。

图3. 不同网格上的结果，使用Superbee的迎风格式：(a) 在 t = 2π 时对应于 Φ = 0.05、Φ = 0.5 和 Φ = 0.95 的等值线；(b) 面积守恒。

虽然圆的面积保持在0.5轮廓线内有轻微的变化，但是实际上没有漂移，即从初始面积的最大偏差不随时间增加。在最粗糙的网格上，我们得到的最大偏差小于0.5%，而在最精细的网格上仅为0.035%。不同网格的初始面积之间的差异是由于初始Φ的离散化误差造成的。显然，这个误差是(Δx)2级别的。

从结果中，我们估计了圆位置和误差相对于以下公式测量的准确性顺序：

其中：

这样，我们可以测量由Φ=0.5定义的锐利界面的误差。气泡的位置由质心定义：

计算出的精度顺序在以下表格中显示。

在之前的计算中，我们使用了迎风格式和由Superbee限制器定义的分段线性重构。完全相同的计算也使用了中心格式，即公式(4)配合公式(5)。所获得的结果显示在图 4(a) 和 (b)中。

图 4. 不同网格上的结果，中心格式：(a) 在 t = 2π 时对应于 Φ = 0.05、Φ = 0.5 和 Φ = 0.95 的等值线；(b) 面积保持。

该结果与使用Superbee格式获得的结果非常相似。对于每种方法，图 5(a) 和 (b) 显示了对应于 Φ 的 0.5 等值集的轮廓。

图 5. 使用不同网格和不同数值方法的 Φ 的 0.5 等值线：(a) 迎风配Superbee；(b) 中心格式

获得了以下精度顺序：

因此，尽管在最粗的网格上迎风格式表现更好，我们仍然为中心格式获得了更高的精度顺序。

5.2. 涡流测试

为了在更复杂的流动中测试我们的输运方案，我们在单位正方形上使用了以下速度场：

以（0.5,0.75）为中心，半径为0.15的圆形被用作初始条件。在t=T时，流场被反转，以便在t=2T时精确解应与初始条件一致。为了获得收敛，我们必须选择ε=(Δx)0.9/2。在四个网格（322、642、1282、2562）上进行了T=1和T=0.5的计算。T=1时在t=0、t=0.5、t=1和t=2的解显示在图6中。T=1时由Φ的0.5等值线界定的区域的面积保持在图7(a)中给出，而t=2时Φ的0.5等值线在图7(b)中显示。T=0.5的相应结果显示在图8(a)和(b)中。

图6. 四种不同网格上的涡流测试，在t=0，t=0.5，t=1，t=2，T=1

图7. 四种不同网格上使用我们的方法进行的涡流测试，T=1：(a)面积保持；(b)t=2时Φ的0.5等值线。

图8. 四种不同网格上使用我们的方法进行的涡流测试，T=0.5：(a)面积保持；(b)t=1时Φ的0.5等值线。

这些结果可以与使用标准level-set方法并重新初始化的结果进行比较，参见图9，图10。对于输运和重新初始化，均使用了二阶ENO方案。

图9. 四种不同网格上使用标准level-set方法的涡流测试，T=1：(a)面积保持；(b)t=2时Φ的0等值线。

图10. 四种不同网格上使用标准level-set方法的涡流测试，T=0.5：(a)面积保持；(b)t=2时Φ的0等值线。

在图6中，我们注意到当被拉伸的圆的厚度接近界面的厚度时会发生断裂。这是一个数值效应，只有当界面的厚度小于两个界面之间的距离时才能避免。我们还从图7(a)中注意到这种断裂导致了小的暂时性质量损失。然而，与标准水平集方法相比，这种质量损失很小，并且在t=2T时质量得到恢复。对于T=0.5，界面保持良好解析，质量保持非常好（图8(a)）。与标准方法相比，我们的方法在质量保持方面明显更好。（注意图7、图9、图8、图10(a)中的比例尺差异。）无论界面是否被良好解析，这一点都成立。

就像上一节中的旋转气泡一样，我们估计了在t=2T时相对于由公式（14）定义的误差的精度顺序：

最后，我们在最细的网格上进行了一个非反转模拟，直到t=4。在t=1、t=2、t=3和t=4的解决方案显示在图11中。显然，网格不够细致以解析界面。

图11. 长时间模拟

5.3 上升气泡

研究了一开始静止在水中的气泡。参考密度和粘度被设定为水的密度和粘度：ρref=1.0×103kg/m3，μref=1.0×10−3Ns/m2，使得ρ1=1，ρ2=0.0013，μ1=1以及μ2=0.016。设定σ=7.3×10−2N/m，lref=5.0×10−3m和uref=0.1m/s，我们得到Re=500，Fr=0.45和We=0.68。界面厚度由ε=(Δx)0.9/2决定。我们选择时间步长以稳定性视角来考虑粘性和对流项。

在t=0.5时不同网格的结果，Δx=2/25，Δx=2/50，Δx=2/100和 Δx=2/200在图12中展示。由

的0.5等高线围成的区域的时间演变在图13(a)中展示。该区域的保持非常好。即使在最粗的网格上，面积的波动也仅约0.1%。在其他测试中，初始面积的差异是由初始

的离散化误差引起的。质心的速度在图13(b)中展示，而

的0.5等高线在图14中展示。实线对应 Δx=2/200，虚线对应Δx=2/100，点划线对应Δx=2/50，点线对应Δx=2/25。我们观察到等高线和速度的收敛，尽管收敛速度相对较慢。可能的原因是表面张力、粘度和密度的模糊化。

图12. 在四种不同网格上，上升气泡在t=0.5的结果

图13. 气泡质心的速度和由 Φ=0.5围成的区域面积：(a)面积保持；(b)速度

图14. 在四种不同网格上，上升气泡在t=0.5的结果。

5.4. 下落的水滴

最后，研究了涉及界面拓扑变化的问题：一个小水滴穿过空气直到撞击水面。初始状态如图15所示。

图15. 下落水滴的初始状态

与前一节相同，我们设定ρref=1.0×103kg/m3，μref=1.0×10−3Ns/m2，ρ1=1，ρ2=0.0013，μ1=1，μ2=0.016以及σ=7.3×10−2N/m。通过选择lref=1.0×10−3m和uref=1.0×10−2m/s，我们得到Re=10，Fr=0.10和We=0.0014。离散化参数被选择为Δx=0.06，Δt通过稳定性决定，如公式(17)。

下落水滴的结果展示在图16和图17中。我们看到在空气中水滴保持相当圆润。这是符合预期的，因为在这种小尺度上，表面张力很大。水滴接近时，水平面保持直线。人们可能期望水滴下的空气会在水平面上形成一个小凸起。然而，在这种小尺度上，这些效应太小而在结果中无法明显看到。我们还注意到，当泡沫与表面的距离接近于界面的厚度时，水滴会轻微地吸引表面。这是由界面的模糊表达所引起的数值效应。当水滴撞击水面时，会产生朝向墙壁传播的波浪。可以看到我们的方法能够很好地处理界面的拓扑变化。在图18中，展示了对应于Φ=0.05，Φ=0.5和Φ=0.95的三条等高线。我们注意到，即使水滴撞击表面时，界面也保持了其厚度。

图16. 从t=0到t=3.5，下落的水滴，Φ=0.5的等高线

图17. 接近表面的下落水滴

图18. 接近表面的下落水滴，Φ=0.05，Φ=0.5和Φ=0.95的等高线

六、结论

我们构建了一种在无散度速度场中传播界面的数值方法。该方法具有守恒性，且扩散界面的厚度保持恒定。我们的方法易于实现，且扩展到三维是直接的。关于拓扑变化，无需特别注意，因为这已自动融入方法中。我们使用了二阶近似，数值测试也显示实际精度约为二阶。与标准水平集方法相比，质量守恒显著提高。

由于我们的方法基于某种平滑水平集函数与一组微分方程，其他数值方法也可以轻松应用。例如，可以使用更高阶精确的和/或有限元离散化。通过这种方式，应该可以构造一个大于二阶的守恒方法。

翻译转载自：《Journal of Computational Physics》：“A conservative level set method for two phase flow”

来源：多相流在线

断裂非线性多相流化学 UM 控制

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2024-06-16