非共軸特性對巖土材料應(yīng)變局部化行為影響的數(shù)值研究

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

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

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

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

1 基于微極理論的非共軸本構(gòu)模型

(1)

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

(2)

共軸塑性應(yīng)變率與傳統(tǒng)塑性流動理論相同，用張量指標(biāo)表示為

(3)

(4)

(5)

非共軸塑性應(yīng)變率[6-7]可表示為

(6)

(7)

式中：hnc為非共軸塑性模量；sij為偏應(yīng)力張量。

最終，非共軸塑性本構(gòu)關(guān)系可以表示為

(8)

(9)

式中：hnc為非共軸塑性模量；Nijkl/hnc為非共軸項應(yīng)變率的柔度矩陣；(De)ijkl為彈性本構(gòu)張量。在經(jīng)典連續(xù)體理論中，彈性本構(gòu)關(guān)系表示為

(10)

(11)

微極理論中引入了物質(zhì)點的微轉(zhuǎn)動自由度ω，其幾何方程除了位移和應(yīng)變的關(guān)系外，還存在微轉(zhuǎn)動位移和微曲率的關(guān)系[22]為

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

(12)

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

同樣，微極理論中的物理方程除了應(yīng)力與應(yīng)變關(guān)系之外，還存在與偶應(yīng)力與微曲率的關(guān)系，其彈性本構(gòu)關(guān)系表示為

(13)

(14)

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

在應(yīng)微極理論中，與非共軸相關(guān)的式(5)需要變?yōu)閷?yīng)柯西應(yīng)力和偶應(yīng)力的兩部分，即

(15)

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

2 應(yīng)變局部化分叉理論

應(yīng)變局部化條件建立在弱不連續(xù)性分叉理論的基礎(chǔ)之上，認(rèn)為聲學(xué)張量行列式為零意味著分叉的開始[5]，即

det[Q]=det[nDn]=0

(16)

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

在經(jīng)典連續(xù)體理論中，對于平面應(yīng)變問題，局部化條件可表示為

(17)

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

在微極連續(xù)體中，三維條件下，經(jīng)推導(dǎo)Q擴展為一個6×6的矩陣，即

(18)

3 數(shù)值模擬

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

3.1 程序驗證

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

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

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

圖3 數(shù)值結(jié)果與試驗對比Fig.3 Comparison of the numerical results and the experiment results

3.2 平面應(yīng)變壓縮試驗

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

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

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

表1 材料參數(shù)Table 1 Material parameters

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

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

3.2.1 應(yīng)變局部化分叉預(yù)測

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

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

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

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

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

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

3.2.2 剪切帶寬度預(yù)測

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

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

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

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

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

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

圖16 網(wǎng)格敏感性調(diào)查(經(jīng)典連續(xù)體理論)Fig.16 Investigation of the mesh sensitivity(classical continuum theory)

4 結(jié)論

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