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

2023-08-02黄佳佳蒋明镜王华宁

同济大学学报（自然科学版） 2023年7期

黄佳佳，蒋明镜，2，王华宁

（1.同济大学土木工程学院，上海 200092；2.苏州科技大学土木工程学院，江苏苏州 215009；3.同济大学航空航天与力学学院，上海 200092）

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

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

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

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

1 力学模型

考虑水合物分解及非稳态渗流-力场全耦合作用下的深海水合物钻井过程中的井壁稳定问题。通常含水合物地层位于上下不透水层之间［19］，流体流动和热量传递可认为仅发生在水平面上。同时，钻井轴向尺寸远大于径向尺寸，储层所处地应力较大且其面内两方向应力近似相等［20］，该问题可简化为水平面内轴对称平面应变问题。钻井活动将改变地层温度、孔压并导致水合物分解，并引起地层力学特性、渗透率、热传导系数等发生显著变化。这些因素之间是相互作用、相互耦合的关系。与以往研究不同，本文将全面考虑渗流场-力场间的相互耦合以及水合物分解作用（三场关系如图1），获得更为接近实际情况的多场耦合井壁稳定解析模型。根据水合物是否分解将地层分为具有不同物理力学性质的2个区域（见图2）：分解区（r0＜r＜r1）和未分解区（r＞r1）［21-22］，其中r0、r1分别为井壁半径、分解半径。分解区内水合物完全分解，弹性模量、渗透率、热传导系数分别为E1、k1、λ1；未分解区内水合物未分解，弹性模量、渗透率、热传导系数分别为E2、k2、λ2。作为等效模型，地层视为均匀、各向同性，渗流、温度、力场边界条件如图2所示，Pwf、Twf分别为钻井液压力、温度，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］。由于温度值和梯度均较小，可以忽略温度应力和温度对渗流的影响。渗流控制方程为

式中：P为孔压；t为时间；εv为地层体积应变；k为地层渗透率；μ为流体黏滞系数；α为Biot系数，α=1-K/Ks，K为含水合物地层体积模量，K=E/［3（1-2ν）］，Ks为地层固体颗粒体积模量；M为Biot模量，，ν、νu分别为含水合物地层的排水、不排水泊松比，νu=，Ku为地层不排水体积模量，［26-27］；E、G分别为含水合物地层的弹性、剪切模量，G=E/［2（1+v）］；Kf为流体体积模量，本文流体指水和气体的混合物；φ为地层孔隙率。

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

式中：Ti为含水合物地层温度；ai为含水合物地层的热扩散系数，ai=λi/（ρichi）；λi、ρi和chi分别为含水合物地层的热传导系数、密度和比热容。

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

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

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

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

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

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

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

式中：Ai、Bi为待定系数， Ei（x）为关于x的指数积分式（15）自动满足初始条件，清晰展示了耦合系数ηi的影响。

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

边界条件

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

水合物的分解取决于压力和温度。通过海水中函数，Ei(x)=由图2可得渗流场的边界条件、连续性条件、初始条件。

边界条件

连续性条件

初始条件

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

连续性条件

初始条件

采用与渗流场相同的求解过程，获得温度的非稳态解答为纯水合物的相平衡试验数据［32］拟合可得水合物的相平衡方程为

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

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

2.3 力场解答

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

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

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

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

图3 增量位移计算的几何模型及其边界条件Fig.3 Geometry and boundary conditions in incremental model of formation for incremental displacement calculation

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

3 模型验证与对比

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

图4 解析结果与数值结果、部分耦合结果的对比Fig.4 Comparison of analytical, numerical, and partial coupling results

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

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

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

为验证本模型在复杂问题中的适用性，将解析解与复杂数值模型［9］结果对比（由于井壁稳定相关的试验和现场数据缺乏，未对比）。为消除不同相平衡方程对结果造成的偏差，解析解采用文献［9］中的分解半径值。数值模型基于南海神狐海域含水合物地区的测井数据，结合水合物模拟器TOUGH+HYDRATE及力场分析软件FLAC3D实现温度-压力-化学-力场的耦合。数值计算考虑了水合物分解对土体力学特性影响，对井壁进行弹塑性分析（Mohr-Coulomb屈服准则），初始时地层孔压为P∞=14.508 MPa，井壁处钻井液压力为Pwf=14.771 MPa，其余关键参数分别为：r0=0.114 3 m，r1=0.12 m，ν1=ν1=0.35，Kf1=10 MPa，Kf2=20 MPa，k2=1×10-14m2，αk=3，μ=1×10-3Pa·s。解析解采用与数值解相同的参数。图5为本文孔压解析结果与文献［9］的数值模型结果及半耦合结果对比。可见解析结果与引文［9］数值模型结果吻合良好，孔压解析结果较引文［9］数值解大4.57%（相对于ΔP=Pwf-P∞，t=2 h，r=21.4r0），误差主要由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］

式中：σ1为最大主应力；σ3为最小主应力；σc为岩石单轴抗压强度；m、s为岩石经验参数，可分别由m=确定，其中GSI为地质强度指数，表征岩石完整程度，m0与组成岩石的矿物成分有关，d为干扰系数。令f(σ1，σ3)=σ1-为HB等效应力，则当f＜0时，井壁处于安全状态；当f＞ 0时，井壁处于失稳状态。借鉴我国南海的水合物地层特点和试采数据，模型参数取值为：σc1=2 MPa、σc2=5 MPa［35］，GSI1=50、GSI2=80［33］，d1=d2=0，m1=15、m2=20［33］。

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

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

图6 不同时间下渗流、温度、力场和分解半径沿极径分布规律Fig.6 Pore pressure, temperature, mechanical field, and the dissociated radius versus polar radius at different times

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

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

钻井液压力对保持井壁稳定非常重要。图7a、7b、7c和7d分别为径向应力、环向应力、增量位移在不同钻井液压力下沿径向分布规律及HB等效应力随钻井液压力变化规律。随着钻井液压力升高，井壁附近的径向应力增大，环向应力减小，而井壁及地层的位移随钻井液压力与地层初始孔压的压差（|Pwf-P∞|）增大而显著增加，如过压钻井时，若钻井液压力Pwf从16 MPa增加到18 MPa，地层位移最大值将增加198.50%（如图7c）。图7d为不同钻井液压力下，基于Hoek-Brown准则的井壁处（r=r0）和分解锋面处（r=r1-，位于分解域内一侧；r=r1+，位于未分解域内一侧）的稳定性分析。无论过压或欠压钻井，井壁处均先于分解锋面处失稳，最危险点位于井壁处。最安全钻井液压力（此时f最小）与地层性质有关，水合物分解引起的地层劣化将降低最安全钻井液压力，如在分解锋面r=r1处，未分解域内的最安全钻井液压力约为17 MPa，而在分解域内为13 MPa。在当前条件下（如图7d所示），钻井液压力Pwf约大于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=E1/E2以表征水合物分解引起的地层刚度劣化程度。图8a、8b、8c和8d分别为径向应力、环向应力、增量径向位移在不同弹性模量比下沿径向分布规律及HB等效应力随弹性模量比αE变化规律。可见随着弹性模量比降低（即分解域的弹性模量降低幅度增大），井壁附近径向应力减小，环向在分解域内减小而在未分解域内增大，径向位移显著增大且欠压钻井时增大更明显。当αE=1/2→1/15时，井壁处位移增大192.85%。图8d为不同弹性模量比αE下，基于Hoek-Brown准则的井壁处（r=r0）和分解锋面处（r=r1-，位于分解域内一侧）的稳定性分析。过压钻井时，f随αE降低呈近似线性减小，当αE约为0.5时，即分解域弹性模量降低50%，井壁发生失稳。欠压钻井时，f并不随αE单调变化，而是随αE的减小先减小后增大，大约在αE=0.3时取得最小值（此时最安全），当αE＜0.15时井壁发生失稳。水和分解引起的地层弹性模量降低会显著影响井壁及地层稳定性，特别是在过压钻井时极易造成井壁失稳。