JCP：全变量笛卡尔网格方法用于不可压缩流和多相流
摘要：本文提出了一种基于笛卡尔网格和多矩概念的不可压缩流体求解器。在所提出的方法中，速度变量位于二维的单元中心、面中心和节点，以及三维的单元中心、面中心、边中心和节点，但压力变量仅位于单元中心。这种笛卡尔网格被称为全变量笛卡尔网格（FVCG），所提出的方法被称为FVCG方法。通过使用FVCG方法，我们可以在不显著增加计算时间的情况下改进速度计算（减少误差和抑制数值振荡）。例如，在三维多相流模拟中，尽管FVCG方法使用的速度分辨率是现有方法的两倍，但计算时间的增加仅为2%-9%。FVCG方法通过各种基准问题进行了验证，并与两个液滴合并和分离的实验以及液滴在超疏水基底上的飞溅进行了比较。这些数值结果表明，FVCG方法能够有效地模拟不可压缩流动和涉及高度复杂表面拓扑变化的自由表面流动，即使在粗网格上也是如此。这些数值结果还表明，FVCG方法是鲁棒的。
不可压缩流体和多相流动的数值模拟在多个研究领域（如工程（机械工程、土木工程、环境工程、航空航天工程、船舶工程等）、科学（物理、化学、生命科学、地球科学等）以及计算机图形学（电影、游戏、艺术等））中发挥着重要作用。在本文中，我们提出了一种基于笛卡尔网格的高效不可压缩流体求解器，并将其应用于一些多相流问题。VSIAM3（基于体积/表面集成平均的多矩方法）[25,28,26]和VSIAM3-TM（带有临时矩的VSIAM3）[39]是基于笛卡尔网格和多矩概念的流体求解器。VSIAM3和VSIAM3-TM已被用于解决各种流体问题，例如牛奶冠现象[34]、海洋流动[31]、液滴撞击干燥表面[41]、液滴在干燥表面上的飞溅[35-37,39]、两个液滴的合并与分离[40,39]等。VSIAM3和VSIAM3-TM是功能强大的流体求解器。据我们所知，其他方法无法像VSIAM3和VSIAM3-TM那样，在[35-37,39]中使用的粗网格上很好地捕捉液滴在干燥表面上的飞溅、两个液滴的合并与分离等物理现象。VSIAM3和VSIAM3-TM的特点是同时使用两种不同类型的变量（单元中心的单元值（单元平均值）和面中心的面值（单元面平均值）），这与标准有限体积方法（仅使用单元平均值）和有限差分方法（仅使用点值），如ENO（本质无振荡）[7]和WENO（加权ENO）方案[15,11]不同。在本文中，我们进一步扩展VSIAM3和VSIAM3-TM的概念，通过添加更多不同类型的变量（在二维中添加节点值，在三维中添加节点值和边缘值），并提出了一种新的流体求解器，称为全变量笛卡尔网格（FVCG）方法。
VSIAM3、VSIAM3-TM和提出的FVCG方法的主要优势之一是，由于附加的自由度（DOF），速度计算得到了改进，因为这些方法使用的速度变量比标准方法多（FVCG使用的速度变量比VSIAM3和VSIAM3-TM多）。缺点是速度的计算成本增加了。然而，在VSIAM3、VSIAM3-TM和FVCG方法中，压力的自由度并没有像标准方法那样增加。压力变量仅位于单元中心，并且压力的计算成本与标准方法相同。众所周知，在使用半隐式方法（压力隐式求解，速度显式求解）的不可压缩流动模拟中，压力计算（线性系统求解器）占主导地位，而速度的计算时间很短。因此，速度的附加自由度不会显著增加总计算时间。这些方法非常适合半隐式方法。
IDO-CF（插值微分算子的守恒形式）方案[10,17]也是一种与提出的方法类似的方法。在IDO-CF方案中，尽管速度的定义类似于FVCG方案，但压力的定义也类似于FVCG方法中的速度。因此，压力的计算成本急剧增加（在二维中，压力变量的数量是提出方案和标准方法的4倍，在三维中是8倍），IDO-CF方案在实际的三维模拟中将过于昂贵而无法使用。所提出的方法不会显著增加计算成本，并且可以用于实际的三维模拟。
FVCG方法（以及VSIAM3和VSIAM3-TM）采用CIP-CSL（约束插值剖面-守恒半拉格朗日）方案[30,16,29,18,1,13,14,38]作为守恒方程求解器。CIP-CSL方案可以有效地同时更新这些速度变量。
在本文中，我们提出了全变量笛卡尔网格（FVCG）方法。第2节给出了FVCG方法的详细内容，并与原始的VSIAM3和VSIAM3-TM进行了一些比较。第3节给出了对流测试（二维正弦波传播和Zalesak问题）、二维不可压缩流动测试（Re=1000、3200、5000和7500的盖驱动腔流）、自由表面流动测试（两个液滴的合并与分离以及液滴飞溅）的数值结果。第4节给出了总结。
在本节中，我们将解释所提方法的细节，以及原始VSIAM3和VSIAM3-TM的一些细节，以便与所提出的方法进行比较。
2.1. 控制方程
在本文中，我们只考虑了不可压缩流动，并使用了以下控制方程
其中 u 是速度，n 是每个控制体积的传出法线Ωe其表面用Γe，ρ 为密度，p 为压力，τ 为粘性应力张量。
2.2. 全变量笛卡尔网格和全变量笛卡尔网格方法
在本节中，我们将详细解释FVCG方法，并阐述FVCG方法与VSIAM3、VSIAM3-TM以及标准方法之间的主要差异。图1（a）显示了二维笛卡尔网格中单元格的索引。在大多数标准方法中，会使用同位网格或交错网格。当使用同位网格时，p、u和v在同一位置定义（例如，在单元格中心pi,j、ui,j和vi,j），其中u和v分别是x和y方向上的速度分量。在交错网格中，p、u和v定义在不同的位置，pi,j位于单元格中心，ui+1/2,j位于一个单元格面上，而vi,j+1/2则位于另一个单元格面上，如图1（b）所示。在VSIAM3以及VSIAM3-TM[39]中，虽然p像标准方法一样仅在单元格中心（pi;j）定义，但u和v不仅在单元格中心定义，还在所有单元格面上定义（ui,j、ui+i/2,j、ui,j+1/2、vi,j、vi+1/2,j、vi,j+1/2），如图1（c）所示。在VSIAM3中，由于为速度分配了比标准方法更多的变量，因此提高了速度计算的准确性。在本文中，我们进一步扩展了这一概念，并增加了速度变量。在二维情况下，我们添加了节点值（u1+1/2,j+1/2、vi+1/2,j+1/2），如图1（d）所示。这样我们可以进一步提高速度计算的准确性。我们将该网格命名为全变量笛卡尔网格（FVCG）。
尽管FVCG方法增加了速度计算的成本，但当使用半隐式方法（该方法隐式求解压力，显式求解速度）时，总计算时间并未显著增加（如后文第3节所示），因为增加的变量不用于压力计算，而压力计算在总计算时间中占主导地位。尽管FVCG方法所需的内存比标准方法（二维中速度变量的内存为4倍，三维中为8倍）和VSIAM3（二维中为4/3倍，三维中为2倍）更多，但这一缺点通过速度计算的改进得到了弥补。
图1. 二维笛卡尔网格示意图。(a) 显示了二维笛卡尔网格中单元格的索引。(b)、(c) 和 (d) 分别展示了交错网格、VSIAM3网格和全变量笛卡尔网格的网格划分和变量分配情况。
图2. 三维全变量笛卡尔网格示意图。(a)显示了三维全变量笛卡尔网格中一个单元的索引。(b)显示了三维全变量笛卡尔网格的网格和变量分配。
在二维全变量笛卡尔网格（2D FVCG）中，对于每个速度分量，使用了四种不同类型的变量：单元值（ui,j）、面值（ui-1/2,j, ui,j-1/2）和节点值（ui-1/2,j-1/2），这些被定义为：
（3）
（4）
（5）
（6）
压力仅在单元中心定义为：
（7）
虽然每个速度分量都使用了四种不同类型的变量，但每个单元都分配了九个变量（1个单元值，4个面值，4个节点值）。
在三维全变量笛卡尔网格（3D FVCG）中，每个速度分量使用了八种不同类型的变量，包括单元值（ui,j,k）、面值（ui-1/2,j,k, ui,j-1/2,k, ui,j,k-1/2）、边值（ui,j-1/2,k-1/2, ui-1/2,j,k-1/2, ui-1/2,j-1/2,k）和节点值（ui-1/2,j-1/2,k-1/2），如图2所示。这里w是速度的z分量。这些变量定义为：（8）-（15）
压力仅在单元中心定义
（16）
虽然每个速度分量使用了八种不同类型的变量，但每个单元分配了27个变量（1个单元值、6个面值、12个边值和8个节点值）。
2.3分步法
分步法的使用方式如下：
（17）
1、对流部分（f A）
（18）
2、非对流部分1（f N A1）
（19）
3、非对流部分2（f N A2）
（20）（21）
这些方程通过FVCG方法求解，其中对流部分由CIP-CSL方法求解。
2.4对流部分（f A）
CIP-CSL方法被用于求解守恒方程
（22）
在这里，Φ是一个标量值（在方程（18）中Φ可以是 u 的每个速度分量）。尽管可以使用任何 CIP-CSL 方案，但本文采用了具有有理函数的 CIP-CSLR（CIP-CSL）方法[29]，它是一种波动较小的 CIP-CSL 方案。
2.4.1. CIP-CSL
这里简要解释一下CIP-CSLR方案。在CIP-CSLR方法[29]中，以下插值函数Φi(x)
（23）
用于在xi-1/2和xi+1/2之间进行插值。系数αi和βi的确定方法如下
（24）（25）
图3. VSIAM3在二维空间中的维度分割方法示意图。
使用以下约束条件
（26）（27）
这里ε是一个很小的数，用于避免零除。本文中所有结果都使用了ε=10-15。通过使用插值函数Φi(x)，可以根据微分形式的守恒方程更新边界值Φi-1/2。
（28）
方程（28）采用以下方法求解
（29）（30）
对于对流方程（29），采用半拉格朗日方法
（31）
方程（30）采用有限差分法求解。单元格值Φi通过有限体积法进行更新
（32）
这里的Fi-1/2是通量
（33）
2.4.2 通过VSIAM3对CSL方案进行多维表述
在这里，我们想回顾一下通过VSIAM3实现的CSL方案的多维方法。在原始的VSIAM3中，为CIP-CSL方案采用了维度分裂方法[28,25]。图3展示了二维表述的示意图。我们现在考虑如何在二维情况下更新ф。如图3所示，ф与u和v位于相同的位置。对于x方向，首先基于фni,j和фni-1/2,j使用一维CIP-CSL求解器（图3中的步骤1）计算ф*i,j和ф*i-1/2,j。但是，由于没有节点值（фi-1/2,,j-1/2和фi+1/2,,j-1/2）用于фi,,j-1/2，因此不能使用一维CIP-CSL求解器更新某些面上的值，如фni,j-1/2。因此，在原始的VSIAM3中，фni,j-1/2通过称为TEC公式的平均过程（图3中的步骤2）进行更新，而无需求解守恒方程，如下所示：
（34）
对于y方向采用了类似的方法。使用一维CIP-CSL方法根据ф*i,j和ф*i-1/2,j计算фn+1i,j和фn+1i,j-1/2（图3中的步骤3）。ф*i-1/2,j通过TEC（步骤4）按以下方式更新：
图4. VSIAM3-TM在二维中使用的维度分裂方法的示意图。
图5. FVCG方法在二维中使用的维度分裂方法的示意图。
（35）
2.4.3. VSIAM3-TM中CSL方案的多维公式
在本节中，我们将回顾VSIAM3-TM[39]。在VSIAM3-TM中，使用了一种不同类型的维度分裂公式，该公式通过创建临时矩（TM）来求解所有矩的方程。在原始的VSIAM3中，当使用一维求解器处理x方向时，虽然我们可以使用一维求解器更新фi,j和фi-1/2,j，但由于没有节点值，因此无法更新фi,j-1/2。因此，在VSIAM3-TM中，我们使用简单的插值来计算фi-1/2,j-1/2和фi+1/2,j-1/2的临时矩：
（36）
如图4所示（步骤1）,一旦计算出临时矩（TM），就可以使用一维CIP-CSL求解器来更新фi,j-1/2（步骤2）。在фi,j-1/2更新后，临时矩将被丢弃。与原始VSIAM3一样，фi,j和фi-1/2,j也在步骤2中更新。对于y方向，过程与x方向几乎相同。对于фi-1/2,j-1/2，我们首先使用以下方式计算临时矩（步骤3）：
（37）
然后使用1D CIP-CSL求解器（步骤4）。该实现方法简单，且易于扩展到3D领域。尽管相较于原始的VSIAM3，VSIAM3-TM的计算成本略有增加，但流体计算的准确性和稳定性得到了显著提升。有关VSIAM3-TM的更多详细信息，请参阅[39]。
2.4.4. 提出的基于全变量笛卡尔网格法的CSL多维公式
在此我们将上述方法进一步扩展，提出了FVCG方法。在VSIAM3和VSIAM3- TM中，采用平均程序(TEC或TM公式)。在本文提出的公式中，我们增加了更多的变量(2D中的节点值，3D中的节点值和边值)，如图1(d)和图2所示。因此，我们可以简单地更新所有变量，而不使用图5所示的任何平均过程，并期望与VSlAM3-TM相比有进一步的改进，因为我们可以在平流部分消除这种平均过程(TEC或TM)，并且与VSIAM3-TM和VSIAM3相比使用了更多的额外变量。
对流部分二维FVCG方法总结:
x方向
1.使用1D CIP-CSL求解器更新单元值（фi,j）和面值（фi-1/2,j, фi+1/2,j）。
2.使用1D CIP-CSL求解器更新面值（фi,j-1/2）和节点值（фi-1/2,j-1/2, фi+1/2,j-1/2）。
y方向
1.使用1D CIP-CSL求解器更新单元值（фi,j）和面值（фi,j-1/2, фi,j+1/2）。
2.使用1D CIP-CSL求解器更新面值（фi-1/2,j）和节点值（фi-1/2,j-1/2, фi-1/2,j+1/2）。
对流部分三维FVCG方法总结:
x方向
1.使用1D CIP-CSL求解器更新单元值（фi,j,k）和面值（фi-1/2,j,k, фi+1/2,j,k）。
2.使用1D CIP-CSL求解器更新面值（фj-1/2,k）和边值（фi-1/2,j-1/2,k, фi+1/2,j-1/2,k）。
3.使用1D CIP-CSL求解器更新面值（фi,j,k-1/2）和边值（фi-1/2,j,k-1/2, фi+1/2,j,k-1/2）。
4.使用1D CIP-CSL求解器更新边值（фj-1/2,k-1/2）和节点值（фi-1/2,j-1/2,k-1/2, фi+1/2,j-1/2,k-1/2）。
y方向
1.使用1D CIP-CSL求解器更新单元值（фi,j,k）和面值（фi,j-1/2,k, фi,j+1/2,k）。
2.使用1D CIP-CSL求解器更新面值（фi-1/2,j,k）和边值（фi-1/2,j-1/2,k, фi-1/2,j+1/2,k）。
3.使用1D CIP-CSL求解器更新面值（фi,j,k-1/2）和边值（фi,j-1/2,k-1/2, фi,j+1/2,k-1/2）。
4.使用1D CIP-CSL求解器更新边值（фi-1/2,j,k-1/2）和节点值（фi-1/2,j-1/2,k-1/2, фi-1/2,j+1/2,k-1/2）。
z方向
1.使用1D CIP-CSL求解器更新单元值（фi,j,k）和面值（фi,j,k-1/2, фi,j,k+1/2）。
2.使用1D CIP-CSL求解器更新面值（фi-1/2,j,k）和边值（фi-1/2,j,k-1/2, ф i-1/2,j,k+1/2）。
3.使用1D CIP-CSL求解器更新面值（фi,j-1/2,k）和边值（фi,j-1/2,k-1/2, ф i,j-1/2,k+1/2）。
4.使用1D CIP-CSL求解器更新边值（фi-1/2,j-1/2,k）和节点值（фi-1/2,j-1/2,k-1/2, фi-1/2,j-1/2,k+1/2）。
在[16]中也提出了类似的方法。尽管在[16]中传输量的分配方式与所提出的方法类似，但速度的分配方式并不与所提出的方法相同（速度采用了一些平均过程）。在所提出的方案中，已经消除了所有平均过程。
2.5. 非对流部分1（f N A1）
粘度项是通过单元格值的标准有限体积公式来计算的。
（38）
我们也对所有其他变量（二维中的面值和节点值，以及三维中的面值、边值和节点值）使用了标准离散化。例如，在二维单相流的情况下，我们使用以下二阶中心差分方案：
这种方法可以直接扩展到三维情况。
2.6. 非对流部分2（f N A2）
我们得到以下泊松方程：
（43）
其中，u*和v*是非对流部分1之后的速度。方程（43）的离散化形式为：
（44）
其中：
（45）
压力泊松方程的离散化与VSIAM3、VSIAM3-TM以及基于交错网格的标准方法中所使用的相同。我们使用了预处理的共轭梯度法（对于二维使用ILU预处理器，对于三维使用多网格预处理器）来求解压力泊松方程。压力泊松方程的收敛容差εp=10-10。通过使用Pi,jn+1，速度的面值更新如下：
（46）（47）
其他速度变量（ui,j, vi,j, ui,j-1/2, vi,j-1/2, ui-1/2,j-1/2, vi-1/2,j-1/2）使用类似TEC公式的平均过程进行更新，如下所示：
（48）-（51）
这种方法可以直接扩展到三维情况。
我们可以通过增加压力变量的数量（如速度）来消除这部分的平均过程，但这肯定会大大增加计算成本，因为压力是隐式求解的。因此，我们对这部分采用了这种平均过程。
2.7. 稳定性条件
本文使用了以下稳定性条件。
（52）
虽然对于FVCG方法、VSIAM3-TM和标准方法来说，该公式基本相同，但在FVCG方法和VSIAM3-TM中，对于un需要额外考虑速度变量。
2.8. 扩展到具有表面张力和重力的多相流
在本文中，我们将FVCG方法应用于一些具有表面力和重力的多相流。在本节中，我们将解释如何使用FVCG方法以及VSIAM3和VSIAM3-TM来模拟多相流。
2.8.1. 具有表面张力和重力的多相流的控制方程
为了处理表面张力和重力，我们在方程（2）中添加了表面张力和重力项。
（53）
其中Fsf是表面张力，g是重力加速度。与第2.3节中所述一样，表面张力和重力项也采用分数步方法（即分别求解每个项）进行求解。
2.8.2. 多相流的分步方法
对于具有表面力和重力的多相流，采用以下分步方法进行求解。
（54）
1、对流部分（f A）
（55）
2、非对流部分1（f N A1）
（56）
3、非对流部分2（f N A2）
（57）
4、非对流部分3（f N A3）
（58）
5、非对流部分4（f N A4）
（59）（60）
2.8.3. 界面捕捉
对于界面捕捉，我们采用了CLSVOF（耦合水平集[19,22]和流体体积[8,12]）方法[21,32]以及THINC/WLIC（双曲线正切界面捕捉/加权线界面计算）方案[27,33,9]。首先，体积分数Ci,j
（61）
通过使用THINC/WLIC方案求解以下对流方程来更新体积分数Ci,j
（62）
其中χ为特征函数，χ的单元平均值即为体积分数Ci,j。然后，根据由Ci,j指示的界面来创建水平集函数ψi,j[21,32]。
密度（颜色）函数фd，用于定义不同材料的物理属性，如密度和粘度，可以通过使用平滑的Heaviside函数来创建。
（63）
（64）
其中2α表示液相和气相之间过渡区域的厚度。在本文中，α= 1.5Δx被用作该值。密度函数被设定为，对于液体фd= 1，对于气体фd= 0。密度ρ和粘度μ计算如下：
（65）（66）
其中ρliquid和ρair分别是液体和空气的密度，μliquid和μair分别是液体和空气的粘度。在本文中，体积分数和密度函数仅在单元中心定义。当需要一些不在单元中心定义的变量（如ρi+1/2,j）时，我们使用简单的平均程序来计算这些变量，如ρi+1/2,j=(ρi,j+ρi+1,j)/2。
2.8.4. 表面张力
对于表面张力，我们使用具有水平集曲率校正的密度缩放平衡连续表面力模型[36,37,40]该公式为：
（67）
其中σ是表面张力系数，фi,jDS是经过密度缩放的密度函数，其计算方式如下：
（68）（69）
ki-1/2,j 是 (xi-1/2, yj) 处的曲率。曲率是通过水平集函数来计算的。
（70）
（71）
首先，我们在单元中心 (ki, j) 计算曲率，并通过从 (xi+1/2, yj) 处最近的液体界面点插值曲率来获得 ki+1/2,j。更多详情，请参阅 [36,37,40]。
通过使用具有水平集曲率校正的密度缩放平衡连续表面力模型，我们可以更新单元面上的某些速度变量。其他速度变量则通过使用在第2.6节中使用的平均程序（48）、（49）、（50）和（51）来更新。
图6. 接触角实现。“虚线代表固体中的虚拟液体界面。接触角由水平集函数表示的虚拟液体界面来考虑。
在一个测试案例中，我们将FVCG方法应用于包含移动接触线的问题。因此，我们想在这里解释如何实现接触角。为了施加接触角，我们使用了Sussman开发的方法[20,41]。该方法的一个重要优点是，我们不需要显式地定位三重点的位置。通过将水平集函数和VOF函数表示的液体界面外推到固体中（如图6所示），来考虑接触角。通过求解扩展方程来外推液体界面。
（72）
（73）
这里θ是接触角。扩展方程通过双线性插值求解。
2.8.5. 重力
重力项的实现很简单，我们只需将所有速度变量加上g△t。
2.8.6. 多相流稳定性条件
本文中多相流采用以下稳定性条件。
（78）
将表面力条件添加到等式（52）中。但是，由于附加条件不包括速度，因此附加条件不受FVCG方法的影响。
我们通过二维正弦波传播测试、Zalesak问题、盖驱动空腔流问题（Re=1000、3200、5000和7500）、无粘性水平剪切层问题、Rayleigh-Taylor不稳定性以及与实验（液滴碰撞与分离[2]、液滴飞溅[24]）的比较，验证了所提出的框架。
3.1 二维正弦波传播
在这个测试中，我们使用VSlAM3-TM和FvCG方法求解了二维守恒方程。初始条件如下：
（79）（80）
定义域为[0,1]x[0,1]，并采用周期性边界条件。使用三种不同的网格尺寸（NxN）=（50x50）、（100x100）和（200x200），其中Δx=Δy=1/N，时间步长Δt=0.2Δx。L1误差定义如下：
（81）
表1 使用CIP-CSLR方法时，二维正弦波在t=1时刻的传播误差。
表1显示了数值结果。FVCG方法的误差始终小于VSIAM3-T'M方法的误差（大约小2/3）。两种方法在此测试中均表现出二阶精度。
3.2 Zalesak问题
Zalesak的带缺口圆旋转测试问题[42]被广泛用作标量平流方法的测试。初始条件如下：
（82）（83）
当使用100x100的笛卡尔网格时，一个旋转周期需要628个时间步长。图7显示了一个旋转周期后VSIAM3-TM和FVCG方法的数值结果（侧视图）。尽管两种结果显然相似，但这些结果表明，与FVCG方法相比，VSIAM3-TM具有轻微的振荡性且更加弥散，如图7（b）、（c）、（e）和（f）所示。如果一维求解器无振荡，则FVCG除数字误差外无振荡。
在此测试中，我们还在计算机（Intel Core i7-9700 3.00 GHz 16 GB内存）上比较了VSIAM3-TM和FVCG方法的计算时间。在二维平流测试中，VSIAM3-T'M和FVCG的计算时间几乎相同（1:1.0038），这可能是因为两种方法都解决了两次一维CSL求解器。
图7. Zalesak问题旋转一周后的数值结果侧视图，分别由原始VlAM3-TM（a-c）和FvcG方法（d-f）得出，使用了100x100的网格。
3.3 盖驱动空腔流
将FVCG方法应用于盖驱动空腔流问题，并将这些数值结果与Ghia等人[61]的数值结果进行了比较。图8显示了FVCG方法在Re=1000、3200、5000和7500下的数值结果（VSIAM3和VSIAM3-T'M的数值结果可以在[39]中找到）。在本次测试中，使用了128x128的笛卡尔网格。FVCG方法能够用相对粗糙的网格很好地模拟Re=5000和7500的空腔流。FVcG方法和VSIAM3-T'M的结果几乎相同，这是因为VSIAM3-TM[39]的结果已经非常接近Ghia等人的解[61]。FVcG方法和VSIAM3-T'M的结果都优于VSiAM3[39]的结果。
图8.通过FVCG方法得到的Re=1000、3200、5000和7500的盖驱动空腔流的数值结果。使用了128x128的笛卡尔网格。实线和点线分别表示在x=0和y=0直线上速度场的x分量和y分量。点表示Ghia等人[6]的数值结果。
3.4 无粘性水平剪切层问题
在本节中，使用VSIAM3-TM和FVCG方法求解欧拉方程。流动由具有小垂直扰动的有限厚度水平剪切层组成[3]。初始条件为：
（84）（85）
采用周期性边界条件。
图9.VSIAM3-TM和 FVCG 方法的无粘性水平剪切层问题的数值结果，参考解使用128x128和256x256网格，t=1.8。VSIAM3-TM和FVCG方法使用128x128的笛卡尔网格进行数值模拟。线条表示士3、土6、士9、士12、士15、18、士21、士24、士27和士30的涡量等值线。
图9展示了当使用128x128的笛卡尔网格时，VSIAM3-TM和FVCG方法的数值结果（涡量等值线图），以及使用128x128和256x256笛卡尔网格在t=1.8时的参考解[3]。两种结果都与参考解合理吻合。看起来VSIAM3-TM和FVCG方法的结果介于使用128x128和256x256网格的参考解之间，并且FVCG方法的数值结果比VSIAM3-T'M更接近于使用256x256网格的参考解。在流动结构简单的情况下（例如，t=0.4和t=0.8），两种方法在流动结构上的差异并不明显，但当流动结构变得更加复杂时（例如，t=1.2和t=1.8），差异变得更加明显。图10展示了VSIAM3-T和FVcG方法的压力分布。与涡量（图9）相比，VSIAM3-TM和FVCG方法的压力分布差异较小，这是因为与速度不同，两种方法中压力变量的数量和压力离散化完全相同。
图10.VSIAM3-TM和 FVCG 方法的无粘性水平剪切层问题的数值结果,使用了128x128的笛卡尔网格。线条表示压力等值线(等值线从-2到1以0.1为增量)
图11. 通过VSIAM3-TM和FVCG方法计算的无粘性水平剪切层问题的数值结果。使用了64x64的笛卡尔网格。线条代表涡量等值线。
图11展示了当使用64x64的笛卡尔网格时，VSIAM3-T'和FVCG方法的数值结果（涡量等值线图）。图11表明，FVCG方法的数值结果比VSIAM3-TM的扩散性更小，且更接近使用128x128网格的数值结果（图9）。图12展示了VSIAM3-TM和FVCG方法的压力分布。图12也表明，FVCG方法的结果更接近使用128x128网格的数值结果（图10）。
图12. 通过VSIAM3-TM和FVcG方法计算的无粘性水平剪切层问题的数值结果。使用了64x64的笛卡尔网格。线条代表压力等值线（等值线以0.1的增量从-2增加到1）。
图13. 无粘性水平剪切层问题的动能随时间变化图。动能已根据初始动能进行了归一化处理。使用了128x128、64x64和32x32的笛卡尔网格。
图13展示了使用三种不同网格分辨率（32x32、64x64和128x128）时，VSIAM3-TM和
FVCG方法的动能时间历程。动能k定义如下：
（86）
在本次测试中，K(t)=1是精确解（能量守恒）。与VSIAM3-TM方法相比，FVCG方法在动能误差方面有所改进。如果网格的分辨率足够高（例如，128x128），则差异几乎可以忽略不计，但是，如果分辨率降低，则差异会变大。随着网格分辨率的降低，误差会增加。这些数值结果表明，在涉及复杂流动模式和有限数值分辨率的场景中，FVCG方法优于VSIAM3-TM方法。与[3]相比，VSIAM3-TM和FVCG方法的误差更大。
图14. 无粘性水平剪切层问题的涡旋强度（涡量平方的积分）随时间变化图。涡旋强度已根据初始涡旋强度进行了归一化处理。使用了128x128、64x64和32x32的笛卡尔网格。
图14展示了使用三种不同网格分辨率（32x32、64x64和128x128）时，VSIAM3-T'M和FVCG方法的涡旋强度（涡量平方的积分）的时间历程。涡旋强度ε定义如下：
（87）
在本次测试中，ε(t)=1是精确解。这些数值结果表明，在所有三种不同的分辨率下，FVCG方法始终优于VSIAM3-TM方法。
3.5. Rayleigh-Taylor不稳定性
计算了具有表面张力的瑞利-泰勒不稳定性[4,5]，以验证所提出的计算具有表面张力的两相流的方法。当一个重流体被一个轻流体支撑以抵抗重力时，会产生Rayleigh-Taylor不稳定性，其中界面的小幅度扰动会随时间呈指数增长，即exp(nt)。增长率n由以下公式给出：
（88）
其中K是扰动的波数，g是重力加速度，A=(ρh-ρl)/(ρl+ρh)，ρl和ρh分别是轻流体和重流体的密度。根据[5]，定义了一个稳定性参数Φ，
（89）
其中σc是n=0时表面张力系数的临界值。对于σc<1，带有扰动的界面会导致不稳定。在本次测试中，ρl=0.25，ρh=1.0，K=1，g=1，μ=0，其中μ是粘度系数。区域宽度和高度分别为2π和6π。通过在界面处的垂直速度分量上施加扰动δvcos(Kx)（本文中使用δv=2×10-3）来引发不稳定性。顶部和底部边界的边界条件是诺伊曼条件，对于侧边界，则使用周期性边界条件。
图15显示了使用三种不同网格分辨率（64x256、32x128和12x64）时Φ=0.25和Φ=0.5的数值结果。实线、虚线和点分别对应从理论获得的梯度、从FVCG获得的振幅和从VSIAM3-TM获得的振幅。随着分辨率的增加，两种方法的解都接近理论的梯度。这一结果证明了两种方法在自由表面流中的有效性。然而，在本次测试问题中，FVCG方法和VSIAM3-TM方法的结果几乎相同。这可以归因于测试问题的简单性，两种方法之间的差异没有显现出来（如第3.4节所示，简单的流动情况几乎没有差异）。
图15. 对于Φ=0.25和Φ=0.5的界面振幅的对数图。实线、虚线和点分别对应从理论获得的梯度、从FVCG获得的振幅和从VSIAM3-TM获得的振幅。
3.6. 液滴的合并与分离
图16. we=23时两个液滴合并与分离的数值结果。(a)、(b)、(c)和(d)分别是使用VSIAM3-TM（Δ=D/8.7）、FVCG方法（Δ=D/8.7）、VSIAM3-TM（Δ=D/17）和FVCG方法（Δ=D/17）得到的数值结果，(e)是相应的实验结果（[2]。时间演变是从右到左。
我们使用VSIAM3-TM和FVCG方法对包含液体界面拓扑变化（液滴的合并与分离）的自由表面流进行了三维数值模拟[2]。图16 (a)、(b)、(c)和(d)分别是使用VSIAM3-TM（网格大小为Δ=D/8.7）、FVCG方法（网格大小为Δ=D/8.7）、VSIAM3-TM（网格大小为Δ=D/17）和FVCG方法（网格大小为Δ=D/17）对韦伯数We=23和雷诺数Re=1820的液滴碰撞进行数值模拟的结果。其中D为初始液滴的直径。(e)为相应的实验结果[2]。在本次测试中，主要使用了64x64x64的网格，除了Δ=D/17的情况（当改变网格大小时，D/10.3情况下的域大小会略有变化，但由于计算边界距离液滴相对较远，因此假设对与域大小相关的现象的影响是有限的）。对于Δ=D/17的情况，使用了128x128x128的网格。在这些数值模拟中，使用了定量参数。密度比为1.25:1000（空气:液体），假设为正面碰撞，韦伯数为23。虽然FVCG方法能够以这种较低的数值分辨率（Δ=D/8.7）捕捉到液滴的分离，但VSIAM3-TM却未能成功。然而，如果将VSIAM3-TM的数值分辨率提高到Δ=D/10.3至Δ=D/8.7，我们就可以在三维空间中减少约40%的网格数。如果压力变量（以及速度变量）也能减少40%，那么差异将不是微不足道的。使用更高数值分辨率（Δ=D/17）的数值结果表明，VSIAM3-TM和FVCG方法在具有表面拓扑变化的自由表面流中表现出合理的收敛性。虽然从高分辨率和低分辨率的数值结果来看，液滴大小似乎略有差异，但这种差异是由于在两次模拟中界面厚度都固定在3个网格上造成的。在低分辨率情况下，相对于初始液滴直径只有8.7个网格，因此3个网格界面平滑的效果更加明显。此外，在低分辨率计算中，平滑的物理长度也会增加。
图17. We=40时，使用VSIAM3-T'M（a）和FVCG方法（b）对中心反射分离现象的数值结果，以及相应的实验结果[2]。时间演变是从右到左。网格大小为Δ=D/7.5。
图17展示了使用VSIAM3-TM（a）和FVCG方法（b）对We=40和Re=240的偏心液滴碰撞进行数值模拟的快照。相应的实验结果可以在[2]中找到。图17的模拟使用了Δ=D/7.5的网格大小。FVCG能够用这样的网格大小捕捉到偏心碰撞，包括分离，而VSIAM3-TM则需要Δ=D/8的网格大小才能再现分离[39]。
图18. 使用PVCG方法对We=56的两个不等大小液滴的反射分离碰撞的数值结果（a）以及相应的实验结果（b） [2]。时间演变是从右到左。网格大小为Δ=D_s/12。液滴的位置与实验略有偏差，但这是因为没有给出液滴的确切速度（只给出了相对速度）。
图18展示了使用FVCG方法对We=56和Re=2840的两个不等大小液滴的反射分离碰撞的数值结果。相应的实验结果可以在[2]中找到。小液滴直径（Ds）与大液滴直径（Dl）的比率为1:2。图18的模拟使用了Δ=Ds/12的网格大小。FVCG还可以以相对较小的网格分辨率模拟We=56的两个不等大小液滴的反射分离碰撞。以We=23为例，我们还比较了在计算机（Intel Core i7-9700，3.00 GHz，16 GB内存）上使用VSIAM3-TM和FVCG方法的计算时间。我们还研究了网格数量对计算时间的影响（所有参数以及Ax、At等网格大小均未改变）。表2和图19显示了结果。关于3D对流部分的计算时间，如表2和图19所示，FVCG比VSIAM3-TM需要更多的计算时间。与二维情况不同，在三维计算中，FVCG的对流部分计算时间比VSIAM3-T'M增加，因为在三维计算中，FVCG需要在每个方向上求解一维CSL求解器4次。另一方面，在二维情况下，VSIAM3-T'M只需要在每个方向上求解一维CSL求解器3次（两种方法都需要求解两次）。然而，如果网格数量增加，对流计算所需的额外计算时间就会减少（表2）。这可能是因为随着数据量的增加，CPU缓存可以被更有效地利用。如果我们比较整个计算时间（包括其他项，如压力泊松方程等），VSIAM3-T'M和FVCG之间的计算时间差异就会变小。如果网格数量增加，差异会进一步缩小，这是因为如图19所示，压力泊松方程占据了总计算时间的主导地位，并且随着网格数量的增加，压力泊松方程需要更多的计算时间。而且，对于这两种方法来说，压力泊松方程的计算成本基本相同（与VSIAM3-TM以及标准方法相比，FVCG中的压力变量并没有增加，如表2和图19所示），而对流部分和其他部分（如粘度、表面张力压力梯度、界面捕捉等）的计算时间相对较小，如图19所示。关于FVCG方法的计算时间，我们的主要结论是，与VSIAM3-T'M相比，FVCG方法的额外时间在整个计算中并不显著。
图19. 对于We=23的两个液滴合并和分离过程，使用三种不同网格分辨率时计算时间的分解。其他项是其他项（如粘性、表面张力、压力梯度、界面捕捉等）计算时间的总和。这些数据已通过VSIAM3-TM对每个网格分辨率的总计算时间进行了归一化处理。
表2：在We=23的两个液滴合并和分离过程中，FVCG方法相对于VSIAM3-TM的额外计算时间。
3.7. 超疏水基底上的液滴飞溅
我们使用FVCG方法对瞬时飞溅进行了数值模拟，相应的实验结果可以在[24]中找到。在比较中，我们使用了定量参数，包括液体密度ρliquid = 1000 kg/m³，空气密度ρair = 1.25 kg/m³，液体粘度μliquid = 1.0 x 10-3 Pa·s，空气粘度μair = 1.82 x 10-5 Pa·s，表面张力σ = 7.2 x 10-2 N/m，重力加速度g = 9.8 m/s²，初始液滴直径Dliquid = 1.86 mm，冲击速度2.98 m/s，以及平衡接触角163°。我们使用了224 x 224 x 48的笛卡尔网格，且Δ = Dsplash/46，其中Dsplash是初始液滴的直径。图20显示了数值模拟结果。FVCG方法能够捕捉到液滴飞溅的物理过程，包括卫星液滴和尖刺。这些数值结果与实验至少达到了定性一致，同时也表明FVCG方法是稳健的。
图20，使用FVCG方法对液滴飞溅的数值模拟结果（a）和相应的实验结果（b）[24]。一个直径为1.86毫米的去离子水液滴冲击在超疏水基底上（平衡接触角为163°）。液滴的冲击速度为2.98米/秒。使用了224 x 224 x 48的笛卡尔网格。
我们提出了全变量笛卡尔网格（FVCG）方法。对流测试的结果表明，FVCG方法在误差、数值扩散和振荡方面优于VSIAM3-TM方法，因为FVCG方法比VSIAM3-TM方法使用了更多的变量（在二维中是节点值，在三维中是节点值和边值），并消除了VSIAM3-TM和VSIAM3对流计算中涉及的一些平均过程。FVCG方法可以模拟各种雷诺数（Re=1000、3200、5000和7500）的腔体流动以及Rayleigh-Taylor不稳定性。无粘性水平剪切层问题的结果表明，在涉及复杂流动模式和有限数值分辨率的场景中，FVCG方法优于VSIAM3-TM方法。与自由表面流动的实验相比，所提出的方法可以模拟各种类型的液滴碰撞和分离（We=23的正面碰撞、We=40的非中心碰撞和We=56的不等大小液滴碰撞）。FVCG方法可以在比VSIAM3-TM更低的数值分辨率下捕捉到液滴分离。同时，也证实了FVCG方法在计算总时间上与VSIAM3-TM相比并没有显著增加。FVCG方法还可以模拟液滴飞溅，并捕捉到液滴飞溅的物理过程，包括卫星液滴和尖刺。这些数值结果也表明，FVCG方法是稳健的。所提出的方法可用于各种不可压缩流动和多相流动的数值模拟。
翻译转载自：《Journal of Computational Physics》：“Full-variable Cartesian grid method for incompressible and multiphase flows”
