红黑着色的相场两相流并行投影算法

摘要：提出一种交错并行的有限体积投影算法求解基于相场法的两相流控制方程组。该算法的实施主要依赖于压力泊松方程的显式推进设计，从而突破投影算法求解不可压Navier-Stokes方程的效率瓶颈。同时，提出了一种交错扫描策略来更新节点上的变量，以实现更紧凑的时空耦合。本算法与相场模型相结合，能够高效、准确地捕获相界面的动态拓扑变化。测试算例结果表明：网格量为131 072，采用8线程CPU并行时，新提出并行算法的效率达到串行标准投影算法的80倍以上。

关键词：两相流；相场法；有限体积法；投影算法；并行计算

中图分类号：U448.213 " " " " "文献标志码：A " " "文章编号：1000-582X（2024）06-075-11

A novel red-black coloring parallel projection algorithm for two-phase flow using the phase field method

WANG Xiaoshuang， ZHANG Liangqi， XIAO Yao， ZENG Zhong

（College of Aerospace Engineering， Chongqing University， Chongqing 400044， P. R. China）

Abstract： In this study， an innovative parallel finite volume projection algorithm with a novel red-black coloring approach is proposed to solve the two-phase flow control equations using the phase field method. This strategy relies on the explicit advancement of the pressure Poisson equation， thus effectively overcoming the efficiency limitation inherent in the traditional projection algorithms for solving the incompressible Navier-Stokes equations. Moreover， we implement an interleaved scanning strategy for updating nodal variables， which significantly enhances the spatiotemporal coupling in a compact manner. The integration of this technique with the phase field method facilitates more accurate capture of interface dynamics and topology at a lower cost. Test results show that， utilizing a grid size of 131 072 and an 8-thread CPU parallelization， the proposed parallel algorithm is more than 80 times more efficient than the serial standard projection algorithm.

Keywords： two-phase flows; phase field; finite volume method; projection algorithm; parallel computing

两相流在自然界和工业界中广泛存在，例如雨滴[1]、喷雾雾化[2]和微流控芯片[3]等。对两相流界面动力学的深入研究，不仅可以帮助人们更好地理解自然现象，还可以为工程技术和生物医学领域的发展提供重要的理论和实践支撑。然而，两相流精确且高效的数值计算仍然是计算流体动力学的重要挑战。当采用基于压力修正的投影算法[4]求解不可压Navier-Stokes方程时，往往将其解耦为3个方程，依次计算得到中间速度 ，当前压力 和当前速度

， （1）

， （2）

， （3）

式中： 是密度；μ是动力黏度； 是速度矢量；P是压力； 是质量通量，即 ； 是外力项（例如重力、表面张力等）。

方程（1）和方程（3）都能显式处理，但泊松方程（2）不含压力瞬态项，只能隐式计算。当采用有限体积法进行离散求解时，隐式格式的每个时间步都需要重新构造系数矩阵并求解，这往往是多相流数值计算效率瓶颈的主要原因。为了提高多相流数值模拟的计算效率，学者们开展了许多研究工作，主要可以分为下述3类方法：第1类是Dong等[5]提出的常系数矩阵时间步进算法求解Navier-Stokes和Cahn-Hilliard方程，避免每个时间步的系数矩阵构造，并结合预处理技术，从而实现大密度比两相流的高效计算；Dodd等[6]也根据常系数矩阵的构造思想和快速傅里叶变换，结合VOF模型加快了大密度比两相流毛细波的计算；Shen等[7]通过引入一组辅助标量，将系数矩阵转变为定常矩阵，实现了高效且无条件能量稳定的梯度流求解。第2类是加速线性方程组本身的收敛过程，文献[8⁃10]采用多重网格法在粗网格和细网格之间反复转换，有效地提高了迭代算法的收敛速度；Fuka[11]发布了基于快速傅里叶变换的开源Poisson方程求解器，可用于基于伪谱法和有限差分法离散线性方程组的高效求解。第3类是并行算法的设计，Dolean等[12]详细介绍了区域分解算法，主要思想是将计算域分解成多个区域并行计算；Djanali等[13]结合区域分解的思想求解泊松方程，加速了投影算法的收敛速度；Adam等[14]通过为泊松方程引入一个瞬态项，在有限差分法框架下实现了不可压Navier-Stokes方程的完全显式并行计算。除了上述的3类方法外，自适应网格[15⁃16]，GPU并行[17⁃18]和人工神经网络[19⁃20]等技术也常被用于改善多相流的计算效率。

总之，隐式或半隐式求解不可压多相流时，线性方程组的存在既占用了大量内存，也严重影响了计算效率，这个难题在三维模型数值计算中将更加突出。文中提出了并行投影的有限体积方法，实现不可压Navier-Stokes方程的全解耦显式计算。同时，结合红黑着色的交替并行策略，提高了算法的守恒性。在气泡上升等测试算例中，结果表明，在保证计算精度的同时，本文算法的效率远高于串行的标准投影算法。

1 界面捕捉方程

选择相场法[21]作为界面捕捉算法，两相界面的运动通过Cahn-Hilliard方程进行描述。

， （4）

， （5）

式中： 为相分数变量，初始化法向分布遵循反正切函数方程（6），范围为[-1， 1]，其中，1表示一相，-1表示另一相，-1～1的分数则表示两相界面。M为迁移率，它定义了相界面区域的扩散率，这里采用1个固定值，M=0.01。 是化学势，由自由能泛函推导得来[22]。 表示相界面的厚度。

。 （6）

文中所有算例中相变量 "和化学势 "的边界条件均为Neumann边界

，

。 （7）

流场物性参数由相变量线性插值计算

，

。 （8）

相场法通过连续表面力模型[22-23]计算表面张力，相较于VOF模型避免了额外的界面重构和曲率计算，σ为表面张力系数

。 （9）

2 求解算法

单一流体模型下不可压缩两相流的计算需要循环求解Navier-Stokes方程和界面捕捉方程。首先，进行初始化，包括计算域、时间步设置、网格划分以及相界面初始化等；然后，随着时间步进，循环求解Cahn-Hilliard方程和Navier-Stokes方程，求解流程如图1所示。

采用同位网格，如图2所示，将所有变量都存储在网格的体心。Cahn-Hilliard方程采用了易于并行化的二阶Runge-Kutta TVD[24]算法显式计算，本文的后续内容主要关注不可压Navier-Stokes方程的并行解耦，描述并行投影算法，再进一步提出红黑着色的并行投影算法。

2.1 并行投影算法

与标准的压力修正投影算法方程（1）～（3）不同，并行投影算法[14]如下。

Function 1：预测步，由式（1）显式计算中间速度 。

Function 2：校正步，考虑到当 "时，，则泊松方程可以修正为

。 （10）

将此方程用一阶欧拉格式显式推进，循环求解K次（K=10～20），此处， ，h为网格宽度， ，因此，泊松方程得以显式并行。

，

。 （11）

Function 3：最终步，更新得到 。

。 （12）

Adam等[14]基于有限差分法和中心差分格式实现了上述算法，有效提高了计算效率，但泊松方程中瞬态项的引入也带来了额外的时间误差，算法的守恒性和精度也有待提高。因此，为了提高并行投影算法的性能，采用守恒性更好的有限体积法，并结合5阶WENO格式[25]离散对流项。此外，还提出了红黑着色的交替并行策略以增强时空耦合。在本文的程序中，瞬态项采用一阶欧拉格式，除了对流项以外的空间算子均采用中心差分格式。

2.2 红黑着色的交错扫描策略

为了说明交错扫描策略，用F1（·）、F2（·）和F3（·）来表示2.1节所述的3个计算步骤，并将网格节点交替染成红色和黑色，如图2所示。

基于差分思想的数值方法（FVM， FDM）用差分代替微分来求解偏微分方程。以中心差分格式为例，每个节点的离散只与相邻节点相关。因此，变量的更新可以交替执行，第一步先由红点的值计算黑点，然后再由上一步计算得到的黑点更新红点的信息，如此反复执行。结合红黑着色的思想，对离散点进行并行交替更新，分为以下6个步骤。

1）预测步（黑点）：黑色节点的 通过周围红色节点的值计算得到，分别用R和B作为下标表示红色节点和黑色节点。

。 （13）

2）校正步（红点）：用上一步计算得到的黑点节点的 ，更新红色节点的压力。

。 （14）

3）最终步（黑点）：最后通过下式计算得到黑色节点的速度 。

。 （15）

至此完成了一半节点的计算，余下节点的计算采用类似的方式。

4）预测步（红点）

。 （16）

5）校正步（黑点）

。 （17）

6）最终步（红点）

。 （18）

红黑着色并行算法的关键是节点信息的交替更新。完全显式的并行投影算法解耦了节点之间的依赖关系，这使得计算易于并行化，文中使用OpenMP实现并行计算。

3 算例测试

通过3个经典的两相流算例测试了所提出算法的精度和效率，所有计算均为二维平面模型。

3.1 静态液滴

初始时刻，一个直径D = 0.4 m的静态液滴放置于1 m×1 m二维平面方腔的中心。不考虑重力的作用，四壁均为无滑移边界，两相的密度和黏度比均为1，采用La数作为该算例的特征值。

。 （19）

通过改变密度来调整La数，而其他参数包括表面张力系数（1 N/m）和黏度（0.1 Pa∙s）保持不变。设置网格量为64×64，时间步长dt=1×10-4。在实际的数值计算中，在相界面附近很难获得精确的数值平衡，会产生所谓的“伪”或“寄生”速度，如图3所示，速度矢量的分布符合Magnini等[26]和Mirjalili等[27]的预期。

为了使伪势速度得到充分发展，计算了10 000步，用毛细数Ca来量化伪势速度的强度。

。 （20）

图4展示了网格细化后伪势速度的演化过程（网格尺寸h分别为1/16、1/32、1/64、1/128），得到本算法的空间收敛阶数在一阶到二阶之间。在本文所提出的算法框架下，相场模型比连续表面力模型[23]（一阶）更快地收敛到尖锐界面极限，这与Mirjalili等[28]的测试是一致的。

3.2 气泡上升

Hysing等[29]发布了一个二维气泡上升的纯数值标准算例，包含有2个测试用例，被广泛用于两相流算法的测试[30⁃31]。计算域和边界条件如图5所示。

该计算域是一个长度为1 m×2 m的全封闭腔体，气泡初始位置为（x， y）=（0.5，0.5），单位m，初始直径D=0.5 m。顶部和底部为无滑移固壁，左右为自由滑移固壁，重力指向定义域的底部。由于流体1（气泡）的密度小于周围流体2（液体），气泡将在浮力作用下缓慢上升，图6和图7为不同时刻，气泡的形状。算例1和算例2的物性参数如表1所示。

为了满足CFL条件和适当的界面宽度，不同分辨率下的计算设置略有不同，网格宽度为h=1/64、h=1/128和h=1/256时，时间步长分别为1×10-4、5×10-5和2.5×10-5，界面宽度 "分别为0.02、0.01、0.005 m。对气泡质心位置和上升速度随时间变化的定量曲线进行比较，与参考文献中的已发表数据吻合[31⁃32]，如图8和图9所示。由于不同的界面宽度下，Case 2中气泡上升的速度略有差别，图9（b）中，Xiao等[31]设置界面宽度为0.02 m时的气泡上升速率曲线，和Aland等[32]设置界面宽度为0.01 m时的速率曲线，与本文所提出算法的测试结果一致。

基于3种不同的网格分辨率来比较本文所提出的并行算法与串行标准投影算法（使用预条件共轭梯度法求解线性方程组）的效率。计算采用的同样的硬件条件（AMD Ryzen 7 4800H，Radeon Graphics 2.90 GHz）和设置。表2中列出了采用8线程OpenMP并行的C++程序在3种网格分辨率情况下的CPU时间。因为完全规避了线性方程组的求解，本文所提出算法的计算效率远高于基于标准投影算法的串行程序，计算规模越大时加速效果也越显著。

3.3 溃坝模型

溃坝算例是海洋与海岸工程水动力学模拟中一个非常经典的基准算例，它与许多复杂的现象有关，比如，地表破裂、高压射流和海洋造浪等。由于其密度比高，重力作用强，且流动会导致剧烈的气液界面拓扑变化，是两相流数值模拟中一个重要且具挑战性的算例。Martin等[33]使用高速摄像机在实验室进行了基准实验，捕捉选定时间间隔的流体运动情况。许多研究人员已经通过数值方法[34⁃36]重复了这个案例，将其作为一个重要的基准算例。文中计算域为0.8 m×0.6 m，液柱长宽为H=W=0.2 m，使用400×300的网格，时间步dt=1×10-7，模拟时间为0.9 s，无量纲时间

。 （21）

ρL=997，ρG=1.185，μL=8.9×10-3，μG=1.789 4×10-5，重力加速度g=9.81 m/s2，忽略表面张力。图10是本文算法与已有实验[37]的比较，结果表明，本文所提出的算法具有较好的精度，具备捕捉实际流动下大密度比两相流界面变形的能力。图11展示了自由液面的前端在t*=5.6时的形状轮廓，数值测试与实验的对比结果吻合。

4 结 "语

提出了求解Navier-Stokes和Cahn-Hilliard方程的红黑着色并行有限体积投影算法，通过引入人工的压力瞬态项，对泊松方程进行不完全迭代，实现了不可压Navier-Stokes方程的全解耦显式计算。此外，结合红黑着色的思想，交替更新离散点，实现了更紧凑的时空耦合。本算法完全避免了线性方程组的求解，并且可以对每个网格点分别进行计算，只需简单地使用几行OpenMP代码，就可以直接实现节点循环的并行化。

采用静态液滴、气泡上升和溃坝模型这3种经典的两相流基准算例，对本文所提出的算法进行了数值验证。结果表明，本文算法具备精确捕捉大密度比两相流复杂界面变化的能力。更重要的是，在相同的设置和硬件环境下，本文算法的效率远高于串行投影算法，在十万网格量时，加速比约80倍，且计算规模越大，加速效果越明显。本文算法为工业流体中不可压缩多相流的高性能计算提供了新的思路。

参考文献

［1］ "Barros A P， Prat O P， Testik F Y. Size distribution of raindrops[J]. Nature Physics， 2010， 6（4）： 232.

［2］ "O’Sullivan J J， Norwood E A， O’Mahony J A， et al. Atomisation technologies used in spray drying in the dairy industry： a review[J]. Journal of Food Engineering， 2019， 243： 57-69.

［3］ "Atencia J， Beebe D J. Controlled microfluidic interfaces[J]. Nature， 2005， 437（7059）： 648-655.

［4］ "Guermond J L， Minev P， Shen J. An overview of projection methods for incompressible flows[J]. Computer Methods in Applied Mechanics and Engineering， 2006， 195（44/45/46/47）： 6011-6045.

［5］ "Dong S， Shen J. A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios[J]. Journal of Computational Physics， 2012， 231（17）： 5788-5804.

［6］ "Dodd M S， Ferrante A. A fast pressure-correction method for incompressible two-fluid flows[J]. Journal of Computational Physics， 2014， 273： 416-434.

［7］ "Shen J， Xu J， Yang J. The scalar auxiliary variable （SAV） approach for gradient flows[J]. Journal of Computational Physics， 2018， 353： 407-416.

［8］ "Lee C S， Hamon F P， Castelletto N， et al. An aggregation-based nonlinear multigrid solver for two-phase flow and transport in porous media[J]. Computers amp; Mathematics With Applications， 2022， 113： 282-299.

［9］ "Liu T. A nonlinear multigrid method for inverse problem in the multiphase porous media flow[J]. Applied Mathematics and Computation， 2018， 320： 271-281.

［10］ "Weymouth G D. Data-driven Multi-Grid solver for accelerated pressure projection[J]. Computers amp; Fluids， 2022， 246： 105620.

［11］ "Fuka V. PoisFFT： a free parallel fast Poisson solver[J]. Applied Mathematics and Computation， 2015， 267： 356-364.

［12］ "Dolean V， Jolivet P， Nataf F. An introduction to domain decomposition methods[M]. Philadelphia， PA： Society for Industrial and Applied Mathematics， 2015.

［13］ "Djanali V， Armfield S W， Kirkpatrick M P， et al. Parallel speed-up of preconditioned fractional step navier-stokes solvers[J]. Applied Mechanics and Materials， 2014， 493： 215-220.

［14］ "Adam N， Franke F， Aland S. A simple parallel solution method for the navier-stokes cahn-hilliard equations[J]. Mathematics， 2020， 8（8）： 1224.

［15］ "Kuhn M B， Deskos G， Sprague M A. A mass-momentum consistent coupling for mesh-adaptive two-phase flow simulations[J]. Computers amp; Fluids， 2023， 252： 105770.

［16］ "Ding L， Shu C， Ding H， et al. Stencil adaptive diffuse interface method for simulation of two-dimensional incompressible multiphase flows[J]. Computers amp; Fluids， 2010， 39（6）： 936-944.

［17］ "Sakane S， Takaki T， Rojas R， et al. Multi-GPUs parallel computation of dendrite growth in forced convection using the phase-field-lattice Boltzmann model[J]. Journal of Crystal Growth， 2017， 474： 154-159.

［18］ "Ha S， Park J， You D. A GPU-accelerated semi-implicit fractional-step method for numerical solutions of incompressible Navier-Stokes equations[J]. Journal of Computational Physics， 2018， 352： 246-264.

［19］ "Qiu R D， Huang R F， Xiao Y， et al. Physics-informed neural networks for phase-field method in two-phase flow[J]. Physics of Fluids， 2022， 34（5）： 052109.

［20］ "Raissi M， Perdikaris P， Karniadakis G E. Physics-informed neural networks： a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[J]. Journal of Computational Physics， 2019， 378： 686-707.

［21］ "Ding H， Spelt P D M， Shu C. Diffuse interface model for incompressible two-phase flows with large density ratios[J]. Journal of Computational Physics， 2007， 226（2）： 2078-2095.

［22］ "Liu H H， Valocchi A J， Zhang Y H， et al. Lattice Boltzmann phase-field modeling of thermocapillary flows in a confined microchannel[J]. Journal of Computational Physics， 2014， 256： 334-356.

［23］ "Kim J. A continuous surface tension force formulation for diffuse-interface models[J]. Journal of Computational Physics， 2005， 204（2）： 784-804.

［24］ "Gottlieb S， Shu C W. Total variation diminishing Runge-Kutta schemes[J]. Mathematics of Computation， 1998， 67（221）： 73-85.

［25］ "Friedrich O. Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[J]. Journal of Computational Physics， 1998， 144（1）： 194-212.

［26］ "Magnini M， Pulvirenti B， Thome J R. Characterization of the velocity fields generated by flow initialization in the CFD simulation of multiphase flows[J]. Applied Mathematical Modelling， 2016， 40（15/16）： 6811-6830.

［27］ "Mirjalili S， Ivey C B， Mani A L. Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[J]. International Journal of Multiphase Flow， 2019， 116： 221-238.

［28］ "Mirjalili S， Khanwale M A， Mani A L. Assessment of an energy-based surface tension model for simulation of two-phase flows using second-order phase field methods[J]. Journal of Computational Physics， 2023， 474： 111795.

［29］ "Hysing S， Turek S， Kuzmin D， et al. Quantitative benchmark computations of two-dimensional bubble dynamics[J]. International Journal for Numerical Methods in Fluids， 2009， 60（11）： 1259-1288.

［30］ "Klostermann J， Schaake K， Schwarze R. Numerical simulation of a single rising bubble by VOF with surface compression[J]. International Journal for Numerical Methods in Fluids， 2013， 71（8）： 960-982.

［31］ "Xiao Y， Zeng Z， Zhang L Q， et al. A spectral element-based phase field method for incompressible two-phase flows[J]. Physics of Fluids， 2022， 34（2）： 022114.

［32］ "Aland S， Voigt A. Benchmark computations of diffuse interface models for two-dimensional bubble dynamics[J]. International Journal for Numerical Methods in Fluids， 2012， 69（3）： 747-761.

［33］ "Martin J C， Moyce W J， Penney W G， et al. Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane[J]. Philosophical Transactions of the Royal Society of London Series A， Mathematical and Physical Sciences， 1952， 244（882）： 312-324.

［34］ "Yu C H， Sheu T W H. Development of a coupled level set and immersed boundary method for predicting dam break flows[J]. Computer Physics Communications， 2017， 221： 1-18.

［35］ "Meng W K， Liao L， Chen M， et al. An enhanced CLSVOF method with an algebraic second-reconstruction step for simulating incompressible two-phase flows[J]. International Journal of Multiphase Flow， 2022， 154： 104151.

［36］ "Li Y L， Wan L， Wang Y H， et al. Numerical investigation of interface capturing method by the Rayleigh-Taylor instability， dambreak and solitary wave problems[J]. Ocean Engineering， 2019， 194： 106583.

［37］ "Kamra M M， Mohd N， Liu C， et al. Numerical and experimental investigation of three-dimensionality in the dam-break flow against a vertical wall[J]. Journal of Hydrodynamics， 2018， 30（4）： 682-693.

（编辑 "郑洁）

doi： 10.11835/j.issn.1000-582X.2023.259

收稿日期：2023-03-02

网络出版日期：2023-07-04

基金项目：国家自然科学基金资助项目（12102071，12172070）；中央高校基本科研业务费资助项目（2021CDJQY-055）。

Foundation：Supported by National Natural Science Foundation of China （12102071， 12172070）， and the Fundamental Research Funds for the Central Universities （2021CDJQY-055）.

作者简介：王小双（1995—），男，硕士研究生，主要从事不可压多相流算法研究与程序开发方向研究，（E-mail）wangxs1995s@cqu.edu.cn。

通信作者：张良奇，男，研究员，博士生导师，主要从事计算流体力学特色数值方法研究及其在前沿流体力学问题中的应用研究，（E-mail）zhangliangqi@cqu.edu.cn。