Cosserat連續體的本構參數對重力壩深層抗滑穩定的影響

涂承義，馬 剛

(1. 中國電建集團 華東勘測設計研究院有限公司，杭州 311122；2. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072)

0 引 言

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

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

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

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

本文計算分析Cosserat剪切模量Gc、軟化模量hp及內部長度參數l對重力壩壩基應變局部化漸進破壞的影響。

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

Cosserat連續體理論在考慮經典連續力學中的Cauchy應力以及應變的基礎上，引入了轉動自由度以及與之相對應的微曲率，與微曲率能量共軛的偶應力，見圖1。在平面應變問題中，每個點不僅具有平動的兩個自由度，還具有轉動自由度，見式(1)。

(1)

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

(3)

式中：l為Cosserat連續介質材料內部特征長度；kzx、kzy為微曲率；mzx、mzy為相應的偶應力。

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

幾何方程可表示為：

ε=Lu

(5)

線彈性的本構關系表示為：

σe=Deεe

(6)

(7)

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

當α→0，1/L→0 時，Cosserat連續體便退化為經典連續體，由此可知，經典連續體可視為Cosserat連續體的一個特例。

在Cosserat連續體理論中，應力張量不變量可表示為：

I1=σ11+σ22+σ33

(8)

(9)

(10)

J2=a1sijsij+a2sijsji+a3mijmij/l2

(11)

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

(12)

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

此外基于Cosserat連續體理論的Drucker-Prager屈服準則可表示為：

式中：c為材料黏聚力；φ為材料內摩擦角。

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

(15)

(16)

(17)

(18)

Drucker-Parager塑性勢函數G用來控制材料的塑性流動方向，其應力空間中可表示為：

(20)

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

2 Cosserat連續體模型中本構參數的研究

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

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

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

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

表1 各材料的物理力學參數Tab.1 Physico-mechanical parameters of materials

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

2.1 Cosserat剪切模量Gc對計算結果的影響

選取軟弱夾層為一層單元的模型計算分析，比較軟弱夾層的不同剪切模量Gc對結果的影響。軟弱夾層的軟化模量hp=-10 MPa，內部長度參數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，計算得到的塑性屈服區貫通的強度儲備系數Kf及有限元迭代不收斂的強度儲備系數Ku如表2所示，相應的Kf與Ku所對應的塑性屈服區分布見圖3～圖5。

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

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

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

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

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

2.2 軟化模量hp對計算結果的影響

選取軟弱夾層為一層單元的模型計算分析，比較軟弱夾層不同的軟化模量hp對計算結果的影響。軟弱夾層的剪切模量Gc=0.50G，內部長度參數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，計算得到的塑性屈服區貫通的強度儲備系數Kf及有限元迭代不收斂的強度儲備系數Ku如表3所示，相應的Kf與Ku所對應的塑性屈服區分布見圖6～圖8。

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

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

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

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

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

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

2.3 內部長度參數l對計算結果的影響

選取軟弱夾層為一層單元的模型計算分析，比較軟弱夾層不同的內部長度參數l對計算結果的影響。軟弱夾層的剪切模量Gc=0.50 G，軟化模量hp=-10MPa；軟弱夾層的厚度B=2m，內部長度參數 分別取為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。

對于不同的內部長度參數l，計算得到的塑性屈服區貫通的強度儲備系數Kf及有限元迭代不收斂的強度儲備系數Ku如表4所示，相應的Kf與Ku所對應的塑性屈服區分布見圖9～圖11。

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

表4 不同的內部長度參數對應的Kf與KuTab.4 Kf and Ku of different internal lengthes

圖9 內部長度參數為0.001 B時塑性屈服區分布Fig.9 Plastic yield zone when internal length is 0.001 B

3 結 語

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

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

圖10 內部長度參數為0.05 B時塑性屈服區分布Fig.10 Plastic yield zone when internal length is 0.05 B

圖11 內部長度參數為0.50 B時塑性屈服區分布Fig.11 Plastic yield zone when internal length is 0.50 B

(3)軟弱夾層的內部長度參數l的取值對Kf、Ku及漸進破壞過程塑性屈服區的發展影響不大。

□

[1] 王小平，徐衛亞．應變局部化問題中的正則化機制探討[J]．巖土力學，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] 李錫夔，唐洪祥．壓力相關彈塑性Cosserat連續體模型與應變局部化有限元模擬[J]．巖石力學與工程報，2005，24(9)：1 497-1 505．

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

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

[10] 蔣春艷，常曉林，周 偉．用于重力壩抗滑穩定分析的分項系數有限元方法[J]．水力發電學報，2006，25(2)：16-20．

[11] 周 偉，常曉林，徐建強．基于分項系數的重力壩深層抗滑穩定分析[J]．巖土力學，2007，28(2)：315-320．

[12] 常曉林，陸述遠．重力壩均質壩基的失穩機理研究[J]．武漢水利電力學院學報，1989，6(1)：1-6．