韦志龙 蒋 勤

(河海大学港口海岸与近海工程学院,南京 210024)

引言

二相流是自然界中常见的流体运动形式.其中,水气二相流与水利工程、海岸工程和海洋工程等实际工程问题密切相关.水气二相流因具有界面变形复杂、密度差大和紊动强烈等特点,需要较高的界面追踪和控制方程对流项的数值求解精度,是计算流体力学的难点和热点问题.相关研究结果表明,当流体速度低于0.3 倍的声速时,流体运动可视为不可压缩流进行处理[1].因此,在对开敞水域的自由表面流工程问题进行数值模拟分析时,为简化模型,一般情况下可以忽略空气的可压缩性,将其视为不可压缩流体.据此,本研究的水气二相流的模型亦采用水和空气的不可压缩假定.

如上所述,水气二相流中界面两侧的介质密度相差较大,而密度的微小变动将会引起界面两侧压强的剧烈变化;此外,水气二相流还具有紊动剧烈的特点;这些均要求所选取的对流项的求解方法在具有较高的精度的同时,能够很好地处理变量的突变问题.加权基本无震荡(weighted essentially nonoscillatory,WENO)格式[2]根据变量变化曲线的光滑程度选取加权值.对于2r−1 阶WENO 格式,在曲线光滑段具有2r−1 阶精度,在间断点处退化为基本无震荡(essentially non-oscillatory,ENO)格式[3],并保持基本无震荡特性,被广泛用于实际工程中流体运动问题,特别是空气动力学问题的数值模拟[4-6].相对而言,WENO 格式在水气二相流中的应用成果较少.鉴于其高精度、自适应、标准化等特点,本文选取WENO 格式离散对流项,同时采用三阶总变差递减(total variation diminishing,TVD)龙格库塔(Runge-Kutta,RK)法[7]对时间项进行离散,求解流体运动的动量方程.

另一方面,对二相流的界面处理以VOF[8](volume of fluid)法应用最为广泛[9-11],VOF 模型理论上具有质量守恒、物理意义明确等特点,可分为几何重构法和代数法[12].在多维情况下,前者的界面几何重构过程过于复杂且计算效率低下,相比而言代数法对界面追踪的计算过程更为简单.Xiao 等[13]提出的双曲正切函数界面捕捉法(tangent of hyperbola for interface capturing,THINC)为典型的代数法.该方法采用分段双曲正切函数重构界面,将流体各物理量沿界面法向方向的变化视为连续函数,通过模型中的界面形态参数控制不同介质界面的厚度,进而确定界面网格的体积分数通量,不仅大大简化了界面追踪的计算流程,而且能够保证流体质量守恒的计算精度.此外,Yokoi[14]提出了将加权线性界面计算(weighted line interface calculation,WLIC)和THINC 法相结合的改进方法,以取代之前将一维THINC 法直接用于多维界面运动问题求解的作法,有效地提高了THINC 法对多维界面追踪的计算精度.THINC/WLIC法的计算过程简单,适用于二维和三维问题的流体界面追踪.除此之外,还有一些对THINC 法的改进方法[12,15-17],但大多包含多维双曲正切函数的重构,涉及的参数过多,程序实现较困难,且计算量亦相对较大.

迄今,尚未见到国内外关于将WENO 格式与THINC 法进行耦合的相关模型研究成果.考虑到WENO 格式的高精度、自适应等特点适用于二相流的对流项的计算,THINC/WLIC 法在自由表面追踪时既简单易行又具有较高的求解精度的特点,适合于复杂界面的追踪计算,本文尝试将两者进行耦合,建立WENO-THINC/WLIC 耦合水气二相流数值模型,并通过对典型流体力学问题的模拟分析与模型验证,探索适用于大多数工程领域中水气二相流运动模拟的高精度数值计算方法.

1 控制方程

本模型中水气两相均视为不可压缩,流体运动的控制方程为NS 方程,其二维形式可表示为

式中,i,j=1,2;ui为速度分量;fi为质量力,在本模型中为重力;τi j为总应力,对于牛顿流体有τi j=−pδi j+2µSij−2µδijSkk/3,其中p为压强,Sij=(∂ui/∂xj+∂uj/∂xi)/2,δi j为Kronecker 算子.

基于VOF 方法追踪水气交界面,水和空气的体积分数αm满足方程

式中,m=1,2;αm满足α1+α2=1.

流体的密度和黏性系数通过下式计算求得

式中,λ 为密度ρ 或黏性系数µ,λ1和λ2分别为对应的水体或气体的特性参数.

2 数值方法

2.1 分步计算法

采用分步计算法[18]对控制方程进行离散求解,计算过程分为3 个阶段:对流项、非对流项(I)、非对流项(II).从时间刻tn到tn+1的具体求解过程如下.

(1)对流项

方程右侧压强梯度∂pn+1/∂xi使用中心差分法离散.

以上为单个时间步的计算流程,单步计算结束后返回步骤(1),根据设定的CFL 条件更新∆t,如此循环可以得到设定时段的计算结果.

2.2 WENO 格式

式(5)的一般形式为

式中,ϕ 为变量u或v.

现以ϕ 在x方向的离散格式为例说明WENO 格式的计算过程

2.3 时间项离散

本文采用三阶TVD RK 格式对时间项进行离散.具体过程如下

式中,ϕ 为基本变量,本模型中为ui;L 为式(1) 和式(2) 中对流项和等号右侧对应的数值解;角标(1)和角标(2)为当前时间步tn和下一时间步tn+1之间的子时间步.

2.4 THINC/WLIC 算法

将水气二相流中水的体积分数输运方程式(3)改写为守恒格式并忽略角标可得

式中,ϕ 为水的体积分数.

其一维形式为

在WLIC 方法中,多维界面重构是根据界面法向量分别对水平界面和竖向界面赋权重来实现的.在主方向使用一维THINC 法计算体积分数通量,在次主方向使用简单线性界面法(simple line interface calculation,SLIC) 计算,从而实现计算量和精度之间的平衡,以较小的计算量获取较高的精度.如图1 所示,设在网格(i,j)内,体积分数为αi,j.根据WLIC 理论,该网格内界面重构表达式为

式中,n为界面法向量,χ 为界面特征表达式,χx和χy分别为竖向界面和水平界面的特征表达式,ωx和ωy为对应的权重.

图1 界面加权分解Fig.1 Weighted decomposition of the interface

以上各量需满足如下条件

ωx和ωy由下式确定

式中,nx,i,j和nx,i,j分别是界面法向量n的x方向和y方向分量,具体计算过程见文献[14].

待权重确定后,通过网格边壁(i+1/2,j)的通量的计算公式如下

其中,式(33) 采用一维THINC 法计算,式(34) 采用SLIC 计算.网格边壁(i,j+1/2)的通量计算过程类似.

3 数值模型验证

3.1 界面处理

采用本研究建立的数值模型对Zalesak’s disk 旋转[21]和shearing flow 问题[22]进行模拟分析,以验证本数值模拟方法对大变形界面的追踪能力.其中,验证案例的计算域均为[0,1]×[0,1],均采用100×100,200×200,400×400 和800×800 四种密度的结构化网格进行模拟计算,并设置CFL=0.25.

(1)Zalesak’s disk 旋转问题

Zalesak’s disk 旋转问题常被用来检验界面处理方法的适用性和模拟精度.如图2 和图3 中黑色实线所示,Zalesak’s disk 是一个带有缺口的圆形,在外加速度场(u,v)=(y−0.5,0.5−x)的作用下绕点(0.5,0.5)匀速旋转,在t=2 时回到原处.Zalesak’s disk 在计算域上的初值定义为

图2 和图 3 分别给出对应于 200×200 和400×400 两种网格的模拟结果.图中a,b,c,d 分别为t=0.5,t=1.0,t=1.5,t=2.0 时刻的模拟结果,黑色实线为理论值.可见,在光滑界面处数值模拟结果与理论值基本吻合,几乎没有形变,在界面转折突变处随着界面运动而逐渐光滑.在网格加密后,这种界面形变得到了很好的抑制.由此可见,模拟得到的Zalesak’s disk 的形状在旋转过程中与理论值基本一致,表明本数值模型具有很好的界面追踪能力.为定量分析界面模拟精度,定义误差计算公式为

图2 Zalesak’s disk 数值模拟结果(网格密度:200×200)Fig.2 Simulation results of Zalesak’s disk with 200×200 mesh

图3 Zalesak’s disk 数值模拟结果(网格密度:400×400)Fig.3 Simulation results of Zalesak’s disk with 400×400 mesh

式中,N为网格数量,H(x)dx为数值模拟得到网格i内的体积分数,为网格i内的体积分数理论参考值.由于网格尺寸倍减,可采用下式计算界面模拟精度

式中,Eq表示qth对应网格的模拟误差.不同网格尺度下t=2.0 时刻的本模型模拟误差见表1.由表可知,采用THINC/WLIC 法本模型对Zalesak’s disk 旋转问题的界面计算精度达到了一阶精度.

表1 Zalesak’s disk 数值模拟误差Table 1 Numerical errors in the Zalesak’s disk test

(2)Shearing flow 问题

Shearing flow 问题是验证大变形情况下数值方法追踪界面能力的另一典型算例.通过模拟初值为圆心在(0.5,0.75)、半径为0.15 的圆,其内部体积分数为1,外部为0 条件下,在随时间变化的外加速度场驱动下圆盘的形状变化来检验数值模型的界面计算精度.其中,外加速度场用如下的流函数表示

分析流函数的时间变化可知,在0

图4 Shearing vortex 数值模拟(网格密度:200×200)Fig.4 Simulation of shearing vortex with 200×200 mesh

图5 Shearing vortex 数值模拟(网格密度:400×400)Fig.5 Simulation of shearing vortex with 400×400 mesh

表2 Shearing flow 数值模拟误差Table 2 Numerical errors in the shearing flow test

此外,对界面追踪模拟过程中计算域内总质量的变化进行了计算分析.设定M=∑αi,M0为初始时刻计算域内总质量.在不同网格尺度条件下,本模型得到的质量变化比|M−M0|/M0均保持在10−13以下,可以认为THINC/WLIC 界面追踪法基本上可以保持质量守恒.

3.2 线性液舱晃荡问题

根据线性势流理论可以得到线性液舱晃荡问题的解析解,该解可用来检验数值计算模型的收敛精度.如图6 所示,假定初始时刻水体和气体都处于静止状态,舱壁为不可渗可滑移边界.在t>0 后舱体以a的加速度水平向左运动,相当于舱内流体在受到竖直向下的重力g还受到水平向右、大小为a的质量力的作用.计算中,取几何尺寸和物理参数与Li等[23]的模拟案例一致.具体地,舱体宽度L=1.00 m,初始水深h1=1.00 m,初始气体高度h2=1.25 m,重力g=9.81 m/s2,水平加速度a=0.01g,水体密度ρ1=1000 kg/m3,气体密度ρ2=1 kg/m3.

根据线性势流理论,水面高程 η(x,t) 的解析解[24]为

其中

图6 线性液舱晃荡问题示意图Fig.6 Setting of the linear sloshing problem

分别选取64 × 144,128 × 288 和256 × 576 三种均匀网格,并设置CFL=0.1 进行模拟.需要指出的是,为了与理论解进行比较,在数值模拟时忽略了控制方程中黏性项的影响.图7 给出了本模型得到的舱体左右边壁处的水面高程(α=0.5) 时间变化的计算结果与解析解的比较.其中,64 × 144 解析度对应的网格尺寸约为0.016 m × 0.016 m,大于水面高程振幅,此时模型模拟结果已较准确地再现了水面高程的变化趋势,说明THINC/WLIC 法能较好地模拟微小的水面变形;256×576 解析度网格对应的网格尺寸约为0.004 m × 0.004 m,由图7 可见,模型模拟结果与解析解基本吻合.

图7 舱体左右边壁处水面高程变化Fig.7 Temporal profile of the interface elevations on the left and right walls

为分析模型的计算精度,增加512×1152 网格算例并将t=4 s 时刻的密度场与理论值进行对比分析.

定义

表3 线性液舱晃荡问题数值计算误差Table 3 Numerical errors in the linear sloshing test

本模型采用有限差分法对控制方程进行求解,而有限差分法在数值守恒性方面偏弱,为确定本模型在质量、动量和能量守恒方面的表现,进行了以下的分析.定义计算域总质量M=∑ρi∆x∆y,M0为初始时刻计算域内总质量,定义质量损失为(|M−M0|)/M0.计算域内总能量为

图8 为对上述4 个特征变量数值守恒性的评估结果.其中,黑色线条代表质量损失,绿色线条代表能量损失,蓝色线条代表水平方向动量损失,红色线条代表竖直方向动量损失.

由图8 可知,数值质量损失保持在10−10量级以下,可认为本案例中数值计算模型满足质量守恒;数值能量损失在10−4量级以下,且随着网格加密能量损失减小.由此可见,本耦合模型在流体的质量守恒方面具有良好表现,在流体的动量和能量守恒方面与采用有限体积法得到模拟结果相比略有不足,但采用较高精度的计算网格会提高模拟结果的能量和动量守恒性.

图8 质量、能量及动量的模拟误差Fig.8 Simulation errors on the conservation of mass,energy and momentum of fluid

3.3 溃坝问题

为验证本研究提出的WENO-THINC/WLIC水气耦合数值模型的模拟精度,针对干床面溃坝问题进行了数值模拟,并与Lobovsk´y 等[25]的物理模型试验结果进行了比较.如图9 所示,数值模拟的计算条件与Lobovsk´y 等的物理模型试验相同.水槽长为1.61 m,宽0.60 m.在距离水槽左边壁600 mm 处设置一竖直挡板,挡板左侧为水体,右侧为空气.水体的高度包括高水位H=600 mm 和低水位H=300 mm两种工况.物理模型试验中,挡板上接定滑轮,连接一重物,在重物的作用下挡板被抽离,左侧水体涌向右侧,形成溃坝流.参照Lobovsk´y 等的模型试验,数值模拟中设置了与物模试验相同的物理量监测点,包括右侧边壁上的两个压力监测点和水槽中部4 个水位监测点.其中,压力监测点P1 距水槽底部30 mm,P2 距水槽底部80 mm,水位监测线具体位置见图9.

图9 溃坝模拟试验示意图(单位:mm)Fig.9 Setting of the dam-breaking simulation test(unit:mm)

首先进行了溃坝过程中水流运动学特征的模拟结果与试验结果的对比,包括溃坝水流形态、前锋位移和沿程的水位变化等.其中,相关数值模拟案例均采用∆x=∆y=0.01 m 的结构化网格并设置CFL=0.1.图10 为不同时刻溃坝水体运动形态的模拟结果与试验结果.由图可见,在低水位H=300 mm情况下,数值模拟的水体运动形态与物理模型试验结果基本一致,溃坝水体前锋撞击右侧边壁后形成的水舌及空气空腔的形状和位置与物模试验结果也基本相同.

图10 不同时刻溃坝水体运动形态对比Fig.10 Comparison of dam-breaking water shapes at different times

图11 为溃坝水体前锋位移随时间变化曲线.将模拟得到的溃坝前锋位移与物理模型试验结果[25-27]进行定量对比可知,在溃坝水体前锋冲击右侧边壁前,数值模拟结果与试验结果吻合良好,说明本模型对溃坝前期水体形状的模拟精度较高.

图11 溃坝前锋线位移计算结果与试验结果比较Fig.11 Comparison of calculated and measured displacements of the wave front

与溃坝水体前锋位移相比,H1 ∼H4 四个监测点处的水位变化能给出更详细的溃坝水体沿程变化过程.考察从无量纲时间t∗=0 到t∗=10 内各监测点水面变化,可以看到溃坝波撞击边壁后形成的水舌和二次波的变化特征.图12 给出了H=300 mm 条件下模拟得到的4 个水位监测点处水位变化与前人物理模型试验研究结果[25]及其他数值模拟结果[28-29]的比较.由图可见,本模型得到的溃坝前锋首次到达各水位监测点的时间与物模试验一致;同时,二次波到达的时间和水位上升趋势与试验结果基本吻合.本模型和其他数值模型结果在溃坝波撞击边壁前均与物理模型试验结果较为吻合;不同点在于,本模型对溃坝水体冲击右侧壁面以后的后期水体运动模拟结果相比其他数值模型更为接近物理模型试验结果.Peltonen 等[28]使用的OpenFOAM 内置的标准VOF求解器interFOAM 所得结果在H3 测点后期出现水位严重偏大现象、H4 测点后期出现水位偏小现象,其改进的GFMFOAM (ghost fluid method FOAM) 方法在H2 测点后期出现水位严重偏小现象;Nguyen 和Yark[29]改进的VOF 法在H2 测点处存在非物理的触顶现象、H3 测点中后期出现水位明显偏大的问题.综合来看,溃坝问题在水体前锋撞击右侧边壁后,水气界面形变复杂,水气紊动更为剧烈,此阶段数值模拟较为困难.相比而言,本模型对溃坝水体整个运动过程的模拟结果与物理模型试验结果更为吻合.可见,本数值计算模型较好地描述了溃坝水体的运动学特征,能够较为准确地再现溃坝水体这种大变形自由表面流的运动过程.

图12 H1 ∼H4 处水位变化Fig.12 Water level elevations at the four locations H1 ∼H4

此外,溃坝水体对下游床面或建筑物的冲击力是溃坝问题研究中的另一个重要物理参数.通过考察高水位H=600 mm 条件下下游压强的计算结果,对本模型再现溃坝水体动力学特征的能力进行分析.其中,设定计算网格分别为∆x=∆y=10 mm和∆x=∆y=5 mm 两种情况,对应的网格数量为60×161 和120×322,并取CFL=0.1,将模拟得到的P1 和P2 两点的压强与试验结果[25]和他人数值结果[30]进行对比.如图13 和图14 所示,当溃坝波未到达右侧边壁时,监测点处的压强为大气压,相对压强为0;当溃坝前锋撞击边壁时,压强急剧增大,之后缓慢变小.沿右边壁,冲击压强峰值随高度上升而下降,即P1 点冲击压强峰值应大于P2 点峰值.本模型得到结果与物理模型试验相比,冲击压强峰值出现时间相同,峰值大小相当,下降趋势也大致相同.考虑到试验数据为典型压力变化,存在一定的波动空间,可认为本模型数值模拟压强结果与物理模型试验结果较为一致.与Sun 等[30]使用MPS (moving particle semiimplicit)法的结果相比,本模型结果的压力波动较小.此外,网格尺寸对压强影响较小,但采用低密度网格的模拟结果出现少许压力局部突变,网格加密后这种现象减少.总体而言,使用网格尺寸∆x=∆y=0.01 m得到的计算结果已与试验结果吻合较好,也就是说本模型通过较少的网格数即可得到较精确的解.

图13 P1 处压强历时曲线Fig.13 Pressure history at P1

图14 P2 处压强历时曲线Fig.14 Pressure history at P2

4 结论

本研究通过采用五阶WENO 格式求解对流项,三阶TVD RK 进行时间离散,THINC/WLIC 法追踪界面,建立了WENO-THINC/WLIC 耦合水气二相流数值计算模型.首先,选取Zalesak’s disk 旋转和shearing flow 问题,验证THINC/WLIC 法对外加速度场下大变形水气界面的追踪能力;其次,以线性液舱晃荡问题为例分析本研究提出的WENO-THINC/WLIC 耦合模型的计算误差以及对质量、动量和能量的数值守恒性;最后,以溃坝水流运动问题为例检验耦合模型模拟不可压缩水气二相流大变形运动问题的能力.主要结论如下:

(1) THINC/WLIC 法可有效捕捉复杂界面变形,模拟收敛精度可达一阶并保持物质质量守恒.

(2) 所建立的WENO-THINC/WLIC 耦合模型对线性液舱晃荡问题的模拟结果与解析解吻合度较高,对微小形变的捕捉较为准确;通过误差分析可知模型模拟精度为一阶,数值分析结果表明本模型保持质量守恒和能量守恒.

(3)耦合模型对溃坝水流运动的计算结果与物理模型试验结果吻合良好,并优于interFOAM、VOF 以及MPS 等数值模型的模拟结果,更够准确地捕捉了溃坝过程中溃坝水体的运动学特征和动力学特征,证明本模型对不可压缩水气二相流问题的数值模拟具有较好的适用性.

将改进的WENO 格式[23,31-36]和THINC 法代入本研究提出的WENO-THINC 耦合模型预期会得到更高精度的模拟结果.