CN102182453B - 井壁坍塌分析方法 - Google Patents

Abstract

井壁坍塌分析方法，包括：采集并输入参数(1)，初应力处理(2)，第一次应力释放(3)，第二次应力释放第一步(4A)，第二次应力释放第二步(4B)，特征值评价(5)，广义力和广义位移(6)，平衡路径曲线(7)，坍塌压力处理(8)。根据得到的最小特征值(μ₁)_m进行评价，当最小特征值(μ₁)_m＞0，时，如果m＜M，则m+1→m，即进入第二次应力释放第二步(4B)；否则，得到的即是总场；当最小特征值(μ₁)_m＜0时，根据插值法计算扰动场；计算坍塌压力q_cr，以及q_cr对应的塑性区分布图，确定防塌钻井液的密度。效果是：能分析非均匀地应力条件下应变软化塑性地层的井壁稳定性；设计的防塌钻井液密度，在现场钻井中起到储层保护效果。

Description

井壁坍塌分析方法

技术领域

本发明涉及石油钻井工程技术领域，特别涉及钻井井壁失稳分析，是一种非均匀地应力条件下应变软化塑性地层井壁稳定性分析的方法。

背景技术

目前，在石油工程中普遍应用弹性力学结合强度准则方法解释井壁坍塌机制和预测坍塌压力。比如，被认为是世界上最好的井壁稳定分析GMI系统，其应用的井壁坍塌模型正是基于弹性力学结合强度准则法。最新的井壁坍塌机理的理论研究表明，井壁坍塌实质上是一种极值点型失稳现象。通过确定井壁压力和井壁位移的平衡路径曲线的极值即可确定力学稳定性意义下的坍塌压力。当然，理论解作了原场地应力是各向同性均匀应力和地层满足Tresca准则的假设，这远离了实际情况。还有重要一点，在分析井壁坍塌问题时仅考虑岩石的理想弹塑性或者弹脆性-塑性特性，而大部分岩石的应力-应变曲线在峰值后下降，称为应变软化，所以不考虑岩石的应变软化特性也是不符合实际的。

国外现有Ansys，Flac，Abaqus，Adina等有限元计算软件，但是缺乏力学意义上的稳定性分析功能，没有将应变软化模型和稳定性分析结合，而且应用领域没有扩展到井壁坍塌分析。国内岩土方面计算系统比较突出的是“理正岩质边坡稳定分析系统”，所采用的也是传统的极限平衡法。未见报道应用应变软化模型和力学稳定性方法开发井壁坍塌分析系统。为了解决非均匀地应力条件下应变软化塑性地层井壁失稳问题，本发明在应变软化理论模型的基础上，应用应变软化模型和力学稳定性方法开发WellCA井壁稳定分析系统。

井壁稳定问题是世界性难题，在世界范围内，每年用于处理井壁失稳的费用高达数亿美元，而我国许多油田也存在井壁失稳复杂情况。近年来，我国在西部塔里木和准噶尔等盆地进行了多口深井、超深井的钻探施工。这些地区山前高陡构造高地应力引起的井眼稳定问题较为突出，特别是深井的井眼稳定受山前构造带地应力影响，问题更为严重。

中国专利授权公告号：CN1239920C，提供了一种“利用地震层速度钻前预测坍塌压力与破裂压力的方法”。包括下列步骤：1)将待钻井与多个相邻已钻井的地震层速度进行数据相关分析，并确定相关系数大于0.75的已钻井为待钻进具有相似构造的已钻井；2)利用声波时差、自然伽玛、密度等测井数据序列，对已钻井全井段进行分层，求出每一层用于表征一定厚度且岩性相似地层的平均声波速度、自然伽玛和地层密度；3)利用每一层的平均声波速度、自然伽玛和地层密度来确定已钻井的坍塌压力与破裂压力；4)根据已钻井的坍塌压力与破裂压力和测井分层速度建立测井模型；5)建立已钻井的层速度钻前预测模型；6)将待钻井的地震层速度代入步骤5)中的模型，获得待钻井的坍塌压力与破裂压力。

发明内容

本发明的目的是：提供一种井壁坍塌分析方法，考虑岩石的应变软化性质，使其更接近客观实际，应用于钻井井壁稳定性理论研究和建立井壁坍塌模型，效解决钻井的井壁坍塌问题。

本发明采用的技术方案是：

井壁坍塌分析方法，包括：采集并输入参数，初应力处理，第一次应力释放，第二次应力释放第一步，第二次应力释放第二步，特征值评价，广义力和广义位移，平衡路径曲线，坍塌压力处理，输出坍塌压力以及塑性区分布图。

(A)采集并输入参数：首先，采集通过相邻已钻井的测井资料和岩石力学实验得到井壁岩石材料参数，包括弹性模量、泊松比、内摩擦角、内聚力、软化系数、屈服函数类型、内变量类型、水平最大主应力、水平最小主应力和深度参数，将井壁岩石材料参数输入计算程序中；输入控制参数，控制参数包括弧长、输出步、切线塑性剪切模量、残余强度。控制参数中的弧长范围10^-5～10^-6m，输出步范围100～200，控制参数中的切线塑性剪切模量、残余强度通过岩石力学实验得到。

(B)初应力处理：建立井壁几何模型，设置井筒半径，内外半径比；划分有限元计算网格单元，采用平面应变等参单元和可描述远场的无限区域单元，在最外一层设置成无限单元它的尺度在规定的方向上延伸至无穷远，起到模拟远处区域的作用，无限区域元的位移插值函数与坐标插值函数有不同的表达形式，位移在无限远处为零，通过衰减函数的引入来实现。

$\stackrel{\‾}{N_i}=N_i{f}_i\left(r\right)---\left(1\right)$

衰减函数采用负指数函数：

$f_i\left(r\right)=e^{-(r-r_i)/L}---\left(2\right)$

剖分节点和单元，施加边界条件，施加情况为2个直边法向约束、径向自由。根据网格节点坐标，由式(3)计算初始应力场各单元高斯点的应力

${\mathrm{\σ}}_r^0=\frac{1}{2}({\mathrm{\σ}}_H+{\mathrm{\σ}}_h)-\frac{1}{2}({\mathrm{\σ}}_H-{\mathrm{\σ}}_h)\cos 2\mathrm{\θ}---\left(3\right)$

${\mathrm{\τ}}_{\mathrm{\θ}}^0=\frac{1}{2}({\mathrm{\σ}}_H-{\mathrm{\σ}}_h)\sin 2\mathrm{\θ}$

(C)第一次应力释放：采用“释放载荷”的方法来模拟实际钻井过程，求扰动场A。首先，根据公式(4)计算第一次释放部分应力分量

公式(5)为保留的应力分量。这样的应力释放方案，等效于式(3)的应力完全释放而同时施加钻井液压力

${\mathrm{\σ}}_r^{\left(1\right)}=-\frac{1}{2}({\mathrm{\σ}}_H-{\mathrm{\σ}}_h)\cos 2\mathrm{\θ}$

(4)

${\mathrm{\τ}}_{\mathrm{\θ}}^{\left(1\right)}={\mathrm{\τ}}_{\mathrm{\θ}}^0=\frac{1}{2}({\mathrm{\σ}}_H-{\mathrm{\σ}}_h)\sin 2\mathrm{\θ}$

${\mathrm{\σ}}_r^{\left(2\right)}=\frac{1}{2}\left({\mathrm{\σ}}_H\text{+}{\mathrm{\σ}}_h\right)---\left(5\right)$

在r＝a，按弹性或弹-理想塑性计算，得扰动场A，σ′，u′，通过第一应力释放，得到总场，这个总场是指初始场加扰动场A：

σ⁰+σ′→σ⁰；a′→a⁰    (6)

(D)第二次应力释放第一步：逐渐释放保留的应力分量，也就是模拟钻井液压力逐渐降低，计算扰动应力场B。引入载荷参数λ，在r＝a，τ′_θ＝0。

第一步用弹性方法确定弹性阶段的载荷参数增量Δλ₁，Δλ₁ →λ_e，求出相应的弹性扰动场Δa′，Δσ′，得到总场，这个总场是指初始场加扰动场A再加弹性扰动场：

σ⁰+Δσ′→σ⁰；a⁰+Δa′→a⁰      (7)

(E)第二次应力释放第二步：由于岩石在高温高压条件下，呈现应变软化塑性性质，所以从第二步开始(m＝1)，指定微弧长增量Δs_m，按弹-软化塑性，用弧长法迭代求解。全过程曲线τ-γ采用三线性形式，下降段斜率是一个负的常数G_T，称为切线剪切模量，而相应的软化曲线则为二线性形式，下降段斜率为G_p＝GG_T/(G-G_T)。由于采用三线性模型，在D-P准则中k就是剪切强度τ_s，它也是初始的屈服应力τ_s(0)。设刚刚进入残余阶段的塑性应变为

那么二线性的软化曲线斜率为

$\frac{\&PartialD;{\mathrm{\&tau;}}_s}{\&PartialD;{\mathrm{\&gamma;}}^p}=\left\{\begin{array}{c}\frac{{\mathrm{GG}}_T}{G-G_T},0\&le;{\mathrm{\&gamma;}}^p\&le;{\mathrm{\&gamma;}}_r^p\\ 0,{\mathrm{\&gamma;}}_r^p<{\mathrm{\&gamma;}}^p\end{array}\right.---\left(8\right)$

剪切强度τ_s和塑性功w^p的关系为：

${\mathrm{\&tau;}}_s\left(w^p\right)={[{\mathrm{\&tau;}}_s^2\left(0\right)+2G_p{w}^p]}^{1/2}---\left(9\right)$

在对岩石等应变软化材料组成的结构应用有限元方法进行非线性分析时，一旦结构接近于极限承载能力或者以后，用一般的计算方法计算分析的稳定性与收敛性相当差以至于无法进行下去。弧长法通过同时约束荷载和位移向量，能较好地计算邻近极值点结构的反应和下降段问题，可以比较好地解决极值点穿越问题。

对于一个给定的迭代步n已知解

待求的下一迭代步n+1的解记为

可有如下的表达式：

${\mathrm{\&psi;}}_{m+1}^n\&equiv;\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)=P\left(a_{m+1}^n\right)-{\mathrm{\&lambda;}}_{m+1}^{n}R---\left(10\right)$

${\mathrm{\&psi;}}_{m+1}^{n+1}=\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)+{\left(K_T\right)}^n(a_{m+1}^{n+1}-a_{m+1}^n)-({\mathrm{\&lambda;}}_{m+1}^{n+1}-{\mathrm{\&lambda;}}_{m+1}^n)R---\left(11\right)$

于是，令

可给出位移解的增量Δaⁿ和载荷参数的增量Δλⁿ之间的关系

其中

将式(12)代入约束方程

${\left(\mathrm{\&Delta;}{a}_m^n\right)}^T\left(\mathrm{\&Delta;}{a}_m^n\right)+c{\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)}^2={\left(\mathrm{\&Delta;s}\right)}^2---\left(14\right)$

可得一个关于Δλⁿ的二次代数方程

$b_1{\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)}^2+b_2\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)+b_3=0---\left(15\right)$

其中

从方程(15)可以得到Δλⁿ的两个根，分别记为(Δλⁿ)₁和(Δλⁿ)₂。为了能够跟踪解的路径得到正确的结果，只选取一组解，它应使解沿平衡曲线往前走，而不是往后返回，即取(Δaⁿ)^T(a_m-a_m-1)为一个最大值的解，得扰动场B，Δλ_m，Δa_m，Δσ_m，Δκ_m，求出总场，这个总场是指初始场加扰动场A再加扰动场B：

λ_m+Δλ_m →λ_m+1

$a_m^0+\mathrm{\&Delta;}{a}_m\&RightArrow;a_{m+1}^0$

(17)

${\mathrm{\&sigma;}}_m^0+\mathrm{\&Delta;}{\mathrm{\&sigma;}}_m\&RightArrow;{\mathrm{\&sigma;}}_{m+1}^0$

${\mathrm{\&kappa;}}_m^0+\mathrm{\&Delta;}{\mathrm{\&kappa;}}_m\&RightArrow;{\mathrm{\&kappa;}}_{m+1}$

(F)特征值评价：结构的切线总刚矩阵不正定才是结构不稳定的充分和必要条件。因此，引入了计算总刚矩阵特征值的算法-子空间迭代法。求的最小特征值(μ₁)_m，根据得到的最小特征值(μ₁)_m进行评价，当最小特征值(μ₁)_m＞0，时，如果m＜M，则m+1→m，即m加上1以后进入第二次应力释放第二步，否则，得到的即为总场；当最小特征值(μ₁)_m＜0时，进行如下计算。

λ_m →λ_N，λ_m-1→λ_N-1

(μ₁)_m →(μ)_N，(μ₁)_m-1→(μ)_N-1                (18)

σ_m→σ_N，σ_m-1→σ_N-1

κ_m→κ_N，κ_m-1→κ_N-1

使用线性插值计算临界扰动应力场：

$\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_N-\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{N-1})$

$\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}\text{-}{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&sigma;}}_N-\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{N-1})---\left(19\right)$

$\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&kappa;}}_N-\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{N-1})$

计算总场，这个总场是指初始场加扰动场A再加临界扰动场：

Δλ_cr+λ_m→λ_cr

Δσ_cr+σ_m→σ_cr                     (20)

Δκ_cr+κ_m→κ_cr

(G)广义力和广义位移：根据弧长法每步骤计算出的载荷因子λ，按照两次解除应力的公式通过积分即可求出当步广义力及广义位移。

在计算总场时，外力功的变分为

$\mathrm{\&delta;W}=\mathrm{F\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}\left[q_r\mathrm{\&delta;}{u}_r\right]\mathrm{ad\&theta;}={\&Integral;}_0^{\mathrm{\&pi;}/2}(1-\mathrm{\&lambda;}){\mathrm{\&sigma;}}_r^{\left(2\right)}\mathrm{\&delta;}{u}_r\mathrm{ad\&theta;}---\left(21\right)$

如果广义力可定义为

$F=\frac{a}{2}(1-\mathrm{\&lambda;})({\mathrm{\&sigma;}}_H+{\mathrm{\&sigma;}}_h)---\left(22\right)$

那么广义位移的变分为

$\mathrm{\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}\mathrm{\&delta;}{u}_r\mathrm{d\&theta;}---\left(23\right)$

第m个载荷步的增量广义位移为

$\mathrm{\&Delta;}{U}_m=\underset{i}{\mathrm{\&Sigma;}}\mathrm{\&Delta;}{u}_r\left({\mathrm{\&theta;}}_i\right)\mathrm{\&Delta;\&theta;}---\left(24\right)$

总广义位移可累加得到：

U_m+1＝U_m+ΔU_m      (25)

(H)平衡路径曲线：平衡路径曲线是由广义力和广义位移为坐标轴绘制的曲线。根据广义力和广义位移的输出结果，绘制平衡路径曲线。

(I)坍塌压力处理：用P和u分别代表系统的广义力和广义位移，那么在稳定分支上有δPδu＞0，在不稳定分支上有δPδu＜0，稳定分支和不稳定分支的交汇点是临界点，在该点δPδu＝0，该点对应于平衡稳定性的临界状态，该点处的广义力称为临界力，它确定了一个结构的承载能力和保持整体性的能力。在平衡路径曲线上找出临界力，根据广义力的公式(22)确定对应的应变软化极限载荷因子λ_cr。

根据得到的应力总场、弹性极限载荷因子λ_e以及平衡路径曲线输出的λ_cr，进入坍塌压力处理，得到弹性坍塌压力和应变软化塑性坍塌压力：

$(1-{\mathrm{\&lambda;}}_e){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_e$

(12)

$(1-{\mathrm{\&lambda;}}_{\mathrm{cr}}){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_{\mathrm{cr}}$

(J)输出坍塌压力q_cr，应变软化极限载荷因子λ_cr对应的塑性区分布图，确定该井的防塌钻井液的密度。

符号说明：

-第i个节点的位移插值函数；N_i-第i个节点的坐标插值函数；L-衰减特征长度，m；r_i-第i个节点位移，m；

-初始场径向应力，MPa；

-初始场切向应力，MPa；σ_H-最大水平主应力，Mpa；σ_h-最小水平主应力，Mpa；

-第一次解除的径向应力，Mpa；

-第一次解除的切向应力，Mpa；θ-与最大水平地应力的夹角，°；

-第一次保留的应力分量，Mpa；λ-载荷因子，无量纲；r-距井眼中心距离，m；a-井眼半径，m；σ⁰-初始场应力，Mpa；σ′-扰动场应力，Mpa；σ′_r-扰动场径向应力，Mpa；τ′_θ-扰动场切向应力，Mpa；Δλ₁-弹性阶段的载荷参数增量，无量纲；λ_e-弹性极限载荷参数，无量纲；G_T为切线剪切模量，MPa；G_p-切线塑性剪切模量，MPa；τ_s(0)-初始的屈服应力，Mpa；τ_s-剪切强度，MPa；w^p-塑性功，MPa；

为残余阶段的塑性应变，％；M-弧长法微弧长增量总量，无量纲；Δs_m微弧长增量，m；Δaⁿ-位移增量向量，m；Δλⁿ-载荷增量系数，无量纲；

-切线刚度矩阵；R-载荷力向量，MPa；ψ_m-不平衡力向量，MPa；κ-内变量，Mpa；(μ₁)_m-第m步的最小特征值，无量纲；q_e弹性坍塌压力，Mpa；F广义力，MPa；U广义位移，m²；q_cr应变软化坍塌压力，MPa；λ_cr应变软化极限载荷因子，无量纲。

本发明的有益效果：本发明井壁坍塌分析方法具有：

(1)能分析非均匀地应力条件下应变软化塑性地层的井壁稳定性；

(2)实现了力学稳定性方法分析井壁稳定的功能，得到平衡路径上的极值点，计算井壁坍塌压力；

(3)现场应用表明，本方法能很好的指导实际钻井作业，根据本方法确定的防塌钻井液密度，在现场钻井中井壁稳定性良好，且能起到储层保护的效果。

附图说明

图1是本发明井壁坍塌分析方法的流程示意图。

图2是实施例中大北101井的平衡路径曲线图。

具体实施方式

实施例1：以一口井的井壁坍塌分析为例，对本发明作进一步详细说明。

以塔里木油田大北101井为例，用本发明井壁坍塌分析方法对该井的1323～1340m井段进行坍塌分析：

井壁坍塌分析方法包括：采集并输入参数1，初应力处理2，第一应力释放3，第二应力释放第一步4A，第二应力释放第二步4B，特征值评价5，广义力和广义位移6，平衡路径曲线7，坍塌压力处理8，输出坍塌压力以及塑性区分布图。参阅图1。

(A)首先，采集通过相邻已钻井的测井资料和岩石力学实验得到井壁岩石材料参数，包括弹性模量、泊松比、内摩擦角、内聚力、软化系数、屈服函数类型、内变量类型、水平最大主应力、水平最小主应力和深度参数，将井壁岩石材料参数输入计算程序中；输入控制参数，控制参数包括弧长、输出步、切线塑性剪切模量、残余强度。控制参数中弧长取为10^-5m，输出步取为100。控制参数中的切线塑性剪切模量、残余强度通过岩石力学实验得到。

(B)初应力处理2：建立井壁几何模型，井筒内半径为0.2m，内半径与外半径比值为1∶10；划分有限元计算网格单元，采用4节点平面应变等参单元和描述远场的4节点无限区域单元，在最外一层设置成无限单元它的尺度在规定的方向上延伸至无穷远，起到模拟远处区域的作用，无限区域元的位移插值函数与坐标插值函数有不同的表达形式，位移在无限远处为零，通过衰减函数的引入来实现。

$\stackrel{\&OverBar;}{N_i}=N_i{f}_i\left(r\right)---\left(1\right)$

衰减函数采用负指数函数：

$f_i\left(r\right)=e^{-(r-r_i)/L}---\left(2\right)$

共剖分为462个节点，420个单元；施加边界条件，施加情况为2个直边法向约束、径向自由。根据网格节点坐标，由式(3)计算初始应力场各单元高斯点的应力

${\mathrm{\&sigma;}}_r^0=\frac{1}{2}({\mathrm{\&sigma;}}_H+{\mathrm{\&sigma;}}_h)-\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\cos 2\mathrm{\&theta;}$

(3)

${\mathrm{\&tau;}}_{\mathrm{\&theta;}}^0=\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\sin 2\mathrm{\&theta;}$

(C)第一次应力释放3：采用“释放载荷”的方法来模拟实际钻井过程，求扰动场A。首先，根据公式(4)计算第一次释放部分应力分量

公式(5)为保留的应力分量。这样的应力释放方案，等效于式(3)的应力完全释放而同时施加钻井液压力

${\mathrm{\&sigma;}}_r^{\left(1\right)}=-\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\cos 2\mathrm{\&theta;}$

(4)

${\mathrm{\&tau;}}_{\mathrm{\&theta;}}^{\left(1\right)}={\mathrm{\&tau;}}_{\mathrm{\&theta;}}^0=\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\sin 2\mathrm{\&theta;}$

${\mathrm{\&sigma;}}_r^{\left(2\right)}=\frac{1}{2}\left({\mathrm{\&sigma;}}_H\text{+}{\mathrm{\&sigma;}}_h\right)---\left(5\right)$

在r＝a，

按弹性或弹-理想塑性计算，得扰动

场A，σ′，u′，通过第一次应力释放，得到总场：

σ⁰+σ′→σ⁰；a′→a⁰    (6)

(D)第二次应力释放第一步4A：释放保留的应力分量，也就是模拟钻井液压力逐渐降低，计算扰动应力场B。引入载荷参数λ，在r＝a，

τ′_θ＝0。

第一步用弹性方法确定弹性阶段的载荷参数增量Δλ₁，Δλ₁→λ_e，得到相应的弹性扰动场Δa′，Δσ′，得到总场：

σ⁰+Δσ′→σ⁰；a⁰+Δa′→a⁰    (7)

(E)第二次应力释放第二步4B：由于岩石在高温高压条件下，呈现应变软化塑性性质，所以从第二步开始(m＝1)，指定微弧长增量Δs_m，按弹-软化塑性，用弧长法迭代求解。全过程曲线τ-γ采用三线性形式，下降段斜率是一个负的常数G_T，称为切线剪切模量，而相应的软化曲线则为二线性形式，下降段斜率为G_p＝GG_T/(G-G_T)。由于采用三线性模型，在D-P准则中k就是剪切强度τ_s，它也是初始的屈服应力τ_s(0)。设刚刚进入残余阶段的塑性应变为

那么二线性的软化曲线斜率为

$\frac{\&PartialD;{\mathrm{\&tau;}}_s}{\&PartialD;{\mathrm{\&gamma;}}^p}=\left\{\begin{array}{c}\frac{{\mathrm{GG}}_T}{G-G_T},0\&le;{\mathrm{\&gamma;}}^p\&le;{\mathrm{\&gamma;}}_r^p\\ 0,{\mathrm{\&gamma;}}_r^p<{\mathrm{\&gamma;}}^p\end{array}\right.---\left(8\right)$

剪切强度τ_s和塑性功w^p的关系为：

${\mathrm{\&tau;}}_s\left(w^p\right)={[{\mathrm{\&tau;}}_s^2\left(0\right)+2G_p{w}^p]}^{1/2}---\left(9\right)$

在对岩石等应变软化材料组成的结构应用有限元方法进行非线性分析时，一旦结构接近于极限承载能力或者以后，用一般的计算方法计算分析的稳定性与收敛性相当差以至于无法进行下去。弧长法通过同时约束荷载和位移向量，能较好地计算邻近极值点结构的反应和下降段问题，可以比较好地解决极值点穿越问题。

对于一个给定的迭代步n已知解

待求的下一迭代步n+1的解记为

可有如下的表达式：

${\mathrm{\&psi;}}_{m+1}^n\&equiv;\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)=P\left(a_{m+1}^n\right)-{\mathrm{\&lambda;}}_{m+1}^{n}R---\left(10\right)$

${\mathrm{\&psi;}}_{m+1}^{n+1}=\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)+{\left(K_T\right)}^n(a_{m+1}^{n+1}-a_{m+1}^n)-({\mathrm{\&lambda;}}_{m+1}^{n+1}-{\mathrm{\&lambda;}}_{m+1}^n)R---\left(11\right)$

于是，令

可给出位移解的增量Δaⁿ和载荷参数的增量Δλⁿ之间的关系

其中

将式(12)代入约束方程

${\left(\mathrm{\&Delta;}{a}_m^n\right)}^T\left(\mathrm{\&Delta;}{a}_m^n\right)+c{\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)}^2={\left(\mathrm{\&Delta;s}\right)}^2---\left(14\right)$

可得一个关于Δλⁿ的二次代数方程

$b_1{\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)}^2+b_2\left(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_m^n\right)+b_3=0---\left(15\right)$

其中

从方程(15)可以得到Δλⁿ的两个根，分别记为(Δλⁿ)₁和(Δλⁿ)₂。为了能够跟踪解的路径得到正确的结果，只选取一组解，它应使解沿平衡曲线往前走，而不是往后返回，即取(Δaⁿ)^T(a_m-a_m-1)为一个最大值的解，得扰动场B，Δλ_m，Δa_m，Δσ_m，Δκ_m，求出总场：

λ_m+Δλ_m→λ_m+1

$a_m^0+\mathrm{\&Delta;}{a}_m\&RightArrow;a_{m+1}^0$

(17)

${\mathrm{\&sigma;}}_m^0+\mathrm{\&Delta;}{\mathrm{\&sigma;}}_m\&RightArrow;{\mathrm{\&sigma;}}_{m+1}^0$

${\mathrm{\&kappa;}}_m^0+\mathrm{\&Delta;}{\mathrm{\&kappa;}}_m\&RightArrow;{\mathrm{\&kappa;}}_{m+1}$

(F)特征值评价5：结构的切线总刚矩阵不正定才是结构不稳定的充分和必要条件。因此，引入了计算总刚矩阵特征值的算法-子空间迭代法。求的最小特征值(μ₁)_m，根据得到的最小特征值(μ₁)_m进行评价，当最小特征值(μ₁)_m＞0，时，如果m＜M，则m+1→m，进入第二次应力释放第二步，否则，得到的即为总场；当最小特征值(μ₁)_m＜0时，进行如下计算。

λ_m→λ_N，λ_m-1→λ_N-1

(μ₁)_m→(μ)_N，(μ₁)_m-1→(μ)_N-1                  (18)

σ_m→σ_N，σ_m-1→σ_N-1

κ_m→κ_N，κ_m-1→κ_N-1

使用线性插值计算临界扰动应力场：

$\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&lambda;}}_N-\mathrm{\&Delta;}{\mathrm{\&lambda;}}_{N-1})$

$\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}\text{-}{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&sigma;}}_N-\mathrm{\&Delta;}{\mathrm{\&sigma;}}_{N-1})---\left(19\right)$

$\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{\mathrm{cr}}=\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}(\mathrm{\&Delta;}{\mathrm{\&kappa;}}_N-\mathrm{\&Delta;}{\mathrm{\&kappa;}}_{N-1})$

计算总场：

Δλ_cr+λ_m→λ_cr

Δσ_cr+σ_m→σ_cr               (20)

Δκ_cr+κ_m→κ_cr

(G)广义力和广义位移6：根据弧长法每步骤计算出的载荷因子λ，按照两次解除应力的公式通过积分即可求出当步广义力及广义位移。

在计算总场时，外力功的变分为

$\mathrm{\&delta;W}=\mathrm{F\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}\left[q_r\mathrm{\&delta;}{u}_r\right]\mathrm{ad\&theta;}={\&Integral;}_0^{\mathrm{\&pi;}/2}(1-\mathrm{\&lambda;}){\mathrm{\&sigma;}}_r^{\left(2\right)}\mathrm{\&delta;}{u}_r\mathrm{ad\&theta;}---\left(21\right)$

如果广义力可定义为

$F=\frac{a}{2}(1-\mathrm{\&lambda;})({\mathrm{\&sigma;}}_H+{\mathrm{\&sigma;}}_h)---\left(22\right)$

那么广义位移的变分为

$\mathrm{\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}\mathrm{\&delta;}{u}_r\mathrm{d\&theta;}---\left(23\right)$

第m个载荷步的增量广义位移为

$\mathrm{\&Delta;}{U}_m=\underset{i}{\mathrm{\&Sigma;}}\mathrm{\&Delta;}{u}_r\left({\mathrm{\&theta;}}_i\right)\mathrm{\&Delta;\&theta;}---\left(24\right)$

总广义位移可累加得到：

U_m+1＝U_m+ΔU_m       (25)

(H)平衡路径曲线7：平衡路径曲线是由广义力和广义位移为坐标轴绘制的曲线。根据广义力和广义位移的输出结果，绘制平衡路径曲线。参阅图2。

(I)坍塌压力处理8：用P和u分别代表系统的广义力和广义位移，那么在稳定分支上有δPδu＞0，在不稳定分支上有δPδu＜0，稳定分支和不稳定分支的交汇点是临界点，在该点δPδu＝0，该点对应于平衡稳定性的临界状态，该点处的广义力称为临界力，它确定了一个结构的承载能力和保持整体性的能力。在平衡路径曲线上找出临界力，根据广义力的公式(22)确定对应的应变软化极限载荷因子λ_cr。

根据得到的应力总场、弹性极限载荷因子λ_e以及平衡路径曲线输出的λ_cr，进入坍塌压力处理，得到弹性坍塌压力和应变软化塑性坍塌压力：

$(1-{\mathrm{\&lambda;}}_e){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_e$

(12)

$(1-{\mathrm{\&lambda;}}_{\mathrm{cr}}){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_{\mathrm{cr}}$

根据弹性强度理论得出的坍塌压力为

防塌钻井液密度应为1.1451.1g/cm³，本系统用应变软化模型力学稳定性分析得出坍塌压力为

依坍塌压力确定防塌钻井液密度为1.053g/cm³。本井段实际用钻井液密度为1.1g/cm³，按照弹性强度理论结果在此井段应该出现坍塌扩径现象，而实际上从井径测井资料来看本井段并无井径坍塌扩大现象，系统计算结果表明在此井段不会发生坍塌，这与实际情况相吻合，表明了本系统的可靠性和实用性。

Claims (2)

1.一种井壁坍塌分析方法，包括：采集并输入参数(1)，初应力处理(2)，第一次应力释放(3)，第二次应力释放第一步(4A)，第二次应力释放第二步(4B)，特征值评价(5)，广义力和广义位移(6)，平衡路径曲线(7)，坍塌压力处理(8)；其特征为：

(A)采集并输入参数(1)：首先，采集通过相邻已钻井的测井资料和通过岩石力学实验得到井壁岩石材料参数，将得到井壁岩石材料参数输入计算程序中；其次输入控制参数，控制参数中的弧长范围10^-5-10^-6m，输出步范围100-200，切线塑性剪切模量、残余强度通过岩石力学实验得到；

(B)初应力处理(2)：建立井壁几何模型，设置井筒半径，内外半径比；划分有限元计算网格单元，采用平面应变等参单元和无限区域单元，最外一层设置成无限单元，它的尺度在规定的方向上延伸至无穷远，位移在无限远处为零，通过衰减函数的引入来实现；

$\stackrel{\&OverBar;}{N_i}=N_i{f}_i\left(r\right)---\left(1\right)$

衰减函数采用负指数函数：

$f_i\left(r\right)=e^{-(r-r_i)/L}---\left(2\right)$

剖分节点和单元，施加边界条件，施加情况为2个直边法向约束、径向自由；根据网格节点坐标，由式（3）计算初始应力场各单元高斯点的应力

${\mathrm{\&sigma;}}_r^0=\frac{1}{2}({\mathrm{\&sigma;}}_H+{\mathrm{\&sigma;}}_h)-\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\cos 2\mathrm{\&theta;}$

（3）

${\mathrm{\&tau;}}_{\mathrm{\&theta;}}^0=\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\sin 2\mathrm{\&theta;}$

(C)第一次应力释放(3)：采用“释放载荷”的方法来模拟实际钻井过程，求扰动场A；首先，根据公式(4)计算第一次释放正应力分量剪应力分量

根据公式(5)计算保留的应力分量

这样的应力释放方案，等效于式(3)的应力完全释放而同时施加钻井液压力；

${\mathrm{\&sigma;}}_r^{\left(1\right)}=-\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\cos 2\mathrm{\&theta;}$

（4）

${\mathrm{\&tau;}}_{\mathrm{\&theta;}}^{\left(1\right)}={\mathrm{\&tau;}}_{\mathrm{\&theta;}}^0=\frac{1}{2}({\mathrm{\&sigma;}}_H-{\mathrm{\&sigma;}}_h)\sin 2\mathrm{\&theta;}$

在r=a，井壁应力为

按弹性或弹-理想塑性计算，得扰动场A；通过第一次应力释放(3)，得到总场，这个总场是指初始场加扰动场A：

σ⁰+σ′→σ⁰；a′→a⁰     （6）

(D)第二次应力释放第一步(4A)：逐渐释放保留的应力分量，计算扰动应力场B；引入载荷参数λ，在r=a，井壁应力为τ′_θ＝0；

第一步用弹性方法确定弹性阶段的载荷参数增量Δλ₁，Δλ₁→λ_e，求出相应的弹性扰动场Δa′，Δσ′，得到总场，这个总场是指初始场加扰动场A再加弹性扰动场：

σ⁰+Δσ′→σ⁰；a⁰+Δa′→a⁰    （7）

(E)第二次应力释放第二步(4B)：从第二步开始，先令m=1，指定微弧长增量Δs_m，按弹-软化塑性，用弧长法迭代求解；全过程曲线τ-γ采用三线性形式，下降段斜率是一个负的常数G_T，称为切线剪切模量，而相应的软化曲线则为二线性形式，下降段斜率为G_p=GG_T/(G-G_T)；由于采用三线性模型，在D-P准则中k就是剪切强度τ_s，它也是初始的屈服应力τ_s(0)；设刚刚进入残余阶段的塑性应变为

那么二线性的软化曲线斜率为：

$\frac{{\&PartialD;\mathrm{\&tau;}}_s}{{\&PartialD;\mathrm{\&gamma;}}^p}=\left\{\begin{array}{c}\frac{{\mathrm{GG}}_T}{G-G_T},0\&le;{\mathrm{\&gamma;}}^p\&le;{\mathrm{\&gamma;}}_r^p\\ 0,{\mathrm{\&gamma;}}_r^p<{\mathrm{\&gamma;}}^p\end{array}\right.---\left(8\right)$

剪切强度τ_s和塑性功w^p的关系为：

${\mathrm{\&tau;}}_s\left(w^p\right)={[{\mathrm{\&tau;}}_s^2\left(0\right)+2G_p{w}^p]}^{1/2}---\left(9\right)$

在对应变软化岩石组成的结构应用有限元方法进行非线性分析时，一旦结构接近于极限承载能力或者以后，用一般的计算方法计算分析的稳定性与收敛性相当差以至于无法进行下去；弧长法通过同时约束荷载和位移向量，能较好地计算邻近极值点结构的反应和下降段问题，可以比较好地解决极值点穿越问题；

对于一个给定的迭代步n已知解

待求的下一迭代步n+1的解记为

可有如下的表达式：

${\mathrm{\&psi;}}_{m+1}^n\&equiv;\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)=P\left(a_{m+1}^n\right)-{\mathrm{\&lambda;}}_{m+1}^{n}R---\left(10\right)$

${\mathrm{\&psi;}}_{m+1}^{n+1}=\mathrm{\&psi;}(a_{m+1}^n,{\mathrm{\&lambda;}}_{m+1}^n)+{\left(K_T\right)}^n(a_{m+1}^{n+1}-a_{m+1}^n)-({\mathrm{\&lambda;}}_{m+1}^{n+1}-{\mathrm{\&lambda;}}_{m+1}^n)R---\left(11\right)$

于是，令

可给出位移解的增量Δaⁿ和载荷参数的增量Δλⁿ之间的关系

其中

将式（12）代入约束方程

${\left({\mathrm{\&Delta;a}}_m^n\right)}^T\left({\mathrm{\&Delta;a}}_m^n\right)+c{\left({\mathrm{\&Delta;\&lambda;}}_m^n\right)}^2={\left(\mathrm{\&Delta;s}\right)}^2---\left(14\right)$

可得一个关于Δλⁿ的二次代数方程

$b_1{\left({\mathrm{\&Delta;\&lambda;}}_m^n\right)}^2+b_2\left({\mathrm{\&Delta;\&lambda;}}_m^n\right)+b_3=0---\left(15\right)$

其中

$b_1={\left(\mathrm{\&Delta;}{\hat{a}}^n\right)}^T\left(\mathrm{\&Delta;}{\hat{a}}^n\right)+c$

从方程（15）可以得到Δλⁿ的两个根，分别记为(Δλⁿ)₁和(Δλⁿ)₂；为了能够跟踪解的路径得到正确的结果，只选取一组解，它应使解沿平衡曲线往前走，而不是往后返回，即取(Δaⁿ)^T(a_m-a_m-1)为一个最大值的解，得扰动场B，求出总场，这个总场是指初始场加扰动场A再加扰动场B：

${\mathrm{\&lambda;}}_m+{\mathrm{\&Delta;\&lambda;}}_m\&RightArrow;{\mathrm{\&lambda;}}_{m+1}$

$a_m^0+{\mathrm{\&Delta;a}}_m\&RightArrow;a_{m+1}^0$

(17)

${\mathrm{\&sigma;}}_m^0+{\mathrm{\&Delta;\&sigma;}}_m\&RightArrow;{\mathrm{\&sigma;}}_{m+1}^0$

${\mathrm{\&kappa;}}_m^0+{\mathrm{\&Delta;\&kappa;}}_m\&RightArrow;{\mathrm{\&kappa;}}_{m+1}$

(F)特征值评价(5)：引入计算总刚矩阵特征值的算法-子空间迭代法，求

的最小特征值(μ₁)_m，根据得到的最小特征值(μ₁)_m进行评价，当最小特征值(μ₁)_m>0，时，如果m<M，M为弧长法微弧长增量总量；则m+1→m，即m加上1以后进入第二次应力释放第二步(4B)，否则，得到的即为总场；当最小特征值(μ₁)_m<0时，进行如下计算；

λ_m→λ_N,λ_m-1→λ_N-1

(γ₁)_m→(μ)_N,(μ₁)_m-1→(μ)_N-1          (18)

σ_m→σ_N,σ_m-1→σ_N-1

κ_m→κ_N,κ_m-1→κ_N-1

使用线性插值计算临界扰动应力场：

${\mathrm{\&Delta;\&lambda;}}_{\mathrm{cr}}={\mathrm{\&Delta;\&lambda;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}({\mathrm{\&Delta;\&lambda;}}_N-{\mathrm{\&Delta;\&lambda;}}_{N-1})$

${\mathrm{\&Delta;\&sigma;}}_{\mathrm{cr}}={\mathrm{\&Delta;\&sigma;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}({\mathrm{\&Delta;\&sigma;}}_N-{\mathrm{\&Delta;\&sigma;}}_{N-1})---\left(19\right)$

${\mathrm{\&Delta;\&kappa;}}_{\mathrm{cr}}={\mathrm{\&Delta;\&kappa;}}_{N-1}+\frac{{\mathrm{\&mu;}}_{N-1}}{{\mathrm{\&mu;}}_{N-1}-{\mathrm{\&mu;}}_N}({\mathrm{\&Delta;\&kappa;}}_N-{\mathrm{\&Delta;\&kappa;}}_{N-1})$

计算总场，这个总场是指初始场加扰动场A再加临界扰动场：

Δλ_cr+λ_m→λ_cr

Δσ_cr+σ_m→σ_cr             (20)

Δκ_cr+κ_m→κ_cr

(G)广义力和广义位移(6)：根据弧长法每步骤计算出的载荷因子λ，按照两次解除应力的公式通过积分即可求出当步广义力及广义位移；

在计算总场时，外力功的变分为

$\mathrm{\&delta;W}=\mathrm{F\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}\left[q_r{\mathrm{\&delta;u}}_r\right]\mathrm{ad\&theta;}={\&Integral;}_0^{\mathrm{\&pi;}/2}(1-\mathrm{\&lambda;}){\mathrm{\&sigma;}}_r^{\left(2\right)}{\mathrm{\&delta;u}}_r\mathrm{ad\&theta;}---\left(21\right)$

如果广义力可定义为

$F=\frac{a}{2}(1-\mathrm{\&lambda;})({\mathrm{\&sigma;}}_H+{\mathrm{\&sigma;}}_h)---\left(22\right)$

那么广义位移的变分为

$\mathrm{\&delta;U}={\&Integral;}_0^{\mathrm{\&pi;}/2}{\mathrm{\&delta;u}}_r\mathrm{d\&theta;}---\left(23\right)$

第m个载荷步的增量广义位移为

${\mathrm{\&Delta;U}}_m=\underset{i}{\mathrm{\&Sigma;}}{\mathrm{\&Delta;u}}_r\left({\mathrm{\&theta;}}_i\right)\mathrm{\&Delta;\&theta;}---\left(24\right)$

总广义位移可累加得到：

U_m+1=U_m+ΔU_m      （25）

(H)平衡路径曲线(7)：平衡路径曲线(7)是由广义力和广义位移为坐标轴绘制的曲线；根据广义力和广义位移(6)的输出结果，绘制平衡路径曲线；

(I)坍塌压力处理(8)：用P和u分别代表系统的广义力和广义位移，那么在稳定分支上有δPδu>0，在不稳定分支上有δPδu<0，稳定分支和不稳定分支的交汇点是临界点，在该点δPδu=0；在平衡路径曲线上找出临界力，根据广义力的公式（22）确定对应的应变软化极限载荷因子λ_cr；

根据得到的应力总场、弹性极限载荷参数λ_e以及平衡路径曲线输出的λ_cr，进入坍塌压力处理(8)，得到弹性坍塌压力和应变软化塑性坍塌压力：

$(1-{\mathrm{\&lambda;}}_e){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_e$

(26)

$(1-{\mathrm{\&lambda;}}_{\mathrm{cr}}){\mathrm{\&sigma;}}_r^{\left(2\right)}\&RightArrow;q_{\mathrm{cr}}$

(J)输出应变软化坍塌压力q_cr，应变软化极限载荷参数λ_cr对应的塑性区分布图，确定该井的防塌钻井液的密度；

符号说明：-第i个节点的位移插值函数；N_i-第i个节点的坐标插值函数；L-衰减特征长度，m；r_i-第i个节点位移，m；

-初始场径向应力，MPa；-初始场切向应力，MPa；σ_H-最大水平主应力，Mpa；σ_h-最小水平主应力，Mpa；-第一次释放正应力分量，Mpa；

-第一次释放剪应力分量，Mpa；θ-与最大水平地应力的夹角，°；-第一次保留的应力分量，Mpa；λ-载荷参数，无量纲；r-距井眼中心距离，m；a-井眼半径，m；σ⁰-初始场应力，Mpa；σ'-扰动场应力，Mpa；σ'_r-扰动场径向应力，Mpa；τ'_q-扰动场切向应力，Mpa；a'-扰动场位移，m；a₀-初始场位移，m；Δλ₁-弹性阶段的载荷参数增量，无量纲；λ_e-弹性极限载荷参数，无量纲；τ-剪应力，MPa；γ-剪应变，%；G-弹性剪切模量，MPa；GT为切线剪切模量，MPa；G_p-切线塑性剪切模量，MPa；τ_s(0)-初始的屈服应力，Mpa；τ_s-剪切强度，MPa；w^p-塑性功，MPa；γ^p为塑性应变，%；

为残余阶段的塑性应变，%；M-弧长法微弧长增量总量，无量纲；Δs_m-微弧长增量，m；m-弧长增量迭代步；Δaⁿ-位移增量向量，m；Δλⁿ-载荷增量系数，无量纲；

-切线刚度矩阵；R-载荷力向量，MPa；ψ_m-不平衡力向量，MPa；c-调节因子，无量纲；κ-内变量，取塑性功w^p为内变量，Mpa；(μ₁)_m-第m步的最小特征值，无量纲；q_e弹性坍塌压力，Mpa；F广义力，MPa；U广义位移，m²；q_cr应变软化坍塌压力，MPa；λ_cr-应变软化极限载荷参数，无量纲；q_r-井壁上的载荷，MPa；δu_r-井壁上的虚位移，m。

2.根据权利要求1所述的一种井壁坍塌分析方法，其特征在于：所述的井壁岩石的材料参数包括：弹性模量、泊松比、内摩擦角、内聚力、软化系数、屈服函数类型、内变量类型、水平最大主应力、水平最小主应力和深度。

Images