非共軸特性對巖土材料應變局部化行為影響的數值研究

常江芳，王偉，牛慶合，聞磊

(1.石家莊鐵道大學力學系，石家莊 050043；2.石家莊鐵道大學土木工程學院，石家莊 050043)

巖土體的應變局部化現象通常表現為宏觀剪切帶，是導致邊坡、堤壩、基坑和擋土墻等巖土工程結構漸進破壞的根本原因[1-2]。巖土顆粒材料一般具有很強的各向異性特征[3]，在主應力方向旋轉的加載條件，將出現主應力與主塑性應變率方向不共軸的現象[4]。這導致傳統共軸理論無法準確預測應變局部化分叉的開始，非共軸理論可以較好地解決這一問題。早期，Rudnicki等[5]建立了具有角點結構屈服面的非共軸彈塑性本構模型，考慮了切向應變率的屈服效應。之后，Papamichos等[6]在此基礎上給出了預測剪切帶萌生和傾角的理論解。王興等[7-8]提出了一種改進的角點理論并將其應用到了砂土狀態相關剪脹模型中，該模型只在主應力方向改變的條件下產生非共軸塑性變形,克服了傳統角點模型的不足。錢建固等[9]推導了有限變形條件下土體變形分叉的非共軸理論。Qian等[10-11]將Rudnicki等[5]的非共軸模型推廣到了三維應力空間，考慮應力第三不變量的影響，在平面應變條件下對剪切帶進行了分析，更好地預測了剪切帶的萌生、傾角以及材料的剪脹性，隨后對平面應變狀態下土體的軟化特性進行了模擬。呂璽琳等[12]采用非共軸理論對真三軸狀態下砂土的分叉行為進行了分析，發現中主應力對分叉特性有顯著影響。Du等[13]建立了黏土的三維非共軸本構模型，很好的描述了單調剪切加載下土體非共軸塑性流動特點。Yang等[14-16]模擬了簡單剪切條件下非共軸特性對巖土顆粒材料的應力應變相應的影響，發現非共軸作用的發揮受硬化定律、流動法則及初始側向壓力等多種因素的影響，并將其非共軸模型應用到了一些結構問題的分析，如淺基礎等，證明如果不考慮非共軸特性的影響，將導致不安全的設計這一共識[17-20]。

研究表明非共軸特性不單影響應變局部化分叉的開始，其對剪切帶的后續發展演化也有影響。這是因為，當剪切帶形成后，其內部發生了劇烈的剪切變形，是一種類似于簡單剪切的變形，存在主應力方向旋轉的情況，進而也存在非共軸現象[21-22]。已經證實采用傳統共軸的彈塑性本構模型不能準確預測土體分叉的開始，建立非共軸本構模型是合適的。但是，基于經典連續體理論的非共軸本構模型，不包含任何內部長度參數，在模擬邊值問題中的應變局部化現象時往往存在網格依賴性而不能很好地描述剪切帶的寬度，更無法體現非共軸特性對剪切帶發展演化的影響。為此，現建立基于微極理論的非共軸彈塑性本構模型，并自主開發相應的用戶自定義單元(user-defined element,UEL)子程序，研究土體非共軸特性對剪切帶的萌生、發展和演化這一漸進失效全過程的影響，克服傳統本構模型的局限性，對進一步完善本土本構模型具有重要意義。

1 基于微極理論的非共軸本構模型

(1)

圖1 非共軸塑性應變率的方向Fig.1 Direction of the non-coaxial plastic strain rate

(2)

共軸塑性應變率與傳統塑性流動理論相同，用張量指標表示為

(3)

(4)

(5)

非共軸塑性應變率[6-7]可表示為

(6)

(7)

式中：hnc為非共軸塑性模量；sij為偏應力張量。

最終，非共軸塑性本構關系可以表示為

(8)

(9)

式中：hnc為非共軸塑性模量；Nijkl/hnc為非共軸項應變率的柔度矩陣；(De)ijkl為彈性本構張量。在經典連續體理論中，彈性本構關系表示為

(10)

(11)

微極理論中引入了物質點的微轉動自由度ω，其幾何方程除了位移和應變的關系外，還存在微轉動位移和微曲率的關系[22]為

εij=uj,i-eijkωk，κij=ωj,i

(12)

式(12)中：κij為由微轉動位移產生的微曲率張量。

同樣，微極理論中的物理方程除了應力與應變關系之外，還存在與偶應力與微曲率的關系，其彈性本構關系表示為

(13)

(14)

式中：μij為偶應力張量；a、b和c為材料參數；G為傳統剪切模量，G與a的乘積為Gc，稱為微極剪切模量；l為內部長度參數,l與b的乘積、l與c的乘積亦可分別表示為lt和lb,分別代表與扭轉和彎曲相關的內部長度尺度，在平面應變條件下參數c不起作用。

在應微極理論中，與非共軸相關的式(5)需要變為對應柯西應力和偶應力的兩部分，即

(15)

最后，在微極理論框架下，基于D-P(Drucker-Prager)屈服準則和非關聯流動法則，引入非共軸塑性應變率，形成了一個基于微極理論的非共軸彈塑性本構模型，對巖土材料應變局部化分叉和剪切帶演化行為進行模擬。

2 應變局部化分叉理論

應變局部化條件建立在弱不連續性分叉理論的基礎之上，認為聲學張量行列式為零意味著分叉的開始[5]，即

det[Q]=det[nDn]=0

(16)

式(16)中：Q為聲學張量；D為本構張量；n為剪切帶的外法線方向矢量，n=(cosθ,sinθ,0)，θ為剪切帶方向角。

在經典連續體理論中，對于平面應變問題，局部化條件可表示為

(17)

式(17)中：Qjk=niDijklnl，Qjk=n1D1jk1n1+n1D1jk2n2+n2D2jk1n1+n2D2jk2n2。

在微極連續體中，三維條件下，經推導Q擴展為一個6×6的矩陣，即

(18)

3 數值模擬

考慮到非共軸流動理論中，存在一個與屈服面相切的塑性應變率，采用完全隱式歐拉算法存在收斂困難的問題，因此采用自帶誤差控制的修正歐拉算法[23]。基于ABAQUS的UEL子程序二次開發，實現了以上本構模型的程序代碼。

3.1 程序驗證

為了驗證程序的有效性，對Leighton Buzzard砂[24]進行了簡單剪切試驗的數值模擬，并與Roscoe的實驗結果[4]進行了對比。簡單剪切試驗模擬如圖2所示，采用了一個20節點減縮單元(C3D20R)，下邊界固定約束，上邊界節點施加水平的位移荷載，豎直向施加垂直壓力134.45 kPa，側向壓力系數0.5，前后表面約束法向位移模擬平面應變加載條件。Leighton Buzzard砂的摩擦角38°，剪切模量100 MPa，泊松比0.3，采用基于D-P準則的理想彈塑性模型。

圖2 簡單剪切試驗Fig.2 Simple shear test

圖3為數值結果與Roscoe的實驗結果的對比，從應力應變響應曲線和大主應力方向、大主塑性應變率方向的對比中可以看出，兩者吻合較好，能夠較好地反映Leighton Buzzard砂在主應力旋轉的加載條件下的非共軸特性。

圖3 數值結果與試驗對比Fig.3 Comparison of the numerical results and the experiment results

3.2 平面應變壓縮試驗

為了研究非共軸特性對應變局部化的影響，對平面應變壓縮試驗進行了數值模擬。為了還原真實試驗條件，建立10 cm×20 cm×2 cm的三維實體模型，采用C3D20R單元剖分網格，上下邊界約束水平向自由度，施加大小相等方向相反的豎向位移荷載uy，前后邊界約束法向位移，實現平面應變加載環境，左右邊界施加圍壓100 kPa，如圖4所示。計算中調用用戶單元子程序(UEL子程序)。

圖4 幾何模型與邊界條件Fig.4 Geometric model and boundary condition

模型采用的材料參數如表1所示[25]。

表1 材料參數Table 1 Material parameters

為了進一步驗證程序的正確性，將Gc、hnc設為零，使程序退化為基于經典連續體的共軸本構模型，并將其與ABAQUS自帶的程序進行了對比驗證，結果如圖5所示，其中圖5(a)左側圖中的PEEQ代表ABAQUS自帶程序計算得到的等效塑性應變，右側圖中的UVARM42是用戶自定義程序得到的等效塑性應變。可見由等效塑性應變表征的剪切帶模式和承載力位移曲線均體現了較好的一致性。

圖5 退化為經典共軸本構模型的UEL子程序與ABAQUS自帶程序的對比Fig.5 Comparison of the results obtained by using an UEL and ABAQUS built-in program

3.2.1 應變局部化分叉預測

如第3節所述，聲學張量行列式為零意味著分叉的開始，圖6和圖7給出了分別基于微極連續體理論建立的共軸和非共軸本構模型所得到的一些典型區域的分叉情況。其中縱軸代表聲學張量行列式與彈性聲學張量行列式的比值，是一個歸一化的結果，橫軸為加載增量步。單元a、c和e位于剪切帶內部，單元b、d位于剪切帶兩側附近，單元f位于遠離剪切帶的位置。可以看出，在加載初期所有位置處歸一化的聲學張量行列式均為1，因為加載初期材料處于彈性階段；隨后縱坐標值出現了驟降，意味著加載進入塑性階段，且數值最終趨于一個定值。共軸模型中曲線縱坐標趨近于一個非零常數，而非共軸模型則趨近于零。原因是非共軸模型中考慮了切向應變率，其彈塑性本構張量中增加了非共軸項，導致聲學張量行列式最終趨于零而出現奇異性，這也是數值計算收斂困難的原因所在。同時觀察到，不同位置處，分叉的時刻不同，最早出現分叉是在剪切帶中心位置單元a，隨后是位于剪切帶內部的單元c和e，單元b和d對應曲線在加載進入塑性階段之后出現了回彈，意味著剪切帶出現后，其內部呈現繼續加載狀態，而外部則出現了卸載，單元f對應曲線一直保持為水平線，說明遠離剪切帶的位置材料一直處于彈性階段。

圖6 共軸模型模擬分叉Fig.6 Bifurcation in the coaxial model

圖7 非共軸模型模擬分叉Fig.7 Bifurcation in the non-coaxial model

圖8給出了由共軸模型和非共軸模型得到的結構承載力位移曲線，結果表明共軸模型得到的承載力峰值要高于非共軸模型，同時，分叉點均出現在峰值點之前，且共軸模型的分叉時刻早于非共軸模型，分叉時所對應的四種情況下剪切帶模式如圖9所示，隨著非共軸塑性模量hnc的降低，非共軸程度加強，承載力峰值進一步降低，分叉時刻更滯后，這與目前的普遍認識一致[12,19-20]。這是因為共軸模型無法考慮由于屈服面切向的塑性應變率所引起的屈服效應，導致共軸模型預測的結果要比非共軸模型預測結果偏危險。

圖8 共軸模型和非共軸模型模擬的分叉時刻對比Fig.8 Bifurcation in the coaxial and the non-coaxial models

圖9 分叉時對應的等效塑性應變云圖Fig.9 Distribution of the equivalent plastic strain when bifurcation occurs

3.2.2 剪切帶寬度預測

圖10～圖12分別為根據共軸塑性應變、非共軸塑性應變和總體塑性應變計算得到的等效塑性應變。可以觀察到，圖10中共軸等效塑性應變隨著非共軸程度的增強，其峰值逐漸降低，圖11中非共軸等效塑性應變則隨著非共軸程度的增強，其峰值逐漸增大，且在共軸模型中為零，但是注意到非共軸等塑性應變要比共軸等效塑性應變的值小一個數量級，最終導致圖12中根據總體塑性應變計算得到的等效塑性應變隨著非共軸程度的增強是逐漸減低的趨勢。

圖11 非共軸等效塑性應變云圖Fig.11Distribution of the non-coaxial part of the equivalent plastic strain

圖13、圖14給出了微極理論中一些物理量的分布云圖，其中圖13為微轉動位移，可以看出，隨著非共軸程度的增強，剪切帶內微轉動位移峰值逐漸降低，轉動變形趨于均勻。圖14為微曲率的分布圖，共軸模型中剪切帶兩側微曲率呈現大小相等方向相反的規律，隨著非共軸程度加強，其峰值也逐漸降低，且剪切帶兩側的微曲率出現不對稱性。這說明非共軸特性從某種程度上削弱了微極轉動效應。為了驗證基于微極理論非共軸本構模型在克服網格依賴性方面的有效性，對10×20、14×28、15×30以及16×32 4種網格密度下的剪切帶寬度進行了比較(圖15、圖16，由等效塑性應變表征的剪切帶)，可以看出隨著網格加密，經典連續體理論模擬得到的剪切帶出現了分支，寬度也隨之改變，而基于微極理論得到的剪切帶寬度幾乎不變。

圖13 微轉動位移分布圖Fig.13 The distribution of the micro rotation displacement in z direction

圖14 微曲率κzx分布圖Fig.14 The distribution of the micro curvature κzx

圖15 網格敏感性調查(微極理論)Fig.15 Investigation of the mesh sensitivity(micropolar theory)

圖16 網格敏感性調查(經典連續體理論)Fig.16 Investigation of the mesh sensitivity(classical continuum theory)

4 結論

本文建立了一個基于微極理論的非共軸彈塑性本構模型，并采用自帶誤差控制的修正的歐拉算法進行了數值實現。通過簡單剪切試驗數值結果與已有室內試驗結果對比，驗證了程序的正確性。重點研究了非共軸特性對應變局部化分叉預測和剪切帶寬度的影響。結果表明非共軸模型預測的局部化分叉點滯后于共軸模型的預測結果，與已有文獻結論一致，且隨非共軸程度的增強，局部化變形趨于均勻，塑性區的范圍擴大，剪切帶寬度有變寬的趨勢。等效塑性應變峰值呈現出隨非共軸程度的增強而減小的趨勢。可以看出本文建立的本構模型既可預測局部化分叉，又可體現非共軸特性對剪切帶寬度的影響。綜合比較下，基于微極理論的非共軸本構模型比傳統的基于經典連續體理論的非共軸模型在模擬巖土材料邊值問題中的應變局部化現象時更具優越性，可為深入探究巖土結構由應變局部化導致的漸進破壞的內在機制提供科學指導。