張 譚，鄭 宏，林 姍

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

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

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

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

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

τ=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準(zhǔn)則的幾何解釋

Mohr-Coulomb準(zhǔn)則可以表示為主應(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)力組合后，就能得到完整的屈服面。在計(jì)算中，一般不直接使用主應(yīng)力表示的形式，而是采用不變量表示的Mohr-Coulomb準(zhǔn)則，即

(5)

其中

(6)

(7)

(8)

I1=σx+σy+σz，

(9)

(10)

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

使用Mohr-Coulomb準(zhǔn)則計(jì)算塑性應(yīng)變?chǔ)舙時(shí)，根據(jù)Koiter法則[1]，可得

(11)

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

(12)

其中

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

(13)

(14)

2 算法及其實(shí)施

2.1 本構(gòu)積分

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

σ=D(ε-εp)

(15)

式中：ε為應(yīng)變；D為彈性矩陣。塑性應(yīng)變?chǔ)舙根據(jù)式(11)中的Koiter法則確定，同時(shí)塑性乘子滿(mǎn)足互補(bǔ)條件

(16)

(a)角點(diǎn)區(qū)域

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

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

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

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

(17)

Mohr-Coulomb準(zhǔn)則可表示為

(18)

然后進(jìn)行時(shí)間離散，令σn、εn、εn,p為初始應(yīng)力、應(yīng)變以及塑性應(yīng)變，給定應(yīng)變?cè)隽喀う挪?yīng)用歐拉向后差分公式，有

(19)

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

(20)

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

(21)

(22)

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

(23)

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

(24)

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

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

(25)

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

(26)

4)更新?tīng)顟B(tài)變量，即

(27)

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

2.2 投影收縮算法

為了方便起見(jiàn)，將互補(bǔ)條件式(26)寫(xiě)成向量的形式，即

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

(28)

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

(29)

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

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

2)判斷是否有k

(30)

(31)

否則，說(shuō)明已達(dá)最大計(jì)算次數(shù)，退出。

(32)

(33)

否則，說(shuō)明已收斂，退出。

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

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

(34)

(35)

(36)

(37)

(38)

否則，進(jìn)入下一步。

5)計(jì)算

(39)

(40)

(41)

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

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

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

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

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

(42)

其中

I1=σx+σy+σz，

(43)

(44)

(45)

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

(46)

整理，得

(47)

即

(48)

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

3 算例

3.1 算例1

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

P—施加的豎向荷載；A、B、C、D—角點(diǎn)。圖3 條形基礎(chǔ)尺寸及單元?jiǎng)澐?/p>

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

(49)

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

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

3.2 算例2

算例2為某寬平臺(tái)二級(jí)路塹式邊坡[23]。一級(jí)邊坡坡率為1∶0.75，坡高為5 m； 二級(jí)邊坡坡率為1∶1，坡高為4 m。邊坡土體自上而下分為3層，幾何尺寸及單元?jiǎng)澐秩鐖D5所示。各土層的土工參數(shù)如表1所示。

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

(b)單元?jiǎng)澐諷1、S2、S3—土層編號(hào)。圖5 邊坡幾何尺寸及單元?jiǎng)澐?/p>

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

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

(50)

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

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

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

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

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

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

4 結(jié)論

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

1)多屈服面的互補(bǔ)方程的投影收縮算法避免了角點(diǎn)光滑化帶來(lái)的誤差以及不易收斂的問(wèn)題。

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

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