Mohr-Coulomb准则的弹塑性积分方法

张 谭，郑 宏，林 姗

(北京工业大学城市与工程安全减灾教育部重点实验室，北京100124)

为了处理实际活跃屈服面未知的问题，已有文献中所采用的方法同样可分为3类: 1)应力指示因子法。de Borst[12]提出用应力指示因子的概念，并根据指示因子正负号来确定哪组屈服面最终被激活；Pankaj等[13]使用并发展了该方法，由于不同的屈服函数对应的应力指示因子不同，因此该方法的通用性受到了限制。2)几何法。Clausen等[14]利用理想Mohr-Coulomb屈服准则在主应力空间中的几何性质，推导了适用于应力返回几何方法，该思想被用于其他塑性模型[15]，但是该方法同样只适用于特定的屈服准则，无法直接推广到其他的屈服函数，需要重新推导。3)互补法。互补理论已经被广泛应用于弹塑性计算[16-17]，Simo等[18]基于Kuhn-Tucker互补条件，采用凸规划的数学思想并结合最近点投影算法，只进行迭代即可确定最终的活跃面，该方法被不断地进行改进[19-20]，并应用于更复杂的弹塑性模型[21]。互补法是目前应用最广泛的一种方法，原因是该方法提供了一个通用的、稳健的框架，排除了对屈服面形状的依赖。

本文中基于Koiter法则和Kuhn-Tucker互补条件，提出一种适用于Mohr-Coulomb准则的本构积分方法。该方法采用Koiter法则决定塑性应变的方向和大小，Kuhn-Tucker条件决定可能的活跃屈服面，然后将Kuhn-Tucker互补方程作为一类特殊的变分不等式使用投影收缩算法求解。在主应力空间中计算应力分量的偏导，然后在一般应力空间中进行应力返回，避免角点处的奇异性，同时不需要进行谱分解和坐标变换。该算法同样能适应各种广泛应用的非光滑多面塑性模型，并可进一步推广到开率各向同性硬化(或软化)的情况。

1 角点问题产生的原因

经典的Mohr-Coulomb准则为

τ=c-σntanφ，

(1)

式中：τ为剪应力；c为黏聚力；σn为正应力；φ为内摩擦角。式(1)的几何解释如图1所示。

τ—剪应力； c—黏聚力；σn—正应力；φ—内摩擦角； σmax、σmin—最大、最小主应力。图1 Mohr-Coulomb准则的几何解释

Mohr-Coulomb准则可以表示为主应力的形式，即

(2)

式中σmax、σmin分别为最大、最小主应力。

整理式(2)，得

(σmax-σmin)+(σmax+σmin)sinφ=2ccosφ。

(3)

令σ1≥σ2≥σ3，其中σ1、σ2、σ3分别为第一、第二和第三主应力，则式(3)为

σ1-σ3=2ccosφ-(σ1+σ3)sinφ。

(4)

在考虑所有能够产生屈服的应力组合后，就能得到完整的屈服面。在计算中，一般不直接使用主应力表示的形式，而是采用不变量表示的Mohr-Coulomb准则，即

(5)

其中

(6)

(7)

(8)

I1=σx+σy+σz，

(9)

(10)

式中：σ为应力张量；σx、σy、σz、τxy、τxz、τyx、τyz、τzx、τzy分别为σ的分量。

使用Mohr-Coulomb准则计算塑性应变εp时，根据Koiter法则[1]，可得

(11)

式中：α为屈服面的序号；λα为对应屈服面fα的塑性乘子。需要计算Mohr-Coulomb屈服函数对应力分量的偏导数，根据链式求导法则，得

(12)

其中

sinφ(tan 3θσ-tanθσ)]，

(13)

(14)

2 算法及其实施

2.1 本构积分

弹塑性本构方程的一般形式为

σ=D(ε-εp)

(15)

式中：ε为应变；D为弹性矩阵。塑性应变εp根据式(11)中的Koiter法则确定，同时塑性乘子满足互补条件

(16)

(a)角点区域

(b)1个最终活跃的屈服面

(c)2个最终活跃的屈服面f1=0、f2=0—屈服面f1、f2的弹性边界； σ—应力； 试应力； Eσ—弹性区域； N1、N2—屈服面f1、f2在角点处法线方向； Δλ1、Δλ2—屈服面f1、f2的塑性乘子增量； C1、C2、C12—角点外不同区域。图2 2个屈服面相交形成的角点区域

对于具有多重屈服面的Mohr-Coulomb塑性模型，在计算时首先需要定义弹性空间Eσ，定义为

Eσ={σ|fα(σ)≤0,α∈{1,2,…,6}} 。

(17)

Mohr-Coulomb准则可表示为

(18)

然后进行时间离散，令σn、εn、εn,p为初始应力、应变以及塑性应变，给定应变增量Δε并应用欧拉向后差分公式，有

(19)

式(16)相应的离散形式为

(20)

在式(20)中，令Δλα=0，即得到试验应力状态为

(21)

(22)

σn+1=σn+Δσ=σn+D(Δε-Δεp)=

(23)

将式(23)代入互补方程(20)中，得

(24)

式(24)可通过投影收缩算法求解，通过迭代执行Kuhn-Tucker互补方程并不断更新约束集合即可，具体计算步骤如下。

1)计算第k步迭代应力，即

(25)

3)确定活跃的屈服面，即

(26)

4)更新状态变量，即

(27)

令k=k+1，回到步骤1)。

2.2 投影收缩算法

为了方便起见，将互补条件式(26)写成向量的形式，即

0≤Δλ⊥F(Δλ)≥0，

(28)

其中Δλ=(Δλ1,Δλ2,…,Δλ6)，F(Δλ)=(-f1(Δλ1),-f2(Δλ2),…,-f6(Δλ6))，⊥表示垂直于，对于每个α∈{1,2,…,6}，有

(29)

式(29)定义的互补方程可利用投影收缩算法进行求解，步骤如下。

1)给定如下参数：塑性乘子初值Δλ、容差T、允许迭代的最大次数Kα，以及控制参数ξ,η∈(0,1)与γ∈[1,2)。令控制参数β=1，迭代次数k=0。

2)判断是否有k

(30)

(31)

否则，说明已达最大计算次数，退出。

(32)

(33)

否则，说明已收敛，退出。

4)判断是否有r>ξ。若有，则计算

β=0.7βmin{1,1/r}，

(34)

(35)

(36)

(37)

(38)

否则，进入下一步。

5)计算

(39)

(40)

(41)

6)判断是否有r≤η。若有，则令β=1.5β；否则，不作调整，进入下一步。

7)令k=k+1，然后进入步骤2)。

2.3 主应力对应力分量的偏导数

在主应力空间中计算Mohr-Coulomb屈服函数对应力分量的偏导数，可以避免在当2个主应力相等时(即角点处)，使用不变量形式表示的屈服函数偏导数不存在的问题

主应力σi满足的特征方程为

(42)

其中

I1=σx+σy+σz，

(43)

(44)

(45)

在式(42)中，分别对9个应力分量σx,σy,…,τxy求微分，得

(46)

整理，得

(47)

即

(48)

类似地，对式(46)再次求微分，即可计算主应力对应力分量的二阶偏导。

3 算例

3.1 算例1

算例1为典型的条形基础[22]，几何尺寸及单元划分如图3所示。地基不考虑重力的影响，采用普通的四边形单元，每个单元取4个高斯积分点。材料参数如下：弹性模量E=206.9 MPa，泊松比ν=0.3，黏聚力c=69 kPa，内摩擦角φ=20°。材料符合关联流动法则的理想弹塑性Mohr-Coulomb模型。

P—施加的竖向荷载；A、B、C、D—角点。图3 条形基础尺寸及单元划分

条形基础在无重量地基上的极限承载力[22]为

(49)

把算例1中的参数代入式(48)，可得基础的极限承载力为1.023 6 MPa。图4所示为使用本文中方法得到的条形基础上角点A的荷载-位移曲线。由图可知，极限荷载为1.043 1 MPa，很接近式(48)给出的参考解，同时也很接近Zienkiewicz等[22]得到的极限荷载1.048 MPa。

图4 条形基础上角点A的荷载-位移曲线

3.2 算例2

算例2为某宽平台二级路堑式边坡[23]。一级边坡坡率为1∶0.75，坡高为5 m； 二级边坡坡率为1∶1，坡高为4 m。边坡土体自上而下分为3层，几何尺寸及单元划分如图5所示。各土层的土工参数如表1所示。

(a)几何尺寸(单位为m)

(b)单元划分S1、S2、S3—土层编号。图5 边坡几何尺寸及单元划分

对于边坡来说，安全系数在传统意义上定义为实际抗剪强度与防止破坏所需的最小抗剪强度之比。正如Duncan[24]所指出，当土体抗剪强度除以某系数F后使得边坡达到破坏边缘，该系数即为安全系数Fs。用有限单元或有限差分方法计算Fs的一个简单方法就是通过不断的降低土体的抗剪强度，直到发生破坏，即

表1 边坡土工参数

(50)

该方法早在1975年就被Zienkiewicz等[25]所使用，而后被国内外学者广泛应用于各种边坡工程的计算[26-29]，积累了大量的经验。

原则上，应当取边坡进入极限平衡状态时所对应的折减系数作为边坡的安全系数。在有限元强度折减技术中，关于边坡进入极限平衡状态的判据有很多种[30-32]。基于不同判据的安全系数差异较小，在这些判据中，塑性区贯通以及计算不收敛是使用最多的2个判据。

采用本文中的算法，利用强度折减法计算边坡的安全系数、折减强度参数的同时，根据文献[33]中的方法，调整弹性模量E和泊松比ν，得到的边坡坡顶的竖向位移与折减系数的关系，如图6所示。从图中可以看出，当折减系数F=1.35时，位移曲线存在一个明显的拐点，此时的塑性区分布如图7所示，此时塑性区刚好贯通。

图6 边坡坡顶竖向位移与折减系数的关系

若取拐点时(即塑性区贯通时)的折减系数为安全系数，则Fs=1.35；若取计算无法收敛时的折减系数为安全系数，则Fs=1.38。二者相差较小，相对误差为2.17%。

图7 边坡塑性区分布

4 结论

在互补理论的基础上，建立了弹塑性互补方程，利用Koiter法则把多曲面的角点问题统一归结为互补方程组，然后利用投影收缩算法进行求解，理论和算例表明：

1)多屈服面的互补方程的投影收缩算法避免了角点光滑化带来的误差以及不易收敛的问题。

2)在一般应力空间中进行应力返回操作，在主应力空间中计算对应力分量的偏导数，既避免了角点处的数值奇异性，又省略了主应力空间应力返回的应力变换。

3)算例得到的基础极限承载力与参考解非常接近，结合强度折减法并根据不同判据得到了边坡的安全系数及破坏机理，表明本文中提出的算法是可靠、有效的。