基于热-流-固耦合模型的石油钻井施工过程数值分析

贾善坡，邹臣颂，王越之，谭贤君，甘新星

（1.长江大学 城市建设学院，湖北 荆州 434023；2.长江大学 油气资源与勘探技术教育部重点实验室，湖北 荆州 434023；3.中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，武汉 430071）

1 引 言

近年来针对岩石多物理场耦合特性的研究已成为国内外岩土工程界关注的热点之一，我国实施与规划中的深部资源开采、放射性废料地质深埋处置、石油天然气地下能源储存等工程，都涉及到多孔介质中的热-水-力（THM）耦合问题。Noorishad等基于扩展的Biot固结理论[1]，首次提出了饱和岩土介质的固-液-热耦合基本方程。在此之后，许多学者对THM耦合理论及数值方法进行了进一步的探讨和研究，特别是从20世纪90代开始，美国、英国、加拿大等国以及欧盟的 12个研究小组开始进行 DECOVALEX的国际合作研究计划，采取不同的方法，对核废料储库围岩体THM耦合行为进行理论和试验研究，并取得了瞩目的成果[2－3]。在80年代末与 90 年代中期，陆续出现过若干 THM模型的专用程序，如THAMES、MOTIF、FRACON、FEMH、FRIP、FRACTURE、GEORACK等，这些程序大都是专用软件，算法的适应性有限，在前后处理、单元库和材料库等方面不及大型商业软件[3]；软件FLAC、UDEC和COMSOL也可用于THM耦合问题计算，但这些程序在求解多种材料组成、复杂的施工过程和三维问题时还存在一定困难。

Rutqvist等[4]将FLAC软件和TOUGH2软件相结合，编制了THM耦合程序，并将该程序应用于核废料地下储库的模拟。陈波等[5]推导了固-液两相介质的温度场、变形场、渗流场三场耦合问题的有限元格式，开发了三场耦合数值分析软件 CDST。张玉军[6]研发出了一个用于分析饱和岩土介质中热-水-应力耦合弹塑性问题的二维有限元程序，并通过2个算例进行了应用尝试。冯夏庭等[7]建立了岩体弹性、弹塑性、黏弹塑性的THMC分析模型，开发了数值分析软件系统，并分析了开挖损伤区温度-水流-应力-化学耦合作用行为。陈益峰等[8]通过选取位移、水压、气压、蒸气压、温度和孔隙率为基本未知量，建立了有限元计算格式，研发三维八自由度多相流THM全耦合有限元程序THYME3D，并采用膨润土THM耦合Mock-up试验对数值模型和计算程序进行验证。傅少君等[9]根据饱和土体热固结问题中温度、渗流及应力之间的复杂耦合关系，系统地探讨了二维有限单元方法的求解过程，并在FORTRAN 语言环境下研制了相应的计算机分析程序 APOSE-T。赵延林等[10]在建立双重介质热-水-力耦合微分控制方程的基础上，开发了双重介质热-水-力耦合三维计算程序DTHM。

笔者在前人研究的基础上，推导出岩土介质温度-渗流-应力耦合的数学模型模型，然后，以MATLAB语言为平台，对ABAQUS软件作为求解器，编制THM耦合计算程序，并通过算例对其正确性进行了试算验证，详细分析了石油钻井施工过程中井壁围岩的热-流-固耦合作用过程。

2 温度-渗流-应力耦合数学模型

2.1 岩土介质热传递模型

取任意微元体，设微元体的孔隙度为n，则单位体积中固相骨架所占的空间为1-n，根据能量守恒原理，结合Fourier定律，给出固相骨架的能量守恒方程[4]：

式中：λs、cs、ρs分别为岩土骨架的导热系数、比热和密度；qs为单位时间内单位体积骨架产生的能量；T为任意时刻的温度；T0为初始温度；βl为线膨胀系数；K为岩土介质体积模量；εv为体积应变。

采用的类似过程可建立流体能量守恒方程。由于流体的渗流速度很小，忽略流体的流动动能，考虑对流和传导的能量方程可表示：

式中：λf、cf、ρf分别为流体的导热系数、比热和密度；qf为单位时间内单位体积流体产生的能量；v为流体的 Darcy渗流速度；(v∇)T称为对流项，表示流体质点移动时引起的温度变化率。

根据混合物理论[3]，按物质组分比例进行叠加，得到岩土介质的总能量方程：

2.2 岩土介质流固耦合渗流模型

结合Darcy定律，给出考虑介质变形的渗流方程[1－2]：

渗透系数、孔隙度与体积应变的关系[1]为

式中：k0、n0分别为初始渗透系数和孔隙度。

2.3 岩土介质热流固耦合力学模型

根据Biot有效应力原理，有效应力可表示为

式中：σij为总应力张量；p为孔隙流体压力，规定以压力为正；α为Biot 系数；δij为Kronecker符号。

为了描述岩土介质的热膨胀及塑性变形等行为，总应力可采用增量形式：

式中：De为弹性刚度张量；ε为总应变张量；εp为塑性应变张量；Tε为温度变化时产生的应变张量。

根据弹塑性理论，式（8）可写为

式中：

ui为岩土骨架的位移分量。

岩土强度准则采用修正的 Mohr-Coulomb准则，屈服函数F和塑性势函数G具体表达式[1]为

3 耦合模型的程序实现

3.1 程序结构

考虑到数值收敛和稳定性的要求，本文以解耦的方法开发软件。热-流-固耦合模型的求解包括温度场、流-固耦合场分别求解以及各物理场之间的解耦关系。

虽然，流-固耦合场和温度场的数值计算有较大差别，但本质上包含了线性化和时步离散（或载荷增量）两个基本内容，即将流-固体场和温度场计算按两个独立的系统分别设计，通过数据通讯的方式，实现在每一时步上的参数耦合；流-固体耦合是以渗透压力变化和温度变化为载荷的，而温度场的计算又是以流-固体耦合变形为基础的，在每一时步上不断修正相关的系数，而这种相互修正又是在一系列时步或载荷上交叉进行的。

以MATLAB为平台，将ABAQUS做为流-固耦合场和温度场的求解器，开发岩土介质多场耦合计算软件，各计算模块之间数据的存储和通讯通过编制ABQMAIN子程序来实现，程序流程见图1。

图1 多场耦合程序系统程序结构Fig.1 Flowchart of the multifield coupling code

另外，在多场耦合的分析中，要求温度场的单元划分是与固体变形计算中的单元划分是一致的，如平面应变问题可统一取为4节点等参单元。

3.2 热-固耦合模型测试

为了验证文中解耦求解方法的正确性和有效性，采用厚壁圆筒热力耦合模型进行对比验证，计算模型如图2所示。厚壁筒弹性模量为300 MPa，泊松比为0.25，热传导系数为2.1 W/(m·ºC)，比热为 800 J/(kg·ºC)，密度为 2 300 kg/m3，热膨胀系数为6.0×10-51/ ºC，屈服强度为0.8 MPa。模型内半径为0.5 cm，外半径为5 cm，内壁温度为100 ℃，外壁温度为0 ℃。

图2 验证模型有限元网格Fig.2 Finite element meshes of verification model

分如下工况进行对比分析，工况1：热-弹塑性力学分析，采用 ABAQUS软件自带耦合方法。工况2：热-弹塑性力学分析，温度对材料力学参数的影响不计，仅考虑塑性应变对热传导系数和比热的影响，采用 ABAQUS软件自带耦合方法，并编制USDFLD子程序来实现热力学参数的动态演化。工况3：热-弹塑性力学分析，温度对材料力学参数的影响不计，仅考虑塑性应变对热传导系数和比热的影响，采用本文推荐的方法。

假设工况2和工况3中材料的热传导系数和比热均按如下公式动态演化：

计算结果如图 3、4所示。较工况 1和工况 2可以看出，变形对热力学参数的影响对计算结果有较大影响，在热力完全耦合模型中，应当考虑材料的变形对热力学参数（如热传导系数、比热等）的影响，但是，从目前的文献调研发现，关于变形对热力学参数影响的研究还不多见。比较工况2和工况3可以发现，采取本文推荐的方法时，在经过5次迭代后计算结果就达到了收敛，并且计算结果与ABAQUS自带的计算模块结果完全一致，说明本文方法的有效性和可靠性。

图3 工况1和工况2计算结果比较Fig.3 Comparison of results between cases 1 and 2

图4 工况2和工况3计算结果比较Fig.4 Comparison of results between cases 2 and 3

3.3 热-流-固耦合模型测试

以一维热结问题为例，进行了有限元数值解与解析解的对比[11]。

岩土体的力学、热学和水力基本参数：弹性模量为600 kPa，泊松比为0.3，土颗粒的体积模量为20 GPa，水的体积模量为5 GPa，土颗粒的热膨胀系数为1.5×10-5/℃，水热膨胀系数为2.0×10-4/℃，孔隙率为0.4，土颗粒的比热为0.8 J/(g⋅℃)，水的比热为4.2 J/(g⋅℃)，水的密度为1.0×106g/m3，土颗粒的密度为 2 600 kg/m3，热传导系数为 0.5 W/(m⋅℃)，Biot系数为1.0，土的渗透系数为1.954 7×10-9m/s。

计算模型如图5所示，土柱高1 m，在土柱的顶面向下施加载荷100 kPa，初始温度为10 ℃，初始孔隙水压力为100 kPa，顶面受60 ℃的温度载荷作用。模型的边界条件为：左右两侧为不透水边界，x方向位移约束；下边界为不透水边界，y方向固定；上表面为渗流自由边界。

图5 一维热固结计算模型Fig.5 One-dimensional thermal consolidation model

为研究的方便，定义时间因数T（无量纲）：

式中：K为岩土介质的热传导系数；t为时间（s）；h为土柱的高度；n为土的孔隙度；ρs、ρw分别为土颗粒和水的密度；cs、cw分别为土的比热。

计算结果如图 6～8所示。对于温度而言，有限元解与解析解完全吻合；对于节点孔隙压力和位移来说，有限元解与解析解大致吻合，而其中存在一定的差别是由于解析解是严格的一维固结问题，而有限元解是二维固结问题。

算例表明，本文采用解耦的数值方法处理热-力耦合问题时比较成功的，本程序成功地和现有的大型有限元软件模块链接，利用现有的成熟的软件资源，实现复杂的多场耦合计算。

4 模型应用

4.1 计算模型

图6 测点B温度变化曲线比较Fig.6 Comparison of temperatures between numerical and analytical solutions of point B

图7 测点A固结位移变化曲线比较Fig.7 Comparison of displacements between numerical and analytical solutions of point A

图8 测点B孔隙压力变化曲线比较Fig.8 Comparison of pore pressures between numerical and analytical solutions of point B

以准噶尔盆地某区块 3 000 m处泥岩地层为例，取准三维模型进行研究。井眼半径为0.108 m，上覆岩体平均密度为2 350 kg/m3，水平最大地应力为75 MPa，水平最小地应力为54 MPa，初始温度为120 ℃，孔隙压力为45 MPa。该地层的平均机械钻速为10 m/h，地层钻开后，井壁位移无约束，液柱压力为50 MPa，钻井液温度为70 ℃，泥饼质量较好，不考虑井壁的渗透性。计算参数见表1。

计算分为如下几个步骤：①初始地应力场平衡计算；②钻孔开挖卸载，将钻孔单元杀死；③热-流-固耦合计算。

表1 计算参数表Table 1 Main parameters for calculation

4.2 井壁失稳系数

基于Mohr-Coulomb准则和最大张应力理论，定义井壁坍塌系数fb和井壁破裂系数ff分别为

式中：c、φ、σt分别为岩层的黏聚力、内摩擦角和抗拉强度；σ1′、σ3′分别为最大和最小有效主应力。

当fb≤1时，井壁发生坍塌破坏，且fb越小井壁破坏越严重；当fb＞1时，井壁处于稳定状况。

5.3 结果分析与讨论

图9为井壁温度分布图。由图可见，井壁温度突然降低使得井壁围岩内温度也逐渐降低，温度沿水平径向的分布规律与垂直径向的分布规律基本类似，这是由于本文给出的地层热力学参数是常数，没有考虑应力对热力学参数的影响。图 10为孔隙压力分布图，钻孔卸载后，孔隙压力在水平向（最大主应力方向）迅速下降，此后孔隙压力梯度使得井眼孔隙压力逐渐增大；井眼钻开后，孔隙压力在垂向（最小主应力方向）迅速增加，此后渗流作用使得孔隙压力逐渐耗散。由于井壁温度降低了 50℃，孔隙压力产生了一个压力降（与恒温相比），24 h后，井壁孔隙压力与远场孔隙压力的压差约为 8 MPa。

图 11为井壁径向应力分布图。对于水平向位置，由于钻孔卸载，井壁附近径向应力迅速降低，此后，温度的降低引起围岩产生收缩，加上孔隙压力梯度引起的渗流作用，使得井壁继续径向应力继续降低；对于垂向位置，卸载后井壁径向应力迅速降低，此后，由于非均匀地应力、井壁围岩的收缩变形以及渗流作用的综合影响，径向应力逐渐增加，由于井壁内岩石的收缩而拉外圈的岩石，使得井壁内壁的应力小于其内部的应力。

图9 THM耦合条件下井壁温度分布示意图Fig.9 The temperature distributions with different times under THM coupling conditions

图10 THM耦合条件下井壁孔隙压力分布示意图Fig.10 The pore pressure distributions with different cases under THM coupling conditions

图11 THM耦合条件下井壁径向应力分布示意图Fig.11 The radial stress distributions with different cases

图12 THM耦合条件下井壁环向应力分布示意图Fig.12 The tangential stress distributions with different cases under THM coupling conditions

图 12为井壁环向应力分布图。对于水平向位置，井眼钻开后，距离井壁0.5rw处应力最大。由于井壁内围岩温度的降低引起，使得围岩产生收缩作用而使得内部应力有增大的趋势；对于垂向位置，卸载后围岩环向应力迅速增加，此后，井壁的收缩变形以及渗流作用使得环向应力继续增大。通过改变钻井液的温度值，研究温度波动对井壁稳定性的影响，具体计算结果见图13、14。

图13 地层受井眼内流体加热作用时井壁失稳系数比较Fig.13 Comparison of coefficients of wellbore stability under fluid heating condition

由图13可见，地层受井眼内流体加热作用时，井壁失稳系数比较图，温差越大，而井壁坍塌系数和破裂系数越小，表明井壁发生剪切破坏和张性破裂的可能性越大。由图14地层受井眼内流体流体冷却时，井壁 失稳系数比较图可以发现，温差越大，井壁坍塌系数和井壁破裂系数也越大，表明井壁发生剪切破坏和张性破裂的可能性越小。

5 结 语

在有关文献的基础上，利用有效应力原理、能量守恒定律和达西定律推导出多孔介质温度场、渗流场和应力场耦合的数学模型及其控制方程，以MATLAB语言为平台，将ABAQUS程序作为一个模块嵌入迭代算法程序中，提出了解决热-水-力耦合问题的解耦计算方法，编制了多场耦合有限元程序，并给出了算例验证，研究了油气钻井过程中的热-流-固耦合作用过程，详细分析了温度波动对井壁稳定性的影响。结果表明，在商业软件ABAQUS高效的内核求解器的辅助下，可以实现大规模多场耦合计算，本文提出的方法能解决热-流-固耦合问题的复杂程度，完全取决于渗流-应力耦合模块，可应用于多场耦合条件下岩体黏-弹塑性分析以及损伤分析。

图14 地层受井眼内流体冷却作用时井壁失稳系数比较Fig.14 Comparison of coefficients of wellbore stability under fluid cooling condition

[1] NOORISHAD J.TSANG C F.WITHERSPOON P A.Coupled thermal-hydraulic-mechanical phenomena in saturated fractured porous rocks: Numerical approach[J].Journal of Geophysical Research, 1984, 89: 10365－10373.

[2] 陈卫忠, 谭贤君, 伍国军, 等.非饱和岩石温度-渗流-应力耦合力学模型研究[J].岩石力学与工程学报, 2007,26(12): 2395－2403.CHEN Wei-zhong, TAN Xian-jun, WU Guo-jun, et al.A thermo-hydro-mechanical coupling model for unsaturated porous rock[J].Chinese Journal of Rock Mechanics and Engineering, 2007, 26(12): 2395－2403.

[3] 周创兵, 陈益峰, 姜清辉, 等.复杂岩体多场广义耦合分析导论[M].北京: 中国水利水电出版社, 2008.

[4] RUTQVIST J, WU Y S, TSANG C F, et al.A modelingapproach for analysis of coupled multiphase fluid flow,heat transfer, and deformation in fractured porous rock[J].International Journal of Rock Mechanics & Mining Sciences, 2002, 39: 429－442.

[5] 陈波, 李宁, 糕瑞花.多孔介质的变形场-渗流场-温度场耦合有限元分析[J].岩石力学与工程学报, 2001,20(4): 467－472.CHEN Bo, LI Ning, ZHUO Rui-hua.Finite element analysis of fully coupled thermo-hydro-mechanic behavior of porous media[J].Chinese Journal of Rock Mechanics and Engineering, 2001, 20(4): 467－472.

[6] 张玉军.热-水-应力耦合弹塑性二维有限元程序的开发与应用尝试[J].岩土力学, 2005, 26(2): 169－174.ZHANG Yu-jun.Development and use try of 2D FEM code for coupled thermo-hydro-mechanical elastoplastic analysis[J].Rock and Soil Mechanics, 2005, 26(2): 169－174.

[7] 冯夏庭, 潘鹏志, 丁梧秀, 等.结晶岩开挖损伤区的温度-水流-应力-化学耦合研究[J].岩石力学与工程学报,2008, 27(4): 656－663.FENG Xia-ting, PAN Peng-zhi, DING Wu-xiu, et al.Thermo-hydro-mechano-chemical coupling studies on excavation damage zone in crystalline rocks[J].Chinese Journal of Rock Mechanics and Engineering, 2008,27(4): 656－663.

[8] 陈益峰, 周创兵, 童富果, 等.多相流传输 THM 全耦合数值模型及程序验证[J].岩石力学与工程学报, 2009,28(4): 649－665.CHEN Yi-feng, ZHOU Chuang-bing, TONG Fu-guo, et al.A numerical model for fully coupled THM processes with multiphase flow and code validation[J].Chinese Journal of Rock Mechanics and Engineering, 2009, 28(4): 649－665.

[9] 傅少君, 吴秋军, 黄振科.考虑温度效应的固结问题有限元分析[J].土木工程学报, 2009, 42(1): 95－100.FU Shao-jun, WU Qiu-jun, HUANG Zhen-ke.Finite element analysis of thermal consolidation[J].China Civil Engineering Journal, 2009, 42(1): 95－100.

[10] 赵延林, 王卫军, 赵阳升, 等.双重介质热-水-力三维耦合模型及应用[J].中国矿业大学学报, 2010, 39(5):709－715.ZHAO Yan-lin, WANG Wei-jun, ZHAO Yang-sheng, et al.3D dual medium model of thermal-hydro-mechanical coupling and its application[J].Journal of China University of Mining & Technology, 2010, 39(5): 709－715.

[11] 白冰.岩土颗粒介质非等温一维热固结特性研究[J].工程力学, 2005, 22(5): 186－191.BAI Bing.One-dimensional thermal consolidation scharacteristics of geotechnical media under nonisothermal condition[J].Engineering Mechanics, 2005,22(5): 186－191.