井壁坍塌分析方法 技术领域 本发明涉及石油钻井工程技术领域， 特别涉及钻井井壁失稳分析， 是一种非均匀 地应力条件下应变软化塑性地层井壁稳定性分析的方法。
背景技术 目前， 在石油工程中普遍应用弹性力学结合强度准则方法解释井壁坍塌机制和预 测坍塌压力。比如， 被认为是世界上最好的井壁稳定分析 GMI 系统， 其应用的井壁坍塌模型 正是基于弹性力学结合强度准则法。最新的井壁坍塌机理的理论研究表明， 井壁坍塌实质 上是一种极值点型失稳现象。 通过确定井壁压力和井壁位移的平衡路径曲线的极值即可确 定力学稳定性意义下的坍塌压力。当然， 理论解作了原场地应力是各向同性均匀应力和地 层满足 Tresca 准则的假设， 这远离了实际情况。还有重要一点， 在分析井壁坍塌问题时仅 考虑岩石的理想弹塑性或者弹脆性 - 塑性特性， 而大部分岩石的应力 - 应变曲线在峰值后 下降， 称为应变软化， 所以不考虑岩石的应变软化特性也是不符合实际的。
国外现有 Ansys， Flac， Abaqus， Adina 等有限元计算软件， 但是缺乏力学意义上的 稳定性分析功能， 没有将应变软化模型和稳定性分析结合， 而且应用领域没有扩展到井壁 坍塌分析。国内岩土方面计算系统比较突出的是 “理正岩质边坡稳定分析系统” ， 所采用的 也是传统的极限平衡法。 未见报道应用应变软化模型和力学稳定性方法开发井壁坍塌分析 系统。为了解决非均匀地应力条件下应变软化塑性地层井壁失稳问题， 本发明在应变软化 理论模型的基础上， 应用应变软化模型和力学稳定性方法开发 WellCA 井壁稳定分析系统。
井壁稳定问题是世界性难题， 在世界范围内， 每年用于处理井壁失稳的费用高达 数亿美元， 而我国许多油田也存在井壁失稳复杂情况。 近年来， 我国在西部塔里木和准噶尔 等盆地进行了多口深井、 超深井的钻探施工。这些地区山前高陡构造高地应力引起的井眼 稳定问题较为突出， 特别是深井的井眼稳定受山前构造带地应力影响， 问题更为严重。
中国专利授权公告号 ： CN1239920C， 提供了一种 “利用地震层速度钻前预测坍塌压 力与破裂压力的方法” 。 包括下列步骤 ： 1) 将待钻井与多个相邻已钻井的地震层速度进行数 据相关分析， 并确定相关系数大于 0.75 的已钻井为待钻进具有相似构造的已钻井 ； 2) 利用 声波时差、 自然伽玛、 密度等测井数据序列， 对已钻井全井段进行分层， 求出每一层用于表 征一定厚度且岩性相似地层的平均声波速度、 自然伽玛和地层密度 ； 3) 利用每一层的平均 声波速度、 自然伽玛和地层密度来确定已钻井的坍塌压力与破裂压力 ； 4) 根据已钻井的坍 塌压力与破裂压力和测井分层速度建立测井模型 ； 5) 建立已钻井的层速度钻前预测模型 ； 6) 将待钻井的地震层速度代入步骤 5) 中的模型， 获得待钻井的坍塌压力与破裂压力。
发明内容
本发明的目的是 ： 提供一种井壁坍塌分析方法， 考虑岩石的应变软化性质， 使其更 接近客观实际， 应用于钻井井壁稳定性理论研究和建立井壁坍塌模型， 效解决钻井的井壁 坍塌问题。本发明采用的技术方案是 ：
井壁坍塌分析方法， 包括 ： 采集并输入参数， 初应力处理， 第一次应力释放， 第二次 应力释放第一步， 第二次应力释放第二步， 特征值评价， 广义力和广义位移， 平衡路径曲线， 坍塌压力处理， 输出坍塌压力以及塑性区分布图。
(A) 采集并输入参数 ： 首先， 采集通过相邻已钻井的测井资料和岩石力学实验得 到井壁岩石材料参数， 包括弹性模量、 泊松比、 内摩擦角、 内聚力、 软化系数、 屈服函数类型、 内变量类型、 水平最大主应力、 水平最小主应力和深度参数， 将井壁岩石材料参数输入计算 程序中 ； 输入控制参数， 控制参数包括弧长、 输出步、 切线塑性剪切模量、 残余强度。控制参 -5 -6 数中的弧长范围 10 ～ 10 m， 输出步范围 100 ～ 200， 控制参数中的切线塑性剪切模量、 残 余强度通过岩石力学实验得到。
(B) 初应力处理 ： 建立井壁几何模型， 设置井筒半径， 内外半径比 ； 划分有限元计 算网格单元， 采用平面应变等参单元和可描述远场的无限区域单元， 在最外一层设置成无 限单元它的尺度在规定的方向上延伸至无穷远， 起到模拟远处区域的作用， 无限区域元的 位移插值函数与坐标插值函数有不同的表达形式， 位移在无限远处为零， 通过衰减函数的 引入来实现。
衰减函数采用负指数函数 ： 剖分节点和单元， 施加边界条件， 施加情况为 2 个直边法向约束、 径向自由。根据网格节点坐标， 由式 (3) 计算初始应力场各单元高斯点的应力
(C) 第一次应力释放 ： 采用 “释放载荷”的方法来模拟实际钻井过程， 求扰动 公式 (5) 为保留的场 A。首先， 根据公式 (4) 计算第一次释放部分应力分量应力分量。这样的应力释放方案， 等效于式 (3) 的应力完全释放而同时施加钻井液压力
(4)在 r ＝ a，按弹性或弹 - 理想塑性计算， 得扰动场 A， σ′，u′， 通过第一应力释放， 得到总场， 这个总场是指初始场加扰动场 A ： 0 0 0
σ +σ′→ σ ； a′→ a (6)
(D) 第二次应力释放第一步 ： 逐渐释放保留的应力分量， 也就是模拟钻井液压力 逐渐降低， 计算扰动应力场 B。引入载荷参数 λ， 在 r ＝ a，
τ′ θ ＝ 0。第一步用弹性方法确定弹性阶段的载荷参数增量 Δλ1， Δλ1 → λe， 求出相应的 弹性扰动场 Δa′， Δσ′， 得到总场， 这个总场是指初始场加扰动场 A 再加弹性扰动场 ：σ0+Δσ′→ σ0 ； a0+Δa′→ a0 (7)
(E) 第二次应力释放第二步 ： 由于岩石在高温高压条件下， 呈现应变软化塑性性 质， 所以从第二步开始 (m ＝ 1)， 指定微弧长增量 Δsm， 按弹 - 软化塑性， 用弧长法迭代求解。 全过程曲线 τ-γ 采用三线性形式， 下降段斜率是一个负的常数 GT， 称为切线剪切模量， 而 相应的软化曲线则为二线性形式， 下降段斜率为 Gp ＝ GGT/(G-GT)。 由于采用三线性模型， 在 D-P 准则中 k 就是剪切强度 τs， 它也是初始的屈服应力 τs(0)。设刚刚进入残余阶段的塑 性应变为 那么二线性的软化曲线斜率为
剪切强度 τs 和塑性功 wp 的关系为 ：在对岩石等应变软化材料组成的结构应用有限元方法进行非线性分析时， 一旦结 构接近于极限承载能力或者以后， 用一般的计算方法计算分析的稳定性与收敛性相当差以 至于无法进行下去。弧长法通过同时约束荷载和位移向量， 能较好地计算邻近极值点结构 的反应和下降段问题， 可以比较好地解决极值点穿越问题。
对于一个给定的迭代步 n 已知解待求的下一迭代步 n+1 的解记为可有如下的表达式 ：
于是， 令可给出位移解的增量 Δan 和载荷参数的增量 Δλn 之间的关系
其中
将式 (12) 代入约束方程 可得一个关于 Δλn 的二次代数方程 其中
从方程 (15) 可以得到 Δλn 的两个根， 分别记为 (Δλn)1 和 (Δλn)2。为了能够 跟踪解的路径得到正确的结果， 只选取一组解， 它应使解沿平衡曲线往前走， 而不是往后返 n T 回， 即取 (Δa ) (am-am-1) 为一个最大值的解， 得扰动场 B， Δλm， Δam， Δσm， Δκm， 求出总 场， 这个总场是指初始场加扰动场 A 再加扰动场 B ：
λm+Δλm → λm+1
(17)(F) 特征值评价 ： 结构的切线总刚矩阵不正定才是结构不稳定的充分和必要条 的最小特征值件。因此， 引入了计算总刚矩阵特征值的算法 - 子空间迭代法。求(μ1)m， 根据得到的最小特征值 (μ1)m 进行评价， 当最小特征值 (μ1)m ＞ 0， 时， 如果 m ＜ M， 则 m+1 → m， 即 m 加上 1 以后进入第二次应力释放第二步， 否则， 得到的即为总场 ； 当最小特 征值 (μ1)m ＜ 0 时， 进行如下计算。
λm → λN， λm-1 → λN-1 (μ1)m → (μ)N， (μ1)m-1 → (μ)N-1 σm → σN， σm-1 → σN-1 κm → κN， κm-1 → κN-1 使用线性插值计算临界扰动应力场 ：(18)计算总场， 这个总场是指初始场加扰动场 A 再加临界扰动场 ： Δλcr+λm → λcr Δσcr+σm → σcr (20)Δκcr+κm → κcr
(G) 广义力和广义位移 ： 根据弧长法每步骤计算出的载荷因子 λ， 按照两次解除 应力的公式通过积分即可求出当步广义力及广义位移。
在计算总场时， 外力功的变分为
如果广义力可定义为那么广义位移的变分为第 m 个载荷步的增量广义位移为总广义位移可累加得到 ：
Um+1 ＝ Um+ΔUm (25)
(H) 平衡路径曲线 ： 平衡路径曲线是由广义力和广义位移为坐标轴绘制的曲线。 根据广义力和广义位移的输出结果， 绘制平衡路径曲线。
(I) 坍塌压力处理 ： 用 P 和 u 分别代表系统的广义力和广义位移， 那么在稳定分支 上有 δPδu ＞ 0， 在不稳定分支上有 δPδu ＜ 0， 稳定分支和不稳定分支的交汇点是临界 点， 在该点 δPδu ＝ 0， 该点对应于平衡稳定性的临界状态， 该点处的广义力称为临界力， 它确定了一个结构的承载能力和保持整体性的能力。在平衡路径曲线上找出临界力， 根据 广义力的公式 (22) 确定对应的应变软化极限载荷因子 λcr。
根据得到的应力总场、 弹性极限载荷因子 λe 以及平衡路径曲线输出的 λcr， 进入 坍塌压力处理， 得到弹性坍塌压力和应变软化塑性坍塌压力 ：
(12)(J) 输出坍塌压力 qcr， 应变软化极限载荷因子 λcr 对应的塑性区分布图， 确定该井 的防塌钻井液的密度。 符号说明 ： - 第 i 个节点的位移插值函数 ； Ni- 第 i 个节点的坐标插值函数 ； - 第一次解除的径向 -第 m ； - 初始场径向应力， MPa ； - 初始场切向 L- 衰减特征长度， m； ri- 第 i 个节点位移， 应力， MPa ； σH- 最大水平主应力， Mpa ； σh- 最小水平主应力， Mpa ； Mpa ； θ- 与最大水平地应力的夹角，° ； 应力， Mpa ； - 第一次解除的切向应力，
一次保留的应力分量， Mpa ； λ- 载荷因子， 无量纲 ； r- 距井眼中心距离， m； a- 井眼半径， m； 0 Mpa ； σ′ - 扰动场应力， Mpa ； σ′ r- 扰动场径向应力， Mpa ； τ′ θ- 扰 σ - 初始场应力， 动场切向应力， Mpa ； Δλ1- 弹性阶段的载荷参数增量， 无量纲 ； λe- 弹性极限载荷参数， 无 量纲 ； GT 为切线剪切模量， MPa ； Gp- 切线塑性剪切模量， MPa ； τs(0)- 初始的屈服应力， Mpa ； p τs- 剪切强度， MPa ； w - 塑性功， MPa ； 为残余阶段的塑性应变，％ ； M- 弧长法微弧长增量 总量， 无量纲 ； Δsm 微弧长增量， m； Δan- 位移增量向量， m； Δλn- 载荷增量系数， 无量纲 ； - 切线刚度矩阵 ； R- 载荷力向量， MPa ； ψm- 不平衡力向量， MPa ； κ- 内变量， Mpa ； (μ1) 无量纲 ； qe 弹性坍塌压力， Mpa ； F 广义力， MPa ； U 广义位移， m2 ； qcr m- 第 m 步的最小特征值， 应变软化坍塌压力， MPa ； λcr 应变软化极限载荷因子， 无量纲。
本发明的有益效果 ： 本发明井壁坍塌分析方法具有 ：
(1) 能分析非均匀地应力条件下应变软化塑性地层的井壁稳定性 ；
(2) 实现了力学稳定性方法分析井壁稳定的功能， 得到平衡路径上的极值点， 计算 井壁坍塌压力 ；
(3) 现场应用表明， 本方法能很好的指导实际钻井作业， 根据本方法确定的防塌钻井液密度， 在现场钻井中井壁稳定性良好， 且能起到储层保护的效果。 附图说明
图 1 是本发明井壁坍塌分析方法的流程示意图。 图 2 是实施例中大北 101 井的平衡路径曲线图。具体实施方式
实施例 1 ： 以一口井的井壁坍塌分析为例， 对本发明作进一步详细说明。
以塔里木油田大北 101 井为例， 用本发明井壁坍塌分析方法对该井的 1323 ～ 1340m 井段进行坍塌分析 ：
井壁坍塌分析方法包括 ： 采集并输入参数 1， 初应力处理 2， 第一应力释放 3， 第二 应力释放第一步 4A， 第二应力释放第二步 4B， 特征值评价 5， 广义力和广义位移 6， 平衡路径 曲线 7， 坍塌压力处理 8， 输出坍塌压力以及塑性区分布图。参阅图 1。
(A) 首先， 采集通过相邻已钻井的测井资料和岩石力学实验得到井壁岩石材料参 数， 包括弹性模量、 泊松比、 内摩擦角、 内聚力、 软化系数、 屈服函数类型、 内变量类型、 水平 最大主应力、 水平最小主应力和深度参数， 将井壁岩石材料参数输入计算程序中 ； 输入控 制参数， 控制参数包括弧长、 输出步、 切线塑性剪切模量、 残余强度。控制参数中弧长取为 -5 10 m， 输出步取为 100。 控制参数中的切线塑性剪切模量、 残余强度通过岩石力学实验得到。 (B) 初应力处理 2 ： 建立井壁几何模型， 井筒内半径为 0.2m， 内半径与外半径比值 为 1 ∶ 10 ； 划分有限元计算网格单元， 采用 4 节点平面应变等参单元和描述远场的 4 节点 无限区域单元， 在最外一层设置成无限单元它的尺度在规定的方向上延伸至无穷远， 起到 模拟远处区域的作用， 无限区域元的位移插值函数与坐标插值函数有不同的表达形式， 位 移在无限远处为零， 通过衰减函数的引入来实现。
衰减函数采用负指数函数 ： 共剖分为 462 个节点， 420 个单元 ； 施加边界条件， 施加情况为 2 个直边法向约束、径向自由。根据网格节点坐标， 由式 (3) 计算初始应力场各单元高斯点的应力
(3) (C) 第一次应力释放 3 ： 采用 “释放载荷”的方法来模拟实际钻井过程， 求扰动 公式 (5) 为保留的场 A。首先， 根据公式 (4) 计算第一次释放部分应力分量应力分量。这样的应力释放方案， 等效于式 (3) 的应力完全释放而同时施加钻井液压力
8(4)102182453 A CN 102182457
说明书7/9 页在 r ＝ a，按弹性或弹 - 理想塑性计算， 得扰动场 A， σ′， u′， 通过第一次应力释放， 得到总场 ： 0 0 0 σ +σ′→ σ ； a′→ a (6) (D) 第二次应力释放第一步 4A ： 释放保留的应力分量， 也就是模拟钻井液压力逐 τ′ θ ＝ 0。渐降低， 计算扰动应力场 B。引入载荷参数 λ， 在 r ＝ a，
第一步用弹性方法确定弹性阶段的载荷参数增量 Δλ1， Δλ1 → λe， 得到相应的 弹性扰动场 Δa′， Δσ′， 得到总场 ： 0 0
σ +Δσ′→ σ ； a0+Δa′→ a0 (7)
(E) 第二次应力释放第二步 4B ： 由于岩石在高温高压条件下， 呈现应变软化塑性 性质， 所以从第二步开始 (m ＝ 1)， 指定微弧长增量 Δsm， 按弹 - 软化塑性， 用弧长法迭代求 解。全过程曲线 τ-γ 采用三线性形式， 下降段斜率是一个负的常数 GT， 称为切线剪切模 量， 而相应的软化曲线则为二线性形式， 下降段斜率为 Gp ＝ GGT/(G-GT)。 由于采用三线性模 型， 在 D-P 准则中 k 就是剪切强度 τs， 它也是初始的屈服应力 τs(0)。设刚刚进入残余阶 段的塑性应变为 那么二线性的软化曲线斜率为
剪切强度 τs 和塑性功 wp 的关系为 ：在对岩石等应变软化材料组成的结构应用有限元方法进行非线性分析时， 一旦结 构接近于极限承载能力或者以后， 用一般的计算方法计算分析的稳定性与收敛性相当差以 至于无法进行下去。弧长法通过同时约束荷载和位移向量， 能较好地计算邻近极值点结构 的反应和下降段问题， 可以比较好地解决极值点穿越问题。
对于一个给定的迭代步 n 已知解待求的下一迭代步 n+1 的解记为可有如下的表达式 ：
于是， 令可给出位移解的增量 Δan 和载荷参数的增量 Δλn 之间的关系
其中
将式 (12) 代入约束方程可得一个关于 Δλn 的二次代数方程 其中
从方程 (15) 可以得到 Δλn 的两个根， 分别记为 (Δλn)1 和 (Δλn)2。为了能够 跟踪解的路径得到正确的结果， 只选取一组解， 它应使解沿平衡曲线往前走， 而不是往后返 n T 回， 即取 (Δa ) (am-am-1) 为一个最大值的解， 得扰动场 B， Δλm， Δam， Δσm， Δκm， 求出总 场：
λm+Δλm → λm+1
(17)(F) 特征值评价 5 ： 结构的切线总刚矩阵不正定才是结构不稳定的充分和必要条 的最小特征值件。因此， 引入了计算总刚矩阵特征值的算法 - 子空间迭代法。求根据得到的最小特征值 (μ1)m 进行评价， 当最小特征值 (μ1)m ＞ 0， 时， 如果 m ＜ M， (μ1)m， 则 m+1 → m， 进入第二次应力释放第二步， 否则， 得到的即为总场 ； 当最小特征值 (μ1)m ＜ 0 时， 进行如下计算。
λm → λN， λm-1 → λN-1
(μ1)m → (μ)N， (μ1)m-1 → (μ)N-1 (18)
σm → σN， σm-1 → σN-1 κm → κN， κm-1 → κN-1 使用线性插值计算临界扰动应力场 ：计算总场 ： Δλcr+λm → λcr Δσcr+σm → σcr Δκcr+κm → κcr(20)(G) 广义力和广义位移 6 ： 根据弧长法每步骤计算出的载荷因子 λ， 按照两次解除 应力的公式通过积分即可求出当步广义力及广义位移。
在计算总场时， 外力功的变分为
如果广义力可定义为那么广义位移的变分为第 m 个载荷步的增量广义位移为总广义位移可累加得到 ：
Um+1 ＝ Um+ΔUm (25)
(H) 平衡路径曲线 7 ： 平衡路径曲线是由广义力和广义位移为坐标轴绘制的曲线。 根据广义力和广义位移的输出结果， 绘制平衡路径曲线。参阅图 2。
(I) 坍塌压力处理 8 ： 用 P 和 u 分别代表系统的广义力和广义位移， 那么在稳定分 支上有 δPδu ＞ 0， 在不稳定分支上有 δPδu ＜ 0， 稳定分支和不稳定分支的交汇点是临 界点， 在该点 δPδu ＝ 0， 该点对应于平衡稳定性的临界状态， 该点处的广义力称为临界 力， 它确定了一个结构的承载能力和保持整体性的能力。 在平衡路径曲线上找出临界力， 根 据广义力的公式 (22) 确定对应的应变软化极限载荷因子 λcr。
根据得到的应力总场、 弹性极限载荷因子 λe 以及平衡路径曲线输出的 λcr， 进入 坍塌压力处理， 得到弹性坍塌压力和应变软化塑性坍塌压力 ：
(12) 根据弹性强度理论得出的坍塌压力为3防塌钻井液密度应为 1.1451.1g/cm ， 本系统用应变软化模型力学稳定性分析得出坍塌压力为 依坍塌压力确定防塌钻井液密度为 1.053g/cm3。本 井段实际用钻井液密度为 1.1g/cm3， 按照弹性强度理论结果在此井段应该出现坍塌扩径现 象， 而实际上从井径测井资料来看本井段并无井径坍塌扩大现象， 系统计算结果表明在此 井段不会发生坍塌， 这与实际情况相吻合， 表明了本系统的可靠性和实用性。