三維應力狀態(tài)下砂的臨界狀態(tài)本構模型分叉分析

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

（上海大學 土木工程系，上海 200072）

1 引 言

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

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

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

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

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

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

其中，

式中： δij為克羅內克爾參數；sij為偏應力張量；I1、I2、I3分別為第1、2及第3應力不變量。

該模型屈服函數f及塑性勢函數g表達式為[7]

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

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

3 分叉理論分析

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

3.1 分叉理論

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

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

3.2 分叉判別式

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

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

其中：

式中：ν為泊松比。因屈服函數與塑性勢函數不同，即f≠g，所以該本構模型為非關聯模型。

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

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

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

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

3.4 三維應力狀態(tài)下分叉結果

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

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

表1 砂的模型參數Table 1 Model parameters for sand

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

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

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

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

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

3.5 孔隙比對分叉結果的影響

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

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

4 真三軸數值模擬

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

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

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

4.1 子程序單元測試

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

4.2 數值分析方案

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

圖8 主應力比-應變關系Fig.8 Relationships between stress ratio and strain

圖9 有限元網格Fig.9 Meshes for finite element analysis

4.3 數值預測結果與理論解比較

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

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

圖10 分叉點的理論和數值預測Fig.10 Theoretical and numerical predictions of bifurcation

5 結 論

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

（2）通過數值方法預測到了分叉點對應的應力狀態(tài)，數值預測值和理論解的結果基本一致。

[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), 黃文雄. 三維彈塑性模型在路堤軟基固結分析中應用[J]. 巖土力學, 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). 不同應力路徑下剪切帶的數值模擬[J]. 巖土力學, 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.