涂承義，馬 剛

(1. 中國電建集團(tuán) 華東勘測設(shè)計(jì)研究院有限公司，杭州 311122；2. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

0 引 言

現(xiàn)階段，重力壩深層抗滑穩(wěn)定的分析方法主要是以剛體極限平衡法為主，輔以有限元法進(jìn)行應(yīng)力和穩(wěn)定復(fù)核。由于有限元計(jì)算方法在單元網(wǎng)格的剖分、材料本構(gòu)關(guān)系模型、單元位移模式、物理力學(xué)參數(shù)以及合理評價(jià)標(biāo)準(zhǔn)等方面存在不足，導(dǎo)致其在重力壩抗滑穩(wěn)定計(jì)算領(lǐng)域的應(yīng)用受到制約。另外，重力壩壩基是非連續(xù)性以及非均勻的介質(zhì)，在材料變形接近或者超過峰值時(shí)，均勻的變形模式被狹窄的帶狀區(qū)域內(nèi)的不連續(xù)位移所代替，該變形集中的帶狀區(qū)域即為應(yīng)變局部化帶。當(dāng)壩基材料的應(yīng)變梯度較高時(shí)，在如此小的體積范圍內(nèi)其應(yīng)力與應(yīng)變將呈現(xiàn)高階次的非線性變化，經(jīng)典連續(xù)介質(zhì)理論所代表的統(tǒng)計(jì)平均值將不能如實(shí)反映壩基材料的強(qiáng)度和變形行為。因此，在高重力壩的抗滑穩(wěn)定計(jì)算中，由于水推力以及水的軟化作用，導(dǎo)致壩基中的相對軟弱的結(jié)構(gòu)面、斷層等發(fā)生較大的塑性變形，形成應(yīng)變局部化帶；那么采用經(jīng)典連續(xù)介質(zhì)理論計(jì)算分析所得到結(jié)果的合理性將無法保證。

對于應(yīng)變局部化問題，常規(guī)的方法是采用正則化機(jī)制[1]，包括非局部化理論、Cosserat連續(xù)體理論、梯度模型、黏性效應(yīng)等。Cosserat連續(xù)體理論最早是由Cosserat兄弟于1909年提出，但直到1980年代末期才被De Borst和Muhlhaus等人應(yīng)用到應(yīng)變局部化的計(jì)算分析中[2-6]。Cosserat連續(xù)體理論是最基本的考慮應(yīng)變梯度效應(yīng)的模型[7]。

由于目前還處于發(fā)展階段，與其他的一些局部化理論一樣，Cosserat連續(xù)體理論應(yīng)用于應(yīng)變局部化研究還缺乏相關(guān)的材料參數(shù)取值的資料。Cosserat連續(xù)體理論的材料參數(shù)與局部剪切帶的寬度有一定的關(guān)系；因此在某些情況下，可通過測量得到剪切帶的寬度，以剪切帶寬度與材料參數(shù)之間的關(guān)系來間接地確定Cosserat連續(xù)體理論的材料參數(shù)，這是十分重要的[8]。目前對于Cosserat連續(xù)體理論，還沒有關(guān)于剪切帶寬度的理論解。

運(yùn)用Cosserat連續(xù)體本構(gòu)模型時(shí)，為了能較好地模擬重力壩壩基深層抗滑穩(wěn)定的應(yīng)變局部化問題，Cosserat連續(xù)體本構(gòu)材料參數(shù)選取很重要，它往往是數(shù)值模擬成敗的關(guān)鍵，弄清這些問題也是進(jìn)一步進(jìn)行理論與應(yīng)用研究的基礎(chǔ)。文獻(xiàn)[8]做了關(guān)于不同的Cosserat連續(xù)體本構(gòu)材料參數(shù)對剪切帶寬度及數(shù)值結(jié)果的影響的工作，指出：①在一定范圍內(nèi)，Cosserat剪模 的取值大小對剪切帶的寬度以及計(jì)算結(jié)果幾乎沒有影響，在具體的數(shù)值模擬計(jì)算中，取Gc=0.5G為宜；②軟化模量的取值對剪切帶寬度與計(jì)算結(jié)果有很大的影響，軟化模量的絕對值越大，剪切帶的寬度越窄，荷載-位移曲線在達(dá)到峰值后的軟化段越陡；③在一定的取值范圍內(nèi)，內(nèi)部長度參數(shù) 越大，剪切帶越寬，模擬應(yīng)變軟化問題中較大變形的能力越強(qiáng)，而對峰值極限荷載的影響不大。

本文計(jì)算分析Cosserat剪切模量Gc、軟化模量hp及內(nèi)部長度參數(shù)l對重力壩壩基應(yīng)變局部化漸進(jìn)破壞的影響。

1 基于Cosserat連續(xù)體理論的D-P彈塑性模型

Cosserat連續(xù)體理論在考慮經(jīng)典連續(xù)力學(xué)中的Cauchy應(yīng)力以及應(yīng)變的基礎(chǔ)上，引入了轉(zhuǎn)動(dòng)自由度以及與之相對應(yīng)的微曲率，與微曲率能量共軛的偶應(yīng)力，見圖1。在平面應(yīng)變問題中，每個(gè)點(diǎn)不僅具有平動(dòng)的兩個(gè)自由度，還具有轉(zhuǎn)動(dòng)自由度，見式(1)。

(1)

與式(1)相對應(yīng)的應(yīng)變及應(yīng)力張量定義為：

(3)

式中：l為Cosserat連續(xù)介質(zhì)材料內(nèi)部特征長度；kzx、kzy為微曲率；mzx、mzy為相應(yīng)的偶應(yīng)力。

圖1 平面應(yīng)變問題中Cosserat連續(xù)體理論的應(yīng)力及偶應(yīng)力Fig.1 Stress and couple stress of cosserat continuum theory in plane strain case

幾何方程可表示為：

ε=Lu

(5)

線彈性的本構(gòu)關(guān)系表示為：

σe=Deεe

(6)

(7)

式中：λ=Ev/2 (1+v)為Lame常數(shù)；E為彈模模量；G為剪切模量；v為泊松比；Cc=αG，為Cosserat剪切模量；α為Cosserat剪切模量系數(shù)。

當(dāng)α→0，1/L→0 時(shí)，Cosserat連續(xù)體便退化為經(jīng)典連續(xù)體，由此可知，經(jīng)典連續(xù)體可視為Cosserat連續(xù)體的一個(gè)特例。

在Cosserat連續(xù)體理論中，應(yīng)力張量不變量可表示為：

I1=σ11+σ22+σ33

(8)

(9)

(10)

J2=a1sijsij+a2sijsji+a3mijmij/l2

(11)

J3=s11s22s33-s33(s12+s21)2/4

(12)

根據(jù)De Bost[5,6]建議，計(jì)算應(yīng)力張量不變量中的系數(shù)可取為：a1=a2=1/4,a3=1/2。

此外基于Cosserat連續(xù)體理論的Drucker-Prager屈服準(zhǔn)則可表示為：

式中：c為材料黏聚力；φ為材料內(nèi)摩擦角。

假定材料的黏聚力c服從線性軟化(硬化)規(guī)則，有：

(15)

(16)

(17)

(18)

Drucker-Parager塑性勢函數(shù)G用來控制材料的塑性流動(dòng)方向，其應(yīng)力空間中可表示為：

(20)

式中：ψ為材料的膨脹角。

2 Cosserat連續(xù)體模型中本構(gòu)參數(shù)的研究

本文計(jì)算分析考慮一個(gè)壩高150 m的混凝土重力壩，重力壩壩基下存在一條2 m厚、上游出露、傾向于下游的軟弱夾層。為了研究比較單元的網(wǎng)格密度對重力壩壩基漸進(jìn)破壞過程模擬的影響，對比采用基于Cosserat連續(xù)體理論的Drucker-Prager彈塑性模型與經(jīng)典的彈塑性模型的計(jì)算結(jié)果及收斂性對單元網(wǎng)格密度的依賴，2 m厚的軟弱夾層分別采用1層、2層及3層單元?jiǎng)澐?，具體模型見圖2。

圖2 重力壩壩基抗滑穩(wěn)定模型Fig.2 The anti-sliding stability model of gravity dam

本文計(jì)算分析采用基于有限元的強(qiáng)度儲備系數(shù)法對重力壩壩基的漸進(jìn)破壞過程進(jìn)行模擬。本文主要采用壩基等效塑性應(yīng)變的分布范圍進(jìn)行對比分析；采用塑性屈服區(qū)貫通法、有限元迭代不收斂判據(jù)作為壩基整體失穩(wěn)的2個(gè)判據(jù)，因此得到相對應(yīng)的2個(gè)強(qiáng)度儲備系數(shù)：假定對應(yīng)塑性屈服區(qū)貫通的強(qiáng)度儲備系數(shù)為Kf，對應(yīng)有限元迭代不收斂的強(qiáng)度儲備系數(shù)為Ku[10-12]。

計(jì)算中采用的材料物理力學(xué)參數(shù)見表1?；鶐r和軟弱夾層采用基于Cosserat連續(xù)體理論的Drucker-Prager彈塑性模型，壩體混凝土采用線彈性模型。為了能準(zhǔn)確地比較軟弱夾層不同單元網(wǎng)格密度對壩基漸進(jìn)破壞過程模擬的影響，本文分析計(jì)算中不考慮基巖的軟化，只考慮軟弱夾層的軟化；因此基巖的軟化模量hp取為0，Cosserat剪模系數(shù)取α為0.001，即Gc=0.001G，內(nèi)部特征長度l取為0.000 1 m。軟弱夾層的Cosserat剪模系數(shù)α、軟化模量hp、內(nèi)部特征長度l則根據(jù)不同的計(jì)算條件選取不同的值。

表1 各材料的物理力學(xué)參數(shù)Tab.1 Physico-mechanical parameters of materials

文獻(xiàn)[9]列出了基于Cosserat連續(xù)體理論Drucker- Prager彈塑性模型在ABAQUS中實(shí)現(xiàn)的UEL自定義單元接口二次開發(fā)；采用強(qiáng)度儲備系數(shù)法，模擬重力壩壩基應(yīng)變局部化漸進(jìn)破壞過程，并對比Cosserat連續(xù)體理論與經(jīng)典連續(xù)體理論的結(jié)果。本文主要分析Cosserat剪切模量Gc、軟化模量hp及內(nèi)部長度參數(shù)l對重力壩壩基漸進(jìn)破壞的影響。

2.1 Cosserat剪切模量Gc對計(jì)算結(jié)果的影響

選取軟弱夾層為一層單元的模型計(jì)算分析，比較軟弱夾層的不同剪切模量Gc對結(jié)果的影響。軟弱夾層的軟化模量hp=-10 MPa，內(nèi)部長度參數(shù)l=0.01 m，彈性模量E=1 GPa，剪切模量G=0.37 GPa；Cosserat剪切模量Gc分別取為0.05G、0.10G、0.25G、0.50G、1.0G、2.5G、10.0G、50.0G。

對于不同的Cosserat剪切模量Gc，計(jì)算得到的塑性屈服區(qū)貫通的強(qiáng)度儲備系數(shù)Kf及有限元迭代不收斂的強(qiáng)度儲備系數(shù)Ku如表2所示，相應(yīng)的Kf與Ku所對應(yīng)的塑性屈服區(qū)分布見圖3～圖5。

表2 不同的Cosserat剪切模量對應(yīng)的Kf與KuTab.2 Kf and Ku of different Cosserat shear modulus

圖3 Cosserat剪切模量為0.05 G時(shí)塑性屈服區(qū)分布Fig.3 Plastic yield zone when Cosserat shear modulus is 0.05 G

圖4 Cosserat剪切模量為0.5 G時(shí)塑性屈服區(qū)分布Fig.4 Plastic yield zone when Cosserat shear modulus is 0.5 G

圖5 Cosserat剪切模量為50 G時(shí)塑性屈服區(qū)分布Fig.5 Plastic yield zone when Cosserat shear modulus is 50 G

由表2及圖3～圖5可知，當(dāng)Cosserat剪切模量Gc為0.05G時(shí)，Cosserat連續(xù)體接近退化為經(jīng)典連續(xù)體，此時(shí)Kf為2.7；較 取其他值時(shí)的Kf=2.6時(shí)大0.1。Cosserat剪切模量Gc小于0.5G時(shí)，Ku=3.1；Cosserat剪切模量Gc大于0.5G時(shí)，Ku=3.2。對于同一強(qiáng)度儲備系數(shù)，不同的Cosserat剪切模量 所對應(yīng)的塑性屈服區(qū)的分布基本相同。因此，在一定的取值范圍內(nèi)，Cosserat剪切模量Gc對Kf與Ku及漸進(jìn)破壞過程影響不大。De Borst，Steinmann，李錫夔等在將Cosserat連續(xù)體理論應(yīng)用于應(yīng)變局部化的研究中，基本按Gc=0.5G，本文亦推薦Gc取0.5G。

2.2 軟化模量hp對計(jì)算結(jié)果的影響

選取軟弱夾層為一層單元的模型計(jì)算分析，比較軟弱夾層不同的軟化模量hp對計(jì)算結(jié)果的影響。軟弱夾層的剪切模量Gc=0.50G，內(nèi)部長度參數(shù)l=0.01 m；令h0=-10 MPa，軟化模量hp分別取為0.01h0、0.1h0、0.2h0、0.5h0、1.0h0、2.0h0、5.0h0、10.0h0、100.0h0，即-0.1、-1.0、-2.0、-5.0、-10.0、-20.0、-50.0、-100.0、-1 000.0 MPa。

對于不同的軟化模量hp，計(jì)算得到的塑性屈服區(qū)貫通的強(qiáng)度儲備系數(shù)Kf及有限元迭代不收斂的強(qiáng)度儲備系數(shù)Ku如表3所示，相應(yīng)的Kf與Ku所對應(yīng)的塑性屈服區(qū)分布見圖6～圖8。

由表3及圖6～圖8可知，當(dāng)軟化模量取值 在-0.1與-10.0 MPa之間時(shí)，各軟化模量對應(yīng)的塑性屈服區(qū)貫通對應(yīng)的Kf及有限元迭代不收斂對應(yīng)的Ku均相同：Kf=2.6，Ku=3.2。當(dāng)hp達(dá)到-50.0 MPa時(shí)，強(qiáng)度儲備系數(shù)為2.0時(shí)，計(jì)算由于迭代不收斂而退出，壩基塑性屈服區(qū)未貫通，即Ku=2.0，Kf不存在；當(dāng)hp=-100.0 MPa時(shí)，Ku=1.3。當(dāng)hp=-1 000.0 MPa時(shí)，當(dāng)強(qiáng)度儲備系數(shù)降低時(shí)，壩基中完好基巖出現(xiàn)較大的塑性屈服區(qū)，而軟弱夾層卻沒有屈服，與實(shí)際認(rèn)知不相符，計(jì)算結(jié)果失去意義。

因此，軟化模量的取值對壩基應(yīng)變局部化的漸進(jìn)破壞模擬有較大的影響。軟化模量的取值在一定的范圍內(nèi)，Kf及Ku相同，漸進(jìn)破壞過程基本相同；當(dāng)軟化模量的絕對值超過一定值時(shí)，計(jì)算在壩基未貫通時(shí)即迭代不收斂；當(dāng)軟化模量的絕對值進(jìn)一步增大時(shí)，降低壩基抗剪斷強(qiáng)度過程中軟弱夾層未屈服，與實(shí)際不符。

表3 不同的軟化模量對應(yīng)的Kf與KuTab.3 Kf and Ku of different softening moduluses

圖6 軟化模量為-10.0 MPa時(shí)塑性屈服區(qū)分布Fig.6 Plastic yield zone when softening modulu is -10.0 MPa

圖7 軟化模量為-50.0 MPa時(shí)塑性屈服區(qū)分布Fig.7 Plastic yield zone when softening modulu is -50.0 MPa

圖8 軟化模量為-1 000.0 MPa時(shí)塑性屈服區(qū)分布Fig.8 Plastic yield zone when softening modulu is -1 000.0 MPa

2.3 內(nèi)部長度參數(shù)l對計(jì)算結(jié)果的影響

選取軟弱夾層為一層單元的模型計(jì)算分析，比較軟弱夾層不同的內(nèi)部長度參數(shù)l對計(jì)算結(jié)果的影響。軟弱夾層的剪切模量Gc=0.50 G，軟化模量hp=-10MPa；軟弱夾層的厚度B=2m，內(nèi)部長度參數(shù) 分別取為0.001 B、0.005 B、0.01 B、0.02 B、0.03 B、0.04 B、0.05 B、0.10 B、0.20 B、0.50 B，即0.002、0.01、0.02、0.04、0.06、0.08、0.10、0.20、0.40、1.00m。

對于不同的內(nèi)部長度參數(shù)l，計(jì)算得到的塑性屈服區(qū)貫通的強(qiáng)度儲備系數(shù)Kf及有限元迭代不收斂的強(qiáng)度儲備系數(shù)Ku如表4所示，相應(yīng)的Kf與Ku所對應(yīng)的塑性屈服區(qū)分布見圖9～圖11。

文獻(xiàn)[4-8]在將Cosserat連續(xù)體理論應(yīng)用于應(yīng)變局部化的研究中，內(nèi)部長度參數(shù)基本取在0.01 B～0.1 B范圍內(nèi)，均能產(chǎn)生一定的正則化效果，并對剪切帶寬度影響大。而本文計(jì)算所得到的軟弱夾層內(nèi)部長度參數(shù)的取值對Kf、Ku及漸進(jìn)破壞過程塑性屈服區(qū)的發(fā)展基本沒影響，原因在于軟弱夾層雖然由于最先屈服，在計(jì)算迭代收斂中起決定作用，但在整個(gè)系統(tǒng)中所占的比例較小，對于基巖的屈服破壞范圍的影響較小，因此內(nèi)部長度參數(shù)的取值不同，對漸進(jìn)破壞過程塑性屈服區(qū)的發(fā)展影響不大。

表4 不同的內(nèi)部長度參數(shù)對應(yīng)的Kf與KuTab.4 Kf and Ku of different internal lengthes

圖9 內(nèi)部長度參數(shù)為0.001 B時(shí)塑性屈服區(qū)分布Fig.9 Plastic yield zone when internal length is 0.001 B

3 結(jié) 語

(1)在一定的取值范圍內(nèi)，Cosserat剪切模量Gc對Kf、Ku及漸進(jìn)破壞過程影響不大。

(2)軟化模量hp的取值對壩基應(yīng)變局部化的漸進(jìn)破壞模擬有較大的影響。hp的取值在一定范圍內(nèi)，Kf、Ku及漸進(jìn)破壞過程基本相同；當(dāng)hp的絕對值超過某一定值時(shí)，降低壩基抗剪斷強(qiáng)度過程中軟弱夾層未屈服，與實(shí)際不符。

圖10 內(nèi)部長度參數(shù)為0.05 B時(shí)塑性屈服區(qū)分布Fig.10 Plastic yield zone when internal length is 0.05 B

圖11 內(nèi)部長度參數(shù)為0.50 B時(shí)塑性屈服區(qū)分布Fig.11 Plastic yield zone when internal length is 0.50 B

(3)軟弱夾層的內(nèi)部長度參數(shù)l的取值對Kf、Ku及漸進(jìn)破壞過程塑性屈服區(qū)的發(fā)展影響不大。

□

[1] 王小平，徐衛(wèi)亞．應(yīng)變局部化問題中的正則化機(jī)制探討[J]．巖土力學(xué)，2008，29(11)：2 997-3 002．

[2] Muhlhaus H B. Application of Cosserat theory in numerical solutions of limit load problems [J]. Ing. Arch., 1989,59(2):124-137.

[3] Muhlhaus H B, Vardoulakis I. The thickness of shear bands in granular materials [J]. Geotechnique, 1987,37(3):271-283.

[4] de Borst R, Sluys L J. Localization in a Cosserat continuum under static and dynamic loading conditions [J]. Comput. Methods Appl. Mech. Eng., 1991,90(2):805-827.

[5] de Borst R. A generalization of J2-flow theory for polar continua[J]. Comput. Methods Appl. Mech. Eng., 1993,103(3):347-362.

[6] de Borst R. Simulation of strain localization：a reappraisal of the Cosserat continuum[J]. Eng. Comput., 1991,8(4):317-332.

[7] 李錫夔，唐洪祥．壓力相關(guān)彈塑性Cosserat連續(xù)體模型與應(yīng)變局部化有限元模擬[J]．巖石力學(xué)與工程報(bào)，2005，24(9)：1 497-1 505．

[8] 唐洪祥．基于Cosserat連續(xù)體模型的應(yīng)變局部化有限元模擬[D]．遼寧大連：大連理工大學(xué)，2007.

[9] 馬 剛，涂承義，常曉林，等．基于 Cosserat 連續(xù)體理論的重力壩深層抗滑穩(wěn)定分析[J]．四川大學(xué)學(xué)報(bào)(工程科學(xué)版)，2012，44(5)：57-63．

[10] 蔣春艷，常曉林，周 偉．用于重力壩抗滑穩(wěn)定分析的分項(xiàng)系數(shù)有限元方法[J]．水力發(fā)電學(xué)報(bào)，2006，25(2)：16-20．

[11] 周 偉，常曉林，徐建強(qiáng)．基于分項(xiàng)系數(shù)的重力壩深層抗滑穩(wěn)定分析[J]．巖土力學(xué)，2007，28(2)：315-320．

[12] 常曉林，陸述遠(yuǎn)．重力壩均質(zhì)壩基的失穩(wěn)機(jī)理研究[J]．武漢水利電力學(xué)院學(xué)報(bào)，1989，6(1)：1-6．