基于流固耦合的含水合物地层井壁稳定非稳态解析模型

摘要

针对深海水合物地层钻井过程中的井壁稳定问题，考虑水合物分解、热传导、力场-渗流场全耦合作用，建立了过压和欠压钻井下渗流、温度、力场随时间和空间变化的非稳态解析模型。解析结果与相同条件下的数值结果吻合良好，且与力场-渗流半耦合解析结果进行了对比。基于解析模型对井壁稳定的关键参数如钻井液压力、水合物分解引起的地层弹性模量劣化程度等进行了分析，结果表明：①与半耦合分析结果相比，考虑体变对渗流的影响后，过（欠）压钻井时孔压减小（增大）、应力增大（减小），增量径向位移减小；②最危险位置在井壁处，过高或过低的钻井液压力均会导致井壁失稳，水合物分解引起的地层劣化将降低最安全钻井液压力；③水合物分解引起的地层刚度降低极易诱发井壁失稳。在通常条件下，过压钻井时分解域弹性模量降低50%即可导致井壁失稳。

关键词

井壁稳定; 流固耦合; 水合物; 非稳态分析; 解析解

天然气水合物是在自然界中由甲烷为主的烃类气体与水在高压、低温条件下形成的似冰状固态结晶物质^［1-2］，因可燃烧，被俗称为可燃冰。因其储量丰富（天然气水合物的含碳量大约为化石能源含碳量的2倍）^［3］、分布广泛、能量密度高、清洁无污染等优势，被认为是未来的理想替代能源。现探明的绝大部分（大于98%的储量）天然气水合物分布于海洋地层环境中^［4］。由于天然气水合物储层环境的特殊性和水合物本身易分解的特点，开采过程涉及相态转变（固体水合物分解为气体甲烷和液体水）、多相渗流和多物理场耦合（渗流压力、温度、化学和力场）^［5-6］，因而其开采难度远大于常规油气。在钻井和开采过程中，外界扰动引发的水合物分解使水合物沉积物胶结弱化，降低其强度和刚度，极易诱发井壁失稳等安全事故^［7］。

井壁失稳是水合物钻井和开采面临的主要问题之一^［8］。由于问题的复杂性，目前大多数学者采用数值模拟方法建立多场耦合模型研究海洋水合物开采的井壁稳定性。这些数值模型大致可以分为2类：①通过TOUGH+HYDRATE软件模拟水合物分解涉及的温度、渗流压力和化学场分布规律，并将分析结果代入力场分析软件（如FLAC^3D）实现与力场的耦合作用，再将力场的作用结果（如渗透率、孔隙率的改变）重新代入TOUGH+HYDRATE软件中，可实现一次渗流-温度-化学-力场的全耦合作用。循环迭代以上过程，可得受水合物分解影响的相对精确的多场耦合分析结果^［9-10］。②首先建立海洋水合物开采的多场耦合数学模型，包括渗流方程、热量传递方程、化学场方程、力场变形方程及辅助方程（如分解动力学方程、水合物相平衡方程、孔隙率方程、渗透率方程等），并利用有限元^［11］或有限差分^［12］等数值方法求解该数学模型，从而实现多场耦合作用。

作为数值方法的验证和补充以及理论分析手段之一，解析模型对实际工程进行合理等效和简化（保留核心影响因素），通过严格数学推导获得满足全部数学方程的封闭或解析形式的解答，能直接显示核心参数的函数形式影响关系和力学控制机理，实现多场耦合条件下井壁力学状态的快速和稳定计算。有关井壁稳定的多场耦合解析模型多针对常规油气开采，根据流固耦合程度可分为半耦合和全耦合模型。半耦合模型通过有效应力原理考虑孔隙压力对力场的单向影响，并结合不同的屈服准则进行稳定性的判断^［13］；全耦合模型则综合考虑了渗流场和力场的相互影响，其中渗流场中引入固体体积变形的影响^［14-15］。由于含水合物地层井壁稳定问题的特殊性和复杂性，目前建立的深海水合物钻井井壁稳定解析模型介于半耦合和全耦合之间，考虑了水合物分解对地层弹塑性力学特性和渗透性的影响^［16-17］以及渗流场对力场作用，尚未考虑流固全耦合作用。Cheng等^［18］建立了半解析渗流-温度-应力耦合模型，结合解析求解应力场和数值求解温度、渗流场，通过迭代计算，研究水合物降压开采过程中的孔压、温度、应力分布。

通过对实际问题合理简化，建立考虑水合物分解、热传导、力场-渗流场全耦合作用的非稳态解析模型，并分析水合物地层钻井井壁稳定性，为实际工程提供参考。

1 力学模型

考虑水合物分解及非稳态渗流-力场全耦合作用下的深海水合物钻井过程中的井壁稳定问题。通常含水合物地层位于上下不透水层之间^［19］，流体流动和热量传递可认为仅发生在水平面上。同时，钻井轴向尺寸远大于径向尺寸，储层所处地应力较大且其面内两方向应力近似相等^［20］，该问题可简化为水平面内轴对称平面应变问题。钻井活动将改变地层温度、孔压并导致水合物分解，并引起地层力学特性、渗透率、热传导系数等发生显著变化。这些因素之间是相互作用、相互耦合的关系。与以往研究不同，本文将全面考虑渗流场-力场间的相互耦合以及水合物分解作用（三场关系如图1），获得更为接近实际情况的多场耦合井壁稳定解析模型。根据水合物是否分解将地层分为具有不同物理力学性质的2个区域（见图2）：分解区（r₀<r< r₁）和未分解区（r> r₁）^［21-22］，其中r₀、r₁分别为井壁半径、分解半径。分解区内水合物完全分解，弹性模量、渗透率、热传导系数分别为E₁、k₁、λ₁；未分解区内水合物未分解，弹性模量、渗透率、热传导系数分别为E₂、k₂、λ₂。作为等效模型，地层视为均匀、各向同性，渗流、温度、力场边界条件如图2所示，P_wf、T_wf分别为钻井液压力、温度，P_∞、T_∞分别为初始地层压力、温度，σ_∞为初始地应力。推导采用极坐标（r， θ），压应力为正，径向位移远离井壁中心为正。

图1  渗流场、温度场与力场耦合关系

Fig.1  Seepage and temperature versus mechanical fields in analytical model

图2  含水合物地层井壁分析的几何模型及其边界条件

Fig.2  Geometry and boundary conditions for wellbore stability in methane hydrate-bearing sediments

2 渗流场、温度场与力场解答

2.1　渗流、温度和力场控制方程及其解耦

由于在钻井过程中水合物分解较少，因此忽略分解产生的水气及热量变化的影响^［23］。采用Coussy提出的温度-渗流-力场耦合模型^［24-25］。由于温度值和梯度均较小，可以忽略温度应力和温度对渗流的影响。渗流控制方程为

$$\frac{k}{\mu }{\nabla }^{\mathrm{2}}P=\frac{\mathrm{1}}{M}\frac{\partial P}{\partial t}+\alpha \frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}$$

（1）

式中：P为孔压；t为时间；ε_v为地层体积应变；k为地层渗透率；μ为流体黏滞系数；α为Biot系数，α=1–K/K_s，K为含水合物地层体积模量，K=E/［3（1–2ν）］，K_s为地层固体颗粒体积模量；M为Biot模量，$M=\frac{\mathrm{2}G\left({\nu }_{\mathrm{u}}-\nu \right)}{{\alpha }^{\mathrm{2}}\left(\mathrm{1}-\mathrm{2}{\nu }_{\mathrm{u}}\right)\left(\mathrm{1}-\mathrm{2}\nu \right)}$，ν、ν_u分别为含水合物地层的排水、不排水泊松比，${\nu }_{\mathrm{u}}=\frac{\mathrm{3}{K}_{\mathrm{u}}-\mathrm{2}G}{\mathrm{2}\left(\mathrm{3}{K}_{\mathrm{u}}+G\right)}$，K_u为地层不排水体积模量，$K_{\mathrm{u}}=K+\frac{{\alpha }^{\mathrm{2}}}{\left(\alpha -\varphi \right)/K_{\mathrm{s}}+\varphi /K_{\mathrm{f}}}$ ^［26-27］；E、G分别为含水合物地层的弹性、剪切模量，G=E/［2（1+v）］；K_f为流体体积模量，本文流体指水和气体的混合物；ϕ为地层孔隙率。

温度场的热量传递方式仅考虑热传导，热对流通过增大分解域内的热传导系数近似考虑^［28］。轴对称下温度T_i （i=1， 2分别代表分解区和未分解区）控制方程为

$${\nabla }^{\mathrm{2}}{T}_i{}^{}=\frac{\mathrm{1}}{a_i}\frac{\partial T_i}{\partial t}$$

（2）

式中：T_i为含水合物地层温度；a_i为含水合物地层的热扩散系数，a_i =λ_i/（ρ_ic_hi）；λ_i、ρ_i和c_hi分别为含水合物地层的热传导系数、密度和比热容。

考虑渗流影响的力场平衡方程为

$$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}^{}}{\mathrm{d}r}+\frac{{\sigma }_{\mathrm{r}}^{}-{\sigma }_{\mathrm{\theta }}^{}}{r}+\alpha \frac{\mathrm{d}P}{\mathrm{d}r}=\mathrm{0}$$

（3）

式中：σ_r 、σ_θ分别为径向、环向正应力（有效应力）。几何方程为

$${\epsilon }_{\mathrm{r}}^{}=\frac{\mathrm{d}{u}_{\mathrm{r}}^{}}{\mathrm{d}r},{\epsilon }_{\mathrm{\theta }}^{}=\frac{u_{\mathrm{r}}^{}}{r}$$

（4）

式中：ε_r 、ε_θ分别为径向、环向正应变；u_r为径向位移。由轴对称条件，在研究平面上其余应力和应变为零。本模型首先按弹性阶段进行计算，然后根据强度准则和应力分布进行破坏分析。本构方程为

$${\sigma }_{\mathrm{r}}^{}=\frac{\left(\mathrm{1}-\nu \right)E}{\left(\mathrm{1}-\mathrm{2}\nu \right)\left(\mathrm{1}+\nu \right)}\left({\epsilon }_{\mathrm{r}}^{}+\frac{\nu }{\mathrm{1}-\nu }{\epsilon }_{\mathrm{\theta }}^{}\right)$$，$${\sigma }_{\mathrm{\theta }}^{}=\frac{\left(\mathrm{1}-\nu \right)E}{\left(\mathrm{1}-\mathrm{2}\nu \right)\left(\mathrm{1}+\nu \right)}\left(\frac{\nu }{\mathrm{1}-\nu }{\epsilon }_{\mathrm{r}}^{}+{\epsilon }_{\mathrm{\theta }}^{}\right)$$

（5）

由平面应变条件得，z向正应力表达式为^［29］

$${\sigma }_{\mathrm{z}}^{}=\frac{\mathrm{1}}{\mathrm{2}}\left({\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{\theta }}\right)$$

（6）

可见，渗流控制方程与力场控制方程是耦合的。为进行解耦，将式（4）代入式（5），并将结果代入式（3），可得u_r的解析表达式为

$$u_{\mathrm{r}}^{}=\frac{(\mathrm{2}\nu -\mathrm{1})(\mathrm{1}+\nu )}{E(\mathrm{1}-\nu )r}{\int }_{r_{\mathrm{0}}}^r\alpha P(r)r\mathrm{d}r+Cr+\frac{D}{r}$$

（7）

式中：C、D为待定系数，可由边界和连续条件确定。对平面应变问题，体积应变可写为${\epsilon }_{\mathrm{v}}={\epsilon }_{\mathrm{r}}+{\epsilon }_{\mathrm{\theta }}$。由式（7）和式（4）可得

$$\frac{\partial {\epsilon }_{\mathrm{v}}^{}}{\partial t}=\alpha \frac{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}{E(\mathrm{1}-\nu )}\frac{\partial P}{\partial t}$$

（8）

$\frac{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}{E(\mathrm{1}-\nu )}$则为一维Biot固结系数^［30］。最后，将式（8）代入式（1），可得考虑渗流-力场全耦合的渗流方程为

$$\frac{k}{\mu }{\nabla }^{\mathrm{2}}P=\left[\frac{\mathrm{1}}{M}+{\alpha }^{\mathrm{2}}\frac{(\mathrm{1}-\mathrm{2}\nu )(\mathrm{1}+\nu )}{E(\mathrm{1}-\nu )}\right]\frac{\partial P}{\partial t}$$

（9）

2.2　渗流场和温度场解答与分解域半径的确定

令$\frac{\mathrm{1}}{{\eta }_i}=\frac{\mu }{k_i}\left[\frac{\mathrm{1}}{M_i}+{\alpha }_i^{\mathrm{2}}\frac{(\mathrm{1}-\mathrm{2}{\nu }_i)(\mathrm{1}+{\nu }_i)}{E_i(\mathrm{1}-{\nu }_i)}\right]$，式（9）可化简为

$${\nabla }^{\mathrm{2}}{P}_i=\frac{\mathrm{1}}{{\eta }_i}\frac{\partial P_i}{\partial t}$$

（10）

式中：η_i为常量，与流体、固体变形特性相关。偏微分方程式（10）可利用分离变量法求得解析解^［31］，为

$$P_i\left(r,t\right)=A_i\mathrm{E}\mathrm{i}\left(-\frac{r^{\mathrm{2}}}{\mathrm{4}{\eta }_{i}t}\right)+B_i$$

（11）

式中：A_i、B_i为待定系数， Ei（x）为关于x的指数积分函数，$\mathrm{E}\mathrm{i}\left(x\right)=-{\int }_{-x}^{\mathrm{\infty }}\frac{{\mathrm{e}}^{-y}}{y}\mathrm{d}y$。由图2可得渗流场的边界条件、连续性条件、初始条件。

边界条件

$${\left.P_{\mathrm{1}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{0}}\\ t=t\end{array}}=P_{\mathrm{w}\mathrm{f}},{\left.P_{\mathrm{1}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}={\left.P_{\mathrm{2}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}},{\left.P_{\mathrm{2}}\left(r,t\right)\right|}_{\begin{array}{l}r=\mathrm{\infty }\\ t=t\end{array}}=P_{\mathrm{\infty }}$$

（12）

连续性条件

$$\frac{\mathrm{2}\mathrm{\pi }{k}_{\mathrm{1}}}{\mu }{\left.r\frac{\mathrm{d}{P}_{\mathrm{1}}\left(r,t\right)}{\mathrm{d}r}\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}=\frac{\mathrm{2}\mathrm{\pi }{k}_{\mathrm{2}}}{\mu }{\left.r\frac{\mathrm{d}{P}_{\mathrm{2}}\left(r,t\right)}{\mathrm{d}r}\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}$$

（13）

初始条件

$${\left.P_{\mathrm{1}}\left(r,\mathrm{0}\right)\right|}_{t=\mathrm{0}}=P_{\mathrm{\infty }}$$ 且r1= r0

（14）

将定解条件式（12）、（13）代入式（11），可确定非稳态孔压的分布为

$$P\left(r,t\right)=\left\{\begin{array}{l}\frac{\left(P_{\mathrm{w}\mathrm{f}}-P_{\mathrm{\infty }}\right)\left[k_{\mathrm{2}}{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{{\eta }_{\mathrm{1}}}-\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}\right)}\mathrm{E}\mathrm{i}\left(-\frac{r^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)+k_{\mathrm{1}}\right]}{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{{\eta }_{\mathrm{1}}}-\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}\right)}{k}_{\mathrm{2}}\left(\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)-\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)\right)+k_{\mathrm{1}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{2}}t}\right)},r_{\mathrm{0}}\le r\le r_{\mathrm{1}}\\ \frac{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{{\eta }_{\mathrm{1}}}-\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}\right)}{k}_{\mathrm{2}}\left(P_{\mathrm{\infty }}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)-P_{\mathrm{w}\mathrm{f}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)\right)+k_{\mathrm{1}}{P}_{\mathrm{w}\mathrm{f}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{2}}t}\right)}{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{{\eta }_{\mathrm{1}}}-\frac{\mathrm{1}}{{\eta }_{\mathrm{2}}}\right)}{k}_{\mathrm{2}}\left(\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)-\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{1}}t}\right)\right)+k_{\mathrm{1}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{2}}t}\right)}\mathrm{E}\mathrm{i}\left(-\frac{r^{\mathrm{2}}}{\mathrm{4}{\eta }_{\mathrm{2}}t}\right)+P_{\mathrm{\infty }},r>r_{\mathrm{1}}\end{array}\right.$$

（15）

式（15）自动满足初始条件，清晰展示了耦合系数η_i的影响。

比较式（2）和（10）可知，温度的控制方程和渗流形式上一致，可采用类似方法求解。温度场的边界条件、连续性条件、初始条件如下。

边界条件

$${\left.T_{\mathrm{1}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{0}}\\ t=t\end{array}}=T_{\mathrm{w}\mathrm{f}},{\left.T_{\mathrm{1}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}={\left.T_{\mathrm{2}}\left(r,t\right)\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}},{\left.T_{\mathrm{2}}\left(r,t\right)\right|}_{\begin{array}{l}r=\mathrm{\infty }\\ t=t\end{array}}=T_{\mathrm{\infty }}$$

（16）

连续性条件

$$\mathrm{2}\mathrm{\pi }{\lambda }_{\mathrm{1}}{\left.r\frac{\mathrm{d}{T}_{}\left(r,t\right)}{\mathrm{d}r}\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}=\mathrm{2}\mathrm{\pi }{\lambda }_{\mathrm{2}}{\left.r\frac{\mathrm{d}{T}_{\mathrm{2}}\left(r,t\right)}{\mathrm{d}r}\right|}_{\begin{array}{l}r=r_{\mathrm{1}}\\ t=t\end{array}}$$

（17）

初始条件

$${\left.T_{\mathrm{1}}\left(r,\mathrm{0}\right)\right|}_{\begin{array}{l}r=r_{\mathrm{0}}\\ t=\mathrm{0}\end{array}}=T_{\mathrm{\infty }}$$ 和 r1= r0

（18）

采用与渗流场相同的求解过程，获得温度的非稳态解答为

$$T\left(r,t\right)=\left\{\begin{array}{l}\frac{\left(T_{\mathrm{w}\mathrm{f}}-T_{\mathrm{\infty }}\right)\left[{\lambda }_{\mathrm{2}}{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{a_{\mathrm{1}}}-\frac{\mathrm{1}}{a_{\mathrm{2}}}\right)}\mathrm{E}\mathrm{i}\left(-\frac{r^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)+{\lambda }_{\mathrm{1}}\right]}{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{a_{\mathrm{1}}}-\frac{\mathrm{1}}{a_{\mathrm{2}}}\right)}{\lambda }_{\mathrm{2}}\left(\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)-\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)\right)+{\lambda }_{\mathrm{1}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{2}}t}\right)},r_{\mathrm{0}}\le r\le r_{\mathrm{1}}\\ \frac{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{a_{\mathrm{1}}}-\frac{\mathrm{1}}{a_{\mathrm{2}}}\right)}{\lambda }_{\mathrm{2}}\left(T_{\mathrm{\infty }}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)-T_{\mathrm{w}\mathrm{f}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)\right)+{\lambda }_{\mathrm{1}}{T}_{\mathrm{w}\mathrm{f}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{2}}t}\right)}{{\mathrm{e}}^{\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}t}\left(\frac{\mathrm{1}}{a_{\mathrm{1}}}-\frac{\mathrm{1}}{a_{\mathrm{2}}}\right)}{\lambda }_{\mathrm{2}}\left(\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)-\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}t}\right)\right)+{\lambda }_{\mathrm{1}}\mathrm{E}\mathrm{i}\left(-\frac{r_{\mathrm{1}}^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{2}}t}\right)}\mathrm{E}\mathrm{i}\left(-\frac{r^{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{2}}t}\right)+T_{\mathrm{\infty }},r>r_{\mathrm{1}}\end{array}\right.$$

（19）

式（19）自动满足初始条件。

水合物的分解取决于压力和温度。通过海水中纯水合物的相平衡试验数据^［32］拟合可得水合物的相平衡方程为

$$P=\frac{\mathrm{1}}{\mathrm{1}\mathrm{000}}{\mathrm{e}}^{\left(\mathrm{38.98}-\frac{\mathrm{8}\mathrm{533.8}}{\mathrm{T}+\mathrm{273.15}}\right)}$$

（20）

式中：P为压力，MPa；T为温度，℃。在水合物分解半径r= r₁处温度和压力满足水合物相平衡方程。将式（15）和式（19）代入式（20）中，即可求得分解域大小r₁。

2.3　力场解答

水合物分解将导致含水合物地层强度和刚度的显著降低，从而影响井壁安全。求解力场的控制方程（3）—（5），可得有效应力及位移通解为

$${\sigma }_{\mathrm{r}-i}^{}=\frac{(\mathrm{1}-\mathrm{2}{\nu }_i)}{(\mathrm{1}-{\nu }_i)r^{\mathrm{2}}}{\int }_{r_{i-\mathrm{1}}}^r{\alpha }_i{P}_i(r,t)r\mathrm{d}r+\frac{E_i}{(\mathrm{1}-\mathrm{2}{\nu }_i)(\mathrm{1}+{\nu }_i)}\left[C_i-(\mathrm{1}-\mathrm{2}{\nu }_i)\frac{D_i}{r^{\mathrm{2}}}\right]-{\alpha }_i{P}_i(r,t)$$

（21）

$${\sigma }_{\mathrm{\theta }-i}^{}=-\frac{(\mathrm{1}-\mathrm{2}{\nu }_i)}{(\mathrm{1}-{\nu }_i)r^{\mathrm{2}}}{\int }_{r_{i-\mathrm{1}}}^r{\alpha }_i{P}_i(r,t)r\mathrm{d}r+\frac{E_i}{(\mathrm{1}-\mathrm{2}{\nu }_i)(\mathrm{1}+{\nu }_i)}\left[C_i+(\mathrm{1}-\mathrm{2}{\nu }_i)\frac{D_i}{r^{\mathrm{2}}}\right]-\frac{{\nu }_i}{\mathrm{1}-{\nu }_i}{\alpha }_i{P}_i(r,t)$$

（22）

$$\begin{array}{l}{u}_{\mathrm{r}-i}^{}=\frac{(\mathrm{2}\nu -\mathrm{1})(\mathrm{1}+\nu )}{E_i(\mathrm{1}-\nu )r}{\int }_{r_{i-\mathrm{1}}}^r{\alpha }_i{P}_i(r,t)r\mathrm{d}r+\\ C_{i}r+\frac{D_i}{r}\end{array}$$

（23）

式中：C_i、D_i为待定系数。可见，孔隙压力通过应力中第1、3项的形式影响应力场，第2项是纯力场问题的通解。定解条件为

$${\left.{\sigma }_{\mathrm{r}-\mathrm{1}}^{}\right|}_{r=r_{\mathrm{0}}}={\alpha }_{\mathrm{1}}{P}_{\mathrm{w}\mathrm{f}}$$，$${\left.{\sigma }_{\mathrm{r}-\mathrm{1}}^{}\right|}_{r=r_{\mathrm{1}}}={\left.{\sigma }_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r=r_{\mathrm{1}}}$$，$${\left.u_{\mathrm{r}-\mathrm{1}}^{}\right|}_{r=r_{\mathrm{1}}}={\left.u_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r=r_{\mathrm{1}}}$$，$${\left.{\sigma }_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r\to \mathrm{\infty }}={\sigma }_{\mathrm{\infty }}-{\alpha }_{\mathrm{2}}{P}_{\mathrm{\infty }}$$

（24）

将式（24）代入式（21）—（23）可确定待定系数C_i、D_i，得到含水合物地层应力和位移。

实际工程中更关心钻井活动引起的增量位移，增量模型（图3a）边界条件可以由全量模型（图3b）减去初始状态模型（图3c）。增量模型边界和连续性条件可表示为

$${\left.{\sigma }_{\mathrm{r}}^{}\right|}_{r=r_{\mathrm{0}}}={\alpha }_{\mathrm{1}}{P}_{\mathrm{w}\mathrm{f}}-\left({\sigma }_{\mathrm{\infty }}-{\alpha }_{\mathrm{2}}{P}_{\mathrm{\infty }}\right)$$，$${\left.{\sigma }_{\mathrm{r}-\mathrm{1}}^{}\right|}_{r=r_{\mathrm{1}}}={\left.{\sigma }_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r=r_{\mathrm{1}}}$$，$${\left.u_{\mathrm{r}-\mathrm{1}}^{}\right|}_{r=r_{\mathrm{1}}}={\left.u_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r=r_{\mathrm{1}}}$$，$${\left.{\sigma }_{\mathrm{r}-\mathrm{2}}^{}\right|}_{r\to \mathrm{\infty }}=\mathrm{0}$$

（25）

将式（25）代入式（21）—（23）可确增量定解条件下的系数C_i、D_i，并得到增量位移分布规律。篇幅所限，在此不展示具体解答。

图3  增量位移计算的几何模型及其边界条件

Fig.3  Geometry and boundary conditions in incremental model of formation for incremental displacement calculation

3 模型验证与对比

为了验证模型数学推导的正确性，将本文得到的解析解与数值解进行对比。数值模型本构、各类参数、考虑因素、边界条件等均和解析模型一致。因温度场和渗流场具有相同形式的控制方程和边界条件，仅对渗流场和力场进行验证。数值模型采用商用软件COMSOL中的偏微分方程模拟渗流场，固体力学模块模拟力场，数值模型最大半径设置为100 m（为1 000r₀）以模拟无穷远边界，其余参数在数值与解析模型中均相同。参数依据实际工程参数范围进行选取，井壁半径r₀=0.1 m，时间t=24 h，Biot系数α₁=α₂=0.7，初始地应力σ_∞=15.6 MPa，初始孔压P_∞=15 MPa，过压、欠压钻井液压力P_wfo=17 MPa、P_wfu=13 MPa，地层弹性模量E₁=40 MPa、E₂=120 MPa，地层渗透率k₁=3×10^-15、k₂=1×10^-15，地层泊松比ν₁=ν₂=0.3，流体体积模量K_f1=10 MPa、K_f2=30 MPa，地层温度T_∞=8 ℃，钻井液温度T_wf=25℃，热传导系数λ₁=6 W·（m·℃）^-1、λ₂=3 W·（m·℃）^-1，地层密度ρ₁=ρ₂=2 000 kg∙m^-3，地层比热容c₁=c₂=2 500 J∙kg∙℃^-1，孔隙率ϕ₁=0.6、ϕ₂=0.3，水的黏滞系数μ=1×10^-3Pa∙s。图4a、4b、4c和4d分别为孔压、径向应力、环向应力和增量径向位移在t=24h时空间分布的解析与数值结果对比，如在r=21.4r₀处，孔压（相对于∆P=P_wf–P_∞）、径向应力、环向应力和增量径向位移的解析结果较数值解的误差分别为-1.99%、0.27%、0.33%和-5.75%，其误差原因主要为数值模型的网格问题。解析和数值结果均吻合良好，验证了解析模型的正确性。

图4  解析结果与数值结果、部分耦合结果的对比

Fig.4  Comparison of analytical, numerical, and partial coupling results

图4同时也给出了考虑渗流-力场全耦合结果与半耦合（仅考虑渗流对力场影响）结果的对比。考虑半耦合的渗流控制方程可仅从流体连续方程推导得到，极坐标下不考虑力场影响的非稳态孔压控制方程可表示为^［31］

$$\frac{k}{\mu }{\nabla }^{\mathrm{2}}P=\varphi C_{\mathrm{t}}\frac{\partial P}{\partial t}$$

（26）

式中：ϕ为含水合物地层孔隙度；C_t为地层压缩系数，C_t=1/K。

图4显示，与半耦合分析结果相比，考虑体变对渗流的影响后，过（欠）压钻井时孔压降低（升高）、应力升高（降低）、增量径向位移减小。因为考虑体变对渗流影响后，阻缓了钻井液对地层的侵入速率，导致孔压较初始孔压的变化幅度降低，地层变形减小。如过压钻井时，全耦合孔压结果与半耦合相比最多降低12.96%（相对于压力差ΔP = P_wf–P_∞，即变化幅度与ΔP的比值）；增量位移最大值减小41.13%，且全耦合位移最大值位置与半耦合相比更接近井壁。

为验证本模型在复杂问题中的适用性，将解析解与复杂数值模型^［9］结果对比（由于井壁稳定相关的试验和现场数据缺乏，未对比）。为消除不同相平衡方程对结果造成的偏差，解析解采用文献［9］中的分解半径值。数值模型基于南海神狐海域含水合物地区的测井数据，结合水合物模拟器TOUGH+HYDRATE及力场分析软件FLAC^3D实现温度-压力-化学-力场的耦合。数值计算考虑了水合物分解对土体力学特性影响，对井壁进行弹塑性分析（Mohr-Coulomb屈服准则），初始时地层孔压为P_∞=14.508 MPa，井壁处钻井液压力为P_wf=14.771 MPa，其余关键参数分别为：r₀=0.114 3 m，r₁=0.12 m，ν₁=ν₁=0.35，K_f1=10 MPa，K_f2=20 MPa，k₂=1×10^-14 m²，α_k=3，μ=1×10^-3 Pa∙s。解析解采用与数值解相同的参数。图5为本文孔压解析结果与文献［9］的数值模型结果及半耦合结果对比。可见解析结果与引文［9］数值模型结果吻合良好，孔压解析结果较引文［9］数值解大4.57%（相对于∆P=P_wf-P_∞，t=2 h，r=21.4r₀），误差主要由2个模型本身的差异（包括文献［9］数值模型考虑了化学场的影响、水合物分解产水及吸热等）引起。相较于传统的半耦合分析结果，全耦合解析模型更接近于符合实际情况的复杂数值模型结果，验证了本文提出的解析模型的适用性。

图5  解析模型的孔压与复杂数值模型^［9］结果对比

Fig.5  Comparison of pore pressure between analytical solution and numerical solution under complex conditions^[9]

4 参数分析

为定量分析水合物地层钻井过程中各因素对井壁稳定的影响，分析三场的时间和空间分布，并对钻井液压力、水合物分解引起的弹性模量劣化程度等参数对井壁稳定的影响进行讨论，并根据强度准则和应力分布进行稳定性分析。井壁失稳即为地层应力状态超过强度准则规定的应力状态。强度准则可以根据实际地层状态选用。由于深海含水合物土体刚度较大，地层性质接近岩石^［33］，强度准则采用Hoek-Brown （HB）准则，其表达式为^［34］

$${\sigma }_{\mathrm{1}}={\sigma }_{\mathrm{3}}+\sqrt[]{m{\sigma }_{\mathrm{c}}{\sigma }_{\mathrm{3}}+s{\sigma }_{\mathrm{c}}^{\mathrm{2}}}$$

（27）

式中：σ₁为最大主应力；σ₃为最小主应力；σ_c为岩石单轴抗压强度；m、s为岩石经验参数，可分别由$m=m_{\mathrm{0}}{\mathrm{e}}^{\left[\left({\mathrm{G}}_{\mathrm{S}\mathrm{I}}-\mathrm{100}\right)/\left(\mathrm{28}-\mathrm{14}\mathrm{d}\right)\right]}$， $s={\mathrm{e}}^{\left[\left({\mathrm{G}}_{\mathrm{S}\mathrm{I}}-\mathrm{100}\right)/\left(\mathrm{9}-\mathrm{3}\mathrm{d}\right)\right]}$确定， 其中G_SI为地质强度指数，表征岩石完整程度，m₀与组成岩石的矿物成分有关，d为干扰系数。令$f\left({\sigma }_{\mathrm{1}},{\sigma }_{\mathrm{3}}\right)={\sigma }_{\mathrm{1}}-\left({\sigma }_{\mathrm{3}}+\sqrt[]{m{\sigma }_{\mathrm{c}}{\sigma }_{\mathrm{3}}+s{\sigma }_{\mathrm{c}}^{\mathrm{2}}}\right)$为HB等效应力，则当f<0时，井壁处于安全状态；当f > 0时，井壁处于失稳状态。借鉴我国南海的水合物地层特点和试采数据，模型参数取值为：σ_c1=2 MPa、σ_c2=5 MPa^［35］，G_SI1=50、G_SI2=80^［33］，d₁=d₂=0，m₁=15、m₂=20^［33］ 。

4.1　渗流、温度、力场的时间和空间分布

分别针对过压钻井（P_wfo=17 MPa）和欠压钻井（P_wfu=13 MPa）工况进行研究，图6a、6b、6c、6d和6e 为时孔压、温度、径向应力、环向应力和增量径向位移在不同时刻沿径向分布规律。井壁压力随时间向井周扩散，压力变化区域（该区域孔压相对于原始地层压力发生变化）随之扩大。温度分布具有与孔压类似的变化规律，但温度的传播速度远小于孔压。如在时间t=1 d时（当前模型参数条件下），温度场的影响半径大约为10r₀而渗流场为70 r₀。

图6  不同时间下渗流、温度、力场和分解半径沿极径分布规律

Fig.6  Pore pressure, temperature, mechanical field, and the dissociated radius versus polar radius at different times

径向应力空间上呈现先减小后增大最后趋于初始值的规律，其最小值发生在分解锋面附近（r=r₁），说明分解造成的地层软化减小了2个域的相互作用。过压钻井时，径向应力最大值在井壁处；而欠压钻井时，最大值发生在邻近井壁的未分解域内。过压钻井时，径向应力随时间减小；欠压钻井时，径向应力随时间增大，井壁附近的径向和环向应力基本不随时间变化。由于分解域弹性模量降低，环向应力在分解锋面两端表现出不同的变化趋势，并在分解锋面处有跳跃增大点（如图6d）。增量位移呈现先增大后减小最后趋于零的趋势，其最大值出现在离井壁一定距离处，数值随时间增加明显，如t=1h到1 d时，过压钻井的井壁处位移和最大值位移分别增长22.73%和94.44%。增量位移在过压和降压时方向相反，且欠压时的井壁处位移（t=24h）比过压时大143.39%。图6f为过压钻井和欠压钻井时分解半径随时间变化规律。可见在钻井液的影响下，水合物分解半径随时间非线性增大，分解速度随时间减小^［36］，且过压时分解速度小于欠压时的分解速度。

4.2　钻井液压力对力场的影响

钻井液压力对保持井壁稳定非常重要。图7a、7b、7c和7d分别为径向应力、环向应力、增量位移在不同钻井液压力下沿径向分布规律及HB等效应力随钻井液压力变化规律。随着钻井液压力升高，井壁附近的径向应力增大，环向应力减小，而井壁及地层的位移随钻井液压力与地层初始孔压的压差（|P_wf–P_∞|）增大而显著增加，如过压钻井时，若钻井液压力P_wf从16 MPa增加到18 MPa，地层位移最大值将增加198.50%（如图7c）。图7d为不同钻井液压力下，基于Hoek-Brown准则的井壁处（r=r₀）和分解锋面处（r=r_1-，位于分解域内一侧；r=r₁₊，位于未分解域内一侧）的稳定性分析。无论过压或欠压钻井，井壁处均先于分解锋面处失稳，最危险点位于井壁处。最安全钻井液压力（此时f最小）与地层性质有关，水合物分解引起的地层劣化将降低最安全钻井液压力，如在分解锋面r=r₁处，未分解域内的最安全钻井液压力约为17 MPa，而在分解域内为13 MPa。在当前条件下（如图7d所示），钻井液压力P_wf约大于16.2 MPa或小于9.8MPa时，井壁分别发生拉伸或剪切破坏而失稳，与地层初始孔压差分别为5.2MPa和1.2MPa，即安全钻井液压力范围为9.8~16.2MPa，且相同压差下欠压钻井更有利于井壁稳定。因此，过高或过低的钻井液压力均会导致井壁失稳，为保持井壁稳定，需使钻井液压力与地层初始孔压的差值限定在一定范围内。

图7  不同钻井液压力下应力、位移和等效应力沿极径分布规律

Fig.7  Stress, displacement, and HB equivalent stress versus polar radius at different drilling fluid pressures

4.3　弹性模量劣化对力场的影响

水合物分解会导致地层刚度和强度的显著降低^［37-38］，降低其承载能力从而引发井壁失稳。如张旭辉等^［37］以粉细砂土和蒙古砂土作为沉积物骨架生成含水合物土体，试验表明水合物分解将使其强度降低至原来的1/7~1/9，Liu等^［38］的试验结果表明水合物分解可使强度和弹性模量分别最大降低89.8%和96.4%。定义分解域与未分解域弹性模量比α_E=E₁/E₂以表征水合物分解引起的地层刚度劣化程度。图8a、8b、8c和8d分别为径向应力、环向应力、增量径向位移在不同弹性模量比下沿径向分布规律及HB等效应力随弹性模量比α_E变化规律。可见随着弹性模量比降低（即分解域的弹性模量降低幅度增大），井壁附近径向应力减小，环向在分解域内减小而在未分解域内增大，径向位移显著增大且欠压钻井时增大更明显。当α_E=1/2→1/15时，井壁处位移增大192.85%。图8d为不同弹性模量比α_E下，基于Hoek-Brown准则的井壁处（r=r₀）和分解锋面处（r=r_1-，位于分解域内一侧）的稳定性分析。过压钻井时，f随α_E降低呈近似线性减小，当α_E约为0.5时，即分解域弹性模量降低50%，井壁发生失稳。欠压钻井时，f并不随α_E单调变化，而是随α_E的减小先减小后增大，大约在α_E=0.3时取得最小值（此时最安全），当α_E<0.15时井壁发生失稳。水和分解引起的地层弹性模量降低会显著影响井壁及地层稳定性，特别是在过压钻井时极易造成井壁失稳。

图8  不同地层弹性模量比下应力、位移和等效应力沿极径分布规律

Fig.8  Stress, displacement, and HB equivalent stress versus polar radius at different elastic modulus ratios

5 结论

针对深海水合物钻井在过压和欠压条件下的井壁稳定性，建立虑水合物分解、热传导、力场-渗流场全耦合作用下解析模型，并分析关键参数对井壁稳定性的影响，得出以下结论：

（1）与半耦合分析结果相比，考虑体变对渗流的影响后，过（欠）压钻井时孔压降低（升高）、应力升高（降低）、增量径向位移减小，如增量位移最多可减小41.13%。

（2）井壁失稳风险随时间升高，最危险位置在井壁处。井壁及地层的位移随钻井液压力与地层初始孔压的压差（| P_wf–P_∞|）增大而显著增加，过高或过低的钻井液压力均会导致井壁失稳。最安全钻井液压力与地层性质有关，水合物分解引起的地层劣化将降低最安全钻井液压力。

（3）水合物分解引起的地层弹性模量劣化会显著增大地层变形，在过压钻井时极易诱发井壁失稳，弹性模量降低50%即可导致井壁失稳。

推导的解析解答可以反映力场对渗流场的部分影响和深海水合物钻井时井壁稳定的力学机理。由于数学解析理论的限制，仅推导了弹性解答，且未考虑渗透率、孔隙率的变化，后续研究将考虑地层塑性变形及渗透率或孔隙率变化对渗流的影响。

作者贡献声明

黄佳佳：公式推导、参数分析、论文撰写。

蒋明镜：项目负责人，论文修改。

王华宁：项目负责人，研究思路指导、论文修改。