当前在深层、超深层油气资源勘探开发中面临着诸多理论和技术难题^[1]，主要表现在：高温高压环境下，环空温度、压力预测难度大，井底压力控制困难；气侵时井筒中往往呈现气液固多相流，流动规律复杂；安全密度窗口窄，易出现漏、溢及卡钻等井下故障。归根结底，核心问题在于对高温高压条件下环空多相流的流动特性认识不清，无法准确掌握环空温度和压力的变化规律，因此，开展深井、超深井多相流水力学模拟是必要的。井筒温度场和压力场相互影响，循环温度、压力的变化均会引起井筒流体热物性参数和流动参数改变，模型求解需要耦合迭代，其涉及井筒水力学、多相流以及传热学等多学科交融。目前，国内外学者针对井筒温度压力场耦合进行了相关研究^{[2,3,4,5,6,7]}，但是已有研究均未基于瞬态温度、压力，同时忽略了钻井流体物性参数是随温度压力的变化而变化的，无法真实模拟高温高压条件下井筒温度压力场耦合变化规律。为此，笔者从质量、动量和能量守恒方程出发，建立了井筒多相流全瞬态温度压力场耦合模型，分析了环空多相流瞬态温度和压力耦合变化规律，以期为深井、超深井水力参数设计提供参考。

1 井筒物理模型

在正常钻井作业中，钻井液以温度T_i由钻柱内进入，钻柱内充满钻井液，后经过钻头喷嘴进入环空，在井底混合天然气、岩屑上返，环空内为气液固多相流体，此时井底温度低于地层温度，多相流体不断从井壁岩石吸收热量，当多相流体达到某一井深时，其温度将大于地层温度，多相流体在上返过程中加热上部井壁岩石，最终多相流体以温度T_o由环空出口返出，如图 1所示。在这一过程中，钻柱内的流体、钻柱壁、环空内的流体、井壁、地层均会在轴向和径向上进行热能传递，钻柱内为单相传热，环空内为两相传热。随着井筒内热量的交换，循环流体的温度实时变化，从而影响井筒压力剖面，反过来井筒压力剖面的改变进一步影响井筒温度场，这是一个相互耦合的过程。

2 井筒瞬态温度压力数学模型 2.1 假设条件

在建立井筒瞬态温度压力数学模型时，进行以下假设：1)钻井液在井筒内作一维流动，忽略径向变化；2)钻井液为水基钻井液，忽略天然气在钻井液中的溶解；3)钻柱、环空和地层均为非稳态传热，其中钻柱壁和井壁边界加入轴向热传导项；4)井筒内流体的轴向热传导与轴向热对流相比可以忽略不计；5)考虑岩屑对于井筒温度和压力的影响。

2.2 质量守恒方程

气相、液相和固相的质量守恒方程分别为：

式中：A为流道面积，m²；ρ_g，ρ_l和ρ_s分别为气相、钻井液、岩屑的密度，kg/m³；α_g，α_l和α_s分别为气相、钻井液、岩屑的体积分数；v_g，v_l和v_s分别为气相、钻井液、岩屑的实际流速，m/s；q_s为单位厚度岩屑产生速度，kg/(s·m)。

2.3 动量守恒方程

动量守恒方程为：

式中：g为重力加速度，m/s²；p_fa为环空内多相流摩阻压降，Pa。

根据文献^[8,9,10]的研究成果,垂直环空气液两相流的流型可划分为泡状流、分散泡状流、段塞流、搅动流和环状流5类。2003年，C.Perez-Tellez^[11]主要针对段塞流进行了压力梯度预测模型的修正(记为Perez模型)，并将数值计算结果与路易斯安那大学全尺寸试验井的试验数据进行了对比，井筒压力平均误差不大于3.3%。因此，气液两相流压降计算采用高精度的Perez模型。

2.4 井筒温度数学模型

钻井液正常循环期间，地层与环空钻井液进行热交换，环空钻井液与钻柱内钻井液进行热交换，钻井液在井内的整个循环过程可以看作是一个具有一定边界条件的热交换器。如图 1所示，井筒温度数学模型在径向具体可以分为5部分，分别为钻柱内、钻杆壁、环空、套管和地层。根据能量守恒原理，得到各组分的温度方程。

钻柱内流体的温度方程为：

钻柱壁的温度方程为：

环空内流体的温度方程为：

井壁的温度方程为：

地层的温度方程为：

式中：ρ为密度，kg/m³；q为流体流量，L/s；Q为井筒摩擦力功率，W/m；T为温度，K；h为对流换热系数，W/(m²·K)；k为导热系数，W/(m·K)；c为比热容，J/(kg·K)；z为轴向坐标，m；r为径向坐标，m。下标：l为液相；g为气相；p为钻柱；w为钻柱壁；a为环空；c为套管；f为地层；pi为钻柱内壁；po为钻柱外壁；ci为套管内壁；co为套管外壁。

气液两相流对流换热系数的计算与单相流相比有很大的差异，为了满足现场对精度的要求，选取不同流型下与试验数据相比误差较低的气液两相对流换热经验公式^[12,13]。在垂直流动过程中，对于泡状流，建议采用Aggour模型；对于分散泡状流，建议采用Knott等人的模型；对于段塞流和搅拌流，建议采用Rezkallah-Sims模型；对于环状流，建议采用Ravipudi-Gobold模型。

2.5 辅助参数求解 2.5.1 天然气的黏度、比热容和导热系数

黏度、比热容和导热系数等热物性参数均是温度、压力的函数，温度和压力变化时，它们亦随之改变。为了保证井筒温度、压力计算的准确性，分别采用与试验数据吻合较好的A.L.Lee^[14]、D.G.Friend^[15]、A.K.Lateef^[16]等人提出的方法计算天然气的黏度、比热容和导热系数。

2.5.2 天然气和钻井液的表面张力

表面张力在多相流计算中非常重要，对环空中气相的滑脱速度和对流换热系数均有影响。与现有模型相比，Sutton模型^[17]的计算精度最高，因此，天然气和水基钻井液的表面张力采用该模型计算。

2.5.3 钻井液物性参数

深井井筒处于高温、高压环境下，此时钻井液的密度/塑性黏度不再是常数。笔者采用解析法推导钻井液物性参数(密度、塑性黏度)与温度、压力的关系式。根据W.C.McMordie Jr.等人^[18]通过试验测出的不同温度、压力下水基钻井液的密度，拟合出水基钻井液密度与温度、压力关系式的系数b_p，b_pp，b_T，b_TT和b_pT分别为4.922 4×10^-10，-9.687 7×10^-19，-3.219 6×10^-4，-1.743 2×10^-6和4.918 6×10^-13。同理可拟合出水基钻井液的塑性黏度表达式。

式中：p₀为地面压力，MPa；T₀为地面温度，K；ρ(p,T)为压力p和温度T下钻井液的密度，kg/m³；μ(p,T)为压力p和温度T下钻井液的塑性黏度，mPa·s。

3 多相流瞬态温度压力场耦合求解 3.1 边界条件

钻井液注入温度和地表温度已知，作为井筒温度边界条件。同时，整个钻井过程中，井口回压维持不变，作为井筒压力边界条件。

3.2 模型离散

井筒多相流瞬态温度压力场耦合求解的先决条件是实现井筒瞬态温度、压力的独立求解。井筒瞬态温度和压力模型的求解分别采用Crank-Nicolson全隐式格式^[19]、以及四点差分显式格式^[20]，网格划分时选用定轴向、定时间及变径向网格。

针对井筒瞬态温度模型的全隐式处理，控制方程中的一阶空间导数采用一阶迎风格式，一阶时间导数采用两点向后差分，二阶空间导数采用三点中心差分。以此原则将井筒温度模型离散化，给出相应的有限差分格式。

钻柱内温度模型的有限差分格式为：

钻柱壁处温度模型的有限差分格式为：

环空内温度模型的有限差分格式为：

井壁处温度模型的有限差分格式为：

地层温度模型的有限差分格式为：

其中

针对井筒瞬态压力模型的显式处理，控制方程中的一阶空间导数采用一阶迎风格式，一阶时间导数采用四点中心差分。以此原则将井筒压力模型离散化，给出相应的有限差分格式。

液相压力模型的有限差分格式为：

气相压力模型的有限差分格式为：

固相压力模型的有限差分格式为：

动量守恒方程的有限差分格式为：

其中

3.3 求解算法

通过迭代方式对井筒多相流瞬态温度压力场耦合模型进行求解，具体步骤如下：1)针对井筒划分网格，在时间和空间上进行离散化处理；2)已知的地层温度作为井筒初始温度，代入井筒压力计算中，将稳态压力计算结果作为初始压力，同时给出井筒各组分流体的热物性参数以及体积分数；3)根据初始时刻的计算值，预估t=t₀+Δt时刻的井筒温度T_j⁰、压力p_j⁰以及各组分流体热物性参数和体积分数，代入瞬态温度模型计算t=t₀+Δt时刻的井筒温度T_j；4)将t=t₀+Δt时刻的井筒循环温度代入瞬态压力模型；5)求解环空内单点参数时，需要先后估算节点压力和气相体积分数；6)节点循环至井口，判断计算的井口回压是否满足精度要求，不满足时返回4)进行校正计算，直至得到满足精度要求的井筒压力p_j；7)判断井筒环空压力和温度是否满足精度要求，即‖p_j-p_j⁰‖<ε_p且‖T_j-T_j⁽⁰⁾‖<ε_T，不满足则返回3)重新进行校正计算；8)t=t₀+Δt时刻计算完毕，耦合迭代计算进入下一时刻。

4 算例分析

利用上述井筒多相流瞬态温度压力场耦合模型，采用迭代求解算法，对塔里木油田某深井的井筒瞬态温度压力耦合变化规律进行了分析。由于已探明储层中不含有H₂S等酸性气体，因此，求解过程中将侵入气体看作纯甲烷气体。计算时用到的基础参数为：井深6 500.00 m；套管内径157.0 mm，套管下深5 500.00 mm；钻头直径152.4 mm；钻杆外径88.9 mm，钻杆内径70.2 mm；喷嘴面积233 mm²；循环排量12 L/s；气侵量0.6 m³/s；机械钻速2.7 m/h；井口回压1 MPa；地表温度25 ℃；地温梯度0.03 ℃/m；钻井液的地面密度1 500 kg/m³，比热容3 900 J/(kg·K)，导热系数1.73 W/(m·K)，入口温度40 ℃；钻柱的密度7 800 kg/m³，比热容400 J/(kg·K)，导热系数43.75 W/(m·K)；水泥的密度1 900 kg/m³，比热容2 000 J/(kg·K)，导热系数1 W/(m·K)；岩石的密度2 640 kg/m³，比热容837 J/(kg·K)，导热系数2.25 W/(m·K)。

4.1 环空温度剖面

图 2为计算出的不同循环时间下的环空温度剖面，图中t=0时刻的环空温度即为初始地层温度。

由图 2可知，环空内多相流体温度的大体规律为：环空井底温度低于地层温度，随着井深减小逐渐降低，径向强制对流换热能量由地层指向环空；当井深达到临界井深(约1 500 m)时，环空温度开始高于地层温度，环空-地层径向能量交换反向；环空多相流体循环最高温度处距离井底的距离为井深的十分之一(井深6 000 m附近)，而环空单相钻井液循环最高温度处距离井底的距离为井深的1/7～1/6，这主要是由气液两相和单相对流换热系数的差异所导致的。

图 3为计算出的不同循环时间下的井口和井底温度。

由图 3可知：随着循环时间增长，井底温度逐渐降低，且降低幅度逐渐变小；井口温度随循环时间增长而升高，且初期升高幅度大，但快速趋于稳定。原因是，近井口段气体膨胀最为明显，两相流由泡状流转换为段塞流，环空钻井液与钻杆外壁及井壁之间的对流换热系数快速增大，导致近井口段微元控制体径向换热充分，井口温度变化迅速趋于平衡。

4.2 岩屑运移速度剖面

图 4为计算出的不同循环时间下的岩屑运移速度剖面。

由图 4可知，岩屑运移速度分别在井深5 500.00和1 200.00 m处变化最为明显，发生阶跃现象。井深5 500.00 m为套管鞋处，此处环空流道面积突然增加，岩屑运移流速突然降低；此后随着井深减小，岩屑流速逐渐升高，原因是在该井段气液两相流流型为泡状流，而泡状流有利于携岩^[21]，因而岩屑运移速度逐渐升高；直至井深1 200.00 m处，岩屑运移速度突然降低，且随着井深减小，岩屑流速继续降低。这是因为此处为泡状流-段塞流转换点，在这里小气泡膨胀聚集形成大的泰勒气泡和液相段塞，而段塞流不利于携岩，因而岩屑流速降低。此外，从图 4还可以看出，随着循环时间增长，岩屑运移速度剖面变化较小，这也说明岩屑运移速度基本不受井筒循环温度的影响。

4.3 物性参数剖面

图 5为计算出的不同循环时间下的钻井液密度剖面。

从图 5可以看出：随着循环时间增长，临界井深1 500.00 m以浅钻井液的密度逐渐增大，这是因为井筒环空循环温度逐渐降低，钻井液发生压缩，密度随之增大；当t=0时刻，即环空温度为地层温度，钻井液的密度随着井深增加逐渐降低。在循环温度影响下，随着井深增加，钻井液的密度先降低后升高，这是温度和压力对钻井液双重作用的结果；环空的最高温度在井深6 000.00m处，随着井深进一步增加环空温度降低，受温度降低和压力升高产生压缩效应的影响，钻井液的密度反而升高。

图 6为计算出的不同循环时间下的钻井液塑性黏度剖面。

从图 6可以看出：随着循环时间增长，临界井深1 500.00 m以浅钻井液的塑性黏度逐渐增大；当t=0时刻，即环空温度为地层温度，钻井液的塑性黏度随着井深增加逐渐降低；在循环温度的影响下，随着井深增加，钻井液的塑性黏度先降低后升高。

从图 5和图 6可以看出，循环8 h后环空井底钻井液的密度由1 360 kg/m³增至1 460 kg/m³，塑性黏度由8.6 mPa·s增至13.8 mPa·s。由以上分析可知，钻井液的物性参数(密度、塑性黏度)受井筒循环温度的影响较大，在建立高温高压井水力学模型时，需要充分考虑钻井液密度和塑性黏度的变化，进而准确计算井筒压力和温度，为现场钻井作业提供可靠的参数。

4.4 井底压力动态变化特征

图 7为计算出的不同循环时间下的井底压力。

由图 7可知：随着循环时间增长，井底压力先迅速降低后逐渐升高，且升高幅度逐渐变小，井底压力渐趋于稳定；在循环时间为0～0.2 h时，井底压力呈线性降低，而在0.2 h以后井底压力呈对数升高。这是因为在循环初期，相比初始地层温度，整个环空温度大多升高且变化幅值较大，仅近井底段环空温度降低，井底压力受井筒环空温度升高的影响逐渐降低。当循环时间达到0.2 h后，环空近井底段温度继续下降，而近井口段温度升高且快速趋于稳定，这样井底压力受近井底段环空温度降低的影响逐渐增大。

4.5 井底压力变化敏感性参数分析

图 8—图 12分别为计算出的不同气侵量、钻井液地面密度、钻井液地面塑性黏度、循环排量、井口回压下井底压力随循环时间的变化曲线。

从图 8—图 12可以看出，井底压力随钻井液的地面密度和塑性黏度，循环排量、井口回压增大而增大，而井底压力随气侵量增大而减小。对比循环8 h后气侵量、钻井液的地面密度和塑性黏度、循环排量和井口回压对井底压力的影响，可知当气侵量从0.2 m³/s升至0.8 m³/s循环8 h后井底压力变化值升高1.2 MPa，钻井液地面密度从1 300 kg/m³升至1 600 kg/m³时升高0.4 MPa，钻井液地面塑性黏度从20 mPa·s升至50 mPa·s时仅升高0.04 MPa，排量从10 L/s升至16 L/s时升高0.2 MPa，井口回压从1.0 MPa升至4.0 MPa时降低0.2 MPa。由此可知，气侵量对井底压力的影响最大，钻井液地面密度、排量、井口回压次之，钻井液地面塑性黏度最小。根据不同参数对井底压力的影响程度， 在钻井过程中，可以通过优化排量和钻井液地面密度以及施加一定的井口回压，使井底压力位于安全窗口之内的同时，尽可能降低由循环温度变化引起的井底压力波动，降低钻井作业风险。

5 结 论

1) 基于井筒多相流、传热学理论，充分考虑循环流体物性参数随温度压力的变化，建立了适用于深井、超深井的井筒多相流全瞬态温度压力场耦合模型，并提出了迭代求解算法，保证了模型的稳定性和收敛性。

2) 环空多相流体最高循环温度位于距离井底十分之一处。随着循环时间增长，岩屑运移速度基本保持不变，而临界井深1 500.00 m以浅钻井液的密度/塑性黏度逐渐升高。随着井深增加，钻井液的密度、塑性黏度先降低后升高，钻井液密度增幅较大。

3) 井底压力随循环时间增长先呈线性降低后呈对数升高。气侵量对井底压力的影响最大，钻井液地面密度、排量、井口回压次之，钻井液地面塑性黏度最小。

4 ) 钻井过程中，适当降低排量和钻井液密度以及施加一定的井口回压，可以有效降低井底压力波动，确保将井底压力维持在安全窗口之内，从而降低钻井作业风险。