王建军， 曾祥俊， 贾善坡， 隋晓凤， 张赟新， 贺海军

(1.中国石油集团工程材料研究院有限公司/石油管材及装备材料服役行为与结构安全国家重点实验室， 西安 710077; 2.东北石油大学地球科学学院， 大庆 163318; 3.东北石油大学非常规油气研究院， 大庆 163318; 4.东北石油大学土木建筑工程学院， 大庆 163318; 5. 中国石油新疆油田分公司呼图壁储气库作业区， 昌吉 831200; 6.中国石油大庆油田有限责任公司采油工程研究院， 大庆 163453)

地下储气库注采井井筒是保证天然气注入、采出的安全通道，注采井受力与变形的影响因素不仅包括重力、浮力、管内外流体压力、底部封隔器处约束方式、操作顺序等常规外界因素，还需要考虑注采井井筒温度、压力非线性分布、注采交变荷载、温度及腐蚀对井筒结构强度影响。由于一般储气库井设计寿命在30年以上，运行时间较长，因此要确保井筒在设计年限内安全使用，必须对引起井筒结构破坏的因素进行深入研究，否则会出现注采井泄露、水泥环内部产生微裂隙或胶结面密封失效等问题[1]。从国内已建成的大港储气库群、京58储气库群等储气库井筒泄露问题来看，主要是由于因交变荷载变化过大，超过接头额定承载能力，密封性不足而引起的泄漏，注采气温度交替变化是产生过大载荷的主要因素。所有的地下储气库注采井的井筒都将面临注采过程中温度的周期性变化，准确简便地预测井筒内温度分布，对于储气库的安全运行至关重要。

目前，获得井筒温度、压力等数据的方法是通过下入井下压力计和温度计对井筒温度和压力的沿程分布进行动态监测，但是该方法成本较高，并且下放的仪器有使用条件的限制，对于一些高温高压气井来说并不适用。此外，就是通过数值模拟的方法和建模解析的方法对井筒温度和压力进行预测。国内外学者开展了大量关于井筒温度和压力预测的理论研究，迄今运用最为广泛的是Ramey[2]建立的井筒热传导理论模型，分析了单相不可压缩流体和理想气体在注入井内温度分布情况，模型中假设井筒内热量传递为稳态过程，在地层的热量传递为非稳态过程，将井筒内油管、环空、套管和水泥环等结构组成的多个热阻用总传热系数来表示，热量在地层的传递比较复杂，引入了无因次时间函数用来表征地层的非稳态传热过程，极大简化的计算过程，并且证明了该模型的适用性。之后的研究者认为模型中井筒的总传热系数不能当作常数考虑，对环空部分流体的对流传热系数和辐射传热系数进行了讨论，并且根据现场测温资料回归出了模型中松弛距离的统计公式，进一步简化了井筒温度的计算过程[3-6]。为了使井筒温度模型适用范围更广，Hasan等[7]将两相流引入的温度模型之中，分析了Joule-Thomson效应对井筒温度的影响，对无因次时间函数进行的分段处理，提高了热量在地层中非稳态传热的计算精度，对于井筒压力分布，运用最为广泛的是Culldener-Smith采用的分段迭代的方法求解。因此，对温度和压力预测的方法主要分为三类：第一类，利用Ramey经典的传热模型，该模型将压力剥离出来，并且没有考虑流体的摩擦和动能改变项对井筒温度分布的影响，认为井筒内气体为理想气体，与实际生产情况相差较大。但该模型解析表达式计算简单，可以简便地计算出井筒任意深度处流体温度。第二类，通过分段平均的方式对井筒压力分布进行迭代计算，这种计算方法需要对压差进行分段，按长度增量迭代或者对长度分段，按压力增量进行迭代，计算出划分管段的平均压力和平均温度。第三类，利用质量守恒定律、能量守恒定律和动力定理，将温度和压力进行耦合计算[8]。部分学者利用以上方法对井筒注入或采出某单一过程的流体温度和压力进行预测，对于地下储气库而言，在生产周期内，包括了注采过程、采气过程和关井阶段，对符合储气库注采规律的温度和压力预测的研究较少。

根据地下储气库实际生产运行规律，对温度、压力预测的方法进行优化，推导出符合地下储气库注采全过程的温度、压力预测方法。本文利用显式四阶Runge-Kutta数值求解方法[9]，将温度、压力、密度和流速等流动参数进行耦合计算，以井底或者井口处气体流动参数为初值，利用单步递推的计算方式输出注/采气过程中温度、压力等气体流动参数的沿程分布数据。为了系统的分析注采井井筒的温度分布，现根据地下储气库实际生产情况，把注采井井筒运行过程分为开井阶段和关井阶段，在关井阶段，以关井前一刻气体流动参数作为计算关井后温度变化的初值，以得到关井后油管内气体温度的沿程分布情况。

1 储气库注采井井筒温度场预测模型

储气库在生产周期内主要分为开井阶段和关井阶段，在开井阶段注/采气过程中井筒内气体会持续不断的注入或者采出。从而假设出气体在井筒内作一维稳定流动，即油管内气体的流动参数不随时间变化，只在流动方向上发生变化，油管内任意截面处流动参数均匀。

为了更为简便和精确地预测各阶段油管内气体温度，以垂直井为例，做出如下假设：①油管内气体到水泥环外缘的传热过程为稳态传热，热量在地层间为非稳态传热；②热量以导热形式从油管传递到地层，径向传热量远大于纵向传热量；③油管、套管和水泥环为同心圆；④关井后，油管内气体由运动变为静止，忽略续流的影响。计算模型如图1所示。

根据图1所示的井身结构，对取出的单元控制体进行径向传热、开/关井阶段的传热分析。在径向传热分析中，认为热量在油管、环空、套管和水泥环这部分做稳定传热，将这部分考虑为由单一热阻构成的传热介质[10]，引入总传热系数Uto来表征热量在稳态传热区的传递；热量在地层间做非稳态传热，将地层的热阻考虑为与时间有关的函数关系，引入无因次时间函数来表征热阻随时间的变化，简化了热量在地层之间的非稳态传热。对于开井阶段，根据储气库实际生产时的规律，结合流体力学中质量守恒定律、能力守恒定律和动量定理，认为天然气在注采井油管内做一维稳定流动，考虑气体的压缩性，通过Runge-Kutta单步递推的数值求解方法，得到了密度、流速、温度和压力的沿程分布情况。在关井阶段，气体速度的稳定比温度的稳定要快得多，认为关井后气体停止运动，以关井前划分井段的温度和密度为计算初值，得到关井后温度的沿程变化情况。

图1 传热模型示意图Fig.1 Schematic diagram of the heat transfer model

1.1 径向传热

油管内气体到地层的传热过程为气体与油管内表面的对流传热、油管导热、环空部分对流传热和辐射传热、套管导热、水泥环到地层的导热，径向传热模型如图2所示。

根据假设条件①，油管内气体到水泥环外缘认为由单一热阻构成的稳态传热，传递的热量[11]为

dQ1=2πrtiUto(Tf-Th)dz

(1)

式(1)中：Uto为总传热系数，W/(m·℃)；dQ1为井筒单位时间内dz长度上的热损失，W。

水泥环外缘到地层的传热过程为非稳态传热，引入无因次时间函数[11]，则传递的热量为

rti为油管内半径，m；rto为油管外半径，m；rci为套管内半径，m；rco为套管外半径，m；rh为水泥环半径，m；Tf为油管内气体温度，℃；Tan为环空保护液温度，℃；Tci为套管温度，℃；Th为水泥环外缘温度，℃；dz为取出的单元控制体长度，m图2 径向传热示意图Fig.2 Radial heat transfer schematic

(2)

式(2)中：ke为地层导热系数，W/(m·℃)；Te为地层温度，℃；dQ2为地层单位时间内dz长度上的热损失，W；f(td)为无因次时间函数，无量纲。

无因次时间函数的公式为

(3)

(4)

式中：α为地层热扩散系数，m2/s;t为生产时间，s。

由于油管内气体传递的水泥环外缘的热量等于水泥环外缘向地层传递的热量，则水泥环外缘温度[12]为

(5)

套管内表面温度[12]为

(6)

式(6)中：kcem为水泥环导热系数，W/(m·℃)。

由于油管的导热性能良好，使得油管内气体温度与油管外表面温度差异很小，则环空部分温度[12]近似

(7)

总传热系数为

(8)

式(8)中：ktub为油管导热系数，W/(m·℃)；hr为环空内辐射传热系数，W/(s·m2·℃)；hc为环空内对流传热系数，W/(s·m2·℃)；kcas为套管导热系数，W/( m·℃)。

1.2 开井阶段

以垂直井为例，气体在油管内作一维稳定流动，规定坐标轴正向与气体流动方向一致，取dz长度的管段为单元控制体，如图3所示。

ρ为气体真实密度，kg/m3；ν为气体流速，m/s；A为油管截面积，m2；g为9.8 m/s2；Wf为摩阻力，N 图3 开井阶段单元控制体示意图Fig.3 Schematic diagram of the unit control body in the well opening stage

根据流体在油管内做一维流动假设，流动参数只在流动方向上发生变化，油管内任意截面处流体参数均匀，由质量守恒定律得到单位时间内流入单元体的质量等于流出单元体的质量，则微分形式的连续方程[13]为

(9)

根据动量定理，则压力梯度[10]为

(10)

式(10)中：f为管壁摩擦系数，无量纲。

根据能量守恒定律[10]：

(11)

式(11)中：h为气体比焓，J/kg；q为单位质量气体热损失，J/kg。

井筒内为稳态传热，井筒与地层之间为非稳态传热[13]，由dQ1=dQ2，可得

(12)

式(12)中：ω为质量流量，kg/s。

流体比焓为

(13)

式(13)中：CJ为Joule-Thomson系数，K/Pa。

由于油管管径变化很小，故Joule-Thomson系数甚小可以忽略，根据式(11)～式(13)，则温度梯度[14]为

(14)

引入松弛距离L为

(15)

气体密度根据实际气体状态方程得到在深度方向上微分形式[13]为

(16)

式(16)中：M为气体摩尔质量，kg/mol；R为气体常数，R=8.314 J/(mol·K)；Z为压缩因子，无量纲；T为热力学温度，K。

根据式(9)、式(10)、式(14)和式(16)，分别得出注气过程[式(17)]和采气过程[式(18)]中温度梯度、压力梯度、密度梯度、速度梯度的微分方程组[8-9，13]为

(17)

(18)

式(18)中：Cp为天然气定压比热。

Cp=1 243+3.14T+7.931×104T2-

6.881×107T3

(19)

油管内任意温度和压力条件下，气体密度和气体流速[14]为

(20)

(21)

式中：γg为气体相对密度；P为气体压力，MPa；qsc为标准状态下气体流量，m3/d。

天然气与油管内壁的摩擦系数f的计算公式为

(22)

式(22)中：ε为管壁的绝对粗糙度，mm；Re为雷诺数，无量纲。

将密度梯度、速度梯度、压力梯度、温度梯度分别记为Fi(i=1，2，3，4)初始条件记为井口或井底处密度、速度、压力、温度，步长记为x。

(23)

利用经典的四阶Runge-Kutta求解此方程组，计算过程[9，13]如下：

(24)

(25)

以井口或者井底的流动参数为初值，通过显式Runge-Kutta的数值求解方法依次计算，从而输出任一深度处气体的密度、流速、温度、压力等参数。

1.3 关井阶段

关井之后，油管内充满一定压力和温度的气体，气体由运动状态变为静止[15]，取长度为dz的单元控制体，如图4所示。

图4 关井阶段单元控制体示意图Fig.4 Schematic diagram of unit control body in shut-in stage

热传导损失于地层中的热量等于内能的变化量加上井筒结构中吸收的热量再减去单位长度控制体动能、势能和焓变的改变量[16]，则能量平衡方程为

(26)

式(26)中：m为油管内气体的质量，kg；m′为井筒结构的质量；E为气体内能，J/kg；E′为井筒系统的内能，J/kg。

关井后，油管内气体由运动状态变为静止，则质量流量为0，式(26)变为

(27)

井筒向地层吸收或损失的热量为，则有

Q=Cp(Te-Tf)L′

(28)

式(28)中：L′为质量流量为0时的松弛距离[7]。

(29)

由于油管、套管、环空和水泥环吸收的能量和井筒内能量变化量存在比例关系为

CtmE=m′E′

(30)

(31)

式中：Ct为热存储系数[7]，在关井压力恢复时取2，在压力降低时取3。

气体温度关于时间和空间上的微分方程[7]为

(32)

在关井的初始时刻，，此时油管内气体温度Tf0等于关井前一刻油管内气体温度Tf，则关井阶段气体温度关于时间的方程[7，15-16]为

(33)

(34)

1.4 计算流程

注采周期内油管内气体流动参数计算流程[17-36]如图5所示。

图5 流动参数计算流程图Fig.5 Flow parameter calculation flowchart

2 模型验证

以ZS 8井注入CO2为例，CO2气体受温度和压力影响，在不同井深处相态会发生改变，由液态转变为超临界状态，达到超临界态的CO2气体具有液体的密度但又具有气体的特性，由于储气库注采的天然气中甲烷含量相对较高，因此用此案例来验证耦合微分方程组计算的准确性。ZS 8井井筒和地层参数、初始注入条件[17]如表1和表2所示。

表1 ZS 8井井筒和地层参数表Table 1 Wellbore and formation parameters of well ZS 8

表2 CO2初始条件Table 2 CO2 initial conditions

根据图6(a)和图6(b)所示，CO2气体温度和压力随着深度的增加而逐渐升高，计算值与实测值趋势一致，计算值略大于实测值，压力计算值与实测值最大相对误差为4.97%，温度计算值与实测值最大相对误差为3.18%，满足工程实际要求。

3 储气库注采井井筒温度影响因素分析

以西部某储气库X注采井为例，预测该井在注气过程和采气过程油管内气体的温度、压力、流速和密度分布情况，并对不同注气温度、注气压力、注/采气时间、注/采气量等敏感性因素进行单因素分析，并且通过只考虑单元控制体径向传热的假设条件下，井筒内环空温度、套管温度、水泥环温度进行预测，以及关井阶段，油管内气体恢复到地层温度所需时间进行分析，X注采直井井身结构如图7所示。

X注采直井井筒和地层参数取值如表3所示。

图6 ZS 8井气体流动参数对比Fig.6 Comparison of gas flow parameters in well ZS 8

表3 X井井筒和地层参数表Table 3 Wellbore and formation parameters of well X

图7 X注采直井井身结构示意图Fig.7 Schematic diagram of the structure of the X injection-production vertical well

3.1 注气过程

本阶段以注气时间为30 d，注气压力为30 MPa，注气温度为30 ℃时，注气量为30×104、50×104、70×104、90×104、110×104m3/d时油管内气体温度、压力、密度、流速的分布情况如图8所示。

如图8(a)所示，油管内气体温度随深度增加而逐渐增大，这是由于随着深度的增加，地层温度逐渐高于注气温度，此时注入的气体从地层中吸收热量。在接近井口0～1 300 m范围内，注气量大时的气体温度高于注气量小时的气体温度，其原因主要是注气量越大，油管内气体流速越快，如图8(c)所示，油管内气体与外界的热交换就越小，注气量越小，气体流速越慢，与外界的热交换就越大，在此范围内油管内气体温度高于地层温度，气体与外界热交换量相对较小。图8(b)反映了随着深度的增加，油管内气体压力逐渐增大，并且随着注气量的增大，油管内气体压力增大。其原因可根据图8(d)可知，注气量越大，油管内气体密度越大，因而压力增加，由密度方程[式(20)]可知，气体密度又是关于气体温度和压力有关的参数，因而相互影响。

注气时间为30 d，注气压力为30 MPa，注气量为50×104m3/d时，注气温度为20、30、40、50 ℃时油管内气体流动参数分布情况如图9所示。

图8 注气阶段不同注气量下气体流动参数分布Fig.8 Distribution of gas flow parameters under different gas injection volumes during gas injection stage

图9(a)反映了不同注气温度下，油管内气体温度沿程分布情况，在注气温度为20 ℃时，油管内气体温度呈现出随深度增加而逐渐增大的趋势，随着注气温度的增大，油管内气体温度先出现了先降低后逐渐上升的趋势。这是由于在接近井口处，随着注气温度的增加，油管内气体温度与地层的温差变大，导致散失的热量增加。随着深度的增加，不同注气温度之间的差异逐渐减小。图9(b)反映了随着注气温度的增加，油管内气体压力减小，可由图9(c)可知，不同注入温度条件下，油管内气体密度不同，因而导致油管内气体压力出现变化。

图9 注气阶段不同注气温度下气体流动参数分布Fig.9 Distribution of gas flow parameters at different gas injection temperatures during gas injection stage

注气时间为30 d，注气温度为30 ℃，注气量时50×104m3/d注气压力分为20、30、40 MPa情况下，油管内流体温度和压力变化情况如图10所示。

图10 注气阶段不同注气压力下气体流动参数分布Fig.10 Distribution of gas flow parameters under different gas injection pressures during gas injection stage

图10(a)反映出通过改变注气时的压力得到油管内气体温度沿程分布的情况，因此改变注气压力对于油管内气体沿程温度分布影响很小。图10(b)反映了随着注气压力的增加，油管内气体压力增大，可由图10(c)可知，压力增加，注入气体的密度会显著增加，随着深度的增加，受到温度增加的影响，气体密度增加趋势减缓。

注气压力为30 MPa，注气温度为30 ℃，注气量为50×104m3/d时，注气时间为10、30、50、70 d油管内气体温度分布如图11所示。

图11(a)反映出随着注气时间的增加，油管内气体温度降低不明显，因此改变注气量对油管内气体温度变化有显著影响。则通过改变注气量来影响井底温度的效果优于改变注气时间。

根据本文的模型假设，只考虑单元控制体径向上的热量传递，得到了注气压力30 MPa，注气温度30，注气量为50×104m3/d时，注气时间持续30 d，注气阶段井筒内套管、环空、水泥环外缘温度预测。

图11(b)反映出，在接近井口处，由于油管内气体温度高于地层温度，则在0～800 m的范围内油管内气体温度依次大于环空、套管和水泥环温度，气体温度到套管温度之间的差异较小；随着深度的增加，地层温度逐渐增大，在注气30 d以及稳定注气量的情况下，油管内气体温度到水泥环外缘的温度依次升高，水泥环外缘温度与线性变化的地层温度之间存在差异，是由于在此注气条件下，井筒的温度对周围地层存在一定的热影响范围。

图11 不同注气时间及井身结构温度分布Fig.11 Different gas injection time and temperature distribution of well structure

3.2 采气过程

图12(a)反映出，随着采气量的增大，井口处气体温度显著增加，越接近于井口处，不同采气量之间温度差异越大，主要原因是越接近于井口油管内气体与地层的温差越大，因而油管内气体的热量损失就会越大，并且采气量的大小决定油管内气体流速，由图12(c)所示，随着采气量的增大，油管内气体流速增加，导致温度差异明显。图12(b)反映出随着采气量的增大，油管内气体压力逐渐降低，根据图12(d)可知油管内气体密度显著下降，在采气量较小时，气体密度变化幅度较小，随着采气量的增加，密度变化较为明显，加上加速度压降损失和摩擦力损失会增大，因而导致油管内气体压力降低。

采气压力为33.96 MPa，地层温度为94.2 ℃，采气时间为30 d，采气量为30×104、50×104、70×104、90×104、110×104m3/d时油管内气体流动参数情况如图13所示。

图13(a)反映出采气量为50×104m3/d时，采气时间为10、30、50、70 d时油管内气体温度沿程分布清情况，油管内气体温度随时间增加而增大，但增加的幅度较小，因而采气时间的长短对油管内气体温度影响较小。

在采气过程中，井底温度比较高，因而油管内高温气体向井口流动的过程中向周围地层传递热量，因此，图14(b)给出了采气时间为30 d，采气量为50×104m3/d时油管、套管、环空液、水泥环和地层的温度沿程分布情况。根据图13(b)可知，在此采气工况下，井筒温度对周围地层的影响趋于稳定，油管温度>环空液温度>套管温度>水泥环温度>地层温度，并且越接近于井口处，油管温度与套管内表面径向温度差异越明显。在油管内高温气体向周围传热的过程中，环空液温度升高，并受热膨胀，套管轴向变形增加，会存在上抬井口的风险，环空液温度发生变化，在温度效应下又会产生环空压力，因而通过本文的预测模型，计算出环空部分温度，为计算环空压力提供可靠的理论数据。

3.3 关井阶段

计算了在采气30 d，油管内气体温度、压力、密度的分布情况，关井之后，以关井前气体流动参数作为初始条件，此时油管内充满一定温度和压力的气体，气体由运动状态变为静止，只存在径向上的热传递，则关井后，油管内气体恢复到地层温度的传热过程如图14所示。

图12 采气阶段不同采气量下气体流动参数分布Fig.12 Distribution of gas flow parameters under different gas production volumes in gas production stage

图13 不同采气时间及井身结构温度分布Fig.13 Temperature distribution of different gas production time and well structure

图14 关井阶段井筒分度分布Fig.14 Wellbore indexing distribution during shut-in stage

根据图14(a)可知，随着关井时间的持续增加，油管内气体温度会逐渐恢复至地层温度，图14(b)反映出，在关井初期，油管内气体温度下降速率较大，并且越接近于井口，气体温度下降速率越明显，主要是油管内气体温度与地层温差较大所致，随着关井时间的持续增加，温度变化速率逐渐放缓。随着关井时间的持续增加，油管气体温度逐渐恢复至地层温度。

4 实例分析

由于本文的计算模型需要以注气时井口温度或者采气时井底温度为初值对气体流动参数的沿程分布情况进行求解，现场只能监测到井口处开井状态下的温度，因此引用西部某储气库X注采井井筒采气过程井口监测数据，计算出从2019年12月10日—2020年1月31日开井期间井口处气体温度，将预测值与实际监测的井口温度进行比对，实测值与预测值如图15所示，随着采气时间的持续增加，井口监测值呈现出上升趋势，最大监测值为50 ℃，最小值为47 ℃。预测值高于实测值，计算温度最高在56 ℃左右，分析误差的来源可能为测量环境、模型中优化的井身结构与实际井身结构存在的差异、模型中的地层温度为线性变化与实际存在差异、实际生产中产量不稳定、地层导热率和地层热扩散系数等因素引起误差。通过对比可知，在该阶段内，预测值与井口实际监测值吻合程度较高。说明本文所建立的温度预测方法，适合于储气库生产运行中的温度预测。

图15 井口监测温度与预测值对比Fig.15 Comparison of wellhead monitoring temperature and predicted value

5 结论

通过建立开井阶段温度梯度、压力梯度、流速梯度和密度梯度耦合的微分方程组，利用Runge-Kutta数值求解方法，计算出不同深度处气体的流动参数，并且考虑单元体径向传热的假设条件下，得到了油管、环空液、套管、水泥环的温度沿程分布情况，结合关井前气体的流动参数，假设关井后气体停止运动，得到了关井后不同深度处气体恢复至地层温度情况。综上所述，得到以下结论。

(1)在注气过程中，随着注气量的增大，井底处气体温度明显降低，气体压力会随着注气量的增大而增加。注气压力和注气量不变的情况下，改变井口处注气温度，发现接近于井底处，温度之间的差异会逐渐缩小。井底处气体温度受注气量的影响大于受注气温度的影响。

(2)随着注气时间的持续增加，井底处气体温度变化较小。只改变注气压力的大小，对油管内气体沿程温度的差异较小，但油管内气体压力会随着注气压力的增大而增大。

(3)采气阶段油管内气体温度会随着采气量的增加而增大，但油管内气体压力随着采气量的增加而减小，因此在生产过程中要控制井口处压力变化，可以通过改变采气量的大小来控制压力变化。

(4)在采气量不变的条件下，随着采气时间的增加，井口处气体温度逐渐增加，但温度增加幅度较小。

(5)关井阶段，油管内气体温度会随着关井时间的持续增加逐渐恢复至地层温度，关井初期越接近井口处温度变化速率越明显，变化幅度越大。