基于Waha 双流体模型修正的UGS环空带压异常物理特征

2019-11-02闫相祯闫怡飞李鹏宇

中国石油大学学报（自然科学版） 2019年5期

闫相祯, 闫行, 闫怡飞, 李鹏宇

(1.中国石油大学(华东)储运与建筑工程学院, 山东青岛 266580;2.中国石油大学(华东)机电工程学院, 山东青岛 266580; 3.海洋石油工程(青岛)有限公司, 山东青岛 266580)

考虑地下储气库(UGS)注采井采取强采强注的作业方式,井筒内环空带压异常或井下层间窜流会严重影响气井的产量,降低采收率,对气田开发后续作业造成不利影响。因此,研究储气库井生产过程井筒环空压力异常形成机制是储气库井安全生产的关键[1-3]。环空压力异常(APB)的主要成因包括油管和套管泄露、固井时顶替效率低、水泥浆体系选择或配方设计不合理和固井后井筒水泥环封隔失效等。其中在国内外多个气田开发中,油套环空带压(A环空)异常井比例均大于56 %[4]。由于储气库气井处于持续循环注采作业过程,井筒内流体压力波动、井深结构等会产生交变应力,易诱发注采管柱振动导致气体泄漏,而完井设计时在A-环空会填充环空保护液,管柱泄漏导致油套环空中存在动态的混气液柱,同时流体与井筒间产生径向传热,考虑流体在环形约束空间中的低膨胀性,从而导致APB效应持续增长甚至过高。因此,对UGS超深井进行环空带压异常物理特征分析显得尤为必要。Kenneth等[5]基于平均温度和密度增量的考量,模拟分析带压异常和泄压过程中环空压力变化。Zhang等[6]开发了一种新的改进的水平管束两相流流量标准,并研究水平管束垂直向上流动中两相压降的影响。郑杰等[7]采用递推法循环迭代分析井筒温度变化对套管环空温度和压力变化的影响。这些研究成果主要以压力-体积-温度(PVT)为基础预测密闭环空压力的变化,用的大都是经验数据,通用性较差,且未充分考虑管柱泄漏过程A-环空气液两相流动特性,导致UGS环空带压异常的力学研究还不够充分。随着Waha 双流体模型的提出与发展,在工程应用分析方面取得了良好的发展成果。其优点是适用于单相和两相流中的各种快速瞬变,又能够分析得出冲击波并描述压力波,适用于湍流及众多非线性科学等领域。国内外已有的工程应用成果体现了目前Waha 双流体模型方法的发展方向[8-9]。笔者考虑UGS井筒结构和传热特点,建立基于Waha 双流体模型修正的UGS力学分析方法。通过有限元软件模拟研究管柱泄漏后A-环空内气液两相流动特性。采用过冷减压试验和数值计算对UGS力学分析方法进行对比研究,并分析储库产量、径向面积和气液比等因素对环空带压的影响。

1 UGS环空带压异常力学分析方法

1.1 UGS环空带压异常力学模型

基于弹性力学和管柱力学相关理论,建立UGS环空带压异常力学模型[10],如图1所示。考虑注采管柱生产实际工作条件,井身结构和井下工具分布较为复杂,包括表层套管、技术套管、生产套管、环空填充液、注采管柱、永久封隔器、尾管等。为了便于进行数值分析,假定井下管柱为长直杆件,且是各向同性、均匀连续的线弹性体;将井下管柱看成一个整体,平均分为n等分。针对管柱内不稳定注采气产生交变载荷,诱发注采管柱气体泄漏并流入油套环空内,气液两相流成为井筒内的主要流状,同时环空填充液具有良好的缓蚀除氧功能,溶解气含量较少,可认为油套环空中存在动态的混气液柱。图1给出了i段管柱泄漏后A-环空中流体流动概念图。

基于Oudeman-Bacretza模型的环空压力为

pAPB=pl+pg+pgl.

(1)

式中,pAPB为环空压力,MPa;pl为环空填充液压力,MPa;pg为环空填充氮气柱压力,MPa;pgl为管柱泄漏环空增长压力,MPa。

通常情况下环空填充液压力pl在pAPB中占据主导地位,约大于80%。一旦发生管柱泄漏,A-环空内两相流将存在许多不稳定性[11],其中Rayleigh-Taylor不稳定性和Kelvin-Helmholtz不稳定性加剧环空内气泡破裂,径向热传递通道内Ledinegg不稳定性诱发两相发生非周期性的漂移,伴生的密度波振荡均将对环空带压异常的力学分析产生重要影响。

1.2 基于Waha 双流体模型修正的环空带压异常数学模型

1.2.1 Waha 双流体模型

Waha 双流体模型是由欧盟在WAHA Loads项目中开发的热工水力分析模型,在涉及两相流的热工水力及安全分析领域有广泛应用。该数学模型考虑了管道的弹性、气液特性和壁面摩擦等,包含相分离和水平分层流动状态下各相之间的热量、质量和动量传递以及壁摩擦的相关性,能够分析具有变横截面的长管道系统中的两相流,范围从具有大的界面长度的流动(例如分层流动)到具有非常小的界面长度尺度的分散流动,大量试验和工程应用证明其对大部分系统瞬态行为的预测能够达到相应的精度要求。Waha 双流体模型平衡方程[8]为

(2)

(3)

AChi|vf|vf-AΓgvaog+AαfρfFgcos-AFf,wall,

(4)

AFg,wall,

(5)

AQvf-AΓg(hf+(vf-vaog)2/2)+AαfρfvfFgcos,

(6)

AQvg+AΓg(hg+(vg-vaog)2/2)+AαgρgvgFgcos.

(7)

式中,p为压力,MPa;α为蒸汽体积分数;v为速度,m/s;w为管道的轴向速度,m/s;e为管道的比总内能,J/kg;vaog为相间速度,m/s;Fwall为壁面摩擦力,N;下标f表示液相,g表示气相;Fgcos和FCVM分别为体积力和虚拟质量力，N;A为管道横截面积,m2;Γg为气相产生率,kg/m3;ρ为密度,kg/m3;pint为界面压力,Pa;Chi为内摩擦系数;t为时间,s;Qv为传热速率,W/m2;u为比焓，J/kg。

图1 UGS环空带压异常力学模型Fig.1 UGS annulus pressure anomaly mechanical model

界面传热过程中，相间能量交换为

Qvg=Hig(Ts-Tg)，

(8)

Qvf=Hif(Ts-Tf).

(9)

相间质量交换为

(10)

式中，Hig和Hif分别为气、液相传热系数,W/(m2·℃);Ts为壁面温度,℃;Tf和Tg分别为气、液体温度,℃。

双曲性是判断两相流模型预测精度和数值稳定性的重要标准之一,这些方程组当且仅当存在实数根时,作为时间和空间的函数具备双曲性[11]。对非常大的相对速度 (混流速度与声速比值)、热损失和流速测量误差等工况计算过程中,Waha双流体模型方程组通过雅可比迭代法求解，会出现复特征值根，从而不具备双曲性特征。为提高模型在储气库井筒安全分析应用中的可行性,基于储气库井筒结构特点对方程组添加的微分项(虚拟质量力、界面压力项和界面传热项)进行分析并修正,以改变方程组的结构形式,避免虚数特征根的出现。

1.2.2 双流体模型双曲特性改进

虚拟质量力在多相流模型中起重要作用[12],当气液两相作相对加速运动时,相间局部质量位移将产生虚拟质量力,是带动其周围的连续相流体运动所需要的附加力,例如两相水锤、环空气液相流动过程中,瞬间发生的压力波动引起流向和加速度变化。Waha 双流体模型中的虚拟质量力为

(11)

式中,CVM为虚拟质量系数,CVM∈[0,1];ρm为两相混合密度,kg/m3。

在实际工程应用中,若虚拟质量力不存在或相对较小,某些工况下Waha 双流体模型将产生复数特征值,并出现在与声速相当的气相滑脱速度vr。 Iztok 等[13]认为相对系数(气相滑脱速度vr与混合物音速vhy比值)ξ≥0.3时,Waha双流体模型中假定的虚拟质量力无法保证方程组的双曲性。只有ξ在数量级上足够小,即气、液相的速度相对值远小于混合物音速时,方程组才具备双曲函数特性。

同时,虚拟质量项[14]在多相流问题中反映了不同相间加速度而出现的相间界面压力、相间速度等。Waha双流体模型中,同时若气液两相压力与相间界面压力相等，仅在分层流动(虚拟质量项CVM=0)时,存在界面压力项,表示为

pint=Sαgαf(ρf-ρg)gD.

(12)

式中，S为流型分类系数。

界面压力项添加项目的是对虚拟质量力效果起到补充作用,但真正的界面压力非常复杂,其缺乏坚实的物理根据,同时若气液两相压力与相间界面压力相等，即产生压力模型时，会造成守恒方程中某些重要信息的缺失,因此原添加项不足以满足所有流动状态的需求,这也是导致模型不具备双曲性的重要因素之一。为此,引入适用于多种流型的相间界面压力项[14]为

(13)

基于虚拟质量项和压力项相间界面的分析,将Waha双流体模型方程组特征化:

kwaha[C1(MN)2+C2MN3+C3M3N+C4M2]+

kwaha(C5N2C6MN+C7N2+C8M2+ΔWint)=0.

(14)

其中

N=(vf-λi).

式中,ΔWint为界面压力项作用产生的独立项。

与虚拟质量力作用方式不同,界面压力项是以独立项形式对特征方程求解产生影响,几何意义为特征方程沿y轴平移ΔWint个单位,目的是减少数值发散提高方程组的双曲性。但式(8)直接求解比较困难,为提高计算效率,忽略两相的可压缩性,则可简化为

αfρgM2+αgρfN2+CvmρmMN-pint=0.

(15)

经判断满足以下条件时方程有实根:

(16)

对应虚拟质量力修正为

(17)

可满足方程组的双曲性要求。

1.2.3A-环空井筒传热

在储气库生产过程中井筒内流体与地层环境温差较大,导致各层套管环空密闭空间内流体温度和环空压力迅速增加[15]。特别是管柱发生泄漏, A-环空内流体受高温地层环境影响,扰乱井筒温度场稳态,在有限密闭空间内形成额外的径向稳态传热(液相-壁面热流Qf、壁面瞬时非稳态导热Qc和气相-壁面间热流Qg,忽略油管壁、油管-套管环空等阻力),环空带压随温差加大而急剧增加,井筒内径向传热效果如图2所示。

图2 井筒内径向传热效果Fig.2 Radial heat transfer effect in wellbore

基于RPI wall boiling模型[16]将原有相间传热项修正为

(18)

(19)

(20)

其中

式中,τ为传热周期,s;Vhi为i段气泡体积,m3;Nw为活性核心分布密度;f为气泡偏离频率,1/s;Ahi为所附气泡所覆盖的面积与受热总体积的比值;λ为扩散率;kl为导热系数;Dw为气泡直径,m。

基于上述修正,针对i段A-环空的动态混合气液的等压双流体模型平衡方程为

(21)

(22)

(23)

(24)

(25)

(26)

其中

式中,A为油套环空横截面面积,m2。

压力脉冲(p(x,t)根据线性关系产生管径的变化。

(27)

式中,d为管壁厚度,m;tw为管壁厚度,m;E为弹性模量,MPa;K取决于管道外壳的边界条件。

假设管壁可以自由地径向移动[7],将修正的守恒方程组改写为矩阵形式为

(28)

式中,ψ为求解变量;A、B和C为含未知量的矩阵。

设λ为式(28)的特征根,其对应特征方程为:det(B-λA)=0,可通过矩阵的行列式det进行求解。

(29)

其中

W=Cvmαfαg(αgρg+αfρf).

2 模型修正的验证

为验证储气库UGS环空带压异常力学分析方法的准确性和精度,采用Yasushi等[17]提出的温度梯度下垂直管道过冷减压试验对本文中所作的模型修正进行验证。图3为温度梯度下垂直管道过冷减压试验设置。

采用不锈钢管全长约3.2 m,直径53.5 mm。观测口在测试段的顶部部分,同时在测试段设置3个测点(A、B和C),并可根据试验需要改变测点位置。测试部分被电加热器加热,其沿试验段线性温度分布,施加的温度范围为室温至160 ℃。数据采集选用6个压力测量端子(PT)和13个测温端子(TC),并沿管轴设置在一个相同的平面上。试验段内部压力采用3种固态压力传感器对压力进行了测量,避免试验过程中压力的损伤和漂移。

图3 试验装置示意图Fig.3 Experimental device diagram

图4给出了原模型分析值、修正模型分析值与试验结果的对比。由计算结果可知,测点A、B和C处的原模型分析值分别从t=3.8、5.9和8.4 ms处无法获取准确计算值,这是由于试验过程中热作用不平衡,管内产生非线性温度分布,Ledinegg不稳定性加剧两相发生非周期性的漂移,但原模型在数值离散方面采用了稀疏的交错网格,使得数值扩散过大无法收敛。而修正模型分析值与试验结果比较接近,其中测点A处t=5 ms时, 修正模型分析值约为 0.41 MPa,试验值约为0.44 MPa, 差值为6.82%,当t=15 ms时, 修正模型分析值约为0.58 MPa,试验值约为0.62 MPa, 差值为6.45%;测点B处t=5 ms时, 修正模型分析值约为0.40 MPa,试验值约为 0.45 MPa, 差值为11.1%,当t=15 ms时, 修正模型分析值约为0.62 MPa,试验值约为0.59 MPa, 差值为-5.08%;测点C处t=5 ms时, 修正模型分析值约为0.53 MPa,试验值约为 0.56 MPa, 差值为5.36%,当t=15 ms时, 修正模型分析值为0.6 MPa,试验值为0.63 MPa, 差值为4.76%。对比发现测点A和B处差值较大,这是由于测试段A、B处受梯度温度加热,管内流体压力受温差影响较大,但在热不平衡作用下,Waha模型修正方程组仍获取较为准确结果,表明本文中修正方程组具备一定稳定性和精度。

图4 原模型分析值、修正模型分析值与试验结果对比Fig.4 Comparison of analysis value of original model, modified model analysis value and experimental results

3 UGS环空带压异常有限元分析

3.1 A-环空气液两相流动特性

为探究UGS管柱泄露后A-环空存在的动态的混气液柱的流动特性,分析UGS环空带压异常力学分析方法在UGS井筒结构环空带压分析的可行性,通过有限元软件Fluent对管柱泄漏后A-环空管柱内动态气液两相流特性进行研究[17-19]。模拟假定距管柱起点等效长度h=2 562.0 m处发生泄漏,管内压p=25 MPa,天然气温度为313.15 K,管壁摩擦系数取0.1,管内天然气甲烷含量为96.1%,摩尔质量为17.097 kg/kmol,沿井筒温度梯度为3.1 ℃/100 m,油管柱114.3 mm×6.88 mm,套管177.8 mm×10.26 mm,由于储气库井筒结构为细长型构件,因此采用双精度模式。湍流设置选用标准双k-e方程,并设定气体泄漏口呈射流状,模拟储库产量为11×104m3/d。

图5为管柱泄漏处环空内气液两相模拟结果。由图5可以看出:A-环空垂直结构内发生气体环状流动,A-环空内出现泡状流、段塞流、环状流与混状流等多种流动形态。当气体泄漏量高于局部流量时,气体开始流动上升,且井筒内以混状流为主,反之井筒内以泡状流或段塞流为主。泄漏处环状切片模拟图显示,随着管柱出气量增大,气体与液体之间会发生剧烈的紊流作用,逐渐在A-环空内聚积形成大气囊,呈现泡状流和段塞流的过渡流型。在A-环空中,不管是泡状流还是泡状-段塞流,气量越大,形成气团的尺寸也越大。

图5 管柱泄漏处环空内气液两相模拟结果Fig.5 Simulation results of gas-liquid two-phase in annular space of column leakage

图6为模拟周期内泄漏处两相中气体体积分数。由图6可知,泄漏点位置(A、B、C、D)处的气液比随模拟采气周期(30 d)的变化状态。由于APB积聚效应气液比分别为0<α<0.7、0.03<α<0.56、0.08<α<0.79和0.1<α<0.9。环空中的气液流动形态与常规油气井有明显不同的规律。A-环空垂直结构有利于气体流动,呈现泡状流、段塞流、环状流与混状流等多种流动形态,环空内气体占比均随着生产时间增长而不断变化。

图6 模拟周期内泄漏处两相中气体体积分数Fig.6 Gas volume fraction in two phases at leak in simulation cycle

图7为温度梯度下环空带压模拟值与原模型值、修正值对比。

图7 温度梯度下环空带压模拟值与原模型分析值、修正模型分析值对比Fig.7 Comparison of simulated value of annulus pressure, original model analysis value and modified model analysis value

由图7可知,当油管柱长度一定时,Waha模型修正值较原模型计算结果偏大,与有限元分析结果比较接近。当管长小于20 m时,梯度温差对A-环空带压异常作用不大。当管长约45 m时,原模型分析值为 0.84 MPa,有限元分析结果为0.79 MPa, 差值约在6.32%,而修正模型分析值为0.4 MPa;当管长约75 m时,原模型分析值为 1.60 MPa,有限元分析结果为1.58 MPa, 差值约在1.27 %,而修正模型分析值为1.04 MPa。随着管长增加,梯度温差增大对管柱的A-环空带压异常影响增强,因此当管柱足够长时,进行相关管柱安全设计时必须考虑井筒径向传热对环空带压异常的影响。

3.2 环空带压异常影响因素

根据笔者推导建立的环空带压异常力学分析方法，研究了生产过程储库产量、A-环空径向面积和气液比等对环空带压异常的影响。图8为不同储库产量对A-环空带压异常的影响曲线。

图8 管柱泄漏后不同储库产量对A-环空带压异常的影响曲线Fig.8 Effect of different reservoir yields on A-annuluspressure anomalies after string leakage

从图8可以看出,在稳定状态下A-环空中的气/液压力梯度几乎恒定,但发生井筒泄漏的储气库环空带压随着注入流量增加而不断升高,且增长趋势明显,当储库产量为7×104m3/d时,泄漏位置处环空带压异常值为21.4 MPa;储库产量为11×104m3/d时,泄漏点处环空带压异常值为24.1 MPa;储库产量为15×104m3/d时,泄漏点处环空带压异常值为 25.9 MPa。表明当储气库井生产过程发生管柱泄漏,储库产量增加对A-环空带压产生了不利影响,这是由于储库产量增加将导致泄漏到A-环空内的气体产生率增加,加剧井筒环空内APB效应。因此在储气库井生产过程中,应根据井口环空带压的实施监测数据,合理安排相应气井的储库产量,确保储气库井的安全生产运行。

图9为不同径向面积对A-环空带压异常的影响。从图9可以看出,随着径向面积增加,发生管柱泄漏的储气库环空带压不断升高,其中近管柱泄漏处的APB增长较为明显,表明径向面积过大将促进APB效应。当套管为139.70 mm×9.71 mm时,泄漏位置处APB异常值为17.0 MPa;储库产量为219.1 mm×10.16 mm时,泄漏处APB异常值为19.8 MPa,储库产量为273.10 mm×11.43 mm时,泄漏处APB异常值为 22.3 MPa。因此,在管柱设计及建设过程中,应避免因管材选用不当引起A-环空压力异常值过高,加剧APB效应甚至破坏井筒结构。

图9 管柱泄漏后不同径向面积对A-环空带压异常的影响Fig.9 Effect of different radial area on A-annulus pressure anomaly after string leakage

图10为管柱泄漏后不同气液比对A-环空带压异常的影响。从图10可以看出,随着气液比增加,发生管柱泄漏的储气库环空带压逐渐增长,其变化趋势与径向面积对环空带压影响趋势相似。当气液比为α<0.25时,泄漏处环空带压异常值为15.4 MPa;当气液比为0.5<α<0.75时,泄漏点处环空带压异常值为16.2 MPa;当气液比为0.8<α<1时,泄漏处环空带压异常值为 16.8 MPa。随着气液比增加,A-环空内气体组分上升,但环形约束空间的低膨胀性将导致泄漏位置处A-环空压力带压异常值逐渐趋于一定范围内稳定。

图10 管柱泄漏后不同气液比对A-环空带压异常影响Fig.10 Effect of different gas-liquid ratio on abnormality of A-ring pressure after string leakage

4 工程实例

4.1 现场数据

某衰竭油气藏储气库UGS-Q区块井群位于华南地域,由原主力生产井改建而成,均具有天然良好的密封构造。建成后用以满足目前陕京线、陕京二线、西气东输等长输管线对储气库季节及安全调峰气量的迫切需求[18-21],因此储气库群设计和安全运行要求比较高。主要针对储气库井UGS-Q5进行分析,其运行压力为10～40 MPa,下深为3 502 m,渗透率为2.46×10-3μm2,孔隙度为16.2%～22.3%,设计库容量为13.4×108m3,工作气量为(5～12)×108m3,储层厚度为10 m,注气时间为90 d。现场地质工况包括:①储层条件。储层以成层分布的溶蚀孔洞为主,岩性致密;②储库数据。原始地层压力为28.3 MPa,储层厚度为10 m,压力系数为0.16;③场地条件。地震峰值加速度为0.15 g,地震基本烈度为Ⅶ度。

UGS-Q储气库井群在完井后直接投产,先后对10口注采井实施注气。在投产期间6口井出现带压异常现象,其中UGS-Q5注气运行时,单井日注气量最低为20.2×104m3,最高可达54.1×104m3;采气运行时,单井日采气量最低为18.5×104m3,最高可达56.3×104m3。但在储气库井UGS-Q5运行过程中,现场井口检测数据显示其在回采第9 d发生带压异常,并呈不断增长趋势,后经井口泄压仍因环空带压过大,储库产量低于预期而被迫关井。图11为储气库井现场数据采集。

图11 UGS-Q5 现场数据采集Fig.11 UGS-Q5 field data collection

4.2 误差分析

根据现场监测数据分析储气库运行周期内UGS-Q5单井环空带压异常情况,并与现场记录结果进行相对误差分析,如图12所示。从数据对比可以看出,在注采周期运行过程中环空带压异常值不断增长,当生产时间为50 d时,计算结果与监测结果分别为10.93、11.34 MPa,误差为3.58%;生产时间为65 d时,计算结果与监测结果分别为12.56、13.23 MPa,误差为5.03%;生产时间为110 d时,计算结果与监测结果分别为12.88、13.5 MPa,误差为4.59%,表明本文模型计算结果计算精度较高。表1给出了UGS-Q储气库带压异常井群现场监测与模型分析结果对比,最大误差分别为4.62%、5.31%、4.28%和4.76%。通过以上误差分析,并结合环空带压异常形成机制研究结果,认为提及的UGS环空带压异常力学分析方法具有计算结果稳定性和工程适用性,可满足实际工程中UGS环空带压异常分析的计算需求。