一维气液两相漂移模型全隐式AUSMV算法研究*

徐朝阳1,孟英峰2,郭劲松1,李 皋2,邱全锋3

(1. 中国石油集团渤海钻探工程有限公司 定向井技术服务分公司,天津 300280;2. 西南石油大学 油气藏地质及开发工程国家重点实验室,成都 610500;3. 西华师范大学 数学与信息学院,四川 南充 637000)

摘要:气液两相漂移模型显式AUSMV(advection upstream splitting method combined with flux vector splitting method)算法的时间步长受限于CFL(Courant-Friedrichs-Lewy)条件,为了提高计算效率,建立了一种全隐式AUSMV算法求解气液两相漂移模型.采用AUSM格式结合FVS(flux vector splitting)格式构造连续方程和运动方程的对流项数值通量,AUSM格式构造压力项数值通量.离散控制方程是非线性方程组,采用六阶Newton(牛顿)法结合数值Jacobi矩阵求解.计算经典算例Zuber-Findlay激波管问题和复杂漂移关系变质量流动问题,结果分析表明:全隐式AUSMV算法,色散效应小,无数值震荡,计算精度高.在压力波波速高的条件下,可以显著提高计算效率,耗散效应小.

关 键 词:气液两相漂移模型; AUSMV格式; 全隐式算法; 六阶Newton法; 数值Jacobi矩阵

引 言

气液两相漂移模型(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]

1 气液两相DFM

DFM是一种TFM的退化模型,两者的主要区别在于TFM有2个运动方程(气相运动方程和液相运动方程),考虑了气相与液相之间的相互作用,而DFM有1个运动方程,考虑气液混相流,通过引入经验性的气液速度之间的漂移关系,补充了缺少的两相之间相互作用关系.相比TFM,DFM的双曲性更强,形式更为简单,计算量小.

1.1 控制方程组

一维等温气液两相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.

1.2 辅助方程

为了封闭方程,需要引入以下辅助方程:

气液相流速间本构关系为

vg=C0vm+vt,

(4)

其中

vm=αlvl+αgvg,

(5)

式中,C0为气体分布系数,无量纲;vm为气液混相平均流速,m/s;vt为气体漂移速度,m/s.根据不同工程问题,C0vt可采用相应的经验公式或半经验公式求得.

液相、气相体积系数的归一化关系为

αl+αg=1.

(6)

根据不同工程问题,液相和气相的状态方程不同,为了验证算法,采用简化线性热力学关系状态方程.

液相状态方程:

(7)

式中,ρl,STP为标况下液相密度,1×103 kg/m3pSTP为标况压力,1×105 Pa;cl为液相压力波波速,1×103 m/s.

气相状态方程:

(8)

式中,cg为气相压力波波速,316 m/s.

2 全隐式AUSMV格式及数值算法

2.1 全隐式AUSMV格式

AUSM类格式基于控制方程的守恒形式离散,DFM矢量形式:

(9)

其中,Uw分别为守恒变量和原始变量,

(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分别为ii+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)

2.2 隐式离散方程数值计算方法

守恒变量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是尾数二进制表示的最高位数.

3 数 值 实 验

数值算例为Evje和Fjelde提出的经典算例:Zuber-Findlay激波管问题和复杂漂移关系质量传输问题[8,21],可以代表井筒中的典型流动状态.隐式计算结果与可靠参考值、显式计算结果对比,因为DFM的原始变量为pαgvl,所以只对比压力剖面、气相体积系数剖面和液相流速剖面.数值实验的目的在于检验全隐式AUSMV算法的计算精度、稳定性和效率.为了便于分析,下文中采用缩写形式EM(explicit method)和IM(implicit method)分别表示显式AUSMV算法和全隐式AUSMV算法.

3.1 Zuber-Findlay激波管问题

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

3.2 复杂漂移关系变质量流动问题

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

4 结 论

本文建立了一种一维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.

Research on the Implicit AUSMV Algorithm for the1D Gas-Liquid Two-Phase Drift Flux Model

XU Chaoyang1, MENG Yingfeng2, GUO Jinsong1, LI Gao2, QIU Quanfeng3

(1.Directional Drilling Service Company,CNPC Bohai Drilling Engineering Company Limited,Tianjin 300280,P.R.China;2.State Key Laboratory of Oil &Gas Reservoir Geology and Exploitation,Southwest Petroleum University,Chengdu 610500,P.R.China;3.School of Mathematics &Information,China West Normal University,Nanchong,Sichuan 637000,P.R.China)

Abstract: The time step of the explicit AUSMV (advection upstream splitting method combined with flux vector splitting) algorithm is limited by the CFL (Courant-Friedrichs-Lewy) conditions. To improve computational efficiency, an implicit AUSMV algorithm was proposed for the gas-liquid two-phase drift flux model. The numerical flux of convective terms in the continuity equations and motion equations was set up with the AUSM scheme plus the FVS (flux vector splitting) scheme, while the numerical flux of pressure terms in the motion equations was built with the AUSM scheme. The nonlinear dynamical discrete governing equation system was solved numerically with the 6th-order Newtonian method and the numerical Jacobian matrix. The classical test examples were simulated, which involved the Zuber-Findlay shock tube problem and the variable mass flow problem with complex slip relation. The numerical results show that, the implicit AUSMV algorithm has small dispersion effects, no numerical oscillation and high computational accuracy. Under the condition of high pressure wave velocity, the algorithm has superior calculation efficiency with low dissipation effects.

Key words: gas-liquid two-phase drift flux model; AUSMV; implicit algorithm; 6th-order Newtonian method; numerical Jacobian matrix

Foundation item: The National Science and Technology Major Project of China(2016ZX05021-004);The National Natural Science Foundation of China(51674217)

中图分类号:O359.1; O241.82

文献标志码:A

DOI: 10.21656/1000-0887.390110

文章编号:1000-0887(2019)04-0386-12

应用数学和力学编委会,ISSN 1000-0887

*收稿日期:2018-04-08; 修订日期:2018-05-21

基金项目: 国家科技重大专项(2016ZX05021-004);国家自然科学基金(51674217)

作者简介:

徐朝阳(1985—),男,工程师,博士(通讯作者. E-mail: 04011xzy@sina.com);

孟英峰(1954—),男,教授,博士,博士生导师(E-mail: cwctmyf@vip.sina.com).

引用本文/Cite this paper:

徐朝阳, 孟英峰, 郭劲松, 李皋, 邱全锋. 一维气液两相漂移模型全隐式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.