网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

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

  • 王小双
  • 张良奇
  • 肖姚
  • 曾忠
重庆大学 航空航天学院,重庆 400044

中图分类号: U448.213

最近更新:2024-06-27

DOI:10.11835/j.issn.1000-582X.2023.259

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

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

两相流在自然界和工业界中广泛存在,例如雨[

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

(ρU)tn*=-mUn+μU+μUTn+fext (1)
1ρn+1ρn+1=1dtU* (2)
ρUt*n+1=-Pn+1 (3)

式中:ρ是密度;μ是动力黏度;U是速度矢量;P是压力;m是质量通量,即ρUfext是外力项(例如重力、表面张力等)。

方程(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方程进行描述。

ϕt+Uϕ=Mψ (4)
ψ=ϕ3-ϕ-ε22ϕ (5)

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

22]ε表示相界面的厚度。

ϕ(z)=tanhz2ε (6)

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

nϕ=0
nψ=0 (7)

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

ρ=ρ1+ρ22+ρ1-ρ22ϕ
μ=μ1+μ22+μ1-μ22ϕ (8)

相场法通过连续表面力模[

22-23]计算表面张力,相较于VOF模型避免了额外的界面重构和曲率计算,σ为表面张力系数

fs=324 σϵψϕ (9)

2 求解算法

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

图1  基于相场法的两相流系统求解流程

Fig. 1  The calculation process of the two-phase flow system based on phase field method

采用同位网格,如图2所示,将所有变量都存储在网格的体心。Cahn-Hilliard方程采用了易于并行化的二阶Runge-Kutta TVD[

24]算法显式计算,本文的后续内容主要关注不可压Navier-Stokes方程的并行解耦,描述并行投影算法,再进一步提出红黑着色的并行投影算法。

图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

2.1 并行投影算法

与标准的压力修正投影算法方程(1)~(3)不同,并行投影算[

14]如下。

Function 1:预测步,由式(1)显式计算中间速度U*

Function 2:校正步,考虑到当d˜t<<dt时,dP/d˜t0,则泊松方程可以修正为

Pn+1t˜=1dtU*-1ρnPn (10)

将此方程用一阶欧拉格式显式推进,循环求解K次(K=10~20),此处,d˜t=ρh2/8h为网格宽度,ρ=minρ1, ρ2,因此,泊松方程得以显式并行。

pn+k+1K-pn+kKd˜t-1ρpn+kK=-U*dtk=0, ..., K-1
pn+k+1K=pn+k+1K+d˜t1ρpn+kK-U*dtk=0, ..., K-1 (11)

Function 3:最终步,更新得到Un+1

ρUt*n+1=-Pn+1 (12)

Adam[

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

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

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

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

1)预测步(黑点):黑色节点的U*通过周围红色节点的值计算得到,分别用R和B作为下标表示红色节点和黑色节点。

UB*=F1URn, PRn, UBn (13)

2)校正步(红点):用上一步计算得到的黑点节点的UB*,更新红色节点的压力。

PRn+1=F2UB*, PBn, PRn (14)

3)最终步(黑点):最后通过下式计算得到黑色节点的速度UBn+1

UBn+1=F3UB*, PRn+1 (15)

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

4)预测步(红点)

UR*=F1UBn+1, PBn, URn (16)

5)校正步(黑点)

PBn+1=F2UR*, PRn+1, PBn (17)

6)最终步(红点)

URn+1=F3UR*, PBn+1 (18)

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

3 算例测试

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

3.1 静态液滴

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

La=σρDμ2 (19)

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

26]和Mirjalili[27]的预期。

图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来量化伪势速度的强度。

Ca=Uμσ (20)

图4展示了网格细化后伪势速度的演化过程(网格尺寸h分别为1/16、1/32、1/64、1/128),得到本算法的空间收敛阶数在一阶到二阶之间。在本文所提出的算法框架下,相场模型比连续表面力模[

23](一阶)更快地收敛到尖锐界面极限,这与Mirjalili[28]的测试是一致的。

图4  压力和表面张力平衡的空间收敛阶数

Fig. 4  Spatial convergence study on the accuracy of pressure/surface tension force calculation

3.2 气泡上升

Hysing[

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

图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图7为不同时刻,气泡的形状。算例1和算例2的物性参数如表1所示。

图6  本算法模拟气泡上升(Case 1)的动态图

Fig. 6  Present algorithm for bubble rise (Case 1)

图7  本算法模拟气泡上升(Case 2)的动态图

Fig. 7  Present algorithm for bubble rise (Case 2)

表1  气泡上升算例的物性参数设置
Table 1  The physical parameters in the bubble rising test
物性参数ρ1ρ2μ1μ2gσ
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×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时的速率曲线,与本文所提出算法的测试结果一致。

图8  气泡质心位置随时间变化的曲线,与参考文献[

31⁃32]的比较

Fig. 8  Comparison of the mass center position between present algorithm and references[31-32]

图9  气泡上升速度随时间变化的曲线,与参考文献[

31⁃32]的比较

Fig. 9  Comparison of the rise velocity between present algorithm and references[31-32]

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

表2  本文算法与串行版本在计算时的效率比较
Table 2  Average CPU time consumption of present algorithm and previous serial version
网格分辨率Case 1Case 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

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,无量纲时间

t*=t(g·H-1)1/2 (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时的形状轮廓,数值测试与实验的对比结果吻合。

图10  溃坝塌陷动态过程,本文算法的模拟结果(左)与Kamra[

37]的实验结果(右)的对比

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 al[

37] (right)

图11  自由液面形状轮廓(t*=5.6

Fig. 11  Free surface shape profile(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, 64): 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, 4377059): 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, 19544/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, 23117): 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 & 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 & 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, 88): 1224. [百度学术] 

15

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. [百度学术] 

16

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, 396): 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, 345): 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, 2262): 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, 2042): 784-804. [百度学术] 

24

Gottlieb S, Shu C W. Total variation diminishing Runge-Kutta schemes[J]. Mathematics of Computation, 1998, 67221): 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, 1441): 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, 4015/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, 6011): 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, 718): 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, 342): 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, 693): 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, 244882): 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, 304): 682-693. [百度学术]