橫觀各向同性巖土材料應變局部化現(xiàn)象的有限元模擬*

常江芳 徐遠杰 楚錫華

(1.武漢大學 工程力學系， 湖北 武漢 430072； 2.石家莊鐵道大學 工程力學系, 河北 石家莊 050043)

橫觀各向同性巖土材料應變局部化現(xiàn)象的有限元模擬*

常江芳1，2徐遠杰1?楚錫華1

(1.武漢大學 工程力學系， 湖北 武漢 430072； 2.石家莊鐵道大學 工程力學系, 河北 石家莊 050043)

巖土材料的各向異性性質與其微結構緊密相關，文中考慮到內摩擦角是應力狀態(tài)、組構張量及材料主方向的函數，發(fā)展了適用于橫觀各向同性巖土材料的修正的Drucker-Prager屈服準則.基于Cosserat連續(xù)體理論，推導了該準則的一致性映射返回算法，形成了一致性切線模量矩陣，并利用有限元軟件Abaqus的用戶單元子程序(UEL)進行了數值實現(xiàn).通過將積分點材料強度隨材料主方向及各向異性程度的變化關系與理論結果進行比較，驗證了程序開發(fā)的正確性.數值算例重點分析了材料主方向和各向異性參數對結構極限承載力及破壞模式的影響，且與經典連續(xù)體的結果進行了比較.結果表明，文中方法能夠較好地模擬具有橫觀各向同性的巖土材料的應變局部化現(xiàn)象.

巖土材料；巖土力學；橫觀各向同性；Cosserat連續(xù)體；屈服準則；應變局部化

自然界的巖土材料一般均呈現(xiàn)很強的各向異性.各向異性又分為固有各向異性和誘導各向異性.由于其沉積過程、孔隙、裂紋等固有屬性而使其在原位狀態(tài)下表現(xiàn)出的各向異性稱為固有各向異性；巖土材料在非比例加載或塑性變形的影響下，由于主應力方向不斷變化，導致其細觀結構發(fā)生演化而表現(xiàn)出的各向異性稱為誘導各向異性.關于此兩種各向異性，國內外學者做了相關的理論及實驗研究[1- 3]，提出了多種各向異性屈服準則[4- 7].姚仰平等[4]以橫觀各向同性土為例，闡述了同時考慮材料外部應力分布與材料內部強度分布下各向異性土的破壞機制.賈乃文[5]以Hill正交異性屈服準則為基礎，討論了橫觀各向同性條件下，巖土中出現(xiàn)塑性屈服的Drucker準則修正表達式.Duveau等[6]對一些強各向異性材料的屈服準則做了一個綜合評價.Gao等[7]在屈服函數的摩擦系數中引入了一個各向異性變量，考慮它為應力不變量及微結構張量的聯(lián)合不變量的函數，提出了一個適用于橫觀各向同性材料的屈服準則.然而基于唯象方法提出的各向異性屈服準則一般包含較多的材料參數，且這些參數與材料微觀結構之間的關系不明確.Pietruszczak等[8- 9]提出了一個各向異性屈服準則，基于物質點給出了各向異性程度及加載方向對材料強度影響的理論分析，其材料參數是應力和微結構張量聯(lián)合不變量的顯式函數，使得實驗驗證更易實現(xiàn).余天堂等[10- 11]基于Pietruszczak等[8- 9]提出的各向異性屈服準則，考慮各向異性參數和單軸抗壓強度是一個由微結構張量和加載方向表示的分布函數，給出了沉積巖的一種各向異性模型.Zhong等[12]對Pietruszczak-Morz屈服準則中的材料參數進行了分析.在Pietruszczak等[8- 9]的研究基礎上，Lade[13- 15]推導了一個基于 Lade 強度準則的各向異性強度準則，并展開了相應的理論及實驗研究.

然而上述研究多關注的是對各向異性破壞準則的描述，是各向異性對一個物質點破壞的影響.關于各向異性對結構承載力及破壞模式影響的數值模擬仍相對較少.Chu等[16]在Pietruszczak等[8- 9]提出的各向異性屈服準則的基礎上，考慮內摩擦角為組構張量、應力狀態(tài)及材料主方向的函數，對Drucker-Prager屈服準則進行了修正，以橫觀各向同性材料為例研究了各向異性對應變局部化的影響，然而不足之處是其基于經典連續(xù)體理論，在伴隨應變軟化出現(xiàn)應變局部化現(xiàn)象的同時會產生網格依賴性.為克服這一問題，文中基于Cosserat連續(xù)體理論[17- 19]，推導了修正的Drucker-Prager屈服準則的一致性映射返回算法及切線模量矩陣，重點分析了材料主方向和各向異性程度對橫觀各向同性材料極限承載力及破壞模式的影響，最后通過與經典連續(xù)體計算結果進行比較來驗證解決網格依賴性問題的有效性.

1 基于Cosserat連續(xù)體的橫觀各向同性彈塑性模型

1.1 Cosserat連續(xù)體的基本方程

與經典連續(xù)體相比，Cosserat連續(xù)體中引入了獨立的轉動自由度，平面應變情況下，其位移矢量u可表示為[17]

u=[uxuyωz]T

(1)

其中，ux、uy為平動位移，ωz為獨立的轉動位移.

應變矢量ε除了與經典意義上所對應的正應變和剪應變有關以外，還包含兩個微曲率，如下所示：

ε=[εxεyεzεxyεyxlcκxzlcκyz]T

(2)

其中：

(3)

可以看到這里引入了一個內部長度參數lc，保證了Cosserat連續(xù)體的控制方程在計算過程中的正定性，可有效克服網格依賴性問題.

幾何方程可表示如下：

ε=Lu

(4)

其中，

(5)

Cosserat連續(xù)體的應力矢量σ引入了兩個偶應力mxz和myz(如圖1所示)：

(6)

圖1 平面應變條件下應力及偶應力示意圖

Fig.1 Sketches of stress and couple stress under strain plane condition

彈性應力-應變關系可以表示為

σ=Deε

(7)

其中本構矩陣為

(8)

準靜態(tài)條件下，若忽略體積力及體力偶的影響，則平衡方程可表示為

LTσ=0

(9)

1.2 修正的Drucker-Prager準則

適用于各向同性材料的Drucker-Prager準則以應力不變量的形式表示為

F=q+Aφp+Bφ=0

(10)

其中,p和q可分別用應力第一不變量及偏應力第二不變量表示如下:

(11)

(12)

其中，P矩陣可表示為

(13)

式(10)中參數Aφ、Bφ為內摩擦角φ和黏聚力c的函數，若Drucker-Prager準則為Mohr-Coulomb準則的內接圓，則可將其分別表示為

(14)

通常對黏聚力考慮線性硬化或軟化，即

(15)

基于應力不變量與組構張量，Pietruszczak與Mroz給出了一個能夠描述各向異性材料的屈服準則[8]：

η=η0(1+Ωijlilj)

(16)

其中：η為描述屈服函數的材料參數，與應力狀態(tài)有關；η0為η的平均值；Ωij為描述η偏離η0程度的偏量.對正交各向異性材料Ωij有兩個不等的特征值，對橫觀各向同性材料，其可用一個標量描述，對各向同性材料，Ωij為零，此時η=η0.li、lj為相對于材料主軸的加載方向.式(16)也可用Ωij的主值表示為

(17)

(18)

圖2 材料主方向示意圖

加載向量按如下方式計算：

(19)

為了簡化，僅考慮內摩擦角的各向異性，即令

(20)

則有

(21)

式(10)、式(3)、式(1)構成了描述橫觀各向同性材料的屈服函數.

對于巖土材料，通常采用非關聯(lián)流動法則，塑性勢可取為

G=q+Aψp+Bψ

(22)

類似地：

(23)

其中，ψ為膨脹角.

由一致性算法最終得到更新后的應力[17- 18]，

(24)

其中：

(25)

(26)

(27)

(28)

其中FE=qE+AφpE+BE.

(29)

以及一致性彈塑性切線模量矩陣：

(30)

其中:

(31)

(32)

可以驗證，當采用關聯(lián)流動法則時，該矩陣具有主對稱性.

2 數值算例

以Abaqus軟件的UEL子程序接口為平臺，對文中發(fā)展的單元類型和彈塑性本構模型進行了數值實現(xiàn)[20].以平板壓縮模型為例，單元類型為平面應變8節(jié)點減縮單元，模型尺寸為0.6m×0.8m，上下邊界約束x方向自由度，左右邊界自由，上下邊界施加對稱均布位移荷載uy，如圖3(a)所示.計算中所采用的材料參數如表1所示.

圖3 有限元模型

參數數值參數數值E/Pa30×106c0/Pa0．06×106υ0．3hcp/Pa0．01×106β/(°)0～180Ω10、-0．05、-0．1、-0．15φ/(°)20Gc/Pa10×106ψ/(°)15lc/m0．02

1)E為彈性模量；υ為泊松比；β為材料主方向角.

Pietruszczak等[8]基于他們提出的各向異性屈服準則對橫觀各向同性材料物質點的理論分析顯示，對于不同的加載方向，材料軸向強度隨Ω1的演化規(guī)律不同，存在一個轉折點，該點材料的軸向強度將獨立于各向異性程度，且Zhong等[12]證明了此轉折點所對應的材料主方向與材料固有屬性有關.這里，為與已有理論結果進行對比，選取了結構的兩個代表性單元，標號為Ⅰ和Ⅱ，分別位于板上邊緣中間位置和剪切帶附近，如圖4所示，給出了兩個單元1號積分點(見圖3(b))上等效應力隨加載位移的變化曲線，其中參數取β=30°，Ω1=-0.1.易知峰值qmax代表了該點的強度，圖5(a)和5(b)則給出了積分點(物質點)上歸一化的材料強度隨各向異性程度及材料主方向的變化趨勢，可見其趨勢與已有理論結果[8]較為吻合(見圖5(c),其中法向強度ζ=Δσ1/(2σ0))，從而驗證了本程序開發(fā)的正確性.下面主要就材料主方向β及各向異性程度Ω1對結構極限承載力及破壞模式的影響做詳細分析.

圖4 積分點等效應力隨加載位移變化曲線

Fig.4 Curves of the equivalent stress versus the loading displacement in the integration points

2.1 材料主方向β的影響

由圖2可知，β為各向同性面相對x軸逆時針轉動的角度，描述了材料主軸相對整體坐標系的偏移，β的改變也可理解為加載方向的改變.圖6為Ω1=-0.10，加載位移為0.05 m時等效塑性應變分布云圖.值得注意的是，Ω1取負值意味著各向同性面法向(e(2)方向)有較大的摩擦強度[14].可以看出剪切帶呈現(xiàn)了X型分布，當β=0°和β=90°時，由于加載方向和材料主方向保持一致或垂直，剪切帶呈對稱分布，如圖6(a)和6(e)所示.當β≠0°和β≠90°時，剪切帶兩個分支呈現(xiàn)不對稱性，沿各向同性面方向的剪切帶較寬，等效塑性應變值較高，文中稱之為強剪切帶，與之共軛的另一條剪切帶稱為弱剪切帶，如圖6(b)和6(h)所示，與姚仰平等[4]的理論分析相比，其破壞面的位置可能就是強剪切帶出現(xiàn)的位置，數值模擬顯示還存在一個弱剪切帶，只是在破壞時強剪切帶占據了主導地位.由圖6(b)-6(d)可見，當0°<β<90°時剪切帶呈“/”型分布，且隨著β的增大“/”方向剪切帶逐漸變寬.反之，當90°<β<180°時，剪切帶呈“”型分布，且隨著β的增大“”方向剪切帶逐漸變窄.對于一對互補的β角，如β=30°和β=150°，等效塑性應變峰值相等，但剪切帶分布模式相反.隨著β從0°演化到180°，強剪切帶上等效塑性應變峰值出現(xiàn)了增大-降低-增大-降低的對稱過程，弱剪切帶則具有互補的變化過程，如圖7所示.

圖5 積分點材料強度隨材料主方向β及各向異性程度Ω1的變化曲線

Fig.5 Variation curves of material strength withβandΩ1in integration points

圖8給出了隨β變化材料的承載力-位移曲線,可以清楚地看到當各向同性面位于水平位置，即β=0°時，結構承載力最大，隨著β的增大承載力逐漸降低，90°時最低，這與橫觀各向同性材料屈服準則中摩擦強度參數各向異性的分布規(guī)律一致.

2.2 各向異性參數Ω1的影響

如前所述，Ωij為描述η偏離η0程度的偏張量，對于橫觀各向同性材料可以用一個標量Ω1來描述.需指出文中各向異性參數取值來自于文獻[8]，依次取Ω1=0，-0.05，-0.1,-0.15,β=30°,加載位移為0.05 m，其他參數同表1.

圖9給出了等效塑性應變隨各向異性程度變化的結果.可以看出，各向異性程度越大，X型剪切帶的非對稱性越強，強剪切帶等效塑性應變的峰值越大，弱剪切帶峰值則越小.圖10為承載力-位移曲線，隨著Ω1絕對值的增大結構承載力逐漸增大.這里同時給出了歸一化的承載力隨各向異性程度及材料主方向的綜合變化規(guī)律，如圖11所示，其中，縱軸表示各承載力峰值除以Ω1=0時的承載力峰值，可見結構整體響應亦呈現(xiàn)出了與圖5類似的規(guī)律，即在β約為47°的位置出現(xiàn)轉折點，β小于該值時，結構承載力隨Ω1絕對值的增大而增大，大于該值時則相反.

圖7 強剪切帶和弱剪切帶上等效塑性應變峰值變化規(guī)律

Fig.7 Variation of the peak value of equivalent plastic strain in the intense shear band and the weak shear band

圖8 不同材料主方向承載力-位移曲線

Fig.8 Curves of load-displacement with different material direction

圖10 不同各向異性程度下承載力-位移曲線

Fig.10 Curves of load-displacement with different anisotropic degrees

2.3 網格依賴性調查

基于經典連續(xù)體理論，采用有限單元法計算彈塑性問題，當應變軟化和局部化現(xiàn)象發(fā)生時，不可避免地會出現(xiàn)網格依賴性問題.本節(jié)通過對經典連續(xù)體和Cosserat連續(xù)體兩種理論框架下的計算結果進行比較，調查了有限元計算的網格依賴性問題.

圖9 等效塑性應變分布(β=30°)

圖11 結構承載力隨材料主方向β及各向異性程度Ω1的變化

Fig.11 Variation of bearing capacity of structure withβandΩ1

取材料主方向β=90°，Ω1=-0.1，其他材料參數同表1，網格密度分別取6×8，12×16，18×24，24×32，30×40，60×80進行比較.圖12給出了基于經典連續(xù)體等效塑性應變分布圖，可以看到隨著網格加密，等效塑性應變的峰值逐漸增大，且剪切帶寬度明顯變窄.圖13顯示結構極限承載力隨網格密度增大而逐漸降低，且模擬軟化段的能力亦逐漸降低.圖14則為基于Cosserat連續(xù)體得到的等效塑性應變分布圖，隨著網格加密，等效塑性應變峰值有所增大但剪切帶寬度基本保持不變.圖15顯示雖然結構承載力有些許降低，但其變化越來越小，計算結果最終收斂于較為接近的值.需說明的是，圖中網格為6×8時計算結果與其他結果相比差異較大，這是由于網格稀疏到一定程度，有限元計算本身的誤差所致.經過比較可見引入Cosserat連續(xù)體后，網格依賴性問題得到了有效的解決.

圖13 經典連續(xù)體不同網格密度承載力-位移曲線

Fig.13 Curves of load-displacement with different mesh density based on classical continuum theory

圖15 Cosserat連續(xù)體不同網格密度承載力-位移曲線

Fig.15 Curves of load-displacement with different mesh density based on Cosserat continuum theory

3 結論

基于Pietruszczak 和 Mroz提出的各向異性屈服準則，考慮摩擦角為材料主方向及應力張量與微結構張量聯(lián)合不變量的函數，提出了一個適用于橫觀各向同性巖土材料的修正的Drucker-Prager屈服準則，并基于Cosserat連續(xù)體推導了其一致性映射返回算法和切線模量矩陣.數值算例重點分析了材料主方向及各向異性程度對結構極限承載力及破壞模式的影響，結果表明：

(1)材料主方向β=0°或β=90°時剪切帶呈對稱分布，前者結構承載力最大，后者最低.0°<β<90°時，隨β增大，剪切帶逐漸變寬，結構承載力逐漸降低，90°<β<180°則相反.β從0°到180°，強剪切帶等效塑性應變峰值經歷了增大-降低-增大-降低的對稱過程，弱剪切帶則相反.β互補時剪切帶呈現(xiàn)相反的分布模式.

(2)Ω1絕對值越大，材料各向同性面的強度越低；對于β≠0°、β≠90°的非對稱情況，剪切帶呈現(xiàn)的X型分布的非對稱性更強.承載力變化趨勢存在一個轉折點，約為β=47°，β小于該值時，結構承載力隨Ω1絕對值的增大而增大，β大于該值時則相反.

(3)與經典連續(xù)體的模擬結果比較表明，文中基于Cosserat連續(xù)體的數值方法較好地解決了網格依賴性問題.

[1] ARTHUR J R F,Menzies B.Inherent anisotropy in a sand [J].Geotechnique,1972,22(1):115- 128.

[2] ARTHUR J R F,CHUA K S,DUNSTAN T.Induced anisotropy in a sand [J].Geotechnique,1979,27(1):13- 30.

[3] ODA M.Inherent and induced anisotropy in plasticity theory of granular soils [J].Mechanics of Materials,1993,16(1/2):35- 45.

[4] 姚仰平,祝恩陽.橫觀各向同性土的簡明破壞機制解釋 [J].巖土力學,2014,35(2):328- 333. YAO Yang-ping,ZHU En-yang.Concise interpretation of damage mechanism for cross-anisotropic soil [J].Journal of Geotechnical Mechanics,2014,35(2):328- 333.

[5] 賈乃文.巖土工程中橫觀各向同性Drucker屈服準則 [J].華南理工大學學報(自然科學版)，1993,21(2):49- 54. JIA Nai-wen.Drucker yield criterion of horizontally isotropic rock and soil [J].Journal of South China University of Technology(Natural Science Edition),1993,21(2):49- 54.

[6] DUVEAU G,SHAO J F,HENRY J P.Assessment of some failure criteria for strongly anisotropic materials [J].Mechanics of Cohesive-Frictional Materials,1998,3(1):1- 26.[7] GAO Z W,ZHAO J D,YAO Y P.A generalized anisotropic failure criterion for geomaterials [J].International Journal of Solids and Structures,2010,47(22/23):3166- 3185.

[8] PIETRUSZCZAK S,MROZ Z.Formulation of anisotropic failure criteria incorporating a microstructure tensor [J].Computer and Geotechnics,2000,26(2):105- 112.

[9] PIETRUSZCZAK S,MROZ Z.On failure criteria for anisotropic cohesive-frictional materials [J].International Journal for Numerical and Analytical Methods in Geomechanics,2001,25(5):509- 524.

[10] 余天堂,盧應發(fā),SHAO J F,等.沉積巖的一種各向異性模型 [J].巖土力學,2002,23(1):47- 50. YU Tian-tang,LU Ying-fa,SHAO J F,et al.An anisotropic model of sedimentary rocks [J].Journal of Geotechnical Mechanics,2002,23(1):47- 50.

[11] 余天堂.巖土材料固有各向異性的模擬 [J].巖石力學與工程學報,2004,23(10):1604- 1607. YU Tian-tang.Modeling of inherent anisotropy for geotechnical material [J].Chinese Journal of Rock Mechanics and Engineering,2004,23(10):1604- 1607.

[12] ZHONG S Y,XU W Y,LING D S.Influence of the parameters in the Pietruszczak-Mroz anisotropic failure criterion [J].International Journal of Rock Mechanics and Mining Sciences,2011,48(6):1034- 1037.

[13] LADE P V.Modeling failure in cross-anisotropic frictional materials [J].International Journal of Solids and Structures,2007,44(16):5146- 5162.

[14] LADE P V.Failure criterion for cross-anisotropic soils [J].Journal of Geotechnical and Geoenvironmental Engineering,2008,134(1):117- 124.

[15] LADE P V.Three-dimensional failure in geomaterials:experimentation and modeling [M]∥Constitutive Modeling of Geomaterials.[S.l.]:Springer Berlin Heidelberg,2013:47- 58.

[16] CHU X H,CHANG J F,JIANG Q H,et al.Numerical simulation of progressive failure in transversely isotropic geomaterials[C]∥Advance in Heterogeneous Material Mechanics.Lancaster:Destech Publication Inc.,2011:965- 968.

[17] DE BORST R.Simulation of strain localization:a reappraisal of the Cosserat continuum [J].Engineering Computations,1991,8(4):317- 332.

[18] LI X K,TANG H X.A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modeling of strain localization [J].Computers & Structures,2005,83(1):1- 10.

[19] 余村,楚錫華,唐洪祥,等.基于Cosserat連續(xù)體的顆粒破碎影響研究 [J].巖土力學,2013,34(增刊1):67- 79. YU Cun,CHU Xi-hua,TANG Hong-xiang,et al.Study of effect of particle breakage based on Cosserat continuum [J].Journal of Geotechnical Mechanics,2013,34(Sup 1):67- 79.

[20] ARSLAN H,STURE S.Finite element analysis of localization and micro-macro structure relation in granular materials(Part I):Formulation [J].Acta Mechanica,2008,197(3/4):135- 152.

Finite Element Simulation of Strain Localization in Transversely Isotropic Geomaterials

CHANGJiang-fang1，2XUYuan-jie1CHUXi-hua1

(1.Engineering Mechanics Department,Wuhan University,Wuhan 430072,Hubei,China; 2. Mechanics Engineering Department,Shijiazhuang Tiedao University,Shijiazhuang 050043,Hebei,China)

The anisotropic properties of geomaterials are significantly related to their inherent microstructures. In this paper,a modified Drucker-Prager yield criterion for transversely isotropic materials is developed by evaluating the internal friction angle with the stress state,the microstructure tensor and the material principal direction.Then,based on the Cosserat continuum theory,a consistent return mapping algorithm for the modified criterion is formulated,and a consistent tangent modulus matrix is achieved.Moreover,the codes are implemented through the user defined element subroutine(UEL) in the finite element software Abaqus,and the correctness of the program is verified by comparing the theoretical results with the relationships of the material strength to the principal direction and anisotropic degree of the material in the integration points.Finally,the influences of the principal direction and anisotropic degree of the material on the bearing capacity and the failure mode of the structure are emphatically analyzed by numerical examples,which are then compared with the results based on the classical continuum theory.It is found that the above-mentioned method is effective in simulating the strain localization of transversely isotropic geomaterials.

geomaterials;geomechanics;transverse isotropy;Cosserat continuum;yield criterion; strain localization

2015- 01- 19

國家自然科學基金資助項目(11172216)；國家重點基礎研究發(fā)展計劃項目(2010CB731502)；湖北省自然科學基金資助項目(2013CFB287)；河北省杰出青年科學基金資助項目(E2015210040) Foundation items: Supported by the National Natural Science Foundation of China(11172216),the National Program on Key Basic Research Project of China(2010CB731502)and the Natural Science Foundation of Hubei Province(2013CFB287)

常江芳(1988-)，女，博士，講師，主要從事巖土材料強度與破壞的數值研究.E-mail:cjf881024@163.com

? 通信作者: 徐遠杰(1956-)，男，博士，教授，主要從事工程結構破壞數值仿真研究.E-mail:xyj9781@163.com

1000- 565X(2016)10- 0070- 11

O 34

10.3969/j.issn.1000-565X.2016.10.011