知识点 弹性模型与开挖模型 MC 与 DP 模型 应变硬化与软化模型 双屈服、剑桥、霍克布朗模型 梁、衬砌、锚索和桩单元 二维单元 三维单元 本章导读 本章主要介绍了有限差分数值模拟原理，先介绍基本原理，包括空间导数的有限差分近 似、运动平衡方程、应变、应力及节点不平衡力和阻尼力的差分形式；再介绍本构方程和有限 差分方程；介绍有限差分数值模拟分析的求解步骤；最后介绍了基于 MORH-COULOMB 塑性 模型的增量弹性理论、屈服准则、流动法则、塑性应力调整的有限差分形式。
3.1 弹性模型和开挖模型
本节主要介绍弹性模型和开挖模型，其中弹性模型包括各向同性和横观各向异性两种 情况。 3.1.1 弹性模型 1．各向同性弹性模型 各向同性弹性模型提供了材料性质最简单的表达方式，这种模型适用于应力应变特性呈 线性关系的无卸载滞后现象的均质、各向同性、连续性材料。 在此种本构模型中，应力应变关系依胡克定律，用应变增量的形式表达。其表面应变的 表达式为： 11 1 11 2 22 22 2 12 1 22 12 12 21 12 33 2 11 22 2 ( ) ( ) e e e e G e e e e e （3-1） 式中，1K(4 / 3)G；2K(2 / 3)G；K 为体积模量；G 为剪切模量。1 2 j i ij j i u u e t x x （3-2） 式中，eij为应力张量增量；u 为体积模量； t 为时步。 在平面应力下，这些等式为： 11 1 11 2 22 22 2 11 1 22 12 12 21 12 33 2 ( ) 0 e e e e G e （3-3） 式中， 2 1 1 ( 2/ 1) ；2 2(22/1)。 对于轴对称几何体： 11 1 11 2 22 33 22 1 22 2 11 33 12 12 21 12 33 1 33 2 11 22 ( ) ( ) 2 ( ) ( ) e e e e e e G e e e e （3-4） 2．横观各向异性弹性模型 各向异性弹性模型适用于模拟在各层的法向方向和切向方向的弹性模量有明显差异的层 状弹性介质。 在 FLAC 中，各向同性平面位于这个模型的 xz 平面，各弹性模量的定义如下： E1（Ex）各向同性面的弹性模量； E2（Ey）与各向同性面垂直平面的弹性模量； G12（Gxy）各向同性面垂直面的交叉剪切弹性模量； G13（Gxz）各向同性面的剪切模量； v12（vyx）由 y 方向的单轴应力引起的各向同性面的 x 方向上的法向应变相关垂直面上的 y 方向上的法向应变的泊松比； v31（vzx）由 z 方向的单轴应力引起的各向同性面的 x 方向上的法向应变相关垂直面上的 z 方向上的法向应变的泊松比； 弹性性质的变化是受限制的（Amadei 1982），因此应用了以下限制条件： 2 2 2 0 0 0 1 1 2 (1 ) 0 x y xy xy xz x yz xz y E E G E E ≤ ≤ ≥ （3-5） 对于常见的正交各向异性弹性体，应力—应变关系式由 Lekhnitskii 给出：
11 11 11 12 22 13 33 16 12 22 12 11 22 22 23 33 26 12 33 13 11 23 22 33 33 36 12 23 44 23 45 13 13 45 23 55 13 12 16 11 26 22 36 33 6 1 ( ) ( ) 2 1 ( ) ( ) 2 1 ( ) ( ) ( ) ( 2 e S S S S e S S S S e S S S S e S S e S S e S S S S
6 12)
（3-6） 其中 Sij（i，j=1，2，3，4，5，6）是与 Ei（i=1，2，3），G12，G13，G23和有关的变量， 其中是从 x 轴逆时针算起的各向异性角。 将 xy 平面进行如下设置，可以获得平面应力状态： 33 13 23 0 根据以上设置，等式（3-6）可以写成如下的矩阵形式： 11 11 12 16 11 22 12 22 26 22 33 16 15 66 12 e s s s e s s s e s s s （3-7） 通过矩阵转置可以很容易地得出大批应力应变关系。 将 xy 平面进行如下设置，可以获得平面应变状态： 33 13 23 0 e e e 这种假设可以获得应力应变的如下关系式： 11 12 13 16 11 11 12 22 23 26 22 22 33 13 23 33 36 12 16 26 36 66 12 0 2 s s s s e s s s s e s s s s e s s s s （3-8） 通过矩阵转置也可以得出到应力应变关系。 3.1.2 开挖模型Null Model 模型用于模拟隧道工程中围岩土体的开挖，即围岩材料一旦被赋予 Null Model 后，表示该部分土体被移去或者开挖掉了。同时，在 Null Model 模型区材料的应力自动设置 为 0。
3.2 M-C 和 D-P 塑性模型
在有限元模拟分析中，采用 M-C 和 D-P 塑性模型来模拟隧道结构周围的土体、加固圈和 注浆材料，以下分别介绍这两种模型。
3.2.1 Morh-Coulomb 塑性模型 1．增量弹性理论 在 FLAC 中，摩尔－库仑模型采用主应力1、2、3和平面外应力zz表示。主应力和 主方向由应力张量分量计算（以压应力为负）。 1 2 3 ≤ ≤ （3-9） 相应的主应变增量 、e1 e2、e3分解为ei eie eip。这里上标 e 和 p 分别表示弹性 和塑性部分，塑性分量只在塑性流动阶段不为零。胡克定律的主应力和主应变的增量表达式为：
1 1 1 2 2 3 2 1 2 2 1 3 3 1 3 2 1 2 e e e e e e e e e e e e e e e e e e （3-10） 式中，1K(4 / 3)G，2K(2/3)G。 2．屈服准则和流动法则 基于式（3-9）的假设，在应力空间和(1，3)平面的破坏准则可表示为图 3-1 的形式。 由摩尔－库仑屈服函数定义的从 A 点到 B 点的破坏包络线为： 1 3 2 s f N c N （3-11） （a） （b） 图 3-1 岩土材料 Mohr-Coulomb 模型及破坏准则由 B 点到 C 点拉应力屈服函数的定义为 ft t 3， 1 sin 1 sin N ， max tan c ，为 摩擦角， c 为凝聚力，t为抗拉强度，材料强度不超过 max 。 剪切势函数g 对应于非关联的流动法则，即：s 1 3 s g N。 势函数g 对应于拉应力破坏的相关联流动法则，即：t 3 t g 。 对于剪切应力和拉应力处于边界的情况，可由摩尔—库仑流动法则，并通过定义三维应 力空间中边界附近的混合屈服函数进行计算。定义函数h(1，3) 0= 用以表示(1，3)平面中 0 s f 和f t 0所代表曲线的对角线，该函数的表达式为： 3 ( 1 ) t p p h （3-12） 这里，p和p为两个常量， p 1 N2 N ，p tN 2c N 。 弹性假设和破坏准则不一样，分别在(1，3)平面中位于 1 区域和 2 区域（对应于h 0区 域内或+区域），如图 3-2 所示。如果位于 1 区，则属于剪切破坏，应用由势函数g 确定的流s 动准则，应力回归到 f 的曲线上；如果位于s 0 2 区，则属于拉应力破坏，应用由势函数g 确t 定的流动准则，应力点回归到 f t 0的曲线上。 图 3-2 流动准则的区域定义 3．塑性应力调整 首先是剪切破坏，其流动法则为 s p s i i g e （i=1，3），这里 s 为待定参数，通过偏微 分后，此式变为： 1 2 3 0 p s p p s e e e N （3-13） 弹性应变增量可以从总增量中减去塑性增量，进一步利用上式的流动法则，式（3-10）中 的弹性法则变为： 区域 2 区域 1 + h = 0 3 1 0 s f 0 t f + +
1 1 1 2 2 3 1 2 2 1 1 2 1 3 2 3 1 1 2 1 2 1 2 ( ) ( ) ( ) (1 ) ( ) ( ) s s s e e e N e e e N e e e N （3-14） 让新旧的应力状态分别由上标 N 和 O 表示，然后通过定义： N O i i i （3-15） 用此式代替式（3-14），并用上标 I 表示由弹性假设得到的应变和原应变之和，则总应变 计算得到的弹性增量为： 1 1 1 1 2 2 3 2 2 1 2 2 1 3 3 3 1 3 2 1 2 ( ) ( ) ( ) I O I O I O e e e e e e e e e （3-16） 对于拉应力破坏的情况，流动法则为： 1 t p t i g e （i=1，3） （3-17） 这里t是待定的参数，用式（3-17）中的g ，通过偏微分后，此式变为：t 1 2 3 0 0 p p p t e e e （3-18） 重复上面相似的原理，可得到： 1 1 2 2 2 2 3 3 1 N I t N I t N I t （3-19） 其中： 3 1 ( ) t I t f （3-20） 3.2.2 Drucker-Prager 塑性模型 德鲁克—普拉格（Drucker-Prager）塑性模型非常适用于模拟摩擦角较小的软粘土，但此 模型在岩土工程材料中应用较少，将它包括进 FLAC 程序中主要是用来同其他程序进行比较。 这种模型的破坏包络线包括 Drucker-Prager 准则和拉应力路径。与剪切流动法则不相关联 而与拉伸法则相关联。 1．增量弹性法则 德鲁克—普拉格模型由两个应力分量表示——剪切应力 和平均法应力，两者分别定义为： 2 11 22 33 1 ( ) 3 J （3-21）
这里J 是偏应力张量第二不变量，它可以表示为： 2 2 2 2 2 2 11 22 22 33 11 33 12 1 [( ) ( ) ( ) ] 6 J （3-22） 同 和有关的剪应变增量 和体积应变增量 e 可表示为： 2 11 22 33 2 J e e e e （3-23） 2 J 是偏应变增量第二不变量，其表达式为： 2 2 2 2 2 11 22 22 33 11 33 22 1 [( ) ( ) ( ) ] 6 J e e e e e e e （3-24） 应变增量可分解为： e p e p e e e （3-25） 上标 e 和 p 分别表示弹性和塑性部分，并且塑性分量仅用在塑性流动阶段不为零的情况。 用通常的应力应变表示的胡克定律的增量表达式为 e e G K e （3-26） G 和 K 分别表示剪切模量和体积模量。 2．屈服函数和势函数 主应力空间中 , 平面内的破坏准则如图 3-3 所示。 图 3-3 德鲁克—普拉格模型破坏准则 由 Drucker-Prager 屈服函数确定的从点 A 到点 B 的破坏包络线为： s f q k （3-27） 由拉应力屈服函数确定从点 B 到点 C 的破坏包络线为： t t f （3-28） 式中， q和 k为有关材料特性的常数。 拉应力强度不能超过 maxt k q （3-29）
剪切势函数g 通常对于不相关联的流动法则，其表达式如下：s s g q （3-30） 式中， q为材料常数，当采用相关流动法则时，其值等于 q。 拉应力破坏对于相关联的流动法则，其表达式来源于势函数g ：t t g （3-31） 对于边界附近，Drucker-Prager 的流动法则有如下所示的表达方法，定义函数 h( , )=0 用以表示（ , ）平面，见图 3-4 中f 和s 0 f t 0所代表曲线的对角线，此函数的表达式为： ( ) p p t h （3-32） p 和p两个常量的定义如下： 2 1 p t p k q q q （3-33） 图 3-4 流动准则的区域定义 弹性假设和破坏准则不一致，分别在 ( , ) 平面中位于 1 区和 2 区（分别对应于 h=0 区域 内正或负区域）。如果位于 1 区，说明是剪切破坏，应用由势函数 s g 确定的流动法则，应力点 回归到f 的曲线上。如果位于s 0 2 区，说明是拉应力破坏，应用由势函数g 确定的流动法t 则，应力点回归到 f t 0的曲线上。 3．塑性修正 首先考虑塑性破坏，流动法则如下： s p s s p s g g e （3-34） 这里s是待定义的参数，用式（3-29）的g ，通过偏微分后，此式变为：s
p s p s e q （3-35） 弹性应变增量可以从式（3-25）表示的总增量减去塑性增量，进一步用式（3-35）的流动 法则，式（3-25）中的弹性法则变为： s s G G k e kq （3-36） 让新旧的应力状态分别由上标 N 和 O 表示，然后通过定义： N O N O （3-37） 上标 I 表示由弹性假设得到的应变和原来应变之和，由总应变计算得到的弹性增量即为： I O I O G K e （3-38） 现在可以定义： ( ) N N I I N ij ij I ij （3-39） 其中ij为克罗内克（Kronecker）符号。 现在考虑拉应力破坏，流动法则的形式为： t p t t p t g g e （3-40） 这里t是待定参数，用式（3-30）中的g ，通过偏微分后，此式变为：t 0 p p t e （3-41） 重复上面相似的推理，可得到： N I N I K t （3-42） 其中： I t t K （3-43） 在该模型的破坏模式中，对应于弹性假设的新的应力偏量可以写成： ( ) N I t I ij ij ij （3-44）
3.3 应变硬化—软化模型
本节将主要介绍应变硬化与软化模型、节理化模型和双线性应变硬化与软化节理模型。3.3.1 应变硬化与软化模型 在 FLAC 中，这种模型是基于剪切流动法则不相关联而与拉力流动法则相关联的摩尔 —库仑模型，差别在于塑性屈服开始后，粘聚力、摩擦角、剪胀扩容和抗拉强度可能发生 变化。在摩尔—库仑模型中，这种性质都假定保持为常量。用户可以自己定义粘聚力、摩 擦角和剪胀为硬化参数的分段线性函数，这些硬化参数量测塑性剪切应变。抗拉强度软化 法则也可被设定为其他量测塑性拉应变的硬化参数的分段性形函数。这种模型通过在每个 时步增加硬化参数以计算总的塑性剪切应变和拉应变，并以此促使材料性质同用户定义的 函数保持一致。 对于屈服函数和势函数，塑性流动法则和应力修正同摩尔—库仑模型完全一致。 1．硬化—软化参数 塑性剪切应变由剪切硬化参数 ps e 计算，eps的增量形式的定义如下（Vermeer.deBorst 1984）： 1 2 2 2 2 1 3 1 1 1 ( ) ( ) ( ) 2 2 2 ps ps ps ps ps ps m m m e e e e e e （3-45） 其中： 1 3 1 ( ) 3 ps ps ps m e e e emps（m=1，2，3）是塑性剪切应变主增量。 抗拉硬化参数ept用于计算累积的张拉塑性应变，它的增量定义为： 3 pt pt e e （3-46） 这里e3pt是主应力方向上的张拉塑性应变增量（拉应力为正）。 以上用到的以及以后在相似的表达式中将用到的一些符号值得进行一些说明。 ps m e 同式 （3-34）中定义的 p m e 是一致的，这里 m=1，2，3。增加上标 s 表示塑性应变同剪切屈服面（而 不是拉伸屈服面 0）相关。注意 ps i e 是塑性主应变增量而不是剪切应变增量。相似地，e3pt同 式（3-39）中定义的e3p是一致的，上标 t 表示塑性应变是同拉伸屈服面相关的。 2．用户自定义函数的材料模型 考虑一维的应力—应变曲线 ，它在达到屈服时开始软化但仍保留一定的残余强度。 e 达到屈服点之前，曲线是线性的，在此阶段，只产生弹性应变，eee，材料屈服后，总 应变由弹性应变和塑性应变两部分组成，eeeep。在软化/硬化模型中，用户定义粘聚力、 摩擦角、剪胀角和抗拉强度这些变量作为总应变中塑性应变部分e 的函数，在 FLAC 中，这p 些函数也可以由各线性线段近似表示。 粘聚力、摩擦角和剪胀角的硬化和软化特性，由用户以表格的形式通过剪切参数e （式ps 3-45）表示，每个表格包含一对值：一个表示参数，另一个表示对应的特性值。它假定特性变 量在表格中的两个连续的参数目之间是线性的。抗拉强度的软化通过同样的方法由参数ept （式 3-46）来表示。
3.3.2 节理化模型 节理化模型是各向异性塑性模型，它包含在摩尔—库仑体内特殊方向上的弱面。根据应 力状态，弱面走向以及模型体和弱面的材料特性不同，屈服可能发生在模型体内，或者发生在 弱面上，或者在两个部位同时发生。 这种模型在 FLAC 中的实现方法是首先判别总体破坏，同时应用与 FLAC 中摩尔—库仑 模型中相同的相关塑性修正。弱面内的破坏准则存在于包含在拉应力路径的摩尔—库仑屈服条 件的局部形式中，与局部剪切流动法则不相关联而与局部拉应力流动法则相关联。 为了简化这部分符号，我们定义ij，它对应于由各阶段总体破坏的塑性修正的应用引起 的应力分量，这种总应力分解成局部应力后可表示为： 2 2 11 11 12 22 2 22 11 12 22 33 33 2 2 12 11 12 12
cos 2 sin cos sin sin 2 sin cos cos
( ) sin cos (cos sin )
Q Q Q Q Q Q Q Q Q Q Q Q （3-47） 式中，Q 为节理夹角（从 x 坐标轴逆时针方向算起）。 依照约定，由表示弱面上切向引力分量的大小，相应的应变变量为 ，可以得到： 12 12 e （3-48） 有了这两个符号，弹性增量法则的局部表达式可以表示为： 11 1 11 2 22 33 22 1 22 2 11 33 33 1 22 2 11 22 ( ) ( ) ( ) 2 e e e e e e e e e e e e e e e e e e e G （3-49） 式中，1K4 / 3G ，2K2 / 3G ，上标 e 代表“弹性部分”。 将摩尔—库仑破坏准则定义为 f ，从点s 0 A 到点 B 的局部破坏包络线的定义为： 22tan s j j f c （3-50） 将拉应力破坏准则定义为f ，从点 B 到点 C 的局部破坏包络线的定义为： 0 22 t t j f （3-51） 式中，j，c 和j ij分别为弱面的摩擦角、粘聚力和抗拉强度。 注意：对于摩擦角不为零的弱面，抗拉强度的最大值定义如下： ,max tan j t j j c （3-52） 剪切和张拉势函数g 和s g 对应于不相关联的流动法则，剪胀角t j 对应于相关联的流动 法则，它们分别表示为：
22 22 tan s j t g g （3-53） 破坏准则边界附近的流动法则用 FLAC 中摩尔—库仑模型中所述的方法定义。这里，由 函数h(22 , ) 0表示(22 , ) 平面中f 和s 0 f t 0曲线的对角线，此函数形式为： 22 ( ) p p t j j j h （3-54） 其中 p j 和 p j 是两个常量，定义如下： 2 tan 1 tan tan p t j j j j p j j j c （3-55） 弹性假设和破坏准则不一致，分别在(22 , ) 平面中位于 1 区或 2 区（分别对应于 h=0 区 域内或+区域）。如果位于 1 区，说明平面内是剪切破坏，应用势函数g 确定的流动法则，应s 力点回归到 f 的曲线上。如果位于s 0 2 区，说明局部拉应力破坏势函数g 确定的流动法则，t 应力点回归到 f t 0的曲线上。 首先考虑平面内的剪切破坏，流动法则如下： 11 11 22 22 33 33 s p s s p s s p s s p s g e g e g e g （3-56） 其中上标 p 表示弱面上同破坏有关的塑性部分，s是待定参数，应用式（3-51）的g ，s 进行偏微分后，上式可写为： 11 22 33 0 tan 0 p p s j p p s e e e （3-57） 最终，由于局部应力修正而分解到总体坐标轴下得到的弱面上剪切破坏的总应力修正 式为： 2 2 11 12 11 22 2 2 22 12 11 22 33 33 2 2 12 12 11 22
2 (cos sin ) cos sin 2 (cos sin ) sin cos
(cos sin ) ( ) cos sin
（3-58） 这些修正被添加到应力分量ij上，ij包括总体破坏的应力修正，即这些应力修正会计算
出该时步新的应力状态。 现在考虑弱面上的拉应力破坏，这种情况下，流动法则的形式为： 11 11 22 22 22 33 33 t p t t p t t p t t p t g e g e g e g e （3-59） 这里t是待定参数，应用式（3-53）的g ，进行偏微分后，上式可写为：t 11 22 33 0 0 0 p p t p p e e e （3-60） 应用如上推理，可以得到： 11 11 2 22 22 33 33 2 22 1 1 ( ) N t N t N t N t t f （3-61） 分解到整体坐标系轴下后，应力修正变成： 2 2 2 11 22 1 2 2 1 22 22 2 2 33 22 1 2 2 12 22 1 ( ) cos sin ( ) sin cos ( ) ( ) 1 cos sin t t t t （3-62） 在大应变模式下，为考虑由于刚体转动和变形引起的转动，对软弱面的方位角 进行调整， 修正值 由一个区域内所用角度的平均值得出，其表达式为： 12 e （3-63） 其中：
2 2
12 11 12 12
. .
2,1 1,2
( ) sin cos (cos sin ) 1 ( ) 2 e e e e u u （3-64） 以 以弧度形式表示。 3.3.3 双线性应变硬化—软化的节理化模型 双线性应变硬化—软化的节理化模型是本章所描述的节理化模型的推广。在双线性模型 中，岩土介质和节理的破坏包络线是由两个摩尔—库仑准则和依照规定法则将硬化或软化的拉 应力路径合成的。非相关联流动法则适用于剪切塑性流动，相关联流动法则适用于张拉塑性流 动。 对于岩土介质和节理的软化特性，可用四个独立的硬化参数描述，其中两个参数用于介 质基体，另外两个用于节理，它们分别计算塑性剪切应变和塑性张拉应变。这种数值模型中， 首先判断总体破坏并应用相应的塑性修正，然后分析弱面上的应力破坏并对应力进行更新。如 果发生塑性流动，且根据用户输入表格的参数调整岩土介质和节理的粘聚力、摩擦角、剪胀角 和抗拉强度参数值，那么硬化参数是增加的。 1．岩土介质破坏依据和流动法则 破坏包络线由两个摩尔—库仑破坏准则定义：f 对应于2s 0 AB 段，f 对应于1s 0 BC 段， 拉应力破坏准则f t 0对应于 CD 段。剪切破坏依据通常表示为f ，它由两个特征值粘聚s 0 力 C 和摩擦角来描述，其中c 和2 2对应于 AB 段，c 和2 2对应于 BC 段。张拉破坏判据由 拉伸强度t（正值）表示，可以得到： 1 3 2 s f N c N （3-65） 3 t t f （3-66） 上式中： 1 sin 1 sin N （3-67） 3 对应于f 和2s 0 f 的交集，由下式给出：1s 0 2 1 2 1 2 1 3 2 2 I c N c N N N （3-68） 注意剪切包络线上 BC 表现出来的拉应力特征，对于摩擦角10的材料，拉应力强度的 最大值为： 1 max 1 tan t c （3-69） 模型的公式表达中，应力的弹性假设首先由总应变增量来判断。如果相应的应力点在 （ 1I, 3I ），在破坏面的外侧，则判断发生塑性屈服。在这种情况下，应力修正必须应用到弹 性假设中。为恢复到与 AB、BC 或 CD 上的应力点位置相关的 f 、2s 0 f 和1s 0 f t 0的状态， 允许发生塑性流动以确定修正应力。 通常假定总应变增量由弹性应变增量和塑性应变增量两部分组成，塑性屈服的流动法
则是： p i i g e i=1，3 （3-70） 对应于非相关联流动法则的剪切屈服势函数g 可表示为：s 1 3 s g N （3-71） 式中，破坏的 AB 段的剪胀角等于1： 1 sin 1 sin N （3-72） 对应于相关联流动法则的张拉屈服势函数g 可表示为：t 3 t g （3-73） 对应于剪切破坏的塑性应变增量的形式可表示为： 1 2 3 0 ps s ps ps e e e N （3-74） 对于剪切破坏的应力修正： 1 1 2 2 2 ) 3 1 2 ( ) (1 ( ) s s s N N N （3-75） 式中： 1 4 3 K G 2 2 3 K G （3-76） 1 3 1 2 1 2 2 ( ) ( ) I I s N c N N N N 同理，张拉破坏的塑性应变增量可表示为： 1 2 3 0 0 pt pt pt t e e e （3-77） 对于张拉破坏的应力修正： 1 2 2 2 3 1 t t t （3-78） 其中： 3 1 I t t （3-79）
2．弱面的破坏判据和流动法则 对岩土介质中塑性流动校正的应力，可以分解为平行和垂直于弱面的分量，并且可以 测试节理面的破坏。破坏依据可以用弱面上的拉应力切向分量12 和法向分量22 的大 小表示。 剪切判据通常表示为 f ，它由两个变量粘聚力和摩擦角s 0 c 、j j描述，其值等于c 、i2 1 i 。张拉破坏判据由拉应力强度 t i （正值）表示，可得到： 22tan s j j f c （3-80） 22 t t j f （3-81） 注意：对于摩擦角j1不为零的弱面，抗拉强度的最大值为： 1 ,max tan j t j ji c （3-82） 塑性屈服的流动法则为： 22 22 ps ps g e g （3-83） 是同 有关的应变变量。 弱面上对应于相关联流动法则的张拉屈势函数g 为：t 22tan s j g （3-84） 式中，j为剪胀角，在破坏的 AB 段的值为j2，在 BC 段的值为j1。 弱面上对应于相关联流动法则的张拉屈服势函数g 为：t 22 t g （3-85） 剪切破坏的局部塑性应变增量为： 22 tan ps s j ps s e （3-86） 对剪切破坏的应力修正为： 11 2 22 1 33 2 tan tan tan 2 s j s j s j s G （3-87） 其中： 22 1 tan 2 tan tan O O j j s j j c G （3-88） 上标 O 表示检查判断弱面上破坏前获得的值。
弱面上的局部剪切应力分量的塑性修正由下式得出： 12 12 O O （3-89） 张拉破坏的局部塑性应变增量为： 22 0 pt t pt e （3-90） 对张拉破坏的应力修正为： 11 2 22 1 33 2 （3-91） 其中： 22 1 O t j t （3-92） 3．硬化参数 在双线性应变硬化—软化的节理化模型中，塑性一旦发生，部分或全域岩土基质和节理 的屈服性质参数（粘聚力、摩擦角、剪胀角和抗拉强度）的值将在一定范围内依照其与硬化参 数的分段线性准则进行自动调整。对于每个软化参数，在 PROPERTY 命令中必须确定表格号 （如果没有确定表格号，则认为参数为常量）。对应的表格中的数据包含成对的关于硬化/软化 参数和性质参数的值，二者中间用了线性变化假设。当硬化参数的值大于其在表格中的最后一 个值时，性质参数使用其在表格中的最后一个值。 这个模型中用到了 4 个独立的硬化变量： （1）s：用来计算岩土材料的塑性剪切应变，还可以用于更新其粘聚力、摩擦角和剪胀 角。 （2） t ：用来计算岩土材料的塑性体积拉应变，以及更新其抗拉强度。 （3） s j ：用来计算节理的塑性剪切应变，控制粘聚力、摩擦角和剪胀角的更新。 （4） t j ：用来计算节理的塑性体积拉应变，控制节理抗拉强度的更新。 对计算域而言，参数都定义成塑性应变的增量和的形式。计算域的硬化增量由域内的所 有三角形的硬化增量的平均值计算得出。 三角形的剪切硬化增量是塑性主应变偏量张量的第二不变量的平方根，对岩土材料： 2 2 1 3 1 ( ) ( ) 2 s ps ps ps ps m m e e e e （3-93） 其中 ps m e 是塑性剪切的体积应变增量， 1( 1 3 ) 3 ps ps ps m e e e 。 塑性应变增量由式（3-74）给出，s由式（3-76）给出。对于节理，方程式为： 2 2 22 12 1 ( ) ( ) 3 s ps ps j e e （3-94） 其中塑性应变增量由式（3-82）给出，s由式（3-84）给出。
四面体的张拉硬化增量是塑性张拉体积应变增量，其相应于岩土材料基质的表达式为： 3 t ept （3-95） 其中塑性应变增量由式（3-73）给出，t由式（3-75）给出。 相应于节理的表达式为： 22 t pt j e （3-96） 其中塑性应变增量由式（3-86）给出，t由式（3-87）给出。
3.4 双屈服、剑桥和霍克布朗模型
3.4.1 双屈服模型 双屈服模型考虑了各向同性压力造成的永久体积变化，在 FLAC 的应变硬化—软化模型 的剪切和张拉破坏包络线的基础上，又包括一个体积屈服面（又称“帽盖”）。为简单起见，帽 盖面由“帽盖压力”p 来定义，它与剪切应力无关，由垂直于剪应力和平均应力之间关系c 0 曲线的直线组成。帽盖压力的硬化特性在塑性体积应变下是可用的，并且遵循由用户提供的表 格所表述的分段线性规律。切向体积模量和剪切模量依照定义系数 R 的法则，以塑性体积应 变的形式出现，这里假定 R 是表示弹性体积模量和塑性体积模量之比的常量。 除了下述与应变软化模型有关的材料参数外，只需要两个额外的材料参数和一张表格。 （1）p 的初始值，它对应于材料在过去所经受的最大平均压力。 c （2）R 值，它控制体积卸载（在土力学术语中，称为“滞涨”）阶段应力—应变曲线的斜 率，其值大于 1。 （3）“硬化曲线”的表格表示法，它将帽盖压力p 同塑性体积应变c epw联系起来。 因此，任何实验室测定的硬化特性都可以利用双参数模型进行模拟。 1．增量弹性法则 此模型在 FLAC 中的实现用到了主应力1，2，3和平面外应力zz，主应力和主方向 由应力张量分量确定，并按顺序排列（压应力为负值），所以有： 1 2 3 （3-97） 相应的主应变增量 、e1 e2、e3可分解如下： e p i i i e e e i=1，3 （3-98） 上标 e 和 p 分别指弹性和塑性部分，且塑性分量仅在塑性流动阶段不为零。假定塑性应变 由剪切屈服和体积屈服产生，并由其相加，可得： p ps pt pv i i i i e e s e （3-99） 上标 ps，pt 和 pv 分别指弹性剪切应变、塑性拉应变和塑性体应变。根据约定，符号 e 与 它的塑性部分ep和弹性部分ee，表示负的值应变增量（ 1 2 3 e e e ）。符号epv 表示塑性体应变（e1pv e2pv e3pv）的负值。 主应力和主应变形式的 Hooke 定律增量表达式为：1 1 1 2 2 3 2 1 2 2 1 3 3 1 3 2 1 2 ( ) ( ) ( ) e e e e e e e e e e e e e e e e e e （3-100） 式中，1Kc4Gc/ 3，2Kc2Gc 3，K 和c G 分别为基于如下考虑的瞬时切向体积c 模量和瞬时切向剪切模量。 在等向压缩试验中，随着压力p 的增加，材料变得更加密实，其塑性刚度（ dc pc/ depv） 通常会增加，既然材料颗粒由于外力作用变得更加紧密，似乎弹性刚度同时增加也是合理的。 因此，在一般加载条件下，双屈服模型应用了一条简单的准则，即增加的弹性刚度K ，是由c 常量 R 乘以当前增加的塑性刚度。用户提供的体积模量 K 和剪切模量 G 的值，认为是K 和c Gc 的上限，假定之比Kc/G 为常量，使用增量符号，此准则定义为以下关系式： c : min( , ) c c pv c c c c p K K K K e K G G K （3-101） 式中，系数 R 已给定，pc/epv是p 当前值位置的曲线斜率。 c 2．屈服函数和势函数 剪切屈服函数 f 和张拉屈服函数s f 分别为：t 1 3 2 s f N c N （3-102） 3 t t f （3-103） 式中，N (1 sin ) /(1 sin ) ；为摩擦角；c 为粘聚力；t抗拉强度。 体积屈服函数 f 定义为：v 1 2 3 1 ( ) 3 v c f p （3-104） 式中，p 为帽盖压力。 c 剪切屈服势函数g 对应于非相关联的流动法则，张拉势函数s g 和体积势函数t g 对应于相v 关联的流动法则，表达式为： 1 3 3 1 2 3 1 ( ) 3 s t v g N g g （3-105） 式中，N (1 sin )(1 sin ) ，为剪账角。 3．硬化软化参数 依靠查找表格确定的硬化准则，剪切和体积屈服面可以硬化或软化，张拉屈服面可以软 化。查找表中的条目记录了一些累积的塑性应变的硬化参数。在剪切和拉伸中，硬化参数的增 量形式为：
1 2 2 1 3 1 1 1 ( ) ( ) ( ) 2 2 2 ps ps ps ps ps ps m m m e e e e e e （3-106） 其中eps ( e1ps e3ps) / 3，ejps（j 1 3，）和e3ps分别是主方向上的塑性剪切和拉伸应变 增量。 4．塑性修正 基于弹性假设，由上标 I 表示原来的应力 O ij 加上由总应变增量计算所得弹性应力增量以 后得到的应力估算，在主轴方向有： 1 1 1 1 2 2 3 2 2 1 2 2 1 3 3 3 1 3 2 1 2 ( ) ( ) ( ) I O I O I O e e e e e e E e e （3-107） 在 FLAC 的 数 值 实 现 过 程 中 ， 如 果 fs( 1I, 3I)0， 认 为 发 生 剪 切 破 坏 ， 如 果 1 2 3 ( , , ) 0 v I I I f ，认为发生体积屈服，如果认为发生张拉屈服，相应的塑性修正分别由下列 方法计算： 首先考虑没有发现张拉破坏但是剪切和体积屈服状态都超出的情况，由式（3-97）和式 （3-98）总应变增量可以表示为： e ps pv i i i i e e e e （i=1，3） （3-108） 剪切和体积屈服的流动法则为： s ps s i i v pv v i i g e g e （3-109） 其中 i=1，3，经求导后，这些表达式变成： 1 2 3 0 ps s ps ps s e e e N （3-110） 及 1 2 3 1 3 1 3 1 3 pv v pv v pv v e e e （3-111） 得到：
1 1 2 2 3 3 / 3 / 3 / 3 e s v e v e s v e e e e e e N （3-112） 有了这些弹性应变增量的表达式，可以得到 Hooke 增量表达式： 1 1 1 2 2 2 2 3 3 2 1 ( ) (1 ) ( ) N I s v N I s v N I s v N K N K N K （3-113） 式中： I i （i=1，3）为初始试验应力，iN iO i（i 1 3，）是当前时步新的主应力。 为了确定系数s和v，做出如下假定：如果剪切和体积屈服同时发生，两个屈服面上都 产生新的应力，此时 fs( 1I, 3I)0和fv( 1I, 2I, 3I)0。可以得到s： 1 2 2 1 (1 ) (1 )(1 ) sI vI s f f N N N N N K N N （3-114） 因此： (1 ) vI v f s N K （3-115） 如果单元仅仅发生剪切屈服，v0，则： 1 2 2 1 sI s f N N N N （3-116） 如果单元仅仅发生体积屈服，v0，则： vI v f K （3-117） 现在考虑由 ft(3I)0确定的张拉破坏的情况，如果体积破坏没有出现，可采用摩尔—库 仑 模 型 中 描 述 过 的 相 同 方 法 和 应 力 修 正 。 如 果 除 张 拉 破 坏 外 还 发 生 了 体 积 破 坏 ， 则 1 3 ( , ) 0 s I I f ≤ 或fs( 1I, 3I)0。假设三种屈服状态同时出现，则假设认为剪切屈服、体积 屈服和张拉屈服的塑性部分相加，即： 1 3 e ps pv pt i i i i i e e e e e i （ ，） （3-118） 张拉屈服的流动法则如下： 1 2 3 0 0 pt pt pt t e e e （3-119） 用上面同样的推理，可以得到： 1 1 1 2 2 2 2 2 2 3 3 2 1 1 ( ) (1 ) ( ) N I s v t N I s v t N I s v t N K N K N K （3-120）
系数s，v和t通过解三个等式 1 3 ( , ) 0 s N N f 、fv(1N,2N,3N)0和 ft(3N) 来确0 定，这里给出如下： 2 1 2 1 2 1 (1 2 ) 3 2 3(1 2 ) 6 3 (1 2 ) 2 3(1 ) (1 2 ) tI vI sI s tI vI sI vI v tI vI sI t f N f f N f f f f K f N N N N f N f （3-121） 代入式（3-120），屈服为： 1 2 3 2 3 (1 ) 2 N t N t c N t N c N P N c N （3-122） 如果仅发生张拉和体积屈服，则在式（3-120）中s 。 0 常量v和t由条件 1 2 3 ( , , ) 0 v N N N f 和 ft(3N) 确定。经过处理，可得到： 0 1 1 1 ( ) vI tI v vI tI t f Kf K K f f K （3-123） 代入式（3-120），得： 1 1 2 2 3 3 2 3 2 vI tI N I vI tI N I N t f f f f （3-124） 3.4.2 修正的剑桥粘土模型 修正的剑桥粘土模型是增量硬化—软化弹塑性模型，其参数包括一个特殊的非线性弹性 部分和一个由体积塑性应变（“密度”确定）确定的硬化—软化特性。破坏包络线在形状上是 相似的，并且同绕主应力空间的平均应力轴转动的椭圆体对应。在这个模型中，同剪切采用了 流动法则相关联法则，没有提供相应于张拉平均应力的阻力。见 Roscoe and Burland(1968) and Wood(1990)对于修正的剑桥模型的详细描述，为方便起见，后文省略掉限定词“修正”；注意 所有的模型都可以有效应力的形式表达。特别地，本节提到的所有压力均为有效应力。 1．屈服函数 相应于特定的固结压力值的屈服函数的形式为： 2 ( ) c f q Mp pp （3-125）
这里 M 是材料常数。屈服状态 f 0由平面内以p 为横轴，以c Mpc为纵轴的椭圆表示。 注意椭圆过原点，因此，此模型的材料并不支持四周受拉应力的情况。破坏准则由主应力空间 内绕平均应力轴转动的椭圆表示。势函数对应于相关联的流动法则，由此可得： 2 2 ( ) c gq M p pp （3-126） 2．硬化—软化准则 屈服函数的形状取决于固结压力p 的值，此压力是塑性体积变化的函数，且随比容而变c 化。 对应于新的值和 p 固结压力p 可以通过c ( ,ln ) p 平面内固结线和膨胀线的交点得出，这 里给出其表达式： ( )( ) 1 c p p e （3-127） 其中： 1 ln p p （3-128） 3．初始应力状态 FLAC 中的剑桥粘土模型仅可用于平均有效应力为压应力的材料。特别地，材料的初始状 态（只是在应用剑桥粘土模型之前）必须同此要求相协调。初始状态可以通过 INITIAL 指定， 或者是其他一应用的本构模型的运行结果。任何情况下，初始有效压力p 在整个介质中必须0 为正值。 4．确定输入参数 超固结比 R——先期固结压力和初始固结压力之比，即： 0 0 c p R p （3-129） 超固结比在描述剑桥粘土材料的行为时很有用。 摩擦常数 M——在临界状态线上 /q p 的比值，所以可以用一系列三轴试验（排水或不排cr 水，并量测孔隙水压力）来确定此常数。这些试验应在大应变条件下进行，确保最终的p 和cr q 值接近临界状态，qpcr（通过回归分析所得）最好吻合线的斜率，即为参数 M。M 同摩 尔—库仑屈服函数的有效应力摩擦角 有关。然而，因为剑桥粘土临界状态线与中主应力2相 关，而摩尔—库仑模型与中主应力2不相关，在屈服时对不用的2值，M 和 的关系也不同 （这种状态相似于摩尔—库仑和德鲁克—普拉格屈服函数间的关系）。对于三轴压缩试验： 6 sin 3 sin M （3-130） 对于三轴拉伸试验，则： 6 sin 3 sin M （3-131） 正常固结和膨胀曲线（和）的斜率——理想状态下这两个参数应由带几个卸载偏移的 各向同性加载三轴（q ）试验确定。图上正常压缩曲线的斜率即为参数0 ，同一图上的卸 载偏移曲线的斜率即为参数。 这两个参数也可以从固结实验的确定假定获得，让V和H为一个固结试验的垂直和水
平应力。在大多数固结仪器中，是不可能测量水平应力H的，所以平均应力p(V 2H) / 3 是未知的。但是，实验数据表明，正常固结条件下水平和垂直有效应力的比值为一个常数。既 然pV(1 2 K0) / 3沿正常压缩曲线，ln p曲线的斜率等于eln的斜率，这里e1 是孔隙比。 压缩指数C 由曲线c log (10 )的斜率计算而得，所以参数为： / ln(10) c C （3-132） 实验数据表明，在固结实验中，沿着膨胀线，不是常数，所以基于膨胀系数的一个估计 值，仅仅是一个近似值。 /1 (10) s C h 在ln p图上定出正常固结线——为了在ln p图上定出正常固结线，必须确定此曲线 上 的 点(, lnp1) 。 确 定 这 点 的 方 法 是 各 向 同 性 的 三 轴 试 验 ， 或 基 于 不 排 水 剪 切 强 度 （BrittoGunn，1987）确定该点。 正常固结线的表达式为： 1 ln p p （3-133） 在pp1的临界状态线上的比容，由下式给出： ( ) ln(2) （3-134） 在土体中，不排水剪切强度 Cu通过下式与比容积cr唯一相关： 1exp 2 cr u Mp C （3-135） 因此，对于一个给定p ，如果对于一个比容积1 cr的不排水剪切强度，连同参数 M，和 都已知， 的值以及r可以计算出来。 先期固结压力p ——先期固结压力确定屈服面的初始大小，表达式： c0 2 2 0 [ ( c p)] q M p p （3-136） 如果试样服从等向加压路径，P 的值达到最大，超过平均有效应力。如果试样服从其他C0 的非等向路径，p 由式（3-136）用 p 和 q 的先期最大值计算。 c0 最大竖向有效应力可以应用卡萨兰德的方法（BrittoGunn，1987）由固结仪试验计算得到。 最大水平有效应力必须由一些假设得出，通常的假设是 Jaky 关系： max max 1 sin h nc v K （3-137） 式中，Knc为正常固结土的静止最大水平应力与最大竖向应力之比。 如某一土体，有效摩擦角为 20 ，已经经受了最大竖向有效应力，然后应用 Jaky 关系式： 1 sin 20 0.658 nc K （3-138） 最大水平应力为max 0.658MPa。 （3-139） p 和 q 的最大值为： max max max 2 0.772MPa 3 V H p （3-140） 将这两个值代入屈服关系式，可得到初始固结压力：
2 max 0 max 2 max 0.961MPa c q p p M p （3-141） 初始比容值和瞬时体积模量——假定初始压力p ，初始比容0 0一定与参数， ， ， p pc0 的选择相一致。初始比容0由膨胀曲线和正常固结线相交的ppc0之点所对应p 的比容的值0 计算得出，表达式为： 0 0 0 1 1 ln pc ln pc p p （3-142） 体积模量的初始值（bulk_current）可用下式算出： 0p0 K （3-143） 在 FLAC 中，0和 K 的默认值在命令执行的第一步由式（3-142）和式（3-143）计算得出。 弹性参数和大最大值——在剑桥粘土模型中，当前体积变量（bulk_current）的值是作为 比容和平均应力的函数变化：
在 FLAC 中，Kmax(bulk 和) G shear 的输入值用于进行质量缩放计算（mass scaling ( ) calculation），以保证数值计算的稳定性。一旦执行 STEP 命令，这种计算就开始进行。像两次 连续的 STEP 命令所计算出的值一样，选择输入的数值须能给出 (K4 / 3 )G 总和上限。然而， 这些值不应设置得太大，否则模型将收敛得较慢，这些值的选取应该基于问题的应力水平。 G 或者——FLAC 中的修正剑桥模型允许用户指定一个常剪切模量或一个恒定的泊松 比。如果没有指定泊松比，则假定常剪切模量等于输入值。然后，泊松比作为比容和平均应力 的函数而变化： 3 2 6 2 p G p G （3-144） 如果指定了一个非零的泊松比，剪切模量将随体积模量同比变化，以使泊松比保持常数： 3 (1 2 ) 2(1 2 ) p G （3-145） 3.4.3 霍克－布朗模型 霍克－布朗模型准则是一个描述完整岩块或岩体引起破坏的经验关系。它在那些使用基 线平衡法的设计方法中得到了十分成功的应用，但是还没有直接用于数值解法中。在特定应力 水平上，选择了与摩尔－库仑模型中应用的等效的摩擦角和粘聚力，与非线性的霍克－布朗模 型强度包线相匹配的。数值解法需要完整的本构模型，它将应力和应变以一种全面的方式联系 起来，此外，还需要一种破坏（屈服）准则，一种“流动法则”以确定破坏时应变速率分量之 间的联系。已经有几种试图从霍克－布朗准则推出一个完整的本构模型的尝试，比如 Pan and Hudson(1988)、carter et al.(1993)以及 Shah(1992)。这些方程假定流动法则和破坏准则之间具有 某种固定的关系，并且流动法则是各向同性的，但是霍克—布朗准则却不是各向同性的。这里
