摘要
提出一种交错并行的有限体积投影算法求解基于相场法的两相流控制方程组。该算法的实施主要依赖于压力泊松方程的显式推进设计,从而突破投影算法求解不可压Navier-Stokes方程的效率瓶颈。同时,提出了一种交错扫描策略来更新节点上的变量,以实现更紧凑的时空耦合。本算法与相场模型相结合,能够高效、准确地捕获相界面的动态拓扑变化。测试算例结果表明:网格量为131 072,采用8线程CPU并行时,新提出并行算法的效率达到串行标准投影算法的80倍以上。
两相流在自然界和工业界中广泛存在,例如雨
, | (1) |
, | (2) |
, | (3) |
式中:是密度;μ是动力黏度;是速度矢量;P是压力;是质量通量,即;是外力项(例如重力、表面张力等)。
方程(1)和方程(3)都能显式处理,但泊松方程(2)不含压力瞬态项,只能隐式计算。当采用有限体积法进行离散求解时,隐式格式的每个时间步都需要重新构造系数矩阵并求解,这往往是多相流数值计算效率瓶颈的主要原因。为了提高多相流数值模拟的计算效率,学者们开展了许多研究工作,主要可以分为下述3类方法:第1类是Dong
总之,隐式或半隐式求解不可压多相流时,线性方程组的存在既占用了大量内存,也严重影响了计算效率,这个难题在三维模型数值计算中将更加突出。文中提出了并行投影的有限体积方法,实现不可压Navier-Stokes方程的全解耦显式计算。同时,结合红黑着色的交替并行策略,提高了算法的守恒性。在气泡上升等测试算例中,结果表明,在保证计算精度的同时,本文算法的效率远高于串行的标准投影算法。
选择相场
, | (4) |
, | (5) |
式中:为相分数变量,初始化法向分布遵循反正切函数方程(6),范围为[-1, 1],其中,1表示一相,-1表示另一相,-1~1的分数则表示两相界面。M为迁移率,它定义了相界面区域的扩散率,这里采用1个固定值,M=0.01。是化学势,由自由能泛函推导得
。 | (6) |
文中所有算例中相变量和化学势的边界条件均为Neumann边界
, |
。 | (7) |
流场物性参数由相变量线性插值计算
, |
。 | (8) |
相场法通过连续表面力模
。 | (9) |
单一流体模型下不可压缩两相流的计算需要循环求解Navier-Stokes方程和界面捕捉方程。首先,进行初始化,包括计算域、时间步设置、网格划分以及相界面初始化等;然后,随着时间步进,循环求解Cahn-Hilliard方程和Navier-Stokes方程,求解流程如

图1 基于相场法的两相流系统求解流程
Fig. 1 The calculation process of the two-phase flow system based on phase field method
采用同位网格,如

图2 交错扫描策略:先由相邻的黑点计算出红点,再由计算出的红点更新黑点的值
Fig. 2 Interlaced scanning strategy: the red dots are first calculated by the adjacent black dots, and then update the value of the black dots by the calculated red dots
与标准的压力修正投影算法方程(1)~(3)不同,并行投影算
Function 1:预测步,由
Function 2:校正步,考虑到当时,,则泊松方程可以修正为
。 | (10) |
将此方程用一阶欧拉格式显式推进,循环求解K次(K=10~20),此处,,h为网格宽度,,因此,泊松方程得以显式并行。
, |
。 | (11) |
Function 3:最终步,更新得到。
。 | (12) |
Adam
为了说明交错扫描策略,用F1(·)、F2(·)和F3(·)来表示2.1节所述的3个计算步骤,并将网格节点交替染成红色和黑色,如
基于差分思想的数值方法(FVM, FDM)用差分代替微分来求解偏微分方程。以中心差分格式为例,每个节点的离散只与相邻节点相关。因此,变量的更新可以交替执行,第一步先由红点的值计算黑点,然后再由上一步计算得到的黑点更新红点的信息,如此反复执行。结合红黑着色的思想,对离散点进行并行交替更新,分为以下6个步骤。
1)预测步(黑点):黑色节点的通过周围红色节点的值计算得到,分别用R和B作为下标表示红色节点和黑色节点。
。 | (13) |
2)校正步(红点):用上一步计算得到的黑点节点的,更新红色节点的压力。
。 | (14) |
3)最终步(黑点):最后通过下式计算得到黑色节点的速度。
。 | (15) |
至此完成了一半节点的计算,余下节点的计算采用类似的方式。
4)预测步(红点)
。 | (16) |
5)校正步(黑点)
。 | (17) |
6)最终步(红点)
。 | (18) |
红黑着色并行算法的关键是节点信息的交替更新。完全显式的并行投影算法解耦了节点之间的依赖关系,这使得计算易于并行化,文中使用OpenMP实现并行计算。
通过3个经典的两相流算例测试了所提出算法的精度和效率,所有计算均为二维平面模型。
初始时刻,一个直径D = 0.4 m的静态液滴放置于1 m×1 m二维平面方腔的中心。不考虑重力的作用,四壁均为无滑移边界,两相的密度和黏度比均为1,采用La数作为该算例的特征值。
。 | (19) |
通过改变密度来调整La数,而其他参数包括表面张力系数(1 N/m)和黏度(0.1 Pa∙s)保持不变。设置网格量为64×64,时间步长dt=1×1

图3 静态液滴内部和周围的伪势速度局部放大图(网格为64×64)
Fig. 3 Zoom-in on spurious currents in and around a stationary drop after equilibrium for the diffuse interface simulations performed on 64×64 grids
为了使伪势速度得到充分发展,计算了10 000步,用毛细数Ca来量化伪势速度的强度。
。 | (20) |

图4 压力和表面张力平衡的空间收敛阶数
Fig. 4 Spatial convergence study on the accuracy of pressure/surface tension force calculation
Hysing

图5 较轻气泡在较重介质中上升算例的初始化示意图
Fig. 5 Schematic of the rise of a lighter bubble in a heavier medium
该计算域是一个长度为1 m×2 m的全封闭腔体,气泡初始位置为(x, y)=(0.5,0.5),单位m,初始直径D=0.5 m。顶部和底部为无滑移固壁,左右为自由滑移固壁,重力指向定义域的底部。由于流体1(气泡)的密度小于周围流体2(液体),气泡将在浮力作用下缓慢上升,

图6 本算法模拟气泡上升(Case 1)的动态图
Fig. 6 Present algorithm for bubble rise (Case 1)

图7 本算法模拟气泡上升(Case 2)的动态图
Fig. 7 Present algorithm for bubble rise (Case 2)
物性参数 | ||||||
---|---|---|---|---|---|---|
Case 1 | 100 | 1 000 | 1.0 | 10 | 0.98 | 24.50 |
Case 2 | 1 | 1 000 | 0.1 | 10 | 0.98 | 1.96 |
为了满足CFL条件和适当的界面宽度,不同分辨率下的计算设置略有不同,网格宽度为h=1/64、h=1/128和h=1/256时,时间步长分别为1×1

图8 气泡质心位置随时间变化的曲线,与参考文献[
Fig. 8 Comparison of the mass center position between present algorithm and references[31-32]

图9 气泡上升速度随时间变化的曲线,与参考文献[
Fig. 9 Comparison of the rise velocity between present algorithm and references[31-32]
基于3种不同的网格分辨率来比较本文所提出的并行算法与串行标准投影算法(使用预条件共轭梯度法求解线性方程组)的效率。计算采用的同样的硬件条件(AMD Ryzen 7 4800H,Radeon Graphics 2.90 GHz)和设置。
网格分辨率 | Case 1 | Case 2 | ||||
---|---|---|---|---|---|---|
红黑着色并行算法/s | 标准投影算法/s | 加速比 | 红黑着色并行算法/s | 标准投影算法/s | 加速比 | |
64×128 | 117 | 1 380 | 11.8 | 117 | 1 860 | 15.9 |
128×256 | 263 | 6 900 | 26.2 | 268 | 7 080 | 26.4 |
256×512 | 1 054 | 85 200 | 80.8 | 1 113 | 84 780 | 76.2 |
溃坝算例是海洋与海岸工程水动力学模拟中一个非常经典的基准算例,它与许多复杂的现象有关,比如,地表破裂、高压射流和海洋造浪等。由于其密度比高,重力作用强,且流动会导致剧烈的气液界面拓扑变化,是两相流数值模拟中一个重要且具挑战性的算例。Martin
。 | (21) |
ρL=997,ρG=1.185,μL=8.9×1

图10 溃坝塌陷动态过程,本文算法的模拟结果(左)与Kamra
Fig.10 The dynamic process of dam break, the simulation results of the present algorithm (left) are compared with the experimental results of Kamra, et a

图11 自由液面形状轮廓(=5.6)
Fig. 11 Free surface shape profile(=5.6)
提出了求解Navier-Stokes和Cahn-Hilliard方程的红黑着色并行有限体积投影算法,通过引入人工的压力瞬态项,对泊松方程进行不完全迭代,实现了不可压Navier-Stokes方程的全解耦显式计算。此外,结合红黑着色的思想,交替更新离散点,实现了更紧凑的时空耦合。本算法完全避免了线性方程组的求解,并且可以对每个网格点分别进行计算,只需简单地使用几行OpenMP代码,就可以直接实现节点循环的并行化。
采用静态液滴、气泡上升和溃坝模型这3种经典的两相流基准算例,对本文所提出的算法进行了数值验证。结果表明,本文算法具备精确捕捉大密度比两相流复杂界面变化的能力。更重要的是,在相同的设置和硬件环境下,本文算法的效率远高于串行投影算法,在十万网格量时,加速比约80倍,且计算规模越大,加速效果越明显。本文算法为工业流体中不可压缩多相流的高性能计算提供了新的思路。
参考文献
Barros A P, Prat O P, Testik F Y. Size distribution of raindrops[J]. Nature Physics, 2010, 6(4): 232. [百度学术]
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. [百度学术]
Atencia J, Beebe D J. Controlled microfluidic interfaces[J]. Nature, 2005, 437(7059): 648-655. [百度学术]
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. [百度学术]
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. [百度学术]
Dodd M S, Ferrante A. A fast pressure-correction method for incompressible two-fluid flows[J]. Journal of Computational Physics, 2014, 273: 416-434. [百度学术]
Shen J, Xu J, Yang J. The scalar auxiliary variable (SAV) approach for gradient flows[J]. Journal of Computational Physics, 2018, 353: 407-416. [百度学术]
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 & Mathematics With Applications, 2022, 113: 282-299. [百度学术]
Liu T. A nonlinear multigrid method for inverse problem in the multiphase porous media flow[J]. Applied Mathematics and Computation, 2018, 320: 271-281. [百度学术]
Weymouth G D. Data-driven Multi-Grid solver for accelerated pressure projection[J]. Computers & Fluids, 2022, 246: 105620. [百度学术]
Fuka V. PoisFFT: a free parallel fast Poisson solver[J]. Applied Mathematics and Computation, 2015, 267: 356-364. [百度学术]
Dolean V, Jolivet P, Nataf F. An introduction to domain decomposition methods[M]. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2015. [百度学术]
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. [百度学术]
Adam N, Franke F, Aland S. A simple parallel solution method for the navier-stokes cahn-hilliard equations[J]. Mathematics, 2020, 8(8): 1224. [百度学术]
Kuhn M B, Deskos G, Sprague M A. A mass-momentum consistent coupling for mesh-adaptive two-phase flow simulations[J]. Computers & Fluids, 2023, 252: 105770. [百度学术]
Ding L, Shu C, Ding H, et al. Stencil adaptive diffuse interface method for simulation of two-dimensional incompressible multiphase flows[J]. Computers & Fluids, 2010, 39(6): 936-944. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
Kim J. A continuous surface tension force formulation for diffuse-interface models[J]. Journal of Computational Physics, 2005, 204(2): 784-804. [百度学术]
Gottlieb S, Shu C W. Total variation diminishing Runge-Kutta schemes[J]. Mathematics of Computation, 1998, 67(221): 73-85. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]
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. [百度学术]