气液两相漂移模型(drift flux model,DFM)最初由Zuber和Findlay提出[1].该模型忽略气、液间的相互作用,基于经验参数构建气液相流速间的本构关系,具有结构简单的特点.当两相流动强耦合时,该模型可以避免双流体模型(two fluid model,TFM)计算中容易中出现的双曲性退化和非守恒项需要数学修正等难题[2-4].而且,气液相流速间本构关系中的经验参数容易通过实验获取,可以方便应用于各类工程问题.所以,DFM在油气井井筒气液两相流动研究中应用广泛,其数值求解方法也是重要的研究课题.
AUSM系列格式(AUSMV、AUSMD、AUSMDV、AUSM+、AUSMPW、AUSM+-up等)不依赖Riemann求解器或Jacobi矩阵构建数值通量,最初由Liou和Steffen提出,具有结构简单、高计算效率、高间断分辨率、高稳健性的特点[5-6].AUSM系列格式最初用于计算单相可压缩流体,Evje等成功将AUSMV、AUSMD、AUSMDV、AUSM+格式扩展至计算气液两相流动模型(TFM和DFM),研究发现计算TFM强激波问题时稳健性弱于近似Riemann求解器,计算DFM的瞬变流和缓变流问题时精度高,结果基本与近似Riemann求解器相同[3,7-9].Niu等先期研究AUSMDV、AUSM+-up格式计算7方程TFM,随后发展出结合混合原始变量Riemann求解器的AUSMD格式,提高了计算精度,但是该方法需要采用Reo平均法处理变量,推导过程比较复杂[10-11].Kitamura和Nonomura成功将AUSM+-up通量函数耦合Godunov精确Riemann求解器的方法拓展至耦合HLLC(Harten-Lax-Van Leer-contact) Riemann求解器,应用于计算TFM,研究表明新方法计算效率高于原耦合方法,计算精度高于原始AUSM+-up格式[12].
上述研究中AUSM类格式都是被用于建立显式算法,虽然显式格式的离散方程结构简单,易于求解,但是时间步长取值受CFL条件限制,限制了计算效率.Evje等提出半隐式WIMF(weakly implicit mixture flux)格式计算TFM和DFM,时间步长小幅增大,满足CFL数等于1,但是该方法引入了1个压力演化方程(偏微分方程),使模型复杂化[13].随后,他们提出了全隐式的SIMF(strongly implicit mixture flux)格式,可以进一步突破CFL条件的限制,增大时间步长,但是该方法仍然需压力演化方程[14].Colonia等采用隐式AUSM+和AUSM+-up格式计算高Mach(马赫)数可压缩流体,计算高效、结果可靠,但是在离散方程组求解中采用解析Jacobi矩阵,推导过程极为复杂[15].Onur和Eyi研究了Euler方程Newton法求解中数值Jacobi矩阵和解析Jacobi矩阵的影响,表明解析Jacobi矩阵的推导过程复杂,尤其是高阶空间精度格式离散方程的解析Jacobi矩阵推导将更为复杂,而采用数值方法可以大幅简化Jacobi矩阵的推导过程,同时计算精度没有明显降低[16].Zeng等研究了隐式和显式AUSM+、 AUSMDV、 AUSM+-up格式计算TFM,离散方程组求解过程采用数值Jacobi矩阵,简化了模型推导,研究中比较了交错网格和同位网格离散方程的影响,结果表明隐式算法可以大幅提高计算效率,隐式算法的间断捕捉精度与显式方法相同,相比同位网格离散方程,交错网格计算结果稳定性微弱提高[17-18].
上述研究主要关注隐式AUSM类格式求解TFM,隐式AUSM类格式求解DFM的研究较少,国内讨论也较少.本文是文献[19]的进一步研究,提出了一种基于AUSMV格式的DFM全隐式算法.该算法采用全隐式AUSMV格式构建DFM的数值通量,同位网格离散控制方程组,采用六阶Newton法结合数值Jacobi矩阵求解离散方程组[17-18,20].通过经典数值算例验证了全隐式AUSMV算法的计算精度、稳定性和计算效率[8,21].
DFM是一种TFM的退化模型,两者的主要区别在于TFM有2个运动方程(气相运动方程和液相运动方程),考虑了气相与液相之间的相互作用,而DFM有1个运动方程,考虑气液混相流,通过引入经验性的气液速度之间的漂移关系,补充了缺少的两相之间相互作用关系.相比TFM,DFM的双曲性更强,形式更为简单,计算量小.
一维等温气液两相DFM由各相连续方程和混相运动方程3个方程组成:
液相连续方程
(1)
气相连续方程
(2)
混相运动方程
(3)
式中,ρg为气相密度,kg/m3;ρl为液相密度,kg/m3;αg为气相体积系数,无量纲;αl为液相体积系数,无量纲;vg为气相流速,m/s;vl为液相流速,m/s;t为时间,s;x为空间长度,m;p为压力,Pa;qg为重力分量,Pa/m;qf为流动阻力分量,Pa/m.
为了封闭方程,需要引入以下辅助方程:
气液相流速间本构关系为
vg=C0vm+vt,
(4)
其中
vm=αlvl+αgvg,
(5)
式中,C0为气体分布系数,无量纲;vm为气液混相平均流速,m/s;vt为气体漂移速度,m/s.根据不同工程问题,C0和vt可采用相应的经验公式或半经验公式求得.
液相、气相体积系数的归一化关系为
αl+αg=1.
(6)
根据不同工程问题,液相和气相的状态方程不同,为了验证算法,采用简化线性热力学关系状态方程.
液相状态方程:
(7)
式中,ρl,STP为标况下液相密度,1×103 kg/m3;pSTP为标况压力,1×105 Pa;cl为液相压力波波速,1×103 m/s.
气相状态方程:
(8)
式中,cg为气相压力波波速,316 m/s.
AUSM类格式基于控制方程的守恒形式离散,DFM矢量形式:
(9)
其中,U,w分别为守恒变量和原始变量,

(10)
分别处理控制方程各相流体的对流项和混相压力项[5-6]:

(11)
式中,Fl(w)表示液相对流项,Fg(w)表示气相对流项,Fp(w)表示气液混相的压力项.
那么,离散单元界面的数值通量为
(12)
结合FVS方法处理连续方程和运动方程的对流项[8],压力项采用AUSM方法处理[5],同位网格离散,可以表示如下:
连续方程对流项数值通量
(13)
运动方程对流项数值通量
(14)
运动方程压力项数值通量
(15)
式(13)、(14)中ph表示气液各相流体,L和R分别表示离散单元界面的左侧与右侧,j表示时间步,i表示空间步,见图1.一阶空间精度格式,L和R分别为i和i+1.全隐式AUSMV格式中数值通量所用原始变量均为j+1时间步.

图1离散单元示意图
Fig. 1 Schematic diagram of discrete cells
对式(13)、(14)中速度分裂函数V+,V-和式(15)中压力分裂函数P+,P-,Evje等推荐采用下面的表达式:
速度分裂函数

(16)
压力分裂函数

(17)
其中
(18)
式(16)的系数满足
χL=αR,χR=αL,χ∈[0,1].
(19)
混相压力波波速为

(20)
那么,根据上述AUSMV格式构建控制方程组的数值通量,DFM的全隐式矢量离散形式为
(21)
根据上述方法,也可以采用显式AUSMV格式得到离散方程,但其数值通量中原始变量均为j时间步,见式(22),用于算例结果对比:
(22)
守恒变量U、数值通量F和源项S均是关于原始变量w的函数,式(21)离散方程组是关于原始变量w的非线性方程组.R是离散方程组的残差矢量,又可以表示为
(23)
Newton-Raphson(牛顿-拉夫逊)法是计算非线性方程组的经典方法,具有局域二阶收敛性.采用改进的六阶Newton法求解离散方程组,提高收敛效率[20].每一时间步的计算结果为下一时间步的初始预测值,直至达到计算时间为止.
六阶Newton法:
y(wk)=wk-[J(wk)]-1R(wk),
(24)
z(wk)=y(wk)-τ[J(wk)]-1R(y(wk)),
(25)
wk+1=z(wk)-τ[J(wk)]-1R(z(wk)),
(26)
其中
τ=2I-[J(wk)]-1J(y(wk)),
(27)
式中,k表示每个时间步中的迭代步.当wk+1-wk的差的2-范数满足如下条件时,则计算收敛,wk+1就为当前时间步的解:
‖wk+1-wk‖<εtol,
(28)
式中,εtol为绝对误差限,无量纲.
I为单位矩阵.J(w)为Jacobi矩阵,是R(w)关于原始变量w的一阶偏微分,即
(29)
AUSMV格式数值通量的完整表达式复杂,因此解析偏微分推导过程很复杂,通过数值微分可以大幅简化推导工作.R(w)关于原始变量w的数值一阶偏微分可以表示为
(30)
其中,en是含有3×N元素的单位向量,N是计算域划分的单元数.e中第n个元素为1,其他都为0.Ri指第i个残差方程,i为从1至3×N.ε是摄动值,取正值为前向差分,取负值为后向差分.ε的取值大小非常关键,直接关系到数值Jacobi矩阵的计算精度.Onur和Eyi提出ε的取值为[16]
(31)
ΣM是计算机精度,为计算机能够识别的最小数,简单估算方法为
(32)
其中,m是尾数二进制表示的最高位数.
数值算例为Evje和Fjelde提出的经典算例:Zuber-Findlay激波管问题和复杂漂移关系质量传输问题[8,21],可以代表井筒中的典型流动状态.隐式计算结果与可靠参考值、显式计算结果对比,因为DFM的原始变量为p,αg,vl,所以只对比压力剖面、气相体积系数剖面和液相流速剖面.数值实验的目的在于检验全隐式AUSMV算法的计算精度、稳定性和效率.为了便于分析,下文中采用缩写形式EM(explicit method)和IM(implicit method)分别表示显式AUSMV算法和全隐式AUSMV算法.
1) 问题描述
水平管长100 m,管径0.1 m,初始管内压力和各相流速在50 m处存在间断,初始条件见表1,属于无黏可压缩气液两相瞬变流.
表1 Zuber-Findlay激波管问题初始条件
Table 1 The Zuber-Findlay shock tube problem’s initial conditions

p/Paαgvl/(m/s)vg/(m/s)left 50 m80 4500.5510.37012.659right 50 m24 2820.550.5611.181
算例中,Zuber-Findlay漂移关系参数C0=1.07,vt=0.216 m/s.计算中忽略摩阻,计算结束时间为1 s,εtol=1×10-6.参考值为文献[21]采用MC限制器的显式Roe算法的计算结果.
2) 结果分析
空间离散Δx=1 m,最大压力波波速cmax=25 m/s,EM计算中CFL数等于0.25,0.5,即Δt=0.01,0.02 s;IM计算中CFL数等于0.25,0.5,1,5,即Δt=0.01,0.02,0.04,0.2 s,计算结果如图2所示(全隐式AUSMV算法100个单元与显式AUSMV算法100个单元,结合MC限制器的显式Roe算法3 200个单元[21]的对比).
Δt=0.01,0.02,0.04 s时,在压力剖面、气相体积系数剖面和液相流速剖面中,在左、右间断的两侧之外的非间断处,EM和IM的计算结果与参考值非常吻合,计算精度高.压力剖面的左、右间断,Δt=0.01,0.02 s时,EM分别需要6个网格捕捉间断,计算精度高;Δt=0.01,0.02 s时,IM分别需要8个网格捕捉间断,计算精度较高;Δt=0.04 s时,IM分别需要16个网格捕捉间断,间断捕捉能力较好;Δt=0.2 s时,IM不能完全捕捉间断,间断捕捉效果较差.气相体积系数剖面和液相流速剖面的间断处,EM和IM的间断捕捉能力与压力剖面的结果相对应,具有一致的现象.
EM采用Δt=0.02 s时,压力剖面、气相体积系数剖面和液相流速剖面的左侧间断处,有轻微的数值震荡,压力剖面比较明显(线框内),Δt=0.01 s时数值震荡减弱.不同时间步长条件下,IM计算结果均无数值震荡存在,色散效应极微弱,稳定性更强.相同CFL条件下,在间断处,IM的计算精度略微低于EM.随着时间步长的增大,在间断处,IM的计算精度降低,间断处被拉宽.当Δt=0.2 s,即CFL数等于5时,气相体积系数剖面中两个间断面处明显拉宽(线框内),尤其是液相流速剖面的小幅间断完全抹平(线框内),表明IM的耗散效应随着时间步长的增大而增大.所以,计算最大压力波波速很小的瞬变流问题时,若对计算精度要求比较高,IM计算时,CFL数不宜过大.

(a) 压力剖面 (b) 气相体积系数剖面
(a) Pressure profiles (b) Gas volume fraction profiles

(c) 液相流速剖面 (d) 迭代计算残差值记录
(c) Liquid velocity profiles (d) Residual value records during iterative calculation
图2Zuber-Findlay激波管问题计算结果剖面
Fig. 2 Comparison of simulation result profiles for the Zuber-Findlay shock tube problem
迭代计算残差值记录表明IM迭代计算收敛速度快,收敛残差值小.表2所示,当CFL数小于1时,IM总迭代步数是EM方法的2倍;提高CFL数,CFL数等于1,IM的总迭代步数是EM在CFL数等于0.25的约1/2,与EM方法在CFL数等于0.25时基本相同.IM方法CFL数等于5的总迭代步数是EM在CFL数等于0.25的约1/5,是EM在CFL数等于0.5的9/25.IM方法的迭代总步数减少明显.
CPU时间方面,CFL数等于0.25,IM的CPU时间约是EM的28倍;CFL数等于0.5,IM的CPU时间约是EM的19倍.因此,当CFL数小于1时,IM的计算效率远低于EM方法.IM在CFL数等于1的CPU时间是EM在CFL数等于0.25的约5倍,是EM在CFL数等于0.5的约7.6倍;IM在CFL数等于5的CPU时间是EM在CFL数等于0.25的约1.6倍,是EM在CFL数等于0.5的约2.5倍.Zuber-Findlay激波管问题中,CFL数等于5,IM的计算效率仍低于EM.
上述分析表明,Zuber-Findlay激波管问题中,相比EM,IM并未明显提高计算效率.原因在于本算例中,气相体积系数大,最大压力波波速cmax=25 m/s,计算时间短,仅为1 s,EM中时间步长取值较大,而IM时间步长的增大幅度有限.虽然,IM的总迭代步数减少,但是每一迭代步的计算量大(计算数值Jacobi矩阵),小幅增大时间步长取值,不能体现IM的优势.
表2 Zuber-Findlay激波管问题EM和IM计算效率对比
Table 2 Comparison of computational efficiency between EM and IM for the Zuber-Findlay shock tube problem

CFLΔt/sschemeiteration step total mTCPU time TCPU/s0.250.01EM1000.163 033IM2004.680 5250.250.02EM500.106 836IM1002.075 72610.04IM510.812 38550.2IM180.265 952
1) 问题描述
水平管长1 000 m,管径0.1 m,管内初始为纯液相(初始αg=1×10-7).出口压力保持为1×105 Pa,入口质量流量均为线性变化,见表3.
表3 各相入口质量流量变化
Table 3 Change of inlet mass flow for each phase

time T/sliquid mass flow rate Qlm/(kg/s)gas mass flow rate Qgm/(kg/s)00010120.0850120.0870120175120
漂移速度是关于气体体积系数的函数,为非线性速度滑移关系,气体分布系数仍是常数.
(33)
该算例属于黏性可压缩气液两相变质量流动问题,液、气相黏度分别为μl=5×10-2 Pa·s,μg=5×10-6 Pa·s,运动方程中重力分量和动阻力分量分别表示为
qg=(ρlαl+ρgαg)gsin(θ),
(34)
(35)
式中,d为管道内径,m.
计算结束时间为175 s,参考值为文献[21]采用MC限制器的显式Roe算法的计算结果.
2) 结果分析
空间离散Δx=10 m,最大压力波波速cmax=1 000 m/s,EM计算中CFL数等于0.25,0.5,即Δt=0.002 5,0.005 s;IM计算中CFL数等于10,100,即Δt=0.1,1 s,计算结果如图3所示(全隐式AUSMV算法 100个单元与显式AUSMV算法100个单元,结合MC限制器的显式Roe算法3 200个单元[21]的对比).
压力剖面中,没有间断,EM在CFL数等于0.25,0.5的计算结果,以及IM在CFL数等于10的计算结果,都与参考值非常吻合,计算精度高;CFL数等于100时,IM在550 m之后的计算结果与参考值非常吻合,在550 m之前的计算结果与参考值吻合较好,略小于其他计算结果,计算精度较高.气相体积系数剖面中,参考值显示,管道中部有强物理间断,即中部气相体积系数较大,而管道两侧气相体积系数基本为零,为纯液相.在550 m之前,EM和IM的计算结果几乎一致,与参考值吻合良好,没有数值震荡;在550 m之后,IM在CFL数等于100的计算结果略大于其他计算结果,所有计算结果与参考值吻合良好.气相体积系数剖面的左、右间断,EM和IM均分别需要18个网格捕捉间断,间断捕捉精度良好.液相流速剖面中,在550 m之前,EM和IM的计算结果都与参考值吻合良好;在550 m之后EM和IM的计算结果都稍大于参考值.本算例中,EM和IM的计算结果都与参考值吻合良好,尤其是压力剖面,计算精度高.
计算结果表明:黏性流动变质量流动问题,IM计算精度高,即使CFL数等于100时仍可以准确捕捉间断,耗散效应随着时间步长的增大而增长微弱.

(a) 压力剖面 (b) 气相体积系数剖面
(a) Pressure profiles (b) Gas volume fraction profiles

(c) 液相流速剖面 (d) 迭代计算残差值记录
(c) Liquid velocity profiles (d) Residual value records during iterative calculation
图3复杂漂移关系变质量流动问题计算结果剖面
Fig. 3 Comparison of simulation result profiles for the variable mass flow problem in complex drift relation
迭代计算残差值记录表明IM迭代计算收敛速度快,收敛残差值小.表4显示,IM在CFL数等于10时,总迭代步数是EM在CFL数等于0.25的约1/40,是EM在CFL数等于0.5的约1/20;IM在CFL数等于100时,总迭代步数是EM在CFL数等于0.25的约1/250,是EM在CFL数等于0.5的约1/125.显然,IM可以大幅降低总迭代计算步.CPU时间方面,CFL数等于100时,IM是EM在CFL数等于0.25的约1/10,是EM在CFL数等于0.5的约1/5,大幅提高了计算效率.
上述分析表明,复杂漂移关系变质量流动问题,通过大幅增大时间步长,IM的计算效率显著提高,且CFL数等于100时计算精度仍然较高.原因在于本算例中,最大压力波波速cmax=1 000 m/s,EM的时间步长取值很小,IM在CFL数等于100的时间步长是EM在CFL数等于0.25的400倍,体现了IM通过大幅增大时间步长而提高计算效率的优势.
表4 复杂漂移关系变质量流动问题EM和IM计算效率对比
Table 4 Comparison of computational efficiency between EM and IM for the variable mass flow problem in complex drift relation

CFLΔt/sschemeiteration step total mTCPU time TCPU/s0.250.002 5EM70 00062.669 6290.50.005EM35 00031.131 551100.1IM1 77268.708 9351001IM2916.087 007
本文建立了一种一维DFM的一阶全隐式AUSMV算法,该算法采用全隐式AUSMV格式构建控制方程组的数值通量,同位网格离散,并采用六阶Newton法结合数值Jacobi矩阵求解离散方程组.通过经典算例验证,该全隐式AUSMV算法具有以下特点:
1) 色散效应小,没有数值振荡;
2) 非黏性激波流动中,CFL数较小时,计算精度较高,随着CFL数的增大,耗散效应明显增大,若时间步长取值过大,间断会被拉宽或抹平,降低间断捕捉精度;
3) 黏性流体变质量流量流动中,耗散效应随着CFL数的增大而微弱增大,时间步长取值较大,间断没有被拉宽或抹平,能够精确地描述间断,计算精度高;
4) 压力波波速较大的问题中,通过大幅增大时间步长,可以显著增大计算效率.
[1] ZUBER N, FINDLAY J A. Average volumetric concentration in two-phase flow systems[J].Journal of Heat Transfer, 1965,87(4): 453-468.
[2] BOURÉ J A. Wave phenomena and one-dimensional two-phase flow models[J].Multiphase Science and Technology, 1997,9(1): 63-107.
[3] EVJE S, FLÅTTEN T. Hybrid flux-splitting schemes for a common two-fluid model[J].Journal of Computational Physics, 2003,192(1): 175-210.
[4] EVJE S, FLÅTTEN T. On thewave structure of two-phase flow models[J].SIAM Journal on Applied Mathematics, 2006,67(2): 487-511.
[5] LIOU M S, STEFFEN C J. A new flux splitting scheme[J].Journal of Computational Physics, 1993,107(1): 23-29.
[6] KITAMURA K, SHIMA E. Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes[J].Journal of Computational Physics, 2013,245: 62-83.
[7] CHANG C H, LIOU M S. A robust and accurate approach to computing compressible multiphase flow: stratified flow model and AUSM+-up scheme[J].Journal of Computational Physics, 2008,227(1): 840-873.
[8] EVJE S, FJELDE K K. Hybrid flux-splitting schemes for a two-phase flow model[J].Journal of Computational Physics, 2002,175(2): 674-701.
[9] EVJE S, FJELDE K K. On a rough AUSM scheme for a one-dimensional two-phase model[J].Computers &Fluids, 2003,32(10): 1497-1530.
[10] NIU Y Y, LIN Y C, CHANG C H. A further work on multi-phase two-fluid approach for compressible multi-phase flows[J].International Journal for Numerical Method in Fluids, 2008,58(8): 879-896.
[11] NIU Y Y. Computations of two-fluid models based on a simple and robust hybrid primitive variable Riemann solver with AUSMD[J].Journal of Computational Physics, 2016,308: 389-410.
[12] KITAMURA K, NONOMURA T. Simple and robust HLLC extensions of two-fluid AUSM for multiphase flow computations[J].Computers &Fluids, 2014,100: 321-335.
[13] EVJE S, FLÅTTEN T. Weaklyimplicit numerical schemes for a two-fluid model[J].SIAM Journal on Scientific Computing, 2005,26(5): 1449-1484.
[14] EVJE S, FLÅTTEN T. CFL-violating numerical schemes for a two-fluid model[J].Journal of Scientific Computing, 2006,29(1): 83-114.
[15] COLONIA S, STEIJL R, BARAKOS G N. Implicit implementation of the AUSM+ and AUSM+-up schemes[J].International Journal for Numerical Method in Fluids, 2014,75(10): 687-712.
[16] ONUR O, EYI S. Effects of the Jacobian evaluation on Newton’s solution of the Euler equations[J].International Journal for Numerical Method in Fluids, 2005,49(2): 211-231.
[17] ZENG Q L, AYDERMIR N U, LIEN F S, et al. Comparison of implicit and explicit AUSM-family schemes for compressible multiphase flows[J].International Journal for Numerical Method in Fluids, 2015,77(1): 43-61.
[18] ZENG Q L, AYDERMIR N U, LIEN F S, et al. Extension of staggered-grid-based AUSM-family schemes for use in nuclear safety analysis codes[J].International Journal of Multiphase Flow, 2017,93:17-32.
[19] 徐朝阳, 孟英峰, 魏纳, 等. 一维气液两相漂移模型的AUSMV算法研究[J]. 应用数学和力学, 2014,35(12): 1373-1382.(XU Chaoyang, MENG Yingfeng, WEI Na, et al. Research on the AUSMV scheme for 1D gas liquid two phase flow drift flux models[J].Applied Mathematics and Mechanics, 2014,35(12): 1373-1382.(in Chinese))
[20] MADHU K. Sixth order Newton-type method for solving system of nonlinear equations and its applications[J].Applied Mathematics E:Notes, 2017,17: 221-230.
[21] FLÅTTEN T, MUNKEJORD S T. The approximate Riemann solver of Roe applied to a drift-flux two-phase flow model[J].ESAIM:Mathematical Modelling and Numerical Analysis, 2006,40(4): 735-764.
徐朝阳(1985—),男,工程师,博士(通讯作者. E-mail: 04011xzy@sina.com);
孟英峰(1954—),男,教授,博士,博士生导师(E-mail: cwctmyf@vip.sina.com).
徐朝阳, 孟英峰, 郭劲松, 李皋, 邱全锋. 一维气液两相漂移模型全隐式AUSMV算法研究[J]. 应用数学和力学, 2019,40(4): 386-397.
XU Chaoyang, MENG Yingfeng, GUO Jinsong, LI Gao, QIU Quanfeng. Research on the implicit AUSMV algorithm for the 1D gas-liquid two-phase drift flux model[J].Applied Mathematics and Mechanics, 2019,40(4): 386-397.