热流固耦合作用下干热岩地热井套管‒水泥环胶结面应力场分析

Stress Field Analysis of Casing‒Cement Sheath Interface in Dry Hot Rock Geothermal Well Based on Thermal‒Hydraulic‒Mechanical Coupling

摘要

干热岩是一种新型清洁能源，对其开发利用，能够有效减少碳排放和缓解化石能源短缺问题。在所有井筒系统结构中，套管‒水泥环胶结面是较为薄弱的结构，在热流固耦合作用下，套管‒水泥环胶结面破裂机理尚不清楚。本文考虑套地层孔隙水压力、套管内注入水与地层之间温差所引起的热载荷、孔隙流体的热渗效应（T‒O），建立干热岩地热井套管‒胶结面‒水泥环‒地层组合体热流固耦合模型。通过拉普拉斯变换和逆变换，对模型进行求解，并与已公开发表的温度场数据和径向应力场数据进行比对验证。以青海共和盆地恰卜恰地区干热岩井为工程背景，得到了井筒系统内各界面随时间发展的温度曲线和径向应力曲线，发现套管‒水泥环胶结面处的温度在2.3 d左右降低到套管内壁温度，而径向应力开始有短暂的降低后增加，径向应力最小值为6.31 MPa。进一步分析胶结面的热流固耦合参数对干热岩地热井套管‒胶结面界面径向应力的影响规律，结果表明：当胶结面的弹性模量、泊松比、T‒O系数与水泥环的相应属性相近时，或者较低的地层温度和较小的胶结面孔隙率，均能减小套管‒胶结面界面处的径向应力，因而更有利于套管‒胶结面‒水泥环‒地层组合体的完整性。该结论对干热岩地热井保护井筒的完整性具有重要的指导意义。

Abstract

Objective The development and utilization of dry hot rock resources help China achieve its dual carbon goals. The reliability of the wellbore system is crucial for the development of dry hot rock resources. However, the casing cement sheath bonding surface in wellbore system components is a relatively weak structure, and when the wellbore is subjected to thermal‒hydraulic‒mechanical (THM) coupling, the casing-cement sheath interface (CCSI) is prone to breakage, reducing the reliability of the wellbore system. This study aims to understand the fracture mechanism of the casing cement sheath bonding surface. Methods The research method of this study was theoretical research and numerical calculation. Using the superposition principle, the boundary conditions of the model were decomposed into axisymmetric and non-axisymmetric boundary conditions. The model under non-axisymmetric boundary conditions obtained displacement and stress field expressions through the elastic thick-wall theory. Under axisymmetric boundary conditions, the model considered the coupling effect of the temperature field, mechanical field, and fluid flow field, and the expressions of each field were functions of position and time. Using the Laplace transform and inverse transform, the solution of the model was obtained. The inverse Laplace transform employed the Stehfest algorithm, and MATLAB software was used for numerical analysis. Results and Discussions The undetermined coefficients of all fields were time functions. Except for λ_E1 = 0, the temperature field coefficients of casing, CCSI, cement sheath, and formation all changed logarithmically with time. The coefficients λ_C1, λ_C2, and λ_D2 increased with time, whereas the coefficients λ_D1 and λ_E2 decreased with time. The coefficients of the mechanical field and fluid flow field of casing, CCSI, cement sheath, and formation followed a logarithmic change law, except for λ_C3 and λ_D3, which followed a linear change law. λ_F1, λ_F2, λ_D3, λ_D4, λ_D5, λ_D6, λ_E3, λ_E4, λ_E5, and λ_E6 increased with time. The coefficients λ_C3, λ_C4, λ_C5, and λ_C6 decreased with time. Published temperature curves were selected to verify the temperature field of the THM coupling model in this study. The results were found to be consistent by comparing the temperature curves at three moments. Published radial stress and radial displacement curves were also selected to verify the mechanical field of the THM coupling model, and the results were similarly consistent. Therefore, the establishment and solution process of the combinatorial THM coupling model proposed in this study was accurate. Based on the engineering background of a dry hot rock well in the Chabcha area, Gonghe Basin, Qinghai Province, the temperature curve and radial stress curve of all surfaces in the wellbore system over time were studied. Due to the thin thickness of CCSI, the transformation law of the bonding surface between the casing and CCSI was essentially the same as that of the bonding surface between CCSI and the cement sheath. The temperature at CCSI decreased to the temperature of the casing inner wall in 2.3 days; however, the cement-formation interface reached the temperature of the casing inner wall in 6.9 days. The radial stress at CCSI decreases briefly at first and then increases, with the minimum value of 6.31 MPa occurring at 0.2 days. The curve exhibits an exponential change, while the radial stress curve of the cement sheath-formation interface shows a linear change. The curve of radial stress versus radius shows that the radial stress increases rapidly in the casing, and the growth rate slows after reaching CCSI. The radial stress at CCSI is larger, while the failure strength of CCSI is relatively low, making CCSI more prone to failure than other locations. Further analysis was conducted on the influence of CCSI THM coupling parameters on the radial stress at the bonding surface between casing and CCSI. The radial stress decreased initially and then increased with the increase of the ratio of elastic modulus between CCSI and cement sheath. When the ratio was 0.33 and 1.67, the radial stress was maximum, and when the ratio was 1.11, the radial stress was minimum. Changes in CCSI Poisson's ratio had little influence on radial stress during the initial period; however, after approximately 4.6 days, the influence became greater. With an increase in Poisson's ratio, the radial stress of CCSI first decreased and then increased. When the Poisson's ratio was 0.30, the radial stress was minimum, and when the Poisson's ratio was 0.34, the radial stress was maximum. With the increase in porosity, the growth rate of radial stress gradually increased. When the porosity was in the range of 0.1~0.3, the growth rate of radial stress was small, while for porosity greater than 0.4, the growth rate of radial stress was large. With the increase in the T‒O coefficient, the radial stress first decreased and then increased. When the T‒O coefficient was 1.0, the radial stress was minimum, and when the T‒O coefficient was 2.0, the radial stress was maximum. After a certain period, the growth rate of the radial stress remained unchanged. The radial stress increased with the increase in temperature difference. The greater the temperature difference, the greater the growth rate of radial stress. When the temperature difference exceeded 200 ℃, the growth rate of radial stress was significantly higher than when the temperature difference was below 200 ℃. The maximum tensile stress criterion was chosen for tension, and the M‒C criterion was selected for compression to assess CCSI damage. Conclusions When the elastic modulus, Poisson's ratio, and T‒O coefficient of CCSI are similar to the corresponding properties of the cement sheath, they are more effective in reducing the radial stress at the interface between the casing and CCSI and in protecting the integrity of the casing-interface-cement-formation system. Lower formation temperatures and lower interface porosity have the same effect in reducing CCSI damage. This study holds great significance for ensuring the integrity of wellbore protection in dry hot rock geothermal wells.

Graphical abstract

关键词

干热岩地热井 / 热流固耦合 / 套管‒水泥环胶结面 / 应力场 / 理论解

Key words

dry hot rock geothermal wells / thermal‒hydraulic‒mechanical coupling / casing-cement sheath interface / stress field / analytical solution

本刊网刊

干热岩，指温度大于200 ℃，埋深数千米，内部不存在流体或仅有少量地下流体的高温岩体^[1‒2]。干热岩资源开发能够有效减少碳排放和缓解化石能源短缺问题，积极促进中国“双碳”战略目标的完成^[3‒4]。

对干热岩资源的开发是通过水力压裂建立增强型地热系统（EGS）^[5‒7]，利用循环流体提取热量，达到地热发电和供暖综合利用的目标。由于干热岩主要岩性为变质岩和火成岩^[3]，且深层干热岩所处的高应力场环境导致干热岩破裂压力极高，因而水力压裂时泵入压力较高^[7]。井筒内高压状态的压裂液对干热岩地热井套管‒水泥环胶结面带来两方面的危害：套管内压急升引起的胶结面破裂和套管温度骤降带来的热冲击引起胶结面产生裂缝^[8‒9]。套管‒水泥环胶结面质量不仅决定着干热岩地热井井筒系统的稳定性，还影响到后续EGS系统的工程建设，是干热岩开发中的关键问题。

考虑到干热岩高温降低了固井水泥浆的热稳定性和水泥水化放热，干热岩地热井比常规油气井在水泥固化过程中更易产生孔洞、裂隙^[10‒11]等缺陷。由于套管和水泥浆之间存在壁效应（wall effect），导致水泥固化过程中水泥大颗粒无法运移至胶结面处，无法形成最优结构，故而套管‒水泥环胶结面包含一层薄弱带，宽度约10～50 μm，称为界面过渡区（ITZ）^[12]。界面过渡区细分包含3个细小区域：以细小晶体为主的区域宽度约4～20 μm，小颗粒贫化区宽度4～20 μm，过渡区宽度2～10 μm，示意图如图1所示。同时，套管残留钻井液^[13]及水泥水化放热引起的套管膨胀^[9]使得界面过渡区的缺陷数量相对较多^[14]。

对于胶结面强度的研究，罗长吉等^[15]设计了界面胶结强度测试仪和界面压窜试验仪，利用该装置进行了界面的胶结强度和抗压强度的规律研究，并且发现了界面胶结强度、水泥浆的收缩特性及水泥水化的热效应之间存在的内在关联。许永志等^[16]通过在钻井液中加入磺化聚合醇（GSN）物质，有效提高固井界面的胶结性能。齐奉忠等^[17]指出套管水泥环界面的胶结性能受到水泥环体积收缩而出现的微环空影响。Wang等^[18]建立了一个轴对称多孔弹性有限元模型，在模型中，地层和水泥岩石界面及套管水泥界面中嵌入了面内厚度为0的黏性界面单元。Nath等^[19]应用数字图像相关（DIC）技术，研究了套管‒水泥环胶结面的胶结强度；Goodwin等^[20]通过加卸载套管内压，开展了水泥环界面失效试验。Carpenter等^[21]搭建了胶结面测试装置，研究了不同类型水泥在温度变化下的胶结面破裂情况，发现当温度变化幅度超过56 ℃时胶结面将会发生严重的破裂。Bois等^[22]研究发现，温度和压力同时变化会加重套管‒水泥环胶结面的破裂。

胶结面的强度受水泥环材料的影响较大，针对水泥环为主体的强度设计研究，殷有泉^[23]和Shi^[24]等选用弹塑性理论，在考虑不同载荷工况下，推导了井筒系统的理论解答，并研究了套管‒水泥环接触面的应力、位移分布规律。胡志强等^[25]基于套管‒水泥环‒地层热固耦合作用，建立了多层套管环空附加压力矩 阵方程。赵新波等^[26‒27]建立了热固耦合下套管‒水泥环‒地层模型的理论模型。Wu等^[28]研究了高温井中的胶结面微环隙的位移机理。Luo等^[29]建立了一个简单的井筒传热模型，以减少井筒胶结面破裂为目标，给出了注入压力的安全区间。Nygaard等^[30]通过数值模拟发现热冲击能够降低近井区域的应力，导致胶结面拥有较大的拉应力，从而增加了胶结面破裂风险。

在先前的研究中，水泥环仍然被视为弹性连续介质，因此忽略了孔隙压力对水泥环力学行为的影响。在自然界中，水泥环是一种多孔结构，由几种类型的孔隙空间组成，例如层间空隙、凝胶孔隙、毛细孔和气孔，其尺寸范围从几纳米到数百微米不等。然而，水泥环内孔隙压力、温度的变化在水泥环的力学行为中起着至关重要的作用，但这在现有研究中并没有得到充分考虑。此外，热渗效应^[31]在描述液体在毛细管或多孔材料中受温度梯度驱动的流动时变得非常普遍，尤其在孔径较小的多孔材料中更为重要。基于数值模拟，Ghassemi等^[31]考虑了热渗效应来分析高盐度泥浆对井筒稳定性的影响，发现这种效应增加了井筒的潜在拉伸破坏。通过一些实验，可测量各种类型地层的热渗系数^[32]。然而，迄今为止很少有研究关注热渗效应对水泥环热‒孔隙‒力学行为的影响。

随着研究的深入，越来越多研究者意识到胶结面和水泥环一样，应当将其看作具有独立力学性质的薄弱层，即界面过渡区。Sevostianov等^[33]将界面过渡区看成一个弹性体，探讨了界面过渡区对复合材料整体性能的影响。李宗利等^[34]提出了一种预测混凝土弹性模量的多相混合夹杂模型，研究发现，混凝土的弹性模量随着界面过渡区厚度的增大而减小，随界面过渡区力学性能的增大而增大。Duan等^[35]基于能量等效原理，提出了考虑界面影响的多相复合材料弹性模量理论模型。邓方茜等^[36]基于细观力学方法，建立考虑界面作用的混杂纤维混凝土弹性模量预测模型。Lavrov等^[37]通过研究水泥水化作用发现，胶结面孔洞、裂隙较多，胶结面非均质性较强。Lavrov等^[13]还利用数值模拟方法研究了缺陷对胶结面力学性能的影响。Torsæter等^[38]研究了胶结面的材料和几何缺陷对套管‒水泥环胶结面的破裂影响规律。

综合上述研究成果发现，压裂液带来的井筒温压变化引起干热岩地热井套管‒水泥环胶结面破裂涉及流体渗流、材料变形、热量传输，但现有研究尚未综合考虑热流固耦合对胶结面破裂的影响。因此，本文综合考虑热流固耦合作用，从理论角度开展干热岩地热井套管‒水泥环胶结面应力场分析，为胶结面破裂分析提供支撑。

1 热流固耦合模型建立

套管‒水泥环之间的界面过渡区称之为胶结面（图1），为单独一层结构，具有独立的材料参数和几何参数。本模型与以往发表的模型不同，在结构上引入胶结面，建立了套管‒胶结面‒水泥环‒地层组合系统。

假设干热岩地热井套管‒胶结面‒水泥环‒地层组合系统无纵向变形，即该组合系统为平面应变问题。根据套管‒胶结面‒水泥环‒地层组合系统的实际服役条件，建立干热岩地热井套管‒胶结面‒水泥环‒地层组合系统热流固耦合模型，如图2所示。

图2中，{\sigma }_{\mathrm{H}}、{\sigma }_{\mathrm{h}}分别为地层水平最大、最小主应力，r_{\mathrm{0}}为套管内半径，r_{\mathrm{1}}为外半径，r_{\mathrm{2}}为胶结面外半径，r_{\mathrm{3}}为水泥环外半径、r_{\mathrm{4}}为地层外半径，T₀为套管内壁温度，P₀为套管内部压力，T_f为无穷远处地层温度，P_f为孔隙压力。

基于传统套管‒水泥环‒地层模型，结合应力张量第一不变量、胡克定律、热应力和孔隙压力引起体积应变等，建立套管和水泥环的本构方程，分别为式（1）和（2）：

\begin{array}{l}{\sigma }_{ij}^{\mathrm{c}\mathrm{a}\mathrm{s}}=\mathrm{2}{G}^{\mathrm{c}\mathrm{a}\mathrm{s}}[{\epsilon }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}\beta \mathrm{\Delta }T{\delta }_{ij}+\\ \frac{{\mu }^{\mathrm{c}\mathrm{a}\mathrm{s}}}{\mathrm{1}-\mathrm{2}{\mu }^{\mathrm{c}\mathrm{a}\mathrm{s}}}(\Phi -\mathrm{\Delta }T\beta ){\delta }_{ij}]\end{array}（1）

\begin{array}{l}{\sigma }_{ij}^{\mathrm{c}\mathrm{e}\mathrm{m}}=\mathrm{2}{G}^{\mathrm{c}\mathrm{e}\mathrm{m}}\left({\epsilon }_{ij}^{}+\frac{{\mu }^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\mathrm{1}-\mathrm{2}{\mu }^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\delta }_{ij}\Phi \right)-\\ {\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}{\delta }_{ij}p-{\omega }_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\delta }_{ij}\mathrm{\Delta }T\end{array}（2）

式（1）～（2）中：{\sigma }_{ij}^{\mathrm{c}\mathrm{a}\mathrm{s}}为套管应力2阶张量，i，j为循环变量，i,j=1,2,3，上标cas代表套管；{\epsilon }_{ij}为应变张量；G^{\mathrm{c}\mathrm{a}\mathrm{s}}为套管剪切弹性模量；{\mu }^{\mathrm{c}\mathrm{a}\mathrm{s}}为套管泊松比；\beta为体积热膨胀系数；\mathrm{\Delta }T为温度改变量；{\delta }_{ij}为Kronecker符号；\Phi为体积应变；{\sigma }_{ij}^{\mathrm{c}\mathrm{e}\mathrm{m}}、G^{\mathrm{c}\mathrm{e}\mathrm{m}}、{\mu }^{\mathrm{c}\mathrm{e}\mathrm{m}}分别为水泥环应力张量、剪切弹性模量、泊松比，上标cem代表水泥环；{\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}为水泥环Biot系数；{\omega }_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}为水泥环排水热弹性有效应力系数。套管不受渗流场影响且温度传输为各向同性。

胶结面，水泥环，地层均看成多孔材料，因此胶结面、地层本构方程参照水泥环本构方程建立。

由于系统受到温度场和渗流场影响，根据热应力影响和孔隙水体积改变得到状态方程：

\left\{\begin{array}{l}{s}_{\mathrm{h}}=m_{\mathrm{d}}^y\mathrm{\Delta }T-{\beta }_{\mathrm{e}}^{y}p+{\omega }_{\mathrm{d}}^y\Phi ,\\ \zeta ={\alpha }^y\Phi +\frac{p{K}_{\mathrm{b}}^y}{{\alpha }^y-{\varphi }_{\mathrm{0}}^y}-{\beta }_{\mathrm{e}}^y\mathrm{\Delta }T\end{array}\right.（3）

式中：s_{\mathrm{h}}为多孔材料的总熵密度；m_{\mathrm{d}}^y为恒定应变和初始温度下排出比热的熵；{\beta }_{\mathrm{e}}^y为体积热膨胀系数， {\beta }_{\mathrm{e}}^y=\alpha {\beta }_{\mathrm{s}}^y+{\varphi }_{\mathrm{0}}^y({\beta }_{\mathrm{f}}-{\beta }_{\mathrm{s}}^y)，其中，{\beta }_{\mathrm{s}}^y和{\beta }_{\mathrm{f}}分别表示多孔材料固相和孔隙水的体积热膨胀系数，上标y为材料类型，材料类型有水泥环、胶结面、地层，并用cem、pcem、for进行区分；{\varphi }_{\mathrm{0}}^y为孔隙率；{\omega }_{\mathrm{d}}^y为排水热弹性有效应力系数；\zeta为孔隙水体积变化；K_{\mathrm{b}}^y为体积模量。

热流固耦合控制方程为：

\left\{\begin{array}{l}{\sigma }_{ij,j}+f_i=\mathrm{0},\\ \frac{\partial \zeta }{\partial t}+\nabla {\mathit{q}}_{\mathrm{f}}\approx \mathrm{0},\\ T_{\mathrm{0}}\frac{\partial s_{\mathrm{h}}}{\partial t}+\nabla {\mathit{q}}_{\mathrm{h}}=\mathrm{0}\end{array}\right.（4）

式中，f_i为体积力，{\mathit{q}}_{\mathrm{f}}为孔隙水通量向量，{\mathit{q}}_{\mathrm{h}}为均匀热通量，t为时间，\nabla为哈密顿算子。

多孔材料中热通量和流体通量的表达式为：

\left\{\begin{array}{l}{\mathit{q}}_{\mathrm{f}}=-{\kappa }^y\nabla p-k_{\mathrm{p}\mathrm{T}}^y\nabla (\mathrm{\Delta }T),\\ {\mathit{q}}_{\mathrm{h}}=-k_{\mathrm{T}\mathrm{p}}^y\nabla p-k_{\mathrm{T}}^y\nabla (\mathrm{\Delta }T)\end{array}\right.（5）

式中：{\kappa }^y为多孔材料的渗透率；k_{\mathrm{p}\mathrm{T}}^y为温度梯度对热传导影响的热渗透系数（T‒O）流体扩散率，k_{\mathrm{p}\mathrm{T}}={k{\epsilon }_{\mathrm{T}}/\mu }_{\mathrm{p}}，{\mu }_{\mathrm{p}}为孔隙水的动态黏度；k_{\mathrm{T}\mathrm{p}}^y为热力系数，而k_{\mathrm{T}\mathrm{p}}^y\nabla p引起的变化很小可以忽略，p为压力；k_{\mathrm{T}}^y为多孔材料的有效热导率；\mathrm{\Delta }T为温度变化量。

热流固耦合模型（图2）需将直角坐标变换为极坐标进行应力场、温度场、渗流场的求解。选取极坐标下的模型外边界（r₄）足够远，不考虑外边界的温度变化^[38]和孔隙压力的变化，且沿环向1周认为温度和孔隙压力没有变化。利用坐标变换，组合系统的应力内外边界条件可以描述为：

\left\{\begin{array}{l}{\sigma }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}={\sigma }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}+{\sigma }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}=\\ -\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{h}}}{\mathrm{2}}-\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{h}}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\theta \right),\\ {\tau }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}={\tau }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}+{\tau }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}=\\ \mathrm{0}+\frac{{\sigma }_{\mathrm{H}}-{\sigma }_{\mathrm{h}}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\theta \right),\\ {\sigma }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}={\sigma }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}+{\sigma }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}=-p_{\mathrm{0}}+\mathrm{0},\\ {\tau }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}={\tau }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}+{\tau }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}=\mathrm{0}+\mathrm{0}=\mathrm{0}\end{array}\right.（6）

式中，r为径向半径，θ为圆周角，{\sigma }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}、{\tau }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}分别为地层外边界r₄的平均径向应力和切应力，{\sigma }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}、{\tau }_r^{\mathrm{f}\mathrm{o}\mathrm{r}}{|}_{r=r_{\mathrm{4}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}分别为地层外边界r₄的偏差径向应力和切应力，{\sigma }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}、{\tau }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{a}\mathrm{v}\mathrm{e}}分别为套管内边界r₀的平均径向应力和切应力，{\sigma }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}、{\tau }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}{|}_{r=r_{\mathrm{0}}}^{\mathrm{d}\mathrm{e}\mathrm{v}}分别为套管内边界r₀的偏差径向应力和切应力，上标ave、dev分别代表平均和偏差。

2 偏差应力场求解

对于偏差应力场，本文仅研究其对结构场的影响。

根据弹性力学基本理论，假设每一层材料满足应力函数\vartheta \left(r,\theta \right)：

\vartheta \left(r,\theta \right)=\left(K_{\mathrm{A}}{r}^{\mathrm{4}}+K_{\mathrm{B}}{r}^{\mathrm{2}}+\frac{K_{\mathrm{C}}}{r^{\mathrm{2}}}+K_{\mathrm{D}}\right)\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\theta )（7）

式中，K_A、K_B、K_C、K_D为待定常数。代入应力函数表达式，得到：

\left\{\begin{array}{l}{\sigma }_r^{}=-\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\theta )(\mathrm{2}{K}_{\mathrm{B}}+\frac{\mathrm{4}{K}_{\mathrm{C}}}{r^{\mathrm{2}}}+\frac{\mathrm{6}{K}_{\mathrm{D}}}{r^{\mathrm{4}}}),\\ {\sigma }_{\theta }^{}=\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\theta )(\mathrm{12}{K}_{\mathrm{A}}{r}^{\mathrm{2}}+\mathrm{2}{K}_{\mathrm{B}}+\frac{\mathrm{6}{K}_{\mathrm{D}}}{r^{\mathrm{4}}}),\\ {\tau }_{r\theta }^{}=\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\theta )(\mathrm{6}{K}_{\mathrm{A}}{r}^{\mathrm{2}}+\mathrm{2}{K}_{\mathrm{B}}-\frac{\mathrm{2}{K}_{\mathrm{C}}}{r^{\mathrm{2}}}-\frac{\mathrm{6}{K}_{\mathrm{D}}}{r^{\mathrm{4}}})\end{array}\right.（8）

式中，{\sigma }_r^{}、{\sigma }_{\theta }^{}、{\tau }_{r\theta }^{}分别为径向应力、环向应力和切应力。

通过物理方程和几何方式，得到径向位移U_r^{}、环向位移U_{\theta }^{}，表达式为：

U_r^{}=((\mathrm{18}-\mathrm{6}\varsigma )K_{\mathrm{A}}{r}^{\mathrm{3}}-\varsigma K_{\mathrm{B}}r+\frac{(\mathrm{2}+\mathrm{4}\varsigma )K_{\mathrm{C}}}{r}+\frac{\mathrm{9}\varsigma K_{\mathrm{D}}}{r^{\mathrm{3}}})\frac{\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{2}\theta \right)}{\mathrm{5}E}（9）

\begin{array}{l}{U}_{\theta }^{}=((\mathrm{2}+\varsigma )K_{\mathrm{A}}{r}^{\mathrm{3}}+\varsigma K_{\mathrm{B}}r-\frac{(\mathrm{2}-\varsigma )K_{\mathrm{C}}}{r}+\\ \frac{\varsigma K_{\mathrm{D}}}{r^{\mathrm{3}}})\frac{\mathrm{6}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\theta \right)}{\mathrm{5}E}\end{array}（10）

式（9）～（10）中：E为弹性模量；\varsigma为中间变量，\varsigma=1+μ，μ为泊松比。

由于组合系统为多层材料，每一层材料的应力场表达式和位移场表达式形式相同，因此只需要将待定常数、材料参数等加下标ξ加以区别，比如K_A变成K_Aξ，E变成E_ξ，ξ的取值域为{1,2,3,4}，分别代表套管、胶结面、水泥环、地层。

3 热流固耦合求解

3.1 热流固耦合模型拉式域求解

套管材料不受渗流场影响，因此套管控制方程可写成：

\frac{\partial {\epsilon }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}}{\partial r}-\frac{\mathrm{1}+{\mu }^{\mathrm{c}\mathrm{a}\mathrm{s}}}{\mathrm{3}(\mathrm{1}-{\mu }^{\mathrm{c}\mathrm{a}\mathrm{s}})}{\beta }^{\mathrm{c}\mathrm{a}\mathrm{s}}\frac{\partial \mathrm{\Delta }T}{\partial r}=\mathrm{0}（11）

式中，{\epsilon }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}为套管径向应变，{\beta }^{\mathrm{c}\mathrm{a}\mathrm{s}}为套管体积热膨胀系数。考虑到套管为薄壁结构，且金属材料导热率较大，因此，可以假设温度在套管内均匀分布，故将式（11）中第2项省略，并将其进行拉普拉斯变换，得到径向应变在拉式域内的通解为：

{{\epsilon }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}}^{°}=F_{\mathrm{1}}（12）

式中，{{\epsilon }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}}^{°}为{\epsilon }_r^{\mathrm{c}\mathrm{a}\mathrm{s}}在拉氏域内表达形式，F₁为常数。根据轴对称模型的几何条件，径向位移在拉式域的解为：

u_r^{\mathrm{c}\mathrm{a}\mathrm{s}°}=\frac{r}{\mathrm{2}}{F}_{\mathrm{1}}+\frac{F_{\mathrm{2}}}{r}（13）

式中，u_r^{\mathrm{c}\mathrm{a}\mathrm{s}°}为套管径向位移在拉氏域内表达形式，F₂为常数。由于胶结面、地层细观结构和水泥环相似，因此，胶结面和地层的温度场、渗流场、结构场解的分布形式与水泥环类似。水泥环模型为轴对称模型，且受到温度场和渗流场双重影响。

根据式（4）第1式得到水泥环结构场控制方程为：

\mathrm{2}{G}^{\mathrm{c}\mathrm{e}\mathrm{m}}\frac{\mathrm{1}-{\mu }^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\mathrm{1}-\mathrm{2}{\mu }^{\mathrm{c}\mathrm{e}\mathrm{m}}}\nabla \Phi -{\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}\nabla p-{\omega }_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}\nabla (\mathrm{\Delta }T)=\mathrm{0}（14）

将式（3）和（5）带入式（4）第2和3式，可得温度场和渗流场控制方程为：

\left\{\begin{array}{l}\frac{\partial {\left(\mathrm{\Delta }T\right)}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-c_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}{(\mathrm{\Delta }T)}^{\mathrm{c}\mathrm{e}\mathrm{m}}=\\ \frac{{\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}}\frac{\partial p^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-\frac{{\omega }_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}\mathrm{d}{g}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}\mathrm{d}t},\\ \frac{\partial p^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-c_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}p=\frac{{\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}}\frac{\partial {\left(\mathrm{\Delta }T\right)}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}+\\ \frac{k_{\mathrm{p}\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\nabla }^{\mathrm{2}}{\left(\mathrm{\Delta }T\right)}^{\mathrm{c}\mathrm{e}\mathrm{m}}-\frac{{\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}}\frac{\mathrm{d}{g}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\mathrm{d}t}\end{array}\right.（15）

式中：{\left(\mathrm{\Delta }T\right)}^{\mathrm{c}\mathrm{e}\mathrm{m}}为水泥环温度场变量；{\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}代表水泥环恒定框架体积下孔隙水含量变化的体积热膨胀系数；S_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}为中间变量，其表达式为S_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}=\mathrm{3}{\omega }_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/(\mathrm{3}{K}_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}+\mathrm{4}{G}^{\mathrm{c}\mathrm{e}\mathrm{m}})+m_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}，K_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}代表水泥环固有渗透率，K为固有渗透率；c_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}为中间变量，c_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}=k_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/(S_{\mathrm{a}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\mathrm{T}}_{\mathrm{0}})；g^cem为水泥环与应变、温度、应力等相关的量；p^{\mathrm{c}\mathrm{e}\mathrm{m}}为水泥环流场压力；S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}为中间变量，S^{\mathrm{c}\mathrm{e}\mathrm{m}}=\mathrm{3}({\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}{)}^{\mathrm{2}}/(\mathrm{3}{K}_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}+\mathrm{4}{G}^{\mathrm{c}\mathrm{e}\mathrm{m}})+K_{\mathrm{b}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/({\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}-{\varphi }_{\mathrm{0}}^{\mathrm{c}\mathrm{e}\mathrm{m}})；c_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}为中间变量，c_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}=K_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}/(S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\mu }_p)，{\mu }_p为孔隙水的动粘度系数。

根据Zimmerman所提出来的理论，如果当油井水泥的热弹性耦合项K^{\mathrm{c}\mathrm{e}\mathrm{m}}({\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{)}^{\mathrm{2}}{T}_{\mathrm{0}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/c_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}\le \mathrm{1}，c_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}为水泥环恒定应变下排出的比热，式（15）可进一步简化为：

\left\{\begin{array}{l}\frac{\partial {(\mathrm{\Delta }T)}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-{\kappa }_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}{(\mathrm{\Delta }T)}^{\mathrm{c}\mathrm{e}\mathrm{m}}=\mathrm{0},\\ \frac{\partial p^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-c_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}{p}^{\mathrm{c}\mathrm{e}\mathrm{m}}=\frac{{\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}+\left(k_{pT}^{cem}/{\kappa }_T^{cem}\right)}{S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}}·\\ \frac{\partial {(\mathrm{\Delta }T)}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\partial t}-\frac{{\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S_{}^{\mathrm{c}\mathrm{e}\mathrm{m}}}·\frac{\mathrm{d}{g}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{\mathrm{d}t}\end{array}\right.（16）

式中：{\kappa }_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}为热扩散系数，{\kappa }_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}=k_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/c_{\mathrm{d}}^{\mathrm{c}\mathrm{e}\mathrm{m}}。

式（16）是位置和时间耦合的方程，直接求解比较困难，为了便于求解，可将式（16）进行拉普拉斯变换转换，在拉氏域内求解后，再用拉普拉斯逆变得到时间域内的解。式（16）通过拉氏变换得到：

\left\{\begin{array}{l}s\mathrm{\Delta }{T}^{°}-{\kappa }_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}(\mathrm{\Delta }{T}^{\circ })=\mathrm{0},\\ s{p}^{\circ }-c^{\mathrm{c}\mathrm{e}\mathrm{m}}{\nabla }^{\mathrm{2}}{p}^{\circ }=s{c}_{\mathrm{T}}\mathrm{\Delta }{T}^{\circ }-s{c}_{\mathrm{g}}{g}^{\circ }\end{array}\right.（17）

式中：\mathrm{\Delta }{T}^{°}为{\left(\mathrm{\Delta }T\right)}^{\mathrm{c}\mathrm{e}\mathrm{m}}拉氏域下的表达式；p^{\circ }为p^{\mathrm{c}\mathrm{e}\mathrm{m}}拉氏域下的表达式；g^{\circ }是g^cem拉氏域下的表达式；s为时间t转换到拉氏域内的对应变量；c_T、c_g为中间变量，c_{\mathrm{T}}=\frac{{\beta }_{\mathrm{c}}^{\mathrm{c}\mathrm{e}\mathrm{m}}+k_{\mathrm{p}\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}/{\kappa }_{\mathrm{T}}^{\mathrm{c}\mathrm{e}\mathrm{m}}}{S^{\mathrm{c}\mathrm{e}\mathrm{m}}}，c_{\mathrm{g}}={\alpha }^{\mathrm{c}\mathrm{e}\mathrm{m}}/S^{\mathrm{c}\mathrm{e}\mathrm{m}}。

式（17）第1式仅为温度变化的方程，因此可以直接求解。利用数学物理方程中的修正贝塞尔方程，可得到变温\mathrm{\Delta }{T}^{°}的表达式：

\mathrm{\Delta }{T}^{°}={\lambda }_{\mathrm{1}}{I}_{\mathrm{0}}({\zeta }_T)+{\lambda }_{\mathrm{2}}{K}_{\mathrm{0}}({\zeta }_T)（18）

式中：{\zeta }_T为中间过程量，{\zeta }_T=\sqrt[]{s/{\kappa }_T}r；T为温度；{\lambda }_{\mathrm{1}}和{\lambda }_{\mathrm{2}}为待定常数；I_{\mathrm{0}}({\zeta }_T)和K_{\mathrm{0}}({\zeta }_T)为{\zeta }_T的函数，其表达式分别为：

I_{\mathrm{0}}({\zeta }_T)={\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{\infty }}}\frac{{{\zeta }_T}^{\mathrm{2}m}}{{\mathrm{2}}^{\mathrm{2}m}m!\mathrm{\Gamma }(m+\mathrm{1})}（19）

K_{\mathrm{0}}({\zeta }_T)=\underset{\theta \to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{\pi }}{\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \mathrm{\pi }\right)}(I_{-\theta }({\zeta }_T)-I_{\theta }({\zeta }_T))（20）

式（19）～（20）中：m为循环变量；I_{\theta }({\zeta }_T)为m取值范围[θ,∞]；I_{-\theta }({\zeta }_T)为m取值范围[‒∞,‒θ]；\mathrm{\Gamma }(m+\mathrm{1})表示在复数范围内的亚纯函数，\mathrm{\Gamma }(x)表达式为：

\mathrm{\Gamma }(x)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}{y}^{x-\mathrm{1}}{e}^{-y}\mathrm{d}y,x>\mathrm{0}（21）

式（17）第2式为孔隙压力和温度耦合方程，通过方程结构研究，设孔隙压力p^{\circ \mathrm{*}}的解为：p^{\circ \mathrm{*}}=a\mathrm{\Delta }{T}^{\circ }+b，其中，a、b为中间过程量，a=c_T/(\mathrm{1}-c/{\kappa }_T)，b=-c_g{g}^{\circ }，进而得到：

p^{\circ }=p^{\circ \mathrm{*}}+{\lambda }_{\mathrm{3}}{I}_{\mathrm{0}}({\zeta }_p)+{\lambda }_{\mathrm{4}}{K}_{\mathrm{0}}({\zeta }_p)（22）

式中：{\zeta }_p为中间过程量，{\zeta }_p=\sqrt[]{s/c^{\mathrm{c}\mathrm{e}\mathrm{m}}}r，下标p为流体压力；{\lambda }_{\mathrm{3}}和{\lambda }_{\mathrm{4}}为待定常数。

胶结面和地层的各场解可参考水泥环求解过程。

3.2 热流固耦合模型时间域求解

本节把上述拉式域内的解，通过拉普拉斯逆变换，并借助边界条件及连续性条件，转变成时间域内的解。

在进行拉普拉斯逆变换时，需要对模型进行如下假设：1）套管、胶结面和水泥环同心；2）在外载荷作用下，套管‒胶结面，胶结面‒水泥环和水泥环‒地层界面未发生脱黏；3）地层中的初始地应力各向同性，套管‒胶结面‒水泥环‒地层模型为轴对称模型；4）忽略远场地层的温度、压力变化。

热流固耦合模型的边界条件有：

\left\{\begin{array}{l}{\sigma }_{rr}^{\mathrm{c}\mathrm{a}\mathrm{s}}(r_{\mathrm{0}})=-p_{\mathrm{0}},\\ T_{{}_{}}^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}(r_{\mathrm{1}})=\mathrm{\Delta }{T}_{\mathrm{0}},\\ q_r^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}(r_{\mathrm{1}})=\mathrm{0},\\ \mathrm{\Delta }T(r_{\mathrm{4}})=\mathrm{0},\\ {\sigma }_{rr}^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})=-\frac{{\sigma }_{\mathrm{H}}+{\sigma }_{\mathrm{h}}}{\mathrm{2}},\\ p_{}^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})=p_{\mathrm{0}}^{},\\ u_r^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})=\mathrm{0}\end{array}\right.（23）

式中：{\sigma }_{rr}^{\mathrm{c}\mathrm{a}\mathrm{s}}(r_{\mathrm{0}})为r_{\mathrm{0}}位置处的套管径向位移；T_{{}_{}}^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}(r_{\mathrm{1}})和q_r^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}(r_{\mathrm{1}})分别为r_{\mathrm{1}}位置处胶结面的温度和流量；\mathrm{\Delta }{T}_{\mathrm{0}}为注入水温度T_{\mathrm{0}}^{\mathrm{W}}和地下水温度T_{\mathrm{0}}^{\mathrm{f}}之间的温差，\mathrm{\Delta }{T}_{\mathrm{0}}=T_{\mathrm{0}}^{\mathrm{W}}-T_{\mathrm{0}}^{\mathrm{f}}；\mathrm{\Delta }T(r_{\mathrm{4}})、{\sigma }_{rr}^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})、p_{}^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})和u_r^{\mathrm{f}\mathrm{o}\mathrm{r}}(r_{\mathrm{4}})分别为r_{\mathrm{4}}位置处的地层温差、径向应力、流体压力和径向位移。

热流固耦合模型的连续性条件有：在两种结构的胶结面上有径向应力，温度，径向位移，孔隙流体压力，热通量相等。

由连续性条件和边界条件，并结合Stehfest算法进行拉普拉斯逆变换，有：

f(t)\approx \frac{\mathrm{l}\mathrm{n}\mathrm{2}}{t}{\displaystyle \sum _{n=\mathrm{1}}^M}{S}_{n}F(n\frac{\mathrm{l}\mathrm{n}\mathrm{2}}{t})（24）

式中：M为给定系数且为偶数，取值范围为\left[\mathrm{6},\mathrm{20}\right]，本文M取值为8；n为正整数；f(t)和F(·)分别表示时域中的函数及其在拉普拉斯域中的对应函数；S_n为中间变量，具体表达式如下：

S_n={\left(-\mathrm{1}\right)}^{n+M/\mathrm{2}}{\displaystyle \sum _{k=(n+\mathrm{1})/\mathrm{2}}^{\mathrm{m}\mathrm{i}\mathrm{n}(n,M/\mathrm{2})}}\frac{k^{M/\mathrm{2}}(\mathrm{2}k)!}{(M/\mathrm{2}-k)!k!(k-\mathrm{1})!(n-k)!(\mathrm{2}k-n)!}（25）

式中，k为正整数。

胶结面，水泥环，地层温度场表达式均可按照式（18）的分布形式，每种结构有两个待定系数。λ_C1、λ_C2为胶结面温度场待定常数；λ_D1、λ_D2为水泥环温度场待定常数；λ_E1、λ_E2为地层温度场的待定系数。由求解基本假设可知，地层远场温度无变化，通过求解可得λ_E1=0，{\lambda }_{\mathrm{E}\mathrm{2}}=T(r_{\mathrm{4}})/K_{\mathrm{0}}({\zeta }_{r_{\mathrm{4}}}^T)，T(r_{\mathrm{4}})为r_{\mathrm{4}}位置处地层温度；{\zeta }_{r\mathrm{4}}^T为{\zeta }_T中r取固定值为r=r₄。胶结面和水泥环求解温度场待定系数的方程为：

\left[\begin{array}{cccc}{I}_{\mathrm{0}}({\zeta }_{r\mathrm{1}}^T)& K_{\mathrm{0}}({\zeta }_{r\mathrm{1}}^T)& \mathrm{0}& \mathrm{0}\\ I_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)& K_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)& \mathrm{0}& \mathrm{0}\\ I_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)& K_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)& I_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)& K_{\mathrm{0}}({\zeta }_{r\mathrm{2}}^T)\\ \mathrm{0}& \mathrm{0}& I_{\mathrm{0}}({\zeta }_{r\mathrm{3}}^T)& K_{\mathrm{0}}({\zeta }_{r\mathrm{3}}^T)\end{array}\right]\left[\begin{array}{c}{\lambda }_{\mathrm{C}\mathrm{1}}\\ {\lambda }_{\mathrm{C}\mathrm{2}}\\ {\lambda }_{\mathrm{D}\mathrm{1}}\\ {\lambda }_{\mathrm{D}\mathrm{2}}\end{array}\right]=\left[\begin{array}{c}{T}_{\mathrm{0}}^{\mathrm{f}}\\ T_{\mathrm{0}}^{\mathrm{f}\text{'}}\\ \mathrm{0}\\ T_{\mathrm{4}}^{\mathrm{f}}\end{array}\right]（26）

式中，T_{\mathrm{0}}^{\mathrm{f}}和T_{\mathrm{0}}^{\mathrm{f}\text{'}}分别为r_{\mathrm{1}}和r_{\mathrm{2}}位置处的温度，T_{\mathrm{4}}^{\mathrm{f}}为中间变量，T_{\mathrm{4}}^{\mathrm{f}}=T(\mathrm{r}\mathrm{4})/(K_{\mathrm{0}}({\zeta }_{r\mathrm{4}}^T)K_{\mathrm{0}}({\zeta }_{r\mathrm{3}}^T))。通过拉普拉斯逆变换，温度场待定系数随时间的变化曲线如图3所示。

温度场确定后，借助边界条件和连续性条件，建立渗流场和结构场待定系数求解方程，如下所示：

\left[\begin{array}{cccccccccccccc}{L}_{\mathrm{1}}^{\mathrm{1}}& -L_{\mathrm{2}}^{\mathrm{1}}& -L_{\mathrm{3}}^{\mathrm{1}}& -L_{\mathrm{4}}^{\mathrm{1}}& -L_{\mathrm{5}}^{\mathrm{1}}& L_{\mathrm{6}}^{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ L_{\mathrm{1}}^{\mathrm{2}}& L_{\mathrm{2}}^{\mathrm{2}}& -L_{\mathrm{3}}^{\mathrm{2}}& L_{\mathrm{4}}^{\mathrm{2}}& -L_{\mathrm{5}}^{\mathrm{2}}& -L_{\mathrm{6}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& L_{\mathrm{3}}^{\mathrm{3}}& -L_{\mathrm{4}}^{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& L_{\mathrm{3}}^{\mathrm{4}}& L_{\mathrm{4}}^{\mathrm{4}}& \mathrm{0}& \mathrm{0}& -L_{\mathrm{7}}^{\mathrm{4}}& -L_{\mathrm{8}}^{\mathrm{4}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -L_{\mathrm{7}}^{\mathrm{5}}& -L_{\mathrm{8}}^{\mathrm{5}}& \mathrm{0}& \mathrm{0}& L_{\mathrm{11}}^{\mathrm{5}}& L_{\mathrm{12}}^{\mathrm{5}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{7}}^{\mathrm{6}}& -L_{\mathrm{8}}^{\mathrm{6}}& \mathrm{0}& \mathrm{0}& -L_{\mathrm{11}}^{\mathrm{6}}& L_{\mathrm{12}}^{\mathrm{6}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{11}}^{\mathrm{7}}& L_{\mathrm{12}}^{\mathrm{7}}& \mathrm{0}& \mathrm{0}\\ L_{\mathrm{1}}^{\mathrm{8}}& -L_{\mathrm{2}}^{\mathrm{8}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& L_{\mathrm{3}}^{\mathrm{9}}& L_{\mathrm{4}}^{\mathrm{9}}& L_{\mathrm{5}}^{\mathrm{9}}& -L_{\mathrm{6}}^{\mathrm{9}}& -L_{\mathrm{7}}^{\mathrm{9}}& -L_{\mathrm{8}}^{\mathrm{9}}& -L_{\mathrm{9}}^{\mathrm{9}}& L_{\mathrm{10}}^{\mathrm{9}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& L_{\mathrm{3}}^{\mathrm{10}}& -L_{\mathrm{4}}^{\mathrm{10}}& L_{\mathrm{5}}^{\mathrm{10}}& L_{\mathrm{6}}^{\mathrm{10}}& -L_{\mathrm{7}}^{\mathrm{10}}& L_{\mathrm{8}}^{\mathrm{10}}& -L_{\mathrm{9}}^{\mathrm{10}}& -L_{\mathrm{10}}^{\mathrm{10}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{7}}^{\mathrm{11}}& L_{\mathrm{8}}^{\mathrm{11}}& L_{\mathrm{9}}^{\mathrm{11}}& -L_{\mathrm{10}}^{\mathrm{11}}& L_{\mathrm{11}}^{\mathrm{11}}& -L_{\mathrm{12}}^{\mathrm{11}}& -L_{\mathrm{13}}^{\mathrm{11}}& L_{\mathrm{14}}^{\mathrm{11}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{7}}^{\mathrm{12}}& -L_{\mathrm{8}}^{\mathrm{12}}& -L_{\mathrm{9}}^{\mathrm{12}}& L_{\mathrm{10}}^{\mathrm{12}}& -L_{\mathrm{11}}^{\mathrm{12}}& L_{\mathrm{12}}^{\mathrm{12}}& -L_{\mathrm{13}}^{\mathrm{12}}& -L_{\mathrm{14}}^{\mathrm{12}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{11}}^{\mathrm{13}}& L_{\mathrm{12}}^{\mathrm{13}}& L_{\mathrm{13}}^{\mathrm{13}}& -L_{\mathrm{14}}^{\mathrm{13}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& L_{\mathrm{11}}^{\mathrm{14}}& -L_{\mathrm{12}}^{\mathrm{14}}& L_{\mathrm{13}}^{\mathrm{14}}& L_{\mathrm{14}}^{\mathrm{14}}\end{array}\right]\left[\begin{array}{c}{\lambda }_{\mathrm{F}\mathrm{1}}\\ {\lambda }_{\mathrm{F}\mathrm{2}}\\ {\lambda }_{\mathrm{C}\mathrm{3}}\\ {\lambda }_{\mathrm{C}\mathrm{4}}\\ {\lambda }_{\mathrm{C}\mathrm{5}}\\ {\lambda }_{\mathrm{C}\mathrm{6}}\\ {\lambda }_{\mathrm{D}\mathrm{3}}\\ {\lambda }_{\mathrm{D}\mathrm{4}}\\ {\lambda }_{\mathrm{D}\mathrm{5}}\\ {\lambda }_{\mathrm{D}\mathrm{6}}\\ {\lambda }_{\mathrm{E}\mathrm{3}}\\ {\lambda }_{\mathrm{E}\mathrm{4}}\\ {\lambda }_{\mathrm{E}\mathrm{5}}\\ {\lambda }_{\mathrm{E}\mathrm{6}}\end{array}\right]=\left[\begin{array}{c}{\psi }_{\mathrm{1}}\\ {\psi }_{\mathrm{2}}\\ {\psi }_{\mathrm{3}}\\ {\psi }_{\mathrm{4}}\\ {\psi }_{\mathrm{5}}\\ {\psi }_{\mathrm{6}}\\ {\psi }_{\mathrm{7}}\\ {\psi }_{\mathrm{8}}\\ {\psi }_{\mathrm{9}}\\ {\psi }_{\mathrm{10}}\\ {\psi }_{\mathrm{11}}\\ {\psi }_{\mathrm{12}}\\ {\psi }_{\mathrm{13}}\\ {\psi }_{\mathrm{14}}\end{array}\right]（27）

式中：λ_F1、λ_F2为套管结构待定系数；λ_C3、λ_C4、λ_C5、λ_C6为胶结面结构待定系数；λ_D3、λ_D4、λ_D5、λ_D6为水泥环结构待定系数；λ_E3、λ_E4、λ_E5、λ_E6为地层结构待定系数；{\psi }_k、L_l^k均为与材料力学参数、贝塞尔函数有关的常数参数，下标k为行号，下标l为列号。

由于系数矩阵庞大，求解复杂，因此本文借助Matlab软件进行求解。各结构的渗流场和结构场系数随时间变化规律如图4所示。

图4中，各类系数经过一定形式的数据处理才投影到图上，如-λ_F1×10⁷表明曲线显示的是λ_F1乘以10⁷并取相反数的值。

3.3 热流固耦合模型验证

选取文献[39]中的地层径向方向上不同位置温度随时间的变化曲线，验证本文热流固耦合模型的温度场。验证过程选取文献[39]中的套管‒水泥环‒地层的力学参数和几何数据，由于本文在推导过程中考虑了胶结面层，因此在验证过程中，选择胶结面的力学参数与水泥环相同，进而本文模型退化为文献[39]的力学模型。通过计算绘制曲线如图5所示，实线为本文模型计算结果，离散点为参考文献[39]中提取的验证温度。通过对比10、100、600 min 3个时刻下的温度曲线发现，本文所获得的温度曲线与文献[39]得到的曲线基本一致，从而验证了热流固耦合模型在温度场计算上的准确性。

选取文献[40]中的地层径向方向上不同位置随距离变化的径向应力和径向位移曲线，验证热流固耦合模型的结构场。验证过程与上述温度场验证思路一致。通过计算绘制了曲线分别如图6和7所示，实线为本文模型计算结果，离散点为文献[40]中提取的径向应力和径向位移。

由图6、7的结果可知，热流固耦合模型所获得的径向应力和位移与文献[40]得到的数据基本一致。

综上所述，通过已发表的温度场和应力场数据分别验证了热流固耦合模型计算得到温度场和应力场的准确性，进而可以认为本文所建立的热流固耦合模型及编程过程是准确的。

4 干热岩地热井算例分析

干热岩地热井通常为高温高压井，高温和高压均对胶结面受力有很大影响。以青海共和盆地恰卜恰地区干热岩井为例^[41]，在钻井深度达到3 700 m时，井筒系统的几何参数如表1所示。

井底地层温度为200 ℃，地层最大水平主应力58 MPa，最小主应力45 MPa，套管内压为32 MPa，压裂液温度为20 ℃。各层结构相关的力学参数如表2、3所示。

4.1 干热岩地热井温度场、应力场分析

利用上述模型，通过Matlab软件计算，得到了热流固耦合作用下干热岩地热井套管‒胶结面‒水泥环‒地层系统的各界面层温差（井筒温度与实时温度的差值）、径向应力，图8为各界面处温度随时间变化曲线。图9为各界面处径向应力随时间的变化曲线。

由图8可知，套管‒胶结面界面、胶结面‒水泥环界面温度快速降低，在2.3 d（1 d=0.864×10⁵ s）左右趋于稳定。水泥环‒地层界面远离套管，其温度降低的时间稍晚于其他两界面，在6.9 d左右接近套管内壁温度。套管‒胶结面界面处温差变化曲线与胶结面‒水泥环界面温差变化曲线相近，分析原因是胶结面厚度较薄，导致其胶结面内外边界温差较小。

由图9可知，所有界面上的径向应力均随着时间的增加而增加，水泥环‒地层界面处的应力增幅较其他界面大。胶结面沿径向方向的厚度较小，导致套管‒胶结面界面和胶结面‒水泥环界面处的径向应力基本相同。套管‒胶结面界面和胶结面‒水泥环界面径向应力在开始有短暂的降低，0.2 d时达到最低点，此时径向应力为6.31 MPa，分析原因是井筒温度降低所产生的径向收缩，导致径向应力降低。

当时间累计到7.0 d时，绘制径向应力沿井眼半径方向的变化曲线，如图10所示。

由图10可知，径向应力在套管内快速增大，到达胶结面后，增长速率减缓。套管‒水泥环胶结面处的径向应力较大，同时胶结面的破坏强度又很小，因此胶结面相对其他位置更容易发生破坏。

4.2 胶结面影响因素分析

设置套管、胶结面、水泥环、地层的计算参数如表2和3所示，研究胶结面处的弹性模量、泊松比、孔隙率、T‒O系数以及地层与套管的温差等参数对套管‒胶结面界面径向应力的影响规律。

4.2.1 弹性模量的影响

保持其他参数不变，改变胶结面的弹性模量，以胶结面与水泥环的弹性模量比值为变量，分析胶结面弹性模量对套管‒胶结面界面的径向应力的影响。胶结面与水泥环弹性模量比值E^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}/E^{\mathrm{c}\mathrm{e}\mathrm{m}}变化范围为[0.33，1.67]，计算得到径向应力随弹性模量比值变化曲线如图11所示。

由图11可知，初期各径向应力变化较小，随着时间延长，套管‒胶结面界面处的径向应力逐渐增大。随着胶结面与水泥环弹性模量比值的增大，径向应力先降低后升高，在E^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}/E^{\mathrm{c}\mathrm{e}\mathrm{m}}为0.33和1.67时，径向应力大于其他几种情况，表明此时套管‒胶结面界面处更容易发生受压破坏；而当E^{\mathrm{p}\mathrm{c}\mathrm{e}\mathrm{m}}/E^{\mathrm{c}\mathrm{e}\mathrm{m}}的值为1.11时径向应力最小，表明胶结面与水泥环的弹性模量相近时，径向应力随着时间的涨幅较小，反之，径向应力涨幅较大。因此，选择胶结面的弹性模量越接近水泥环的弹性模量，胶结面越不容易发生破坏，进而有利于套管‒胶结面‒水泥环‒地层组合体的完整性。

4.2.2 泊松比的影响

保证其他参数不变，改变胶结面的泊松比，分析胶结面泊松比对套管‒胶结面界面的径向应力的影响。胶结面泊松比变化范围为[0.26，0.34]，计算得到径向应力随胶结面泊松比变化曲线如图12所示。

由图12可知，在初始时间段内，胶结面泊松比的变化对径向应力的影响较小，在较长时间内，胶结面泊松比的变化对径向应力的影响较大。随着泊松比的增加，胶结面的径向应力先降低后增加。当泊松比为0.30时，径向应力最小，而当泊松比为0.34时，径向应力最大。本算例选择的水泥环泊松比为0.30，因此可以认为当胶结面的泊松比接近水泥环泊松比时，其最大径向应力越小，胶结面越不容易发生破坏，反之，当胶结面的泊松比与水泥环泊松比相差越大时，其径向应力增加越快，越不利于组合结构的完整性。

4.2.3 孔隙率的影响

保证其他参数不变，改变胶结面的孔隙率，分析胶结面孔隙率对套管‒胶结面界面径向应力的影响。孔隙率变化范围为[0.1，0.5]，计算得到径向应力随孔隙率变化曲线如图13所示。

由图13可知，随着孔隙率的增加，胶结面径向应力的增长速度逐渐增加，当孔隙率为[0.1，0.3]时，径向应力增长速率较小，当孔隙率大于0.4时，径向应力增长速率变大；在初期时间段内（t<2×10⁵ s），径向应力增长速率较缓慢，在后期时间段内（2×10⁵ s \let\le1×10⁶ s），径向应力增长速率变大，但随着时间的增加基本保持恒定，因而径向应力曲线展现为线性形式。综上分析，胶结面孔隙率会显著影响套管‒胶结面界面处的径向应力。通过文献[42]中不同养护时间胶结面平均孔隙率图可知，胶结面在养护条件较好的情况下最大孔隙率不会超过0.3，因此得出以下结论：孔隙率对胶结面的径向应力的影响较小。

4.2.4 T‒O系数的影响

保证其他参数不变，改变胶结面的T‒O系数，分析T‒O系数对套管‒胶结面界面径向应力的影响。T‒O系数变化范围为[0.5，2.0]，计算得到径向应力随T‒O系数变化曲线如图14所示。

由图14可知，随着T‒O系数的增加，径向应力先降低后增加，当T‒O系数为1.0时，径向应力最小，同时径向应力增长速率也最小；当T‒O系数为0.5和2.0时，径向应力较大，且径向应力曲线接近。由此可以认为，合适的T‒O系数能有效降低胶结面的径向应力，因此，T‒O系数是研究套管‒水泥环胶结面应力场不可忽略的参数，由此可以认为，本文建立的热流固耦合模型考虑T‒O系数具有一定的必要性。本算例中水泥环的T‒O系数为1.0，因此，当胶结面的T‒O系数与水泥环接近的时候，套管‒水泥环胶结面的径向应力最小。

4.2.5 地层温度的影响

保证其他参数不变，改变地层温度，分析地层温度对套管‒胶结面界面径向应力的影响。地层与套管内注入水的温差变化范围为[160，240] ℃，计算得到径向应力随温差变化曲线如图15所示。

由图15可知，随着温差的增加，径向应力随之增加；温差不同，径向应力增长速率差异较大，温差越大，径向应力增长速率越大，且温差\mathrm{\Delta }T>200 ℃时的径向应力增长速率远大于温差\mathrm{\Delta }T\le200 ℃时的增长速率。因此，从整体结构不发生破坏角度分析，尽可能提高套管内温度，能够有效降低胶结面处的径向应力。

4.3 胶结面破裂分析

通过文献^[43]调研发现，胶结面的破坏形式一般被认为有沿着井眼轴线的剪切破裂、垂直井眼轴线的张拉破裂和压缩破坏等3种主要形式，在本文的平面应变模型中，胶结面的破裂形式以张拉破裂为主。由于本文模型相较于已发表的模型在结构和参数上考虑得更为精细，因此，计算得到胶结面应力场的精度也会更高，从而为胶结面破裂分析提供精确的数据基础。

胶结面为脆性材料，其张拉破裂满足最大拉应力准则，即：

（28）

式中：{\sigma }_{\mathrm{t}\mathrm{e}\mathrm{n}}为胶结面的破裂拉应力，可通过套管‒水泥环压入试验得到；σ₁为胶结面上任一点的最大主应力。

胶结面压缩破坏，满足Mohr‒Coulomb准则，即：

（29）

式中，\phi、C分别为胶结面内摩擦角和黏聚力，σ₃为胶结面上任一点的最小主应力。

由于偏差应力的影响，最大主应力方向会偏离径向或环向方向，但在偏差应力较小的情况下，主应力可以用环向应力和径向应力近似代替。

在工程中，胶结面沿着环向张拉破裂更容易引起胶结面密封失效。因此，现场施工应避免径向应力达到破裂拉应力的工况发生。

5 结论

1）考虑套管内水压引起的孔隙水压力、套管内注入水与地下地层之间温差引起的热载荷和T‒O效应，建立了干热岩地热井套管‒胶结面‒水泥环‒地层组合体的热流固耦合模型。

2）通过拉普拉斯变化和逆变换，并借助边界条件及连续性条件，得到了热流固耦合模型中各个结构的应力场、渗流场、温度场表达式。并通过Matlab软件，求解了各表达式的待定系数。

3）通过改变套管‒水泥环胶结面的弹性模量、泊松比、孔隙率、T‒O系数和地层温度等参数，分析了热流固耦合参数对干热岩地热井套管‒胶结面界面径向应力的影响，结果表明：当胶结面与水泥环具有相近的弹性模量、泊松比、T‒O系数时，套管‒水泥环胶结面处的径向应力较小；地层温度和孔隙率越小，套管‒水泥环胶结面处的径向应力就越小。

参考文献

[1] Song Guofeng, Li Gensheng, Song Xianzhi,et al.Multi-objective balanced method of optimizing the heat extraction performance for hot dry rock[J].Natural Gas Industry B,2022,9(6):497‒510. doi:10.1016/j.ngib.2022.11.002

[2] Song Xianzhi, Shi Yu, Li Gensheng,et al.Numerical simulation of heat extraction performance in enhanced geothermal system with multilateral wells[J].Applied Energy,2018,218:325‒337. doi:10.1016/j.apenergy.2018.02.172

[3] Zhou Zhou, Jin Yan, Lu Yunhu,et al.Present challenge and prospects of drilling and hydraulic fracturing technology for hot dry rock geothermal reservoir[J].Scientia Sinica (Physica,Mechanica & Astronomica),2018,48(12):97‒102. doi:10.1360/SSPMA2018-00184

[4] 周舟,金衍,卢运虎,等.干热岩地热储层钻井和水力压裂工程技术难题和攻关建议[J].中国科学:物理学 力学 天文学,2018,48(12):97‒102. doi:10.1360/SSPMA2018-00184

[5] Huang Wenbo, Cen Jiwen, Chen Juanwen,et al.Heat extraction from hot dry rock by super-long gravity heat pipe:A field test[J].Energy,2022,247:123492. doi:10.1016/j.energy.2022.123492

[6] Samin M Y, Faramarzi A, Jefferson I,et al.A hybrid optimisation approach to improve long-term performance of enhanced geothermal system (EGS) reservoirs[J].Renewable Energy,2019,134:379‒389. doi:10.1016/j.renene.2018.11.045

[7] Gong Facheng, Guo Tiankui, Sun Wei,et al.Evaluation of geothermal energy extraction in Enhanced Geothermal System(EGS) with multiple fracturing horizontal wells (MFHW)[J].Renewable Energy,2020,151:1339‒1351. doi:10.1016/j.renene.2019.11.134

[8] Zhang Wei, Wang Chunguang, Guo Tiankui,et al.Study on the cracking mechanism of hydraulic and supercritical CO₂ fracturing in hot dry rock under thermal stress[J].Energy,2021,221:119886. doi:10.1016/j.energy.2021.119886

[9] Feng Chenchen, Wang Huaijiu, Jing Zefeng.Investigation of heat extraction with flowing CO₂ from hot dry rock by numerical study[J].Renewable Energy,2021,169:242‒253. doi:10.1016/j.renene.2021.01.037

[10] Lavrov A, Andrade J D, Torsæter M,et al.Experimental set-up for testing cement sheath integrity in arctic well[C]//Processings of Arctic Technology Conference,Houston:OTC,2014:24587. doi:10.4043/24587-ms

[11] Liu Lingrong.Research and application of thermostable cement slurry system[J].Exploration Engineering (Rock & Soil Drilling and Tunneling),2018,45(9):29‒33.

[12] 刘玲榕.抗高温水泥浆体系研究及应用[J].探矿工程(岩土钻掘工程),2018,45(9):29‒33.

[13] Hailesilassie B W, Partl M N.Adhesive blister propagation under an orthotropic bituminous waterproofing membrane[J].Construction and Building Materials,2013,48:1171‒1178. doi:10.1016/j.conbuildmat.2013.07.047

[14] Torsæter M, Todorovic J, Lavrov A,et al.Avoiding damage of CO₂ injection wells caused by temperature variations[J].Energy Procedia,2017,114:5275‒5286. doi:10.1016/j.egypro.2017.03.1645

[15] Lavrov A, Todorovic J, Torsæter M.Impact of voids on menchanical stability of well cement[J].Engery Procedia, 2016,86:401‒410. doi:10.1016/j.egypro.2016.01.041

[16] Tian Wei, Zhang Pengkun, Xie Yongli,et al.3D distribution characteristics on concrete porous structure under freeze-thaw environment based on CT technique[J].Journal of Chang'an University(Natural Science Edition),2016,36(3):49‒55.

[17] 田威,张鹏坤,谢永利,等.冻融环境下基于CT技术混凝土孔隙结构的三维分布特征[J].长安大学学报(自然科学版),2016,36(3):49‒55.

[18] Luo Changji, Wang Yunliang, Zhang Bin. Experimental study on bond strength at the cement sheath interface in well cementing[J].Oil Drilling and Production Process,1993,15(3):47‒51.

[19] 罗长吉,王允良,张彬.固井水泥环界面胶结强度实验研究[J].石油钻采工艺,1993,15(3):47‒51.

[20] Xu Yongzhi, Zheng Tao, Jing Zengxiu.The laboratory study of slurry and mud additive[J].Oil Drilling & Production Technology,1994,16(2):46‒48.

[21] 许永志,郑韬,敬增秀.固井、钻井液添加剂室内实验研究[J].石油钻采工艺,1994,16(2):46‒48.

[22] Qi Fengzhong, Yang Chengji, Liu Zishuai.Technical research on improving cementing quality for complex oil and gas wells[J]. Petroleum Science and Technology Forum,2013,32(1):19‒22.

[23] 齐奉忠,杨成颉,刘子帅.提高复杂油气井固井质量技术研究[J].石油科技论坛,2013,32(1):19‒22.

[24] Wang W, Dahi‒Taleghani A.Emergence and propagation of delamination cracks along the casing-cement interface[C]//Processdings of 46th U.S. Rock Mechanics/Geomechanics Symposium.Chicago:ARMA,2012:2012449.

[25] Nath F, Kimanzi R J, Mokhtari M,et al.A novel method to investigate cement-casing bonding using digital image correlation[J].Journal of Petroleum Science and Engineering,2018,166:482‒489. doi:10.1016/j.petrol.2018.03.068

[26] Goodwin K J, Crook R J.Cement sheath stress failure[J].SPE Drilling Engineering,1992,7(4):291‒296. doi:10.2118/20453-pa

[27] Carpenter R B, Brady J L, Blount C G.The effects of temperature and cement admixes on bond strength[J].Journal of Petroleum Technology,1992,44(8):936‒941. doi:10.2118/22063-pa

[28] Bois A P, Garnier A, Galdiolo G,et al.Use of a mechanistic model to forecast cement-sheath integrity[J].SPE Drilling & Completion,2012,27(2):303‒314. doi:10.2118/139668-pa

[29] Yin Youquan, Chen Zhaowei, Li Ping'en.Theoretical solutions of stress distribution in casing-cement and stratum system[J].Chinese Journal of Theoretical and Applied Mechanics,2006,38(6):835‒842.

[30] 殷有泉,陈朝伟,李平恩.套管‒水泥环‒地层应力分布的理论解[J].力学学报,2006,38(6):835‒842.

[31] Shi Yucai, Guo Boyun, Guan Zhichuan,et al.Influence of the initial loaded state on the stress distribution of a wellbore system[J].Journal of Petroleum Science and Engineering,2017,157:547‒557. doi:10.1016/j.petrol.2017.07.054

[32] Hu Zhiqiang, Yang Jin, Liu Shujie,et al.Prediction model of multilayer casing annular pressure buildup based on casing-cement sheath-formation with thermo-structural coupling effects[J].Journal of Engineering Thermophysics,2018,39(8):1824‒1832.

[33] 胡志强,杨进,刘书杰,等.基于套管-水泥环-地层热固耦合作用的多层套管环空附加压力预测模型[J].工程热物理学报,2018,39(8):1824‒1832.

[34] Zhao Xinbo, Han Shengchao, Yang Xiujuan,et al.Mechanical characteristics analysis of casing-cement sheath-formation multilayer composite system with thermo-structural coupling effects[J].Journal of Central South University(Science and Technology),2017,48(3):837‒843.

[35] 赵新波,韩生超,杨秀娟,等.热固耦合作用下的套管‒水泥环‒地层多层组合系统应力分析[J].中南大学学报(自然科学版),2017,48(3):837‒843.

[36] Zhao Xinbo, Yang Xiujuan, Li Xiangyang,et al.Integrity analysis of high temperature and high pressure wellbores with thermo-structural coupling effects[J].Chinese Journal of Engineering,2016,38(1):11‒18.

[37] 赵新波,杨秀娟,李向阳,等.考虑热固耦合作用的高温高压井井筒完整性分析[J].工程科学学报,2016,38(1):11‒18.

[38] Wu Jiwei, Wang Xuegang, Song Lin,et al.Microannulus formation mechanism at the cementing interface of a thermal recovery well during cyclic steam injection[J].Advances in Civil Engineering,2020,2020(1):8217013. doi:10.1155/2020/8217013

[39] Luo Zhiyuan, Bryant S.Influence of thermo-elastic stress on fracture initiation during CO₂ injection and storage[J].Energy Procedia,2011,4:3714‒3721. doi:10.1016/j.egypro.2011.02.304

[40] Nygaard R, Salehi S, Weideman B,et al.Effect of dynamic loading on wellbore leakage for the wabamun area CO₂-sequestration project[J].Journal of Canadian Petroleum Technology,2014,53(1):69‒82. doi:10.2118/146640-pa

[41] Ghassemi A, Tao Q, Diek A.Influence of coupled chemo-poro-thermoelastic processes on pore pressure and stress distributions around a wellbore in swelling shale[J].Journal of Petroleum Science and Engineering,2009,67(1/2):57‒64. doi:10.1016/j.petrol.2009.02.015

[42] Trémosa J, Gonçalvès J, Matray J M,et al.Estimating thermo-osmotic coefficients in clay-rocks:Ⅱ.In situ experimental approach[J].Journal of Colloid and Interface Science,2010,342(1):175‒184. doi:10.1016/j.jcis.2009.09.055

[43] Sevostianov I, Kachanov M.Effect of interphase layers on the overall elastic and conductive properties of matrix composites.Applications to nanosize inclusion[J].International Journal of Solids and Structures,2007,44(3/4):1304‒1315. doi:10.1016/j.ijsolstr.2006.06.020

[44] Li Zongli, Deng Chaoli, Zhang Guohui.A predictive model of effective elastic modulus of concrete under influence of aggregate gradation[J].Journal of Hydraulic Engineering,2016,47(4):575‒581.

[45] 李宗利,邓朝莉,张国辉.考虑骨料级配的混凝土有效弹性模量预测模型[J].水利学报,2016,47(4):575‒581.

[46] Duan H L, Yi X, Huang Z P,et al.A unified scheme for prediction of effective moduli of multiphase composites with interface effects.Part Ⅰ:Theoretical framework[J].Mechanics of Materials,2007,39(1):81‒93. doi:10.1016/j.mechmat.2006.02.009

[47] Deng Wenqian, Xu Lihua, Chi Yan,et al.Calculation of effective elastic modulus for hybrid fiber reinforced concrete based on homogenization theory[J].Journal of the Chinese ceramic society,2019,47(2):151‒170.

[48] 邓文茜,徐礼华,池寅,等.基于均匀化理论的混杂纤维混凝土有效弹性模量计算[J].硅酸盐学报, 2019,47(2):151‒170.

[49] Lavrov A, Panduro E A C, Torsæter M.Synchrotron study of cement hydration:Towards computed tomography analysis of interfacial transition zone[J].Energy Procedia,2010,114:5109‒5117.

[50] Torsæter M, Todorovic J, Lavrov A.Structure and debonding at cement‒steel and cement‒rock interfaces:Effect of geometry and materials[J].Construction and Building Materials,2015,96:164‒171. doi:10.1016/j.conbuildmat.2015.08.005

[51] Lv Xinrun.Study on temperature field of CO₂ fractured wellbore in shale gas horizontal well[D].Beijing:China University of Petroleum (Beijing),2019.

[52] 吕欣润.页岩气水平井CO₂压裂井筒温度场研究[D].北京:中国石油大学(北京),2019.

[53] Li Xiaorong, Gu Chenwang, Ding Zechen,et al.THM coupled analysis of cement sheath integrity considering well loading history[J].Petroleum Science,2023,20(1):447‒459. doi:10.1016/j.petsci.2022.09.001

[54] Jin Yuxiang.Study on key technologies of drilling in dry-hot rock in Gonghe basin,Qinghai province[D].Beijing:China University of Petroleum (Beijing),2020.

[55] 金宇翔.青海共和盆地干热岩钻井关键技术研究[D].北京:中国石油大学(北京),2020.

[56] Scrivener K L, Crumbie A K, Laugesen P.The interfacial transition zone (ITZ) between cement paste and aggregate in concrete[J].Interface Science,2004,12(4):411‒421. doi:10.1023/b:ints.0000042339.92990.4c

[57] Restrepo M M, Teodoriu C, Salehi S,et al.A novel way to look at the cement sheath integrity by introducing the existence of empty spaces inside of the cement (voids)[J].Journal of Natural Gas Science and Engineering, 2020,77:103274. doi:10.1016/j.jngse.2020.103274

基金资助

国家自然科学青年基金项目(52104004)

国家自然科学青年基金项目(12302264)