基于水化反应动力学的深水固井井筒温度与压力耦合预测模型

王雪瑞，孙宝江,2，刘书杰，李中，刘正礼，王志远,2，李昊,2，高永海,2

（1.中国石油大学（华东）石油工程学院，山东青岛 266580；2.海洋水下设备试验与检测技术国家工程实验室，山东青岛266580；3.中海油研究总院，北京 100027；4.中海石油（中国）有限公司湛江分公司，广东湛江 524057；5.中海石油（中国）有限公司深圳分公司，广东深圳 518067）

0 引言

近年来，国际海洋油气勘探开发逐渐从浅水区域向深水区域发展。尽管深水海域具有广阔的勘探开发前景，但复杂的深水海洋环境给油气勘探开发带来巨大的挑战[1]，其中典型挑战之一即为深水固井技术。固井作业通过在套管与地层之间建立完整的水泥环封隔，达到封隔油、气、水层的目的，进而避免层间窜流事故的发生。准确预测井筒温度场和压力场是固井水泥浆体系以及固井工艺流程设计的基础，尤其在固井候凝期间，剧烈的水泥浆水化反应与温度、压力之间的相互作用直接影响固井质量与安全，具体表现在：水泥浆水化反应速率受温度、压力影响明显，导致水泥浆水化度、水泥浆胶凝强度等水化性能呈现随时空变化的特性[2-3]；水泥浆水化热的释放改变水泥温度场，诱发水泥环内热应力场，因此形成的水泥环内微裂缝为气窜提供了运移通道[4-5]；水泥浆水化过程中存在“失重现象”，其内部压力随时间不断下降，为气窜提供动力[6]。此外，特殊的深水温压环境导致水泥浆水化反应、温度、压力之间呈现出更为复杂的相互作用。因此，建立准确的固井井筒温压场预测模型是深水安全高效固井的前提。

当前工业界最常用的预测固井井筒温度场的方法是美国石油学会（API）提出的半经验方法[7]，但Wedelich等[8]通过与现场测量结果对比指出该方法总是过高地预测井筒温度。因此，国际上不同学者针对固井过程建立了各种热传导物理模型[9-11]；国内刘洋等[12]、王清顺等[13]针对不同固井工艺建立了相应的温度场预测模型。固井井筒压力场方面，Carter等[14]首先发现了候凝期间压力降低的现象；在此基础上，国内外诸多学者通过实验发表了相关的理论解释以及压力预测公式[15-17]。

目前不同学者针对固井井筒温度压力场的研究仍将温度与压力变化视为两个独立的过程，且未充分考虑到水泥浆水化反应的影响。鉴于此，本文考虑了水泥浆水化反应、温度、压力三者的相互作用，结合深水特殊环境，建立了基于水化反应动力学的深水固井候凝井筒温度压力耦合预测模型，并利用此模型分析固井候凝过程中井筒温度、压力、水化度随时空的演化规律，进行了敏感性分析，为固井施工提供了安全指导。

1 深水固井候凝阶段井筒温度压力耦合模型

深水固井过程中水泥浆被注替到套管与地层之间的环形空间中（见图1）。固井进入候凝阶段时，水泥浆开始发生剧烈的水化反应，在该过程中，水化反应、温度、压力三者之间存在复杂的相互作用：①水泥浆水化反应过程中释放出大量的水化热，这些热量不断改变水泥浆内部的温度。②随着水化反应的深入，水泥浆自身产生一定的胶凝强度，分担了部分原水泥浆静液柱压力，进而导致水泥浆内孔隙压力的下降。③在地温梯度以及水泥浆初始静液柱压力的作用下，井筒内不同位置水泥浆处于不同的温度压力环境，因此，水泥浆的水化反应速率随空间变化，不同深度处水化反应进度也不同。④某一位置处水泥浆的温度、压力受水化反应影响随时间不断变化；同时，瞬态的温度、压力也改变水泥浆的水化反应速率，不同时间点水泥浆水化度、水泥浆胶凝强度等水化性能也不同。

图1 深水固井工艺示意图

1.1 深水固井候凝阶段井筒瞬态温度场计算模型

固井水泥热传导过程中，候凝期间井筒内水泥浆温度变化的影响因素主要包括3个方面：水泥浆水化热导致的水泥浆温度升高；水泥浆与套管内钻井液之间的热量交换；水泥浆与周围地层之间的热量交换。根据热量交换过程，需要分别建立环空内水泥浆热传导模型、套管内流体热传导模型以及地层热传导模型以得到水泥浆温度。考虑整个固井系统以井筒为中心轴对称，建立如图1所示的柱坐标系。

1.1.1 套管内流体热传导模型

如图1所示，取长度为zΔ的单元格为研究对象。根据能量守恒原理，单元格内热量的变化等于流入的热量减去流出的热量：

（1）式左侧为单元格内流体的能量变化量；右侧第1项表示与环空的热交换量，第2项为与下方相邻单元格的热交换量，第3项代表与上方相邻单元格的热交换量。将上式写成差分形式可以得到套管内流体热传导模型：

1.1.2 地层热传导模型

热量在井筒周边地层的传导过程是一个轴对称的扩散过程。根据建立的柱坐标系，地层内热传导模型可以根据轴坐标系下的热扩散方程获得[18-19]：

1.1.3 环空内水泥浆热传导模型

水泥浆与套管以及地层之间的热交换可以根据（2）式和（3）式进行定量描述；水泥浆水化热的释放量根据水化反应动力学进行计算。取长度为Δz的单元格为研究对象，根据能量守恒原理，环空内水泥浆热传导模型可以表达为：

（4）式左侧表示环空单元格内内能的变化量；右侧第1项表示水泥浆水化释放的热量，第2项表示与套管的热交换量，第3项表示与地层的热交换量，第4、5项分别为与下方和上方相邻单元格的热交换量。将上式写成差分形式可以得到环空内水泥浆热传导模型：

1.2 深水固井候凝阶段井筒瞬态压力场计算模型

室内以及现场实验证实[6,15]，水泥环内部的孔隙压力随着水化反应的深入逐渐下降，该现象被称为水泥失重。针对水泥失重现象，水泥浆胶凝强度理论[20-23]认为：在固井的循环阶段，水泥浆呈现出液体的性质并能够传递全部静液柱压力；一旦候凝期水泥浆水化反应开始，水泥浆的胶凝强度逐渐增加并由液态转变为胶凝态，由于胶凝悬挂力分担了部分原静液柱压力，水泥环内部孔隙水的压力逐渐降低。结合该理论，本文给出胶凝强度对水泥环压力的影响关系模型。

1.2.1 水泥环应力状态

针对水泥环内部压力下降的现象，Drecq等[20]，Parcevaux[21]，以及Erik等[22]将水泥浆的水化过程比作沉积土固结成石的过程，并通过实验数据对比发现土力学中的应力状态模型——太沙基定律可以准确地描述水泥环内部应力状态：

1.2.2 水泥浆总应力

假设水泥浆的总应力是恒定不变的，即该位置处的上覆水泥环压力为[20-22]：

1.2.3 水泥浆有效应力

水泥浆的有效应力与水泥浆的静胶凝强度有关，并随时间不断变化。Moore[23]提出应用经典的应力应变方程计算移动水泥环所需要的压力和水泥浆胶凝悬挂力。根据水泥环单元格受力示意图可知（见图2），水泥环单元格主要受到4个作用力。

①上部水泥环单元格施加向下的作用力：

②下部水泥环单元格施加向上的作用力：

③水泥环单元格自身重力：

④井壁以及套管外表面胶凝悬挂力：

根据图2可得水泥环受力平衡方程为：

由（12）式可知，水泥环单元格的有效应力等于单元格自重加上上下单元格压差，即水泥环单位横截面积上的胶凝悬挂力为水泥浆的有效应力。结合（8）式—（11）式，得到水泥浆有效应力表达式：

1.2.4 瞬态压力场计算模型

结合（6）式—（7）式以及（13）式，可以得到水泥浆孔隙压力的表达式为：

水泥浆的孔隙压力随胶凝强度的增加不断减小，该过程与水泥浆水化反应密切相关。因此，水泥浆内孔隙压力的计算应考虑水泥浆水化反应的影响。当水泥浆的有效应力增加到一定程度，水泥浆的强度便能够支撑起自身的重量，该时刻水泥浆内的孔隙压力下降到孔隙水的静液柱压力并维持恒定[21]：

图2 水泥环单元受力示意图

结合水泥浆的水化动力学计算模型，水泥浆孔隙压力可以表达为：

基于（16）式，深水固井过程中不同时刻以及不同深度处水泥浆孔隙压力的表达式为：

1.3 不同养护温度、压力下的水泥浆水化动力学模型

化学反应动力学主要研究化学反应的内在因素（反应物的浓度、状态等）以及外在因素（温度、压力等）对反应速率以及反应过程的影响。借助水泥浆水化动力学能够很好地描述温度、压力与水化反应速率之间的相互作用。水化反应动力学的直接研究对象为水化度，定义为水泥内部已经水化的水泥占原水泥的质量百分比。水化度通常由水化热等易量化的参数来进行描述，水化度可以表达为：

Krstulovic与Dabic[24]提出的模型被认为是最经典的水化反应动力学模型之一，广泛应用于工业领域。基于水泥浆水化热释放曲线，Krstulovic-Dabic模型（K-D模型）给出了水泥浆水化度和水化时间之间的定量关系。根据该模型，水泥浆水化反应包含3个过程：成核结晶与晶体生长过程（NG）、相边界反应过程（I）以及扩散过程（D）。这3个过程可以同时发生，也可以只发生一个或两个，但是整个水化反应速率取决于其中最慢的反应过程。根据K-D模型，初始阶段控制水化反应的是NG过程，而后逐渐被I过程或D过程代替，3个反应阶段的数学模型分别表达为[24]：

此外，温度和压力主要通过影响化学反应速率常数来改变化学反应进程。结合阿伦尼乌斯公式描述化学反应速率与温度、压力的关系[2]，即：

由（19）式—（21）式可知，K-D模型求解的关键是各阶段的反应速率常数以及反应级数。等温量热计可以连续获取水泥浆水化过程中热量的释放曲线，根据（18）式描述的水泥浆水化度与释放热量之间的关系，可以获得水泥浆水化度随时间变化的曲线。在此基础上，对（19）式两边取对数便可得到-ln（1-α）与t的关系曲线，再根据曲线的斜率与截距可以确定n和KNG的值；使用同样方法，根据（20）式—（22）式，通过绘制相应对数曲线可以获得KI、KD以及Ea的值。

2 温度压力耦合数值求解方法

本文利用数值差分方法来耦合求解水泥浆水化动力学模型、温度场模型以及压力场模型，并给出了优化耦合迭代的计算流程。

2.1 初始和边界条件

2.1.1 初始条件

①初始压力等于水泥浆静液柱压力：

②水泥浆初始水化度为零：

③井筒内水泥浆的初始温度可根据Sun等[18-19]给出的井筒循环温度预测模型进行计算：

④初始地层温度可由地温梯度进行计算：

2.1.2 边界条件

①在井筒与地层的边界处，单元格内能的变化量是环空水泥浆传入的热量与传出到外部地层的热量之和：

②地层在无穷远处不受扰动，温度保持初始值不变：

③在泥线井口位置处，水泥浆的孔隙压力为环境压力：

2.2 单元格划分

考虑到井眼直径远小于井筒长度，故将套管与水泥环划分成一维单元格；而地层相对于井筒为轴对称，故可划分成二维单元格。图3为计算过程中单元格划分示意图。

图3 单元格划分示意图

2.3 方程离散

根据有限差分方法的计算原理，将（5）式和（17）式以及（19）式—（21）式转化成有限差分迭代格式：

2.4 耦合迭代流程

为了清晰地表达模型耦合迭代求解过程，结合初始和边界条件、单元格划分以及差分格式，绘制了求解流程图（见图4）。

图4 耦合迭代求解过程示意图

3 模型验证

3.1 水化反应动力学模型验证

为了验证本文模型的精确度，将新模型应用到Dillenbeck等[25]的室内实验中。Dillenbeck等为了探讨水泥浆水化放热对井筒温度的影响，利用等温量热计分别测量了26.6 ℃和52.4 ℃时水泥浆的水化放热曲线（见图5）。

图5 不同温度下水泥浆水化热释放曲线

根据（18）式，利用水泥浆水化热曲线可得到水泥浆水化度随时间的变化曲线，由此可获取K-D水化动力学模型参数，如表1所示。

表1 水泥浆K-D水化动力学模型参数

根据表1不同温度下水化动力学参数，可得到NG、I和D过程的化学反应活化能分别为20 787.4，25 407.5，40 436.3 J/mol。结合（22）式，可以定量描述不同温度下水泥浆水化反应过程。图6所示为本模型计算结果与Dillenbeck实验测量结果对比，结果显示本模型能够很好地描述不同温度下水泥浆水化反应过程。

图6 不同温度下模型计算结果与实验测量结果对比

为了进一步验证本文模型的准确性，笔者还分析了阎培渝等[26]针对硅酸盐水泥浆、掺膨胀剂硅酸盐水泥浆在不同温度下的水化动力学实验，以及Pang等[2]针对某H级水泥浆体系开展的不同压力下的水化动力学实验。如图7所示，通过对比实验测量数据以及模型计算的瞬态水泥浆水化度，发现本文模型在不同水泥浆体系、温度、压力条件下均能准确反映水化反应过程，满足现场工程要求。

3.2 固井井筒温度压力验证

3.2.1 Dillenbeck现场试验验证

Dillenbeck等[25]将电子温度传感器下入现场井筒，进而实时监测固井过程中井筒温度瞬态变化。如图8所示，将Dillenbeck现场实测温度与本文模型计算温度进行了对比。对比结果显示，本文温度预测模型在考虑水化反应影响的前提下最大误差为5.6%，忽略水化反应影响时温度预测误差可高达24.5%，表明水泥浆水化反应对井筒温度的影响不可忽略。

3.2.2 Cooke现场试验验证

图7 不同研究条件下水泥浆水化动力学实验数据与本模型计算结果对比

图8 固井候凝期间水泥浆温度计算结果与现场数据对比

为了进一步验证固井井筒温度压力预测模型的准确性，将新模型计算结果与Cooke等[6]的现场测量结果进行对比。Cooke等在固井过程中将4个传感器下入到井筒内并依附于套管的外壁，传感器分别位于1 108，1 673，2 106，2 668 m深度处。利用一根录井电缆收集各个传感器获取的数据，以此测量水泥环内不同深度的温度和压力。实验井总井深为2 712 m，采用外径为73 mm的套管，井眼尺寸为200 mm，所在区块地温梯度为2.12 ℃/100 m。

对比了本文模型的计算结果与Cooke实验测量结果，发现本文建立的瞬态温度压力计算模型能够较好拟合实验结果，满足工程实践要求（见图9、图10）。

图9 瞬态温度计算结果与现场测量数据对比

由温度变化趋势看出，由于早期水泥浆水化反应释放出大量的热，在该阶段水泥浆温度呈现出明显的增长趋势；随着水化反应速率的降低，水泥浆的温度逐渐接近地层温度。而实测压力和计算压力都随着时间而不断降低，根据计算结果可知，随着水泥浆胶凝强度的提升，水泥浆孔隙压力逐渐降低至孔隙水静液柱压力。

图10 瞬态压力计算结果与现场测量数据对比

4 深水固井候凝阶段井筒温度、压力及水化度瞬态演化规律

为了分析深水固井候凝期间井筒内温度、压力及水化度的演化规律，利用新建立的模型对一口深水模拟井进行数值模拟计算，并分析其演化规律。模拟井采用Dillenbeck实验中的水泥浆体系，井身结构如图1所示，基础参数见表2。

表2 模拟井基础参数

4.1 瞬态温度演化规律

图11 不同时间环空温度剖面演化

图11所示为固井候凝期间水泥浆温度剖面随时间变化曲线，其中0 h表示固井循环结束，即候凝阶段开始的时刻。对于不同位置处的水泥浆，温度都是从初始循环温度逐渐变为地层温度。在该过程中由于水泥浆水化反应释放出大量的热，导致水泥浆温度呈现复杂的演化规律。图12所示为固井候凝期间井筒不同位置处水泥浆瞬态温度变化曲线，水泥浆温度演化过程可以划分为3个阶段：阶段Ⅰ较为短暂，主要为水泥浆与套管内流体之间的热交换过程，在该过程中水泥浆温度略微下降，两者温度差缩小，2 250 m及2 500 m深度处水泥浆与套管内流体温度接近，所以温度降低过程不明显；水泥浆水化反应进入加速期标志了阶段Ⅱ的开始，该阶段大量的水化热被释放，导致水泥浆温度的迅速升高并达到某一峰值；阶段Ⅲ水泥浆水化反应进入了减速期，水化放热速率大大降低，该阶段主要是水泥浆与地层之间发生缓慢的热交换，水泥浆温度逐渐趋向于地层温度。

图12 不同位置处环空内温度演化

4.2 瞬态压力演化规律

如图13和图14所示，随着水泥浆的水化，其孔隙压力随着时间逐渐下降，从不同时间环空压力剖面可以观察到10.5 h后水泥浆孔隙压力开始低于地层压力，此时有发生气窜的风险。水泥浆必须在孔隙压力降至地层压力之前，产生足够的强度来预防气窜。本文模型能够为抗气窜水泥浆体系设计提供理论依据。

图13 不同时间环空压力剖面演化

图14 不同位置处环空压力演化

4.3 水泥浆水化度演化规律

如图15和图16所示，水泥浆的水化度随时间逐渐增加并趋向于一个最终值。在地温梯度与液柱压力的作用下，水泥浆温度、压力随着深度的增加而增加，高温高压环境能够加速水泥浆水化反应速率，因此井筒深处水泥浆具有更高的水化度。模拟结果表明，2 500 m处的水泥浆水化度在14.1 h后达到0.5；而位于1 500 m（泥线附近）处的水泥浆水化度需要43.8 h才能达到0.5。因此，深水固井作业中，泥线附近的低温环境会延长水泥浆的候凝时间从而延长固井工作周期。

图15 不同时间水泥浆水化度剖面演化

图16 不同位置处水泥浆水化度演化

5 结论

考虑固井候凝期间水泥浆温度、压力以及水化反应之间的复杂耦合作用，本文基于水泥浆水化反应动力学建立了固井井筒温度与压力耦合预测模型，并利用新建立的模型对一口深水模拟井进行了数值模拟计算。忽略水泥浆温度、压力与水化反应之间的耦合作用会给井筒温度压力预测带来较大误差，而经过与实验测量数据对比，本文模型计算精度在5.6%以内，能够较好地满足工程要求。水泥浆温度在水化热的作用下迅速升高；水泥浆孔隙压力随着水泥浆胶凝强度的提高逐渐降低，甚至低于地层压力，从而引发气窜风险。瞬态变化的温度压力场又会反作用于水泥浆水化反应速率，井筒深部的水泥浆在高温高压环境下具有更快的水化反应速率。对于深水固井作业来说，要注意泥线附近的低温环境会延长水泥浆的候凝时间，从而延长固井工作周期。本文建立的新模型能够协助固井工程师设计水泥浆体系，及时发现与规避安全隐患。

符号注释：

Aa——环空的横截面积，m2；Ac——套管横截面积，m2；cc——环空内水泥浆的比热容，J/(kg·K)；ce——地层比热容，J/(kg·K)；cf——套管内钻井液的比热容，J/(kg·K)；di——水泥环的内径，m；dw——水泥环的外径，m；Ea——化学反应的活化能，J/mol；F——井壁、套管外表面胶凝悬挂力，N；F1——上部水泥环单元格作用力，N；F2——下部水泥环单元格作用力，N；g——重力加速度，m/s2；G——水泥环单元格自身重力，N；gG——地温梯度，℃/m；i——沿井筒方向单元格编号；j——垂直井筒方向单元格编号；kc——环空内水泥浆的导热系数，W/(m·K)；ke——地层导热系数，W/(m·K)；kf——套管内钻井液的导热系数，W/(m·K)；K——化学反应速率常数，s-1；KD——D过程中的水化反应速率常数，s-1；KI——I过程中的水化反应速率常数，s-1；KNG——NG过程中的水化反应速率常数，s-1；Kr——参考温度压力下化学反应速率常数，s-1；m——时间离散节点；M——垂直井筒方向地层单元格总数；n——化学反应级数；N——沿井筒方向单元格总数；p——水泥浆孔隙压力，Pa；pen——环境压力，Pa；pfinal——水泥浆最终孔隙压力，Pa；pr——化学反应的参考压力，Pa；p0——井口回压，Pa；p1——上部水泥环单元格压力，Pa；p2——下部水泥环单元格压力，Pa；P——任意时刻水泥浆某种特性的特征值，比如胶凝强度；P∞——水泥浆某种特性的最终值；Q——某一时刻水泥浆水化热，J/kg；Qc——水化热量，J；Qca——水泥浆与套管间交换热量，J；Qf——水泥浆与地层间交换热量，J；Q∞——水泥浆最终水化热，J/kg；rci——套管内径，m；rw——井筒的半径，m；R——气体常数，8.13 J/(mol·K)；t——水泥浆水化时间，s；T——化学反应的温度，K；Tc——环空水泥浆的温度，K；Te——地层温度，K；Te,0——井筒与地层界面处温度，K；Te,1——紧邻井筒与地层界面处第一个单元格内温度，K；Tf——套管内钻井液温度，K；Tr——化学反应的参考温度，K；Ua——从环空到地层的总传热系数，W/(m2·K)；Uc——从套管到环空的总传热系数，W/(m2·K)；vc——环空内水泥浆的循环流速，m/s；vf——套管内钻井液的循环流速，m/s；wc——环空水泥浆质量流量，kg/s；wf——套管内钻井液质量流量，kg/s；x——与井筒的水平距离，m；z——沿井筒方向到井口的距离，m；α——水泥浆水化度；α1——NG过程结束进入I过程时的水化度；α2——I过程结束进入D过程时的水化度；α500——胶凝强度达到239 Pa（500 lbf/100 ft2）时的水化度；γ——水泥浆容重，N/m3；ρc——环空内水泥浆的密度，kg/m3；ρe——地层密度，kg/m3；ρf——套管内钻井液的密度，kg/m3；ρw——孔隙水的密度，kg/m3；σ——水泥浆的总应力，Pa；σ′——水泥浆的有效应力，Pa；τ——水泥浆胶凝强度，Pa；τ500——胶凝强度500 lbf/100 ft2对应的标准单位下的数值，即239 Pa；Δt——时间步长，s；ΔV——表观活化体积，m3/mol；Δx——地层单元格垂直井筒方向网格步长，m；Δz——水泥环单元格高度，m。