三維應(yīng)力狀態(tài)下砂的臨界狀態(tài)本構(gòu)模型分叉分析

甄文戰(zhàn)，孫德安，陳瑤瑤

（上海大學(xué) 土木工程系，上海 200072）

1 引 言

在三軸剪切試驗中，常觀察到密砂試樣整體承載能力在達(dá)到一個峰值荷載水平之后，隨著變形的增加而承載力降低。即：巖土材料在應(yīng)力達(dá)到強(qiáng)度極限后，隨變形的繼續(xù)增加，其強(qiáng)度值將迅速降低到一個較低的水平，這種由變形引起的材料性能劣化的現(xiàn)象稱為應(yīng)變軟化現(xiàn)象。

巖土材料從硬化到軟化階段的轉(zhuǎn)變一般認(rèn)為是材料變形模式分叉造成的。對于服從關(guān)聯(lián)流動法則的巖土材料，在材料峰值強(qiáng)度前材料一般是穩(wěn)定的，在峰值點(diǎn)處，材料切向剛度矩陣奇異，材料失穩(wěn)，變形模式開始出現(xiàn)分叉，通常認(rèn)為峰值點(diǎn)即為分叉失穩(wěn)點(diǎn)，從而分叉后引起應(yīng)力-應(yīng)變關(guān)系不再具有惟一性。但對于服從非關(guān)聯(lián)流動法則的巖土材料，由于材料的剛度矩陣不對稱，則分叉點(diǎn)有可能出現(xiàn)應(yīng)力-應(yīng)變關(guān)系的硬化階段。

分叉理論為研究服從非關(guān)聯(lián)流動法則的軟化材料失穩(wěn)提供了理論基礎(chǔ)。在分叉理論研究方面，Hill[1－2]最早用分叉理論建立了預(yù)測應(yīng)變局部化的理論框架。Rudnicki和Rice[3]，Rice[4]在Hill理論框架下繼續(xù)深入研究，并在此基礎(chǔ)上提出應(yīng)變局部化判別準(zhǔn)則的具體表達(dá)試，并成為了目前應(yīng)變局部化研究的主要理論基礎(chǔ)。此后，許多學(xué)者[5－6]仍致力于分叉理論的研究，但主要集中對平面應(yīng)變二維問題的分析，而對不同應(yīng)力路徑下分叉的三維特性研究還非常少。

本文基于姚仰平等[7]提出可以考慮平均應(yīng)力水平及孔隙比影響的砂臨界狀態(tài)三維彈塑性本構(gòu)模型，就服從非關(guān)聯(lián)流動法則的軟化型本構(gòu)模型在不同加載路徑下進(jìn)行了分叉理論分析。討論了不同應(yīng)力路徑下分叉特性，給出了分叉點(diǎn)對應(yīng)的應(yīng)力比、主應(yīng)變及剪切帶傾角隨應(yīng)力洛德角變化的趨勢。最后，編寫本構(gòu)模型材料子程序，實現(xiàn)了該模型與有限元 ABAQUS軟件的結(jié)合，通過對立方體密砂試樣的數(shù)值模擬，證明了分叉現(xiàn)象的存在，并捕捉到了分叉點(diǎn)的應(yīng)力狀態(tài)，驗證了理論解的正確性。

2 臨界狀態(tài)本構(gòu)模型

砂的變形特性受平均應(yīng)力水平及孔隙比的影響。初始狀態(tài)很松的砂呈硬化型應(yīng)力-應(yīng)變關(guān)系，而中密砂受剪時呈現(xiàn)先硬化后軟化的應(yīng)力-應(yīng)變關(guān)系。為了反映砂的上述變形特性，Yao等[7]提出了一個砂的臨界狀態(tài)本構(gòu)模型，該模型不僅可以反映松砂強(qiáng)度低、體縮大的變形特性，還可以反映中密砂及密砂的剪脹、應(yīng)力-應(yīng)變關(guān)系軟化的變形特性。該模型認(rèn)為，存在一條臨界狀態(tài)線(CSL)，斜率為λ，如圖 1所示。參考線(RCL)平行于臨界狀態(tài)線，該線以上的點(diǎn) A表示松砂狀態(tài)。砂在最密時的狀態(tài)線(DCL)對應(yīng)于加荷再卸荷曲線，斜率為κ。砂的臨界狀態(tài)本構(gòu)模型采用 SMP準(zhǔn)則，通過變換應(yīng)力[8]的方法實現(xiàn)了該模型三維化。該方法假定在三軸壓縮條件下變換應(yīng)力與原應(yīng)力相同，使得在原應(yīng)力空間下π平面上為曲邊三角形的SMP準(zhǔn)則，在變換應(yīng)力空間的π平面上為圓形，變換應(yīng)力σ~ij公式為

其中，

式中： δij為克羅內(nèi)克爾參數(shù)；sij為偏應(yīng)力張量；I1、I2、I3分別為第1、2及第3應(yīng)力不變量。

該模型屈服函數(shù)f及塑性勢函數(shù)g表達(dá)式為[7]

式中：上標(biāo)～表示為變換應(yīng)力空間中的值； λ為壓縮指數(shù)；κ為回彈指數(shù)；e0為初始孔隙比；為塑性體積應(yīng)變；分別為初始平均應(yīng)力、平均應(yīng)力；為應(yīng)力比；比容 v0= 1+e0；Mfmax、M分別為最大及臨界狀態(tài)時的應(yīng)力比； H~為硬化參數(shù)；為當(dāng)p =1 kPa時參考線、臨界狀態(tài)線及最密狀態(tài)線的比容；為p時參考線、應(yīng)力比η~、臨界狀態(tài)線及最密狀態(tài)線的比容。

圖1 比容－平均應(yīng)力平面中砂的密實狀態(tài)表示[7]Fig.1 Sand with different initial densities in v-lnp plane[7]

3 分叉理論分析

理論研究認(rèn)為：應(yīng)變軟化并不是發(fā)生分叉現(xiàn)象的惟一原因。Rudnicki和Rice證明了在采用非關(guān)聯(lián)流動法則的情況下，即使在應(yīng)變硬化時仍可能發(fā)生分叉現(xiàn)象[3]。而服從非關(guān)聯(lián)流動法則的軟化模型在不同應(yīng)力路徑下的分叉特性如何？下面就此問題進(jìn)行理論分析。

3.1 分叉理論

Rice等[4]認(rèn)為，速度場在剪切帶內(nèi)保持連續(xù)，而速度梯度通過剪切帶產(chǎn)生跳躍，但速度和速度梯度在平行于局部化變形帶方向仍保持均勻，而變形場在帶外保持均勻。如圖2所示，在S面上速度場保持連續(xù)，而速度梯度場產(chǎn)生跳躍。在這一假設(shè)條件下提出了經(jīng)典的分叉判別準(zhǔn)則：

圖2 三維應(yīng)力狀態(tài)下剪切帶示意圖Fig.2 Shear band under 3D stress state

3.2 分叉判別式

根據(jù)巖土材料的真三軸試驗現(xiàn)象[9]，假定剪切帶的法線方向垂直于中主應(yīng)力方向，即剪切帶沿中主應(yīng)力方向不發(fā)生變化，但中主應(yīng)力的大小影響分叉的產(chǎn)生條件。將n2=0代入行列式(10)中，即得：

砂的臨界狀態(tài)本構(gòu)模型本構(gòu)張量表達(dá)式為

其中：

式中：ν為泊松比。因屈服函數(shù)與塑性勢函數(shù)不同，即f≠g，所以該本構(gòu)模型為非關(guān)聯(lián)模型。

方程(11)有實數(shù)解，出現(xiàn)局部化分叉判別式為

3.3 三維應(yīng)力狀態(tài)的分叉條件

在等平均應(yīng)力(p=200 kPa)下，不同應(yīng)力洛德角下的條件為

當(dāng)θσ為25o、±15o、0o及±30o時，聯(lián)立式(13)～(15)，可求得在等p下初始孔隙比為0.68時，不同應(yīng)力洛德角（如圖3所示）對應(yīng)的局分叉理論解。模型材料參數(shù)取自文獻(xiàn)[7]，見表1。

3.4 三維應(yīng)力狀態(tài)下分叉結(jié)果

表2為不同應(yīng)力路徑下分叉理論解，并給出了分叉時對應(yīng)的應(yīng)力比、主應(yīng)變及剪切帶傾角值。

圖3 π平面上應(yīng)力洛德角Fig.3 Lode’s angle in π-plane

表1 砂的模型參數(shù)Table 1 Model parameters for sand

表2 分叉理論解Table 2 Theoretical solutions to bifurcation

圖4表示了分叉時(不分叉區(qū)域的上圍包線)對應(yīng)的應(yīng)力比或主應(yīng)變隨應(yīng)力洛德角的變化趨勢。由圖可知，分叉受應(yīng)力路徑的影響，即隨著應(yīng)力洛德角的變化表現(xiàn)出不同的分叉特性。羅德角在-25°～15°之間變化時均能產(chǎn)生分叉，且該分叉解隨洛德角變化規(guī)律與文獻(xiàn)[10]試驗結(jié)果基本一致。并且分叉時對應(yīng)的應(yīng)力比或主應(yīng)變均隨洛德角的增加先增加而后減小。圖4還分別給出了不同洛德角時不分叉區(qū)域范圍。即當(dāng)洛德角和對應(yīng)的應(yīng)力比或主應(yīng)變在不分叉區(qū)域內(nèi)就不會出現(xiàn)分叉現(xiàn)象。根據(jù)該圖可以對是否出現(xiàn)分叉進(jìn)行判別。

圖4 不同應(yīng)力路徑下不分叉區(qū)域Fig.4 Ranges of no bifurcation at different stress paths

圖5為剪切帶傾角隨洛德角變化情況，由圖可知，隨著洛德角的增加(-25°～0°)，剪切帶傾角增長較大，而后稍微減小，且分叉時剪切帶傾角(對應(yīng)剪切帶出現(xiàn)最早時刻值)均大于45°。

圖5 不同應(yīng)力路徑下剪切帶傾角變化曲線Fig.5 Relationship between inclination angle of shear band and Lode’s angle

3.5 孔隙比對分叉結(jié)果的影響

通過以上分析可認(rèn)為：是否出現(xiàn)局部化分叉受應(yīng)力路徑的影響。但同一應(yīng)力路徑下的分叉結(jié)果又受孔隙比的影響，即孔隙比的大小影響著分叉點(diǎn)在應(yīng)力-應(yīng)變關(guān)系上出現(xiàn)的位置。為了分析孔隙比對分叉的影響，在平面應(yīng)變-應(yīng)力路徑下分析不同孔隙比下的分叉和峰值應(yīng)力比。圖6表示了圍壓一定時的計算結(jié)果。由圖可知，分叉點(diǎn)對應(yīng)的應(yīng)力比(分叉應(yīng)力比)隨孔隙比的增加而減小，但分叉點(diǎn)對應(yīng)的分叉應(yīng)力比小于峰值應(yīng)力比，則說明分叉出現(xiàn)在應(yīng)力-應(yīng)變關(guān)系硬化階段。

圖6 應(yīng)力比隨初始孔隙比變化曲線Fig.6 Relationships between stress ratio and initial void ratio

4 真三軸數(shù)值模擬

經(jīng)典本構(gòu)關(guān)系理論建立在連續(xù)介質(zhì)力學(xué)基礎(chǔ)之上，雖具有嚴(yán)密的數(shù)學(xué)理論框架，但卻難以描述應(yīng)變局部化問題，在有限元分析中存在很大的局限性。在接近破壞條件出現(xiàn)失穩(wěn)分叉時，高應(yīng)變區(qū)呈現(xiàn)局部化，且該區(qū)域內(nèi)許多單元積分點(diǎn)處的聲張量成為奇異(剛度矩陣失去正定性，特征值出現(xiàn)負(fù)值)，擬靜力問題的偏微分控制方程喪失橢圓型，這使得所考慮的邊值問題變?yōu)椴B(tài)，其解出現(xiàn)非適定性[11]。

為了在不同應(yīng)力路徑下通過數(shù)值方法預(yù)測分叉現(xiàn)象，作者在在文獻(xiàn)[12－13]的基礎(chǔ)上，采用半隱式回映應(yīng)力更新算法，編寫了臨界狀態(tài)本構(gòu)模型的子程序，程序流程圖如圖7所示。把該模型嵌入到有限元程序 ABAQUS中，從而實現(xiàn)了在數(shù)值計算中得以捕捉分叉(剛度矩陣出現(xiàn)負(fù)特征值)時對應(yīng)的應(yīng)力狀態(tài)的方法。

圖7 子程序流程圖Fig.7 Flow chart of the subroutine

4.1 子程序單元測試

為了驗證子程序的精度，取一個單元進(jìn)行子程序特性測試。并與豐浦砂的三軸壓縮試驗結(jié)果進(jìn)行了比較，如圖8所示，初始孔隙比e0=0.68時的應(yīng)力-應(yīng)變的計算曲線與試驗結(jié)果基本一致，驗證了程序的正確性。

4.2 數(shù)值分析方案

為了通過數(shù)值方法對本構(gòu)模型在不同應(yīng)力路徑的分叉特性進(jìn)行預(yù)測，取一個各向同性均勻的立方體，邊長為10 cm。將試樣劃分為14×14×14=2 744個單元，有限元網(wǎng)格劃分如圖9所示，其中1方向為大主應(yīng)力方向。為了在ABAQUS中實現(xiàn)等p時不同應(yīng)力路徑下的數(shù)值模擬，采用多步試算的方法確定圍壓變化幅值曲線，通過圍壓設(shè)定荷載幅值曲線的方法，成功實現(xiàn)了等p時不同應(yīng)力路徑下的數(shù)值預(yù)測。數(shù)值預(yù)測采用的模型參數(shù)同上表1。

圖8 主應(yīng)力比-應(yīng)變關(guān)系Fig.8 Relationships between stress ratio and strain

圖9 有限元網(wǎng)格Fig.9 Meshes for finite element analysis

4.3 數(shù)值預(yù)測結(jié)果與理論解比較

圖10給出了不同應(yīng)力路徑下數(shù)值預(yù)測分叉點(diǎn)對應(yīng)的應(yīng)力-應(yīng)變狀態(tài)。數(shù)值計算時取負(fù)特征值(negative eigenvalue)出現(xiàn)對應(yīng)的應(yīng)力狀態(tài)為分叉點(diǎn)。其原因是從 ABAQUS計算方法上分析，負(fù)特征值表示在應(yīng)力-應(yīng)變關(guān)系曲線的硬化階段出現(xiàn)負(fù)斜率，與原應(yīng)力-應(yīng)變關(guān)系不同，控制偏微分方程由橢圓型轉(zhuǎn)化為雙曲線型，即分叉；從數(shù)學(xué)上分析，出現(xiàn)負(fù)特征值，剛度矩陣不正定，方程組存在多組解，從而可判定數(shù)值計算的分叉點(diǎn)。

由圖可知，理論上只要在該應(yīng)力路徑下存在分叉點(diǎn)時，通過數(shù)值方法也同樣捕捉到了分叉點(diǎn)的狀態(tài)，證明了分叉現(xiàn)象在數(shù)值計算中的存在。同時，通過圖中理論解及數(shù)值預(yù)測分叉點(diǎn)的對比，可以發(fā)現(xiàn)數(shù)值預(yù)測結(jié)果也與理論解基本一致(應(yīng)力比最大誤差為θσ=0°對應(yīng)的值11.4%)。

圖10 分叉點(diǎn)的理論和數(shù)值預(yù)測Fig.10 Theoretical and numerical predictions of bifurcation

5 結(jié) 論

（1）理論分析表明：分叉現(xiàn)象強(qiáng)烈依賴于應(yīng)力路徑，即應(yīng)力洛德角在-25°～15°之間變化時，均會在應(yīng)力-應(yīng)變關(guān)系硬化階段出現(xiàn)分叉現(xiàn)象，而其他應(yīng)力路徑下不會產(chǎn)生分叉。但同一應(yīng)力路徑下的分叉結(jié)果受孔隙比的影響。且分叉時對應(yīng)的應(yīng)力比、主應(yīng)變及剪切帶傾角都隨洛德角而變化，其值均先增加而后減小。

（2）通過數(shù)值方法預(yù)測到了分叉點(diǎn)對應(yīng)的應(yīng)力狀態(tài)，數(shù)值預(yù)測值和理論解的結(jié)果基本一致。

[1]Hill R. A general theory of uniqueness and stability in elastic-plastic solids[J]. Journal of the Mechanics and Physics of Solids, 1958, 5: 236－249.

[2]HILL R A, HUTCHINSON J W. Bifurcation phenomena in the plane tension test[J]. Journal of the Mechanics and Physics of Solids, 1975, 23: 239－264.

[3]RUDNICKI J W, RICE J R. Conditions for the localization of deformation in pressure-sensitive dilatant materials[J]. Journal of the Mechanics and Physics of Solids, 1975, 23: 371－394.

[4]RICE J R. The localization of plastic deformation[C]//Proceedings of 14th International Conference on Theoretical and Applied Mechanics. [S. l.]: [s. n.], 1976:207－220.

[5]VARDOULAKIS I, GOLDSCHEIDER M, GUDEHUS G.Formation of shear bands in sand bodies as a bifurcation problem[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1978, 2(2): 99－128.

[6]OTTOSEN N S, RUNESSON K. Properties of discontinuous bifurcation solutions in elasto-plasticity[J].International Journal of Solids and Structures, 1991,27: 231－254.

[7]YAO Y P, SUN D A, LUO T. A critical state model for sands dependent on stress and density[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2004, 28(4): 323－337.

[8]MATSUOKA H, YAO Y P, SUN D A. The Cam-clay model revised by the SMP criterion[J]. Soils and Foundations, 1999, 39(1): 81－95.

[9]松岡元, 中井照夫. Closure to discussion on “Stressdeformation and strength characteristics of soil under three-different principal stress[M]. 日本: [s. n.], 1976.

[10]WANG Q, LADE P V. Shear banding in true triaxial tests and its effect on failure in sand[J]. Journal of Engineering Mechanics, ASCE, 2001, 127(8): 754－761.

[11]BAZANT Z P, BELYTSCHKO T, CHANG T P.Continuum theory for strain softening[J]. Journal of Engineering Mechanics, ASCE, 1984, 110: 1666－1692.

[12]孫德安, 甄文戰(zhàn), 黃文雄. 三維彈塑性模型在路堤軟基固結(jié)分析中應(yīng)用[J]. 巖土力學(xué), 2009, 30(3): 669－674.SUN De-an, ZHEN Wen-zhan, HUANG Wen-xiong.Application of 3D elastoplastic model to analysis of consolidation behavior of embankment on soft soils[J].Rock and Soil Mechanics, 2009, 30(3): 669－674.

[13]孫德安, 甄文戰(zhàn). 不同應(yīng)力路徑下剪切帶的數(shù)值模擬[J]. 巖土力學(xué), 2010, 31(7): 2253－2258.SUN De-an, ZHEN Wen-zhan. Numerical simulations of shear bands along different stress paths[J]. Rock and Soil Mechanics, 2010, 31(7): 2253－2258.