Mohr-Coulomb準則的彈塑性積分方法

張 譚，鄭 宏，林 姍

(北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點實驗室，北京100124)

為了處理實際活躍屈服面未知的問題，已有文獻中所采用的方法同樣可分為3類: 1)應(yīng)力指示因子法。de Borst[12]提出用應(yīng)力指示因子的概念，并根據(jù)指示因子正負號來確定哪組屈服面最終被激活；Pankaj等[13]使用并發(fā)展了該方法，由于不同的屈服函數(shù)對應(yīng)的應(yīng)力指示因子不同，因此該方法的通用性受到了限制。2)幾何法。Clausen等[14]利用理想Mohr-Coulomb屈服準則在主應(yīng)力空間中的幾何性質(zhì)，推導(dǎo)了適用于應(yīng)力返回幾何方法，該思想被用于其他塑性模型[15]，但是該方法同樣只適用于特定的屈服準則，無法直接推廣到其他的屈服函數(shù)，需要重新推導(dǎo)。3)互補法。互補理論已經(jīng)被廣泛應(yīng)用于彈塑性計算[16-17]，Simo等[18]基于Kuhn-Tucker互補條件，采用凸規(guī)劃的數(shù)學(xué)思想并結(jié)合最近點投影算法，只進行迭代即可確定最終的活躍面，該方法被不斷地進行改進[19-20]，并應(yīng)用于更復(fù)雜的彈塑性模型[21]。互補法是目前應(yīng)用最廣泛的一種方法，原因是該方法提供了一個通用的、穩(wěn)健的框架，排除了對屈服面形狀的依賴。

本文中基于Koiter法則和Kuhn-Tucker互補條件，提出一種適用于Mohr-Coulomb準則的本構(gòu)積分方法。該方法采用Koiter法則決定塑性應(yīng)變的方向和大小，Kuhn-Tucker條件決定可能的活躍屈服面，然后將Kuhn-Tucker互補方程作為一類特殊的變分不等式使用投影收縮算法求解。在主應(yīng)力空間中計算應(yīng)力分量的偏導(dǎo)，然后在一般應(yīng)力空間中進行應(yīng)力返回，避免角點處的奇異性，同時不需要進行譜分解和坐標變換。該算法同樣能適應(yīng)各種廣泛應(yīng)用的非光滑多面塑性模型，并可進一步推廣到開率各向同性硬化(或軟化)的情況。

1 角點問題產(chǎn)生的原因

經(jīng)典的Mohr-Coulomb準則為

τ=c-σntanφ，

(1)

式中：τ為剪應(yīng)力；c為黏聚力；σn為正應(yīng)力；φ為內(nèi)摩擦角。式(1)的幾何解釋如圖1所示。

τ—剪應(yīng)力； c—黏聚力；σn—正應(yīng)力；φ—內(nèi)摩擦角； σmax、σmin—最大、最小主應(yīng)力。圖1 Mohr-Coulomb準則的幾何解釋

Mohr-Coulomb準則可以表示為主應(yīng)力的形式，即

(2)

式中σmax、σmin分別為最大、最小主應(yīng)力。

整理式(2)，得

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

(3)

令σ1≥σ2≥σ3，其中σ1、σ2、σ3分別為第一、第二和第三主應(yīng)力，則式(3)為

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

(4)

在考慮所有能夠產(chǎn)生屈服的應(yīng)力組合后，就能得到完整的屈服面。在計算中，一般不直接使用主應(yīng)力表示的形式，而是采用不變量表示的Mohr-Coulomb準則，即

(5)

其中

(6)

(7)

(8)

I1=σx+σy+σz，

(9)

(10)

式中：σ為應(yīng)力張量；σx、σy、σz、τxy、τxz、τyx、τyz、τzx、τzy分別為σ的分量。

使用Mohr-Coulomb準則計算塑性應(yīng)變εp時，根據(jù)Koiter法則[1]，可得

(11)

式中：α為屈服面的序號；λα為對應(yīng)屈服面fα的塑性乘子。需要計算Mohr-Coulomb屈服函數(shù)對應(yīng)力分量的偏導(dǎo)數(shù)，根據(jù)鏈式求導(dǎo)法則，得

(12)

其中

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

(13)

(14)

2 算法及其實施

2.1 本構(gòu)積分

彈塑性本構(gòu)方程的一般形式為

σ=D(ε-εp)

(15)

式中：ε為應(yīng)變；D為彈性矩陣。塑性應(yīng)變εp根據(jù)式(11)中的Koiter法則確定，同時塑性乘子滿足互補條件

(16)

(a)角點區(qū)域

(b)1個最終活躍的屈服面

(c)2個最終活躍的屈服面f1=0、f2=0—屈服面f1、f2的彈性邊界； σ—應(yīng)力； 試應(yīng)力； Eσ—彈性區(qū)域； N1、N2—屈服面f1、f2在角點處法線方向； Δλ1、Δλ2—屈服面f1、f2的塑性乘子增量； C1、C2、C12—角點外不同區(qū)域。圖2 2個屈服面相交形成的角點區(qū)域

對于具有多重屈服面的Mohr-Coulomb塑性模型，在計算時首先需要定義彈性空間Eσ，定義為

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

(17)

Mohr-Coulomb準則可表示為

(18)

然后進行時間離散，令σn、εn、εn,p為初始應(yīng)力、應(yīng)變以及塑性應(yīng)變，給定應(yīng)變增量Δε并應(yīng)用歐拉向后差分公式，有

(19)

式(16)相應(yīng)的離散形式為

(20)

在式(20)中，令Δλα=0，即得到試驗應(yīng)力狀態(tài)為

(21)

(22)

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

(23)

將式(23)代入互補方程(20)中，得

(24)

式(24)可通過投影收縮算法求解，通過迭代執(zhí)行Kuhn-Tucker互補方程并不斷更新約束集合即可，具體計算步驟如下。

1)計算第k步迭代應(yīng)力，即

(25)

3)確定活躍的屈服面，即

(26)

4)更新狀態(tài)變量，即

(27)

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

2.2 投影收縮算法

為了方便起見，將互補條件式(26)寫成向量的形式，即

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

(28)

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

(29)

式(29)定義的互補方程可利用投影收縮算法進行求解，步驟如下。

1)給定如下參數(shù)：塑性乘子初值Δλ、容差T、允許迭代的最大次數(shù)Kα，以及控制參數(shù)ξ,η∈(0,1)與γ∈[1,2)。令控制參數(shù)β=1，迭代次數(shù)k=0。

2)判斷是否有k

(30)

(31)

否則，說明已達最大計算次數(shù)，退出。

(32)

(33)

否則，說明已收斂，退出。

4)判斷是否有r>ξ。若有，則計算

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

(34)

(35)

(36)

(37)

(38)

否則，進入下一步。

5)計算

(39)

(40)

(41)

6)判斷是否有r≤η。若有，則令β=1.5β；否則，不作調(diào)整，進入下一步。

7)令k=k+1，然后進入步驟2)。

2.3 主應(yīng)力對應(yīng)力分量的偏導(dǎo)數(shù)

在主應(yīng)力空間中計算Mohr-Coulomb屈服函數(shù)對應(yīng)力分量的偏導(dǎo)數(shù)，可以避免在當2個主應(yīng)力相等時(即角點處)，使用不變量形式表示的屈服函數(shù)偏導(dǎo)數(shù)不存在的問題

主應(yīng)力σi滿足的特征方程為

(42)

其中

I1=σx+σy+σz，

(43)

(44)

(45)

在式(42)中，分別對9個應(yīng)力分量σx,σy,…,τxy求微分，得

(46)

整理，得

(47)

即

(48)

類似地，對式(46)再次求微分，即可計算主應(yīng)力對應(yīng)力分量的二階偏導(dǎo)。

3 算例

3.1 算例1

算例1為典型的條形基礎(chǔ)[22]，幾何尺寸及單元劃分如圖3所示。地基不考慮重力的影響，采用普通的四邊形單元，每個單元取4個高斯積分點。材料參數(shù)如下：彈性模量E=206.9 MPa，泊松比ν=0.3，黏聚力c=69 kPa，內(nèi)摩擦角φ=20°。材料符合關(guān)聯(lián)流動法則的理想彈塑性Mohr-Coulomb模型。

P—施加的豎向荷載；A、B、C、D—角點。圖3 條形基礎(chǔ)尺寸及單元劃分

條形基礎(chǔ)在無重量地基上的極限承載力[22]為

(49)

把算例1中的參數(shù)代入式(48)，可得基礎(chǔ)的極限承載力為1.023 6 MPa。圖4所示為使用本文中方法得到的條形基礎(chǔ)上角點A的荷載-位移曲線。由圖可知，極限荷載為1.043 1 MPa，很接近式(48)給出的參考解，同時也很接近Zienkiewicz等[22]得到的極限荷載1.048 MPa。

圖4 條形基礎(chǔ)上角點A的荷載-位移曲線

3.2 算例2

算例2為某寬平臺二級路塹式邊坡[23]。一級邊坡坡率為1∶0.75，坡高為5 m； 二級邊坡坡率為1∶1，坡高為4 m。邊坡土體自上而下分為3層，幾何尺寸及單元劃分如圖5所示。各土層的土工參數(shù)如表1所示。

(a)幾何尺寸(單位為m)

(b)單元劃分S1、S2、S3—土層編號。圖5 邊坡幾何尺寸及單元劃分

對于邊坡來說，安全系數(shù)在傳統(tǒng)意義上定義為實際抗剪強度與防止破壞所需的最小抗剪強度之比。正如Duncan[24]所指出，當土體抗剪強度除以某系數(shù)F后使得邊坡達到破壞邊緣，該系數(shù)即為安全系數(shù)Fs。用有限單元或有限差分方法計算Fs的一個簡單方法就是通過不斷的降低土體的抗剪強度，直到發(fā)生破壞，即

表1 邊坡土工參數(shù)

(50)

該方法早在1975年就被Zienkiewicz等[25]所使用，而后被國內(nèi)外學(xué)者廣泛應(yīng)用于各種邊坡工程的計算[26-29]，積累了大量的經(jīng)驗。

原則上，應(yīng)當取邊坡進入極限平衡狀態(tài)時所對應(yīng)的折減系數(shù)作為邊坡的安全系數(shù)。在有限元強度折減技術(shù)中，關(guān)于邊坡進入極限平衡狀態(tài)的判據(jù)有很多種[30-32]。基于不同判據(jù)的安全系數(shù)差異較小，在這些判據(jù)中，塑性區(qū)貫通以及計算不收斂是使用最多的2個判據(jù)。

采用本文中的算法，利用強度折減法計算邊坡的安全系數(shù)、折減強度參數(shù)的同時，根據(jù)文獻[33]中的方法，調(diào)整彈性模量E和泊松比ν，得到的邊坡坡頂?shù)呢Q向位移與折減系數(shù)的關(guān)系，如圖6所示。從圖中可以看出，當折減系數(shù)F=1.35時，位移曲線存在一個明顯的拐點，此時的塑性區(qū)分布如圖7所示，此時塑性區(qū)剛好貫通。

圖6 邊坡坡頂豎向位移與折減系數(shù)的關(guān)系

若取拐點時(即塑性區(qū)貫通時)的折減系數(shù)為安全系數(shù)，則Fs=1.35；若取計算無法收斂時的折減系數(shù)為安全系數(shù)，則Fs=1.38。二者相差較小，相對誤差為2.17%。

圖7 邊坡塑性區(qū)分布

4 結(jié)論

在互補理論的基礎(chǔ)上，建立了彈塑性互補方程，利用Koiter法則把多曲面的角點問題統(tǒng)一歸結(jié)為互補方程組，然后利用投影收縮算法進行求解，理論和算例表明：

1)多屈服面的互補方程的投影收縮算法避免了角點光滑化帶來的誤差以及不易收斂的問題。

2)在一般應(yīng)力空間中進行應(yīng)力返回操作，在主應(yīng)力空間中計算對應(yīng)力分量的偏導(dǎo)數(shù)，既避免了角點處的數(shù)值奇異性，又省略了主應(yīng)力空間應(yīng)力返回的應(yīng)力變換。

3)算例得到的基礎(chǔ)極限承載力與參考解非常接近，結(jié)合強度折減法并根據(jù)不同判據(jù)得到了邊坡的安全系數(shù)及破壞機理，表明本文中提出的算法是可靠、有效的。