知道岩石何时屈服， 为何还算不出变形？
PETROLEUM PAPER · 高级地质力学
L31 / 流动法则与塑性应变
知道岩石何时屈服，
为何还算不出变形？
从屈服面上的法向箭头，到可以核算的应变增量：讲清关联流动、非关联流动、剪胀角与塑性乘子的作用。
假设已经判断井壁附近一小块岩石达到屈服。接下来，它会沿哪个方向产生不可恢复的变形？会压密还是膨胀？变形量又有多大？
这些问题无法只靠一个强度判据回答。屈服准则确定边界，流动法则确定塑性应变增量的方向；增量的大小，还要由加载条件、弹性关系和一致性条件共同确定。把这几件事分开，才有可能建立完整的弹塑性计算。
本讲的符号与适用范围
采用压应力、压缩应变为正；σ 表示有效应力，主应力按 σ₁ ≥ σ₂ ≥ σ₃ 排序。讨论小应变、率无关的弹塑性响应，不计温度和化学应变。εp 上标 p 表示塑性；d 或 Δ 表示增量，绝不是把已有的全部应变重新计算一次。
01 / 完整模型屈服以后，还缺哪几条关系？
图 01｜计算塑性变形需要同时组织屈服条件、硬化关系、应变分解、流动法则和卸载条件。仅有一条强度包络线，还没有把应变问题闭合。
配图：D. Nicolas Espinoza 对应公开课程讲义，Figure 4.10 局部。
在弹性阶段，已知弹性参数和应力增量，就能求出相应的弹性应变增量。进入塑性阶段后，总应变中还包含不可恢复的部分，必须先进行分解：
dε = dεe + dεp
dσ = C : dεe
C 是弹性刚度张量；冒号表示相应分量的求和运算。本讲把 C 用作刚度；原板书中的 D 在 ε = Dσ 关系中表示柔度，二者不能因字母不同而混用。
因此，即使两种模型采用完全相同的屈服面，只要塑性势或硬化规律不同，也可能给出不同的位移、体积变化和残余变形。强度拟合接近，不代表变形预测也已经被验证。
02 / 流动法则法向箭头决定的是应变分量比例
常用的塑性流动法则写成：
dεijp = dλ · ∂g / ∂σij
dλ ≥ 0
g 是塑性势函数，dλ 是塑性乘子增量。偏导数给出流动方向，dλ 控制本步的幅度；它不是可以任意指定的岩石常数。
把 g 看成应力空间中的一个等值面，它的梯度指向该面的法向。如果 g 与屈服函数 f 相同，就得到关联流动法则；若用另一个塑性势控制方向，则称为非关联流动法则。
这里的“法向”位于应力空间。箭头的分量表示塑性应变增量的比例，不能把它直接当成真实岩体中的裂缝走向或剪切带倾角。图的横纵坐标若采用不同缩放，视觉上的直角也不能代替偏导数计算。
为什么公式里不是单位法向量？
梯度不一定长度为 1。若把 g 乘以一个正常数，方向不会改变，但要相应改变 dλ 才能保持同一个 dεp。因此，报告塑性乘子数值时，必须同时说明塑性势的具体定义。
03 / 关联流动用 Mohr–Coulomb 算出三个应变分量
先考虑主应力严格满足 σ₁ > σ₂ > σ₃、且当前点位于莫尔—库仑曲面的一个光滑面上的情形。选用如下归一化：
f = σ₁ − Nφσ₃ − UCS
Nφ = (1 + sin φ) / (1 − sin φ)
UCS = 2c cos φ / (1 − sin φ)
c 为黏聚力，φ 为内摩擦角。f 与 UCS 的单位均为应力。板书用 q 表示这里的 Nφ；本讲使用 Nφ，避免与后面偏应力强度 q 混淆。
图 02｜红线表示选定主应力平面中的屈服边界，蓝色法向箭头给出关联流动下的应变增量方向。右侧岩样中的破坏面角度属于物理空间，与左侧法向箭头的含义不同。
配图：D. Nicolas Espinoza 对应公开课程讲义，Figure 4.11 局部。
对于关联流动，取 g = f。分别对三个主应力求偏导：
∂f/∂σ₁ = 1
∂f/∂σ₂ = 0
∂f/∂σ₃ = −Nφ
梯度向量 n = (1, 0, −Nφ)。它给出方向，不是本步已经发生的塑性应变。
于是三个主方向的塑性应变增量为：
dε₁p = dλ
dε₂p = 0
dε₃p = −Nφdλ
dεvp = (1 − Nφ)dλ
dεᵥᵖ 是三个塑性正应变增量之和。在本讲压缩为正的约定下，dεᵥᵖ < 0 表示塑性体积膨胀，即剪胀。
若 φ = 30°，则 Nφ = 3。只要 dλ > 0，塑性应变分量比例就是 1 : 0 : −3，塑性体积应变增量为 −2dλ。也就是说，第一方向压缩，第三方向伸长，而且伸长量更大，从而产生剪胀。
图 03｜理想塑性阶段中，塑性轴向应变与侧向应变继续发展，并产生相应的体积变化。曲线用于解释机制；应变正负应以正文约定为准。
配图：D. Nicolas Espinoza 对应公开课程讲义，Figure 4.12 局部。
dε₂ᵖ = 0 不是普遍的平面应变条件
这个结果来自特定光滑面上的塑性势偏导。它既不说明总应变 dε₂ 必然为零，也不能直接套到 σ₂ = σ₃ 的轴对称棱线上；在棱线或顶点，可能需要组合多个活动面的流动方向。
04 / 非关联流动让摩擦角管强度，让剪胀角管体积变化
关联莫尔—库仑流动把强度参数 φ 与体积膨胀趋势绑定在一起。对于一些岩石和土体，这会给出偏大的剪胀。一个常见处理是：屈服函数仍使用 φ，塑性势另引入剪胀角 ψ。
g = σ₁ − Nψσ₃
Nψ = (1 + sin ψ) / (1 − sin ψ)
dεvp = (1 − Nψ)dλ
塑性势中可以加上与应力无关的常数，其应力梯度不变，所以这里省略该常数。g 的具体数值无需等于零；流动方向取决于它在当前点的梯度。
图 04｜剪胀角的正、零、负，分别对应这类塑性势中的剪胀、等体积塑性流动和剪缩。ψ 应由相关变形数据约束，不能仅用强度数据替代。
配图：D. Nicolas Espinoza 对应公开课程讲义，Figure 4.13 局部。
一个容易忽略的比较前提
当 φ = 30° 时，关联情形有 dεᵥᵖ/dλ = −2；若改取 ψ = 10°，则 Nψ ≈ 1.4203，比值变为 −0.4203。
这比较的是同一 dλ 下的方向比例。若对同一外部加载重新求解，dλ 本身也可能改变，不能据此直接宣称实际体积膨胀必然减少某个固定百分比。
05 / 塑性乘子方向知道了，幅度还要靠一致性条件
率无关塑性常用以下互补条件组织加载与卸载：
f ≤ 0，dλ ≥ 0，f dλ = 0
f < 0 时，必须有 dλ = 0；f = 0 时，只说明当前位于边界上，仍要区分继续塑性加载、切向加载或弹性卸载。
因此，岩石曾经达到屈服，并不意味着此后每一步都继续产生塑性应变。卸载进入弹性域时，已有塑性应变可以保留，但本步塑性增量为零。若屈服面发生硬化，塑性加载过程还必须满足随硬化变量一起变化的一致性条件。
先做弹性试算，再校正应力
对给定的总应变增量 Δε，先假设本步全为弹性，得到试算应力 σᵗʳ = σⁿ + C : Δε。如果它越过屈服边界，就扣除塑性部分对应的应力贡献：
σn+1 = σtr − Δλ C : m
m = ∂g/∂σ
n 与 n+1 表示相邻计算步，tr 表示弹性试算。塑性应变沿 m 发展，而应力校正沿 C:m 进行；二者不应混为同一个几何方向。
对于理想塑性、线性屈服面、常量梯度且没有跨越棱线的一步，令 n = ∂f/∂σ，代入 f(σⁿ⁺¹) = 0，可以直接得到：
Δλ = f(σtr) / (n : C : m)
这是本讲算例的闭合关系。若存在硬化、非线性势函数、主方向变化或多个活动面，就需要相应的一致性方程与更新算法，不能原样套用这个单面公式。
06 / 完整算例从总应变增量，算到塑性应变
假设条件 / 全部数值仅用于教学
采用小应变、各向同性线弹性—理想塑性模型，孔压保持不变；本步三个主方向固定，且应力点始终处于同一光滑面。
E = 10 GPa = 10 000 MPa，ν = 0.25；
c = 10 MPa，φ = ψ = 30°，采用关联流动。
因此 Nφ = 3，UCS ≈ 34.6410 MPa。
第 1 步：给定初始应力和总应变增量
σn = (64.6410，30，10) MPa
Δε = (0.001，0，−0.0002)
初始点满足 64.6410 − 3 × 10 − 34.6410 = 0，位于屈服面。第二方向的总应变增量为零，是本算例给定的加载条件。
由 E、ν 可得三个正应变方向之间的刚度分量：Cᵢᵢ = 12 000 MPa；当 i ≠ j 时，Cᵢⱼ = 4 000 MPa。本步没有剪切应变，因此只使用这三个主方向之间的关系。
第 2 步：检查弹性试算点
C : Δε = (11.2，3.2，1.6) MPa
σtr = (75.8410，33.2，11.6) MPa
ftr = 75.8410 − 3 × 11.6 − 34.6410
= 6.4 MPa > 0
试算点越过屈服面，需要进行塑性校正。不能把这个纯弹性外推点直接当作本步最终应力。
第 3 步：求塑性乘子与塑性应变
关联流动下 n = m = (1，0，−3)。计算得到 C:m = (0，−8 000，−32 000) MPa，因此 n:C:m = 96 000 MPa。
Δλ = 6.4 / 96 000
= 6.6667 × 10−5
Δεp = (0.000066667，0，−0.0002)
本算例中 g 的单位为应力，所以其应力梯度无量纲，Δλ 具有应变的量纲。第三方向出现塑性伸长。
第 4 步：检查应变分解和屈服一致性
将总应变、塑性应变与弹性应变放在一起比较。下表统一使用微应变，1 με = 10⁻⁶：
由校正式得到最终应力约为 (75.8410，33.7333，13.7333) MPa。代回 f = σ₁ − 3σ₃ − UCS，得到约 0，剩余微小差值来自小数取舍；三个主应力的排序也仍满足假设。
算例的物理含义
本步塑性体积应变增量为 0.000066667 − 0.0002 = −0.000133333，即约 −0.01333%。按压缩为正的约定，它表示塑性部分发生体积膨胀。
但本步总的体积应变增量为 0.001 − 0.0002 = 0.0008，仍为压缩。塑性剪胀与总体压缩可以同时发生，因为还存在弹性体积变化。
07 / 工程理解同一块岩石，为什么可能先剪胀、再剪缩？
图 05｜在平均有效应力 p′ 与偏应力强度 q 的示意空间中，不同位置的流动方向可分别对应剪胀、近似等体积和剪缩。底部颗粒图用于解释不同约束下的结构变化。
配图：D. Nicolas Espinoza 对应公开课程讲义，Figure 4.13 局部。
这里 p′ = I₁/3，q = √(3J₂)。低约束条件下，颗粒滑动和绕越可能伴随孔隙张开；较高约束下，颗粒重排、破碎或孔隙塌缩可能使体积收缩。因此，变形方向与应力状态、结构和加载历史有关。
图中的弯曲边界表达了状态改变时流动方向可以转变这一思想。它不意味着固定 ψ 的线性莫尔—库仑塑性势会自动产生整条剪胀—剪缩转换曲线。要描述这种演化，可能需要状态相关的剪胀规律，或包含压密机制的塑性势与屈服面。
对于井壁变形、储层压密或剪切带附近的位移预测，仅用峰值强度数据不足以验证流动法则。应同时测量轴向和侧向应变，区分弹性与塑性体积变化，并记录围压、排水条件与加载路径。这样才能判断 ψ 或更一般的塑性势是否合理。
本讲所用弹性试算、塑性势法向及应力校正的基本关系，可进一步核对 Itasca 的 Mohr–Coulomb 模型说明。不同软件可能采用相反的应力符号与排序，移植公式前要先统一约定。
08 / 课堂自检用五个问题检验你的理解
1．只要 f = 0，就一定有新的塑性应变吗？
不一定。还要判断加载方向。处在屈服面上，也可能开始弹性卸载，此时本步 dλ = 0。
2．ψ = 0，是不是意味着材料仍处于弹性？
不是。在这里的塑性势中，它表示塑性体积应变增量为零。塑性剪切变形仍可发展。
3．φ 与 ψ 的区别是什么？
φ 参与定义强度边界，ψ 参与定义塑性流动方向。在关联流动的这一特例中二者相同；在非关联流动中，应分别用相应数据约束。
4．dε₂ᵖ = 0，能否说明整个问题是平面应变？
不能。这只是在某一光滑面上得到的塑性分量结果。总应变还包含弹性部分，且棱线处的处理可能不同。
5．塑性体积应变为负，总体积就一定增大吗？
不一定。按本讲的符号，它只说明塑性部分剪胀。总的体积变化还取决于弹性部分；本讲完整算例中，总体仍然压缩。
掌握流动法则以后，看一幅塑性区图就可以多问一步：这里的塑性应变方向由什么决定，幅度如何求得，体积变化是否有独立数据支持？把这三个问题回答清楚，模型才真正从“判断是否屈服”走向“预测留下多少变形”。
课程对应主题：D. Nicolas Espinoza，Advanced Geomechanics（PGE 383），L31：Determination of plastic strains with the flow rule。本篇依据对应公开课程讲义整理；配图采用作者板书原图局部，公式用可复 制文字重排。完整数值算例为独立教学演算，假设条件已在正文说明。
原图及讲义：Advanced Geomechanics · Inelasticity，第 4.3 节 Plastic strains，Figures 4.10—4.13。
模型关系核对：Itasca · Mohr–Coulomb Model；Plastic Model Group。
