一种井壁坍塌渐进破坏过程的数值模拟方法

摘要

本发明公开一种井壁坍塌渐进破坏过程的数值模拟方法，包括：根据室内实验和测井资料确定目标井地层的岩石基本物性参数、岩石力学参数、地层地应力、地层孔隙压力、井筒压力；建立井周应力分布的流‑固耦合有限元数学模型；计算井周应力分布；计算损伤变量F；确定破坏区域，更新破坏区域的杨氏模量、孔隙度、渗透率；重复步骤进行迭代，直到损伤变量F≤0；绘制井壁失稳区域图，确定失稳区宽度Φb和深度rb。本发明克服了常规有限元模型更新几何模型形状来更新内边界的不足，本发明通过更新破坏区域地层参数杨氏模量、孔隙度和渗透率，令破坏范围内的单元令为等效消失，无法承受应力加载，从而达到更新内边界目的。

权利要求

1.一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，包括以下步骤：步骤S10、根据室内实验和测井资料确定目标井地层的岩石基本物性参数、岩石力学参数、地层地应力、地层孔隙压力、井筒压力；

步骤S20、建立井周应力分布的流‑固耦合有限元数学模型；

步骤S21、建立控制方程：

2

G▽u+(G+λ)▽divu‑α▽p＝0式中：G和λ是拉梅常数；k是多孔介质渗透率；μ是流体的粘度；u和p分别是多孔介质的位移和孔隙压力；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和岩石T

的体积模量；I＝[1,1,1,0,0,0]；D是弹性刚度矩阵；α是Biot系数；

步骤S22、使用伽辽金法有限元法逼近控制方程，可以得到控制方程的有限元求解格式；

其中，

T

M＝∫VBDBdV

B＝LNu

式中：M、H、S、C分别是弹性刚度、流动刚度、流动能力和耦合矩阵；B为应变和位移相关的应变矩阵；上标T为矩阵转置；Nu和Np分别为位移形函数和压力形函数；L微分算子；ρs和ρw分别为岩石和流体的密度；g为重力加速度； 和 分别对应边界Γ中的力边界条件和流量边界条件；qw为作用在边界上的压力流量；u和p分别为未知变量u和p的向量；ut和pt分别u p

为变量u和p的时间导数；f、f分别为节点载荷的向量和流体源汇的向量；

步骤S30、根据井周应力分布的流‑固耦合有限元数学模型计算井周应力分布；

步骤S40、根据井周应力分布计算损伤变量F；

F1＝‑σ3‑fto

其中，

I1＝σ1+σ2+σ3

式中，I1为应力第一不变量；J2为应力偏张量第二不变量；为岩石内摩擦角；C为内聚力；fto是岩石抗拉强度；α0和k0为材料常数，与岩石内聚力和内摩擦角有关；F1为最大拉应力状态函数；F2为Drucker‑Prager状态函数；F为损伤变量；σ1为最大主应力；σ2为中间主应力；

σ3为最小主应力；

步骤S50、根据损伤变量F确定破坏区域，更新破坏区域的杨氏模量、孔隙度、渗透率；

E＝(1‑F)E0

式中：E、E0分别为岩石损伤前后的杨氏模量；φ、φ0分别为岩石损伤前后的孔隙度；K、K0分别为岩石损伤前后的渗透率；F为损伤变量；

步骤S60、重复步骤S30至步骤S50进行迭代，直到损伤变量F≤0；此时无新的损伤区出现，井眼最终趋于稳定，绘制井壁失稳区域图，确定失稳区宽度Φb和深度rb。

2.根据权利要求1所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，所述步骤S10中所述岩石基本物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力；所述地层地应力包括最大水平地应力、最小水平地应力。

3.根据权利要求2所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，所述步骤S30的具体过程为：根据井周应力分布流‑固耦合有限元数学模型建立井壁坍塌渐进破坏过程模拟的有限元模型，在分别对有井壁坍塌渐进破坏过程模拟的有限元模型赋值材料参数、施加边界条件和划分有限元网格，最后再计算井周应力分布，得到最大主应力σ1、中间主应力σ2和最小主应力σ3。

4.根据权利要求3所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，根据井周应力分布流‑固耦合有限元数学模型建立井壁坍塌渐进破坏过程模拟的有限元模型包括：根据平面应变、轴对称条件和井周应力分布流‑固耦合有限元数学模型，建立井壁坍塌渐进破坏过程模拟的有限元模型，其井筒半径为R，有限元模型为50R×50R的几何模型，以排除边界效应的影响，并为了减少计算量。

5.根据权利要求3所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，对有限元模型赋值材料参数包括：对有限元模型求解域赋值材料参数，设定求解域的岩石基本物性参数和岩石力学参数。

6.根据权利要求3所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，所述边界条件包括：内边界施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp。

7.根据权利要求3所述的一种井壁坍塌渐进破坏过程的数值模拟方法，其特征在于，所述步骤30中由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，远离井筒，网格逐渐变粗，以减少要计算的单元数量。

说明书

一种井壁坍塌渐进破坏过程的数值模拟方法

技术领域

[0001] 本发明涉及一种井壁坍塌渐进破坏过程的数值模拟方法，属于石油钻井工程、岩石力学与工程领域。

背景技术

[0002] 井壁失稳是由于应力集中和井眼连续破坏导致最初圆形井眼的扩展。井壁失稳是一个渐进破坏的过程，如果破坏程度较低，则称为井壁崩落；过多的井壁破坏发生后，则导
致井壁垮塌。井壁崩落和井壁垮塌是预测地应力的重要参考指标，经过实验和理论证明，直
井的井壁崩落的方向与最小水平地应力方向一致。利用井壁崩落去预测地应力的状态的优
势是它们广泛存在井眼中，以及可以使用测井技术出检测它们的位置、尺寸和崩落方位。在
深井中，水力压裂很难实现地应力的测量，而井眼崩落是唯一可以得到地层应力状态信息。

[0003] 为了模拟井壁渐进破坏失稳，国内外学者已经提出了大量的模型，例如，实验：通过真三轴实验，观察井眼的崩落或破坏的形状，由于实验装置尺寸较小，结果通常受到边界
效应的影响，并且花费需要大量的人力和物力；解析解：应用弹性理论、弹塑性理论和孔弹
塑性理论计算井周应力分布，结合准则预测地层岩石的破坏，解析解具有输入参数少和便
易计算的优点，但仅适用于预测井壁崩落的开始和简单的几何模型；离散元模型：把岩石离
散为单个独立的颗粒，研究颗粒之间的破坏，其局限性是由于它需要过高的计算机能力，颗
粒之间的接触行为的精确描述。深度学习模型：利用AI算法高度线性/非线性拟合地应力与
井壁稳定崩落形状关系，不足处是需要大量的数据集来训练AI算法，不能预测崩落或垮塌
区的裂纹萌生、扩展、崩落形成；为此有限元是便易实现的方法，有限元模型：可以计算任意
井眼形状的应力分布，根据准则确定破坏区域，更新几何内边界。但对于更新几何边界的方
法往往难以实现。

[0004] 为此，本发明提出一种井壁坍塌渐进破坏过程的数值模拟方法。采用连续损伤理论，更改已破坏区域的力学参数，视为该区域等效挖去，模拟了井壁渐进破坏发展的四个阶
段：裂纹萌生、扩展、崩落形成和稳定。

发明内容

[0005] 本发明主要是克服现有技术中的不足之处，提出一种井壁坍塌渐进破坏过程的数值模拟方法；该方法采用连续损伤理论，更改已破坏区域的力学参数，视为该区域等效挖
去，从而模拟了井壁渐进破坏发展的四个阶段:裂纹萌生、扩展、崩落形成和稳定。

[0006] 本发明解决上述技术问题所提供的技术方案是：一种井壁坍塌渐进破坏过程的数值模拟方法，包括以下步骤：

[0007] 步骤S10、根据室内实验和测井资料确定目标井地层的岩石基本物性参数、岩石力学参数、地层地应力、地层孔隙压力、井筒压力；

[0008] 步骤S20、建立井周应力分布的流‑固耦合有限元数学模型；

[0009] 步骤S30、根据井周应力分布的流‑固耦合有限元数学模型计算井周应力分布；

[0010] 步骤S40、根据井周应力分布计算损伤变量F；

[0011] 步骤S50、根据损伤变量F确定破坏区域，更新破坏区域的杨氏模量、孔隙度、渗透率；

[0012] 步骤S60、重复步骤S30至步骤S50进行迭代，直到损伤变量F≤0；此时无新的损伤区出现，井眼最终趋于稳定，绘制井壁失稳区域图，确定失稳区宽度Φb和深度rb。

[0013] 进一步的技术方案是，所述步骤S10中所述岩石基本物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力；所述地层地应力
包括最大水平地应力、最小水平地应力。

[0014] 进一步的技术方案是，所述步骤S20的具体建立过程为：

[0015] 步骤S21、建立控制方程：

[0016]

[0017]

[0018] 式中：G和λ是拉梅常数；k是多孔介质渗透率；μ是流体的粘度；u和p分别是多孔介质的位移和孔隙压力；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和
T
岩石的体积模量；I＝[1,1,1,0,0,0]；D是弹性刚度矩阵；α是Biot系数；

[0019] 步骤S22、使用伽辽金法有限元法逼近控制方程，可以得到控制方程的有限元求解格式；

[0020]

[0021] 其中，

[0022] M＝∫VBTDBdV

[0023]

[0024]

[0025]

[0026]

[0027] B＝LNu

[0028]

[0029]

[0030] 式中：M、H、S、C分别是弹性刚度、流动刚度、流动能力和耦合矩阵；B为应变和位移相关的应变矩阵；上标T为矩阵转置；Nu和Np分别为位移形函数和压力形函数；L微分算子；ρs
和ρw分别为岩石和流体的密度；g为重力加速度； 和 分别对应边界Γ中的力边界条件
和流量边界条件；qw为作用在边界上的压力流量；u和p分别为未知变量u和p的向量；ut和pt
u p
分别为变量u和p的时间导数；f、f分别为节点载荷的向量和流体源汇的向量。

[0031] 进一步的技术方案是，所述步骤S30的具体过程为：根据井周应力分布流‑固耦合有限元数学模型建立井壁坍塌渐进破坏过程模拟的有限元模型，在分别对有井壁坍塌渐进
破坏过程模拟的有限元模型赋值材料参数、施加边界条件和划分有限元网格，最后再计算
井周应力分布，得到最大主应力σ1、中间主应力σ2和最小主应力σ3。

[0032] 进一步的技术方案是，根据井周应力分布流‑固耦合有限元数学模型建立井壁坍塌渐进破坏过程模拟的有限元模型包括：根据平面应变、轴对称条件和井周应力分布流‑固
耦合有限元数学模型，建立井壁坍塌渐进破坏过程模拟的有限元模型，其井筒半径为R，模
型为50R×50R的几何模型，以排除边界效应的影响，并为了减少计算量。

[0033] 进一步的技术方案是，对井壁坍塌渐进破坏过程模拟的有限元模型赋值材料参数包括：对井壁坍塌渐进破坏过程模拟的有限元模型求解域赋值材料参数，设定求解域的岩
石基本物性参数和岩石力学参数。

[0034] 进一步的技术方案是，所述边界条件包括：内边界施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp。

[0035] 进一步的技术方案是，所述步骤30中由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，远离井筒，网格逐渐变粗，以减少要计算的单元数量。

[0036] 进一步的技术方案是，所述步骤S40中的计算公式如下：

[0037]

[0038] F1＝‑σ3‑fto

[0039]

[0040] 其中，

[0041]

[0042] I1＝σ1+σ2+σ3

[0043]

[0044]

[0045] 式中，I1为应力第一不变量；J2为应力偏张量第二不变量；为岩石内摩擦角；C为内聚力；fto是岩石抗拉强度；α0和k0为材料常数，与岩石内聚力和内摩擦角有关；F1为最大拉
应力状态函数；F2为Drucker‑Prager状态函数；F为损伤变量；σ1为最大主应力；σ2为中间主
应力；σ3为最小主应力。

[0046] 进一步的技术方案是，所述步骤S50中的计算公式为：

[0047] E＝(1‑F)E0

[0048]

[0049]

[0050] 式中：E、E0分别为岩石损伤前后的杨氏模量；φ、φ0分别为岩石损伤前后的孔隙度；K、K0分别为岩石损伤前后的渗透率；F为损伤变量。

[0051] 本发明具有以下有益效果：本发明克服了常规有限元模型更新几何模型形状来更新内边界的不足，本发明通过更新破坏区域地层参数杨氏模量、孔隙度和渗透率，令破坏范
围内的单元令为等效消失，无法承受应力加载，从而达到更新内边界目的。

附图说明

[0052] 图1为实施流程图；

[0053] 图2为四分之一几何模型图；

[0054] 图3为施加边界条件结果图；

[0055] 图4为有限元网格划分图；

[0056] 图5为每次迭代步下蹦落区域图；

[0057] 图6为每次迭代步下崩落的形状数值图。

具体实施方式

[0058] 下面结合实施例和附图对本发明做更进一步的说明。

[0059] 如图1所示，本发明的一种井壁坍塌渐进破坏过程的数值模拟方法，包括以下步骤：

[0060] 步骤S10、根据室内实验和测井资料，确定井壁坍塌渐进破坏过程模拟目标井地层的岩石基本物性参数、岩石力学参数、地层地应力、地层孔隙压力、井筒压力等基础参数；

[0061] 其中所述岩石基本物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力；所述地层地应力包括最大水平地应力、最小水平
地应力；

[0062] 步骤S20、建立井周应力分布的流‑固耦合有限元数学模型；

[0063] 具体过程如下：

[0064] 步骤S21、建立控制方程：

[0065]

[0066]

[0067] 式中，G和λ是拉梅常数；k是多孔介质渗透率；μ是流体的粘度；u和p分别是多孔介质的位移和孔隙压力；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和
T
岩石的体积模量；I＝[1,1,1,0,0,0]；D是弹性刚度矩阵；α是Biot系数；

[0068] 步骤S21、使用伽辽金法有限元法逼近控制方程，可以得到控制方程(1)和(2)的有限元求解格式：

[0069]

[0070] 其中，

[0071] M＝∫VBTDBdV‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑(4)

[0072]

[0073]

[0074]

[0075] 式中，M、H、S、C分别是弹性刚度、流动刚度、流动能力和耦合矩阵；B为应变和位移相关的应变矩阵；上标T为矩阵转置；Nu和Np分别为位移形函数和压力形函数；u和p分别为未
u p
知变量u和p的向量；ut和pt分别为变量u和p的时间导数；f、f 分别为节点载荷的向量和流
体源汇的向量；

[0076]

[0077] B＝LNu‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑(9)

[0078]

[0079]

[0080] 式中，L微分算子；ρs和ρw分别为岩石和流体的密度；g为重力加速度； 和 分别对应边界Γ中的力边界条件和流量边界条件；qw为作用在边界上的压力流量；

[0081] 步骤S30、根据井周应力分布流‑固耦合有限元数学模型，建立井壁坍塌渐进破坏过程模拟的有限元模型，对有限元模型赋值材料参数、施加边界条件和划分有限元网格，计
算井周应力分布；

[0082] 步骤S31、根据平面应变和轴对称条件、井周应力分布流‑固耦合有限元数学模型，建立井壁坍塌渐进破坏过程模拟的有限元模型，其井筒半径为R，整个有限元模型为50R×
50R的几何模型，以排除边界效应的影响，并为了减少计算量；

[0083] 步骤S32、对有限元模型求解域赋值材料参数，设定求解域的基础物性参数和岩石力学参数，包括孔隙度、渗透率和Biot系数、杨氏模量、泊松比、内摩擦角、内聚力等基础参
数；

[0084] 步骤S33、对有限元模型施加边界条件，其中，内边界(井筒)施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp；

[0085] 步骤S34、对有限元模型划分有限元网格，由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，特别是最小水平地应力方向进行网格加密，远离井
筒，网格逐渐变粗，以减少要计算的单元数量；

[0086] 步骤S35、最后对有限元模型的网格模型进行求解，计算井周应力分布，从而得到最大主应力σ1、中间主应力σ2和最小主应力σ3；

[0087] 步骤S40、根据步骤S30得到的应力计算损伤变量F；

[0088] 步骤S41、确定地层岩石破坏时的应力状态函数：最大拉应力状态函数和Drucker‑Prager应力状态函数。

[0089] F1＝‑σ3‑fto‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑(12)

[0090]

[0091] 其中，

[0092]

[0093] I1＝σ1+σ2+σ3‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑(15)

[0094]

[0095]

[0096] 式中，I1为应力第一不变量；J2为应力偏张量第二不变量；为岩石内摩擦角；C为内聚力；fto是岩石抗拉强度；α0和k0为材料常数，与岩石内聚力和内摩擦角有关；F1为最大拉
应力状态函数；F2为Drucker‑Prager状态函数。

[0097] 步骤S42、根据最大拉应力状态函数F1和Drucker‑Prager状态函数F2，计算地层岩石的损伤变量F；

[0098]

[0099] 式中，F为损伤变量；当F＝0对应无损伤单元状态，F＝1对应完全损伤状态，即岩石发生断裂或破坏损伤；F1≥0表示地层岩石发生拉伸破坏损伤；F2≥0表示地层岩石发生剪切
破坏损伤；dF1>0和dF2>0分别为两种损伤后的继续加载状态，可引起损伤变量的增加；当dF1
<0和dF2<0表示卸载状态，不产生新的损伤，损伤变量保持上一个加载步(时间步)的数值；

[0100] 步骤S50、根据损伤变量F确定破坏区域，更新破坏区域地层参数杨氏模量、孔隙度、渗透率，使得破坏区域内的单元丧失承载应力能力；

[0101] 步骤S51、根据各向同性弹性损伤理论，单元的杨氏模量会随着损伤过程逐渐退化，损伤后的杨氏模量可表示为：

[0102] E＝(1‑F)E0‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑(19)

[0103] 式中：E、E0分别为岩石损伤前后的杨氏模量；由于此处采用的是弹性损伤理论，当F＝1时岩石发生破坏性损伤，为了消除数值计算带来的误差问题，假定F＝1的单元体杨氏
模量由万分之一的初始杨氏模量(0.0001E0)代替，即更新F＝1的单元体杨氏模量为
0.0001E0；

[0104] 步骤S52、受应力集中破坏的岩石受钻井液循环而带离原有位置，破坏区域应等效为完全由钻井液支撑，即破坏区域被视为新的内边界；因此，有必要加入孔隙度和渗透率与
岩石单元损伤的函数，其物理场参数如式(20)和(21)所示：

[0105]

[0106]

[0107] 对于破坏区域，渗透率应无限大，对于数值模拟无法输入无限大这个变量，对于渗透率大于4达西对渐进破坏结果无明显影响。因此，后续数值模拟时，发生井壁破坏时，渗透
率至少取值4达西；

[0108] 步骤S60、重复步骤S30至步骤S50进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，实现井壁坍塌渐进破坏过程的数值模拟，并绘制井壁失稳区域
图，确定失稳区宽度Φb和深度rb。

[0109] 实施例

[0110] 步骤S10、根据室内实验和测井资料，确定井壁坍塌渐进破坏过程模拟目标井地层的岩石基本物性参数、岩石力学参数、地层地应力、地层孔隙压力、井筒压力等基础参数。其
参数如表1所示：

[0111] 表1基础参数

[0112]

[0113] 步骤S20、建立流‑固耦合的数学模型；

[0114] 步骤S30、根据步骤S20的井周应力分布流‑固耦合有限元数学模型，建立井壁坍塌渐进破坏过程模拟的有限元模型，对有限元模型赋值材料参数、施加边界条件和划分有限
元网格，计算井周应力分布；

[0115] 步骤S40、根据步骤S30得到的应力计算损伤变量F；

[0116] 步骤S50、根据步骤S40的损伤变量F确定破坏区域，更新破坏区域地层参数杨氏模量、孔隙度、渗透率，使得破坏区域内的单元丧失承载应力能力；

[0117] 步骤S60、重复步骤S30至步骤S50进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，绘制井壁失稳区域图，确定失稳区宽度Φb和深度rb。

[0118] 图5和图6分别为每个迭代步得到的井壁崩落区域和对应的失稳区宽度Φb和深度rb。不难看出：随着迭代步数的增加，失稳区的宽度Φb在第一个迭代步后保持不变，而失稳
区的深度rb随着迭代步的增加而不断增加，迭代6步后达到稳定状态，说明井壁坍塌渐进破
坏过程主要是在失稳区的深度rb上逐渐发展和演化的，这与现场实际情况和认识是吻合
的。

[0119] 以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围
内，当可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未
脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、
等同变化与修饰，均仍属于本发明技术方案的范围内。