土體邊界面模型在ABAQUS軟件中的研發(fā)與驗(yàn)證

欽亞洲，孫 鈞

（1.同濟(jì)大學(xué) 巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海200092；2.南通大學(xué) 交通學(xué)院，江蘇 南通226019；3.杭州豐強(qiáng)土建工程研究院，浙江 杭州310006）

天然土體一般都存在明顯的各向異性.Casagrande［1］認(rèn)為土體各向異性可分為固有各向異性和應(yīng)力誘發(fā)各向異性.前者是指土顆粒在沉積過程中，顆粒長軸總存在垂直于最大主應(yīng)力排列的趨勢；而后者則由于應(yīng)力作用后土體變形而產(chǎn)生.一般的方法是將土體各向異性分為初始各向異性和應(yīng)力誘發(fā)各向異性，而初始各向異性包含了固有各向異性和由于偏壓固結(jié)產(chǎn)生的誘發(fā)各向異性兩個部分.對各向異性的定量描述也存在兩種方法，一種是基于微觀結(jié)構(gòu)，采用描述土顆粒空間排列特征的組構(gòu)張量［2－3］；另一種是由宏觀不等向固結(jié)現(xiàn)象，采用描述宏觀土體各向異性特性的各向異性張量［4］，這兩者之間可以相互轉(zhuǎn)換.

各向異性會對土體剪切強(qiáng)度、應(yīng)力應(yīng)變行為，以及屈服面傾向等產(chǎn)生影響.一些研究結(jié)果表明，對于正常固結(jié)土，考慮土體各向異性后，其力學(xué)行為將趨于超固結(jié)土.

目前有多種能反映土體各向異性的本構(gòu)模型，它們間的區(qū)別在于屈服面選擇、采用的硬化形式等方面.多屈服面本構(gòu)模型一般采用運(yùn)動硬化來反映各向異性，單屈服面本構(gòu)模型則一般采用旋轉(zhuǎn)硬化來反映.

Wheeler［5］針對芬蘭黏土，提出了一個能考慮土體初始各向異性及應(yīng)力誘發(fā)各向異性的彈塑性本構(gòu)模型：S－CLAY1.此前的各向異性模型都采用塑性體應(yīng)變εpv作為各向異性張量αij的硬化因子，而SCLAY1模型則同時采用塑性體應(yīng)變εpv及塑性偏應(yīng)變εps作為旋轉(zhuǎn)硬化因子，綜合考慮兩種塑性變形對各向異性張量的作用，所以更為合理.

在S－CLAY1模型的基礎(chǔ)上，本文通過引入邊界面理論，在S－CLAY1模型的基礎(chǔ)上將其拓展為能同時模擬正常固結(jié)土和中等超固結(jié)土的本構(gòu)模型，在ABAQUS軟件中編程開發(fā)實(shí)現(xiàn)，并對所建議模型進(jìn)行了模擬驗(yàn)證.

1 模型構(gòu)建

彈塑性模型可以通過引入邊界面理論，從而改造為一個邊界面模型，分述如下.

1.1 邊界面方程

在3維應(yīng)力空間，邊界面方程表達(dá)式如下：

1.2 硬化法則

硬化法則包括兩個部分：等向硬化法則和旋轉(zhuǎn)硬化法則.

等向硬化法則由前期固結(jié)壓力p0演化來反映，與修正的劍橋模型一致，以塑性體應(yīng)變εpv作為硬化因子

式中：λ為e—lnp空間正常固結(jié)線的斜率；κ為e—lnp空間回彈曲線的斜率；v0＝1＋e0；e0為初始孔隙比.

各向異性張量αij的硬化法則為旋轉(zhuǎn)硬化法則，以塑性體應(yīng)變εpv及塑性偏應(yīng)變εps作為硬化因子

式中：μ，β為常量參數(shù)，μ控制αij變化量的大小，β為控制塑性體應(yīng)變率與塑性偏應(yīng)變率對αij變化影響比例的因子；〈〉為Macaulay符號，含義為當(dāng)時，當(dāng)時為塑性偏應(yīng)變分量增量.

Wheeler給出了參數(shù)μ，β的確定方法和范圍.根據(jù)Wheeler的建議，μ一般取10／λ～15／λ，β／M取值在0.5～1.0.

1.3 映射法則

采用徑向映射法則［6－7］，如圖1所示.

圖1 邊界面及徑向映射示意圖Fig.1 Schematic illustration of bounding surface and mapping rule

徑向映射為從投影中心引出直線，經(jīng)過當(dāng)前應(yīng)力點(diǎn)（p，q），與邊界面交于點(diǎn)，此點(diǎn)定義為當(dāng)前應(yīng)力點(diǎn)的像點(diǎn).

圖1中，l為當(dāng)前應(yīng)力點(diǎn)與像點(diǎn)之間的距離，r0為投影中心到像點(diǎn)之間的距離，定義比例因子b為

式中b必須滿足b≥1.

則像點(diǎn)處應(yīng)力與當(dāng)前應(yīng)力存在有以下關(guān)系：

1.4 塑性模量

當(dāng)前應(yīng)力點(diǎn)的塑性模量Kp可由及當(dāng)前應(yīng)力點(diǎn)至投影中心距離插值求出.Kp插值函數(shù)的選擇，影響模型的準(zhǔn)確性.

本模型插值函數(shù)表達(dá)式為

式中：pa代表一個大氣壓，取100kPa；Rij為塑性流動方向；h代表加載純彈性域的大小.對于黏土而言，幾乎不存在一個加載純彈性域，所以取h＝0；W為常量參數(shù)，一般可取W＝2.0.

2 數(shù)值計算

2.1 彈性應(yīng)力應(yīng)變關(guān)系

彈性應(yīng)力應(yīng)變關(guān)系與修正劍橋模型一致，為

式中：K為體積模量，隨球應(yīng)力變化；e·eij為彈性偏應(yīng)變增量.

式中：G為剪切模量；ν為泊松比.

2.2 塑性應(yīng)力更新算法

塑性應(yīng)力更新算法采用圖形返回算法，這是一種完全隱式積分算法，無條件穩(wěn)定.它強(qiáng)化在時間步結(jié)束時的一致性，即fn＋1＝0，避免應(yīng)力點(diǎn)離開屈服面而漂移，被證明是強(qiáng)健的和精確的［8］.算法示意圖如圖2.

圖2 圖形返回算法示意圖Fig.2 Schematic illustration of return mapping algorithm

算法分兩步，第一步為彈性預(yù)測，第二步為塑性修正.

在彈性預(yù)測階段，首先假定第n＋1步應(yīng)變增量全部為彈性增量，得到彈性試探應(yīng)力Δε，帶入邊界面方程中判斷：若F＜0，只發(fā)生彈性變形，直接更新應(yīng)力；若F≥0，則發(fā)生塑性變形，此時應(yīng)力點(diǎn)可能位于屈服面之外，需要對應(yīng)變增量、內(nèi)變量、塑性加載因子等進(jìn)行迭代修正，使得應(yīng)力點(diǎn)沿塑性流動方向逐步拉回至新的屈服面，這個過程如圖2所示.

Manzari［9］和費(fèi)康［10］等采用隱式算法實(shí)現(xiàn)了對各向同性邊界面模型的研發(fā).

3 模型驗(yàn)證

借助于大型商業(yè)軟件ABAQUS 所提供的UMAT（user－defined material）子程序接口完成了模型開發(fā)，編程時注意以下幾點(diǎn)：

①ABAQUS中以拉應(yīng)力為正；

②ABAQUS中剪應(yīng)變?yōu)楣こ碳魬?yīng)變；

③ABAQUS中應(yīng)力應(yīng)變分量的順序?yàn)棣?1，σ22，

Stipho［11］及Ling［12］對高嶺土進(jìn)行了一系列三軸不排水剪切試驗(yàn)，可用來驗(yàn)證本模型.此次驗(yàn)證模擬了k0＝0.8的初始偏壓固結(jié)情況，選用OCR＝1的正常固結(jié)土和OCR＝4的超固結(jié)土試樣.模擬采用位移加載，直至軸向應(yīng)變量達(dá)17%后停止加載.

3.1 模型參數(shù)選取

對于修正劍橋模型參數(shù)λ，κ，ν，M，p0按Stipho試驗(yàn)值選取；若土體處于k0正常固結(jié)狀態(tài)，參數(shù)β可由下式確定：

式中，ηk0為土體k0狀態(tài)時應(yīng)力比

因此取β＝0.6；參數(shù)μ的選取，Wheeler只給出大致取值區(qū)間.文獻(xiàn)［13－14］給出一種確定參數(shù)μ的方法，若對k0正常固結(jié)土樣施加大小為2～3 倍前期固結(jié)壓力的等向圍壓后，土體宏觀各向異性基本消失，從而由下式確定參數(shù)μ：

式中，αk0為土體正常k0狀態(tài)時各向異性量，且

此外，在模型其他參數(shù)已確定的情況下，也可采用擬合試驗(yàn)數(shù)據(jù)方法確定參數(shù)μ的值.本文模型參數(shù)選取見表1.

表1 高嶺土模型參數(shù)Tab.1 Model parameters for Kaolin clay

Wheeler文中已給出參數(shù)β隨α／M，η／M的變化規(guī)律，所以本文只考察參數(shù)μ對土體工程行為的影響.在保持其他參數(shù)取值不變的情況下，分別取μ＝0，15，30，60對OCR＝1，k0＝0.8的初始偏壓固結(jié)土不排水剪切試驗(yàn)進(jìn)行了模擬，結(jié)果如圖3和圖4.

由圖3和圖4可見，隨著μ值的增大，土體的剪切強(qiáng)度增大，并且在剪切過程中，土體應(yīng)力路徑外擴(kuò)，這都反映了旋轉(zhuǎn)硬化的影響.

3.2 試驗(yàn)?zāi)M結(jié)果

若將本模型參數(shù)μ設(shè)置為零，則模型退化為各向同性邊界面模型.本文分別采用各向同性邊界面模型和各向異性邊界面模型對試驗(yàn)進(jìn)行了模擬，模擬結(jié)果如圖5，6，7所示.

由模擬結(jié)果可見，總體上各向異性模型模擬的結(jié)果與試驗(yàn)值較吻合.根據(jù)土體臨界狀態(tài)理論，當(dāng)土體應(yīng)力狀態(tài)到達(dá)破壞線（M線）時，此時土體塑性體應(yīng)變εpv不再發(fā)生變化，也就是土體不再發(fā)生等向硬化.所以各向同性模型在達(dá)到臨界狀態(tài)后，強(qiáng)度曲線會出現(xiàn)水平段，如圖5虛線所示.這使得各向同性模型低估了土的抗剪強(qiáng)度，并使得孔壓也與試驗(yàn)值存在些偏差.各向異性模型由于采用了包含塑性剪應(yīng)變εps的旋轉(zhuǎn)硬化法則，所以當(dāng)土體應(yīng)力達(dá)到臨界狀態(tài)后，由于εps的存在，土體仍然會發(fā)生旋轉(zhuǎn)硬化，強(qiáng)度曲線表現(xiàn)為沿M線繼續(xù)向上一小段，如圖6實(shí)線所示.

對于超固結(jié)土的不排水剪切試驗(yàn)的模擬，兩種模型都存在一些偏差.

4 結(jié)論

在Wheeler彈塑性模型的基礎(chǔ)上，結(jié)合邊界面理論，將彈塑性模型拓展為各向異性邊界面模型.模型采用旋轉(zhuǎn)硬化法則來反映土體的初始各向異性及隨后的應(yīng)力誘發(fā)各向異性.利用軟件ABAQUS提供的UMAT 子程序接口，采用隱式積分算法——圖形返回算法編程實(shí)現(xiàn).通過對k0＝0.8的初始偏壓固結(jié)土體三軸不排水剪切試驗(yàn)的模擬，并與試驗(yàn)結(jié)果驗(yàn)證表明，模型能夠合理描述正常固結(jié)土及中等超固結(jié)土的應(yīng)力應(yīng)變行為、孔壓曲線及應(yīng)力路徑等.此外，相比于其他各向異性模型，本模型參數(shù)量少，除修正劍橋模型參數(shù)外，只另增加3個參數(shù)，且參數(shù)易于確定.

［1］ Casagrande A，Carrillo N.Shear failure of anisotropic materials［J］.Journal of Boston Society of Civil Engineers，1944，31（4）：74.

［2］ Anandarajah A，Dafalias Y F.Bounding surface plasticity.Ⅲ：application to anisotropic cohesive soils［J］.Journal of Engineering Mechanics，1986，112（12）：1292.

［3］ Liang R Y，Ma F G.Anisotropic plasticity model for undrained cyclic behavior of clays ［J］.Journal of Geotechnical Engineering，1992，ASCE，118（2）：227.

［4］ Whittle A J. Evaluation of a constitutive model for overconsolidated clays［J］.Geotechnique，1993，43（2）：289.

［5］ Wheeler S J，Naatanen A，Karstunen，et al.An anisotropic elastoplastic model for soft clays［J］.Canadian Geotechnical Journal，2003，40：403.

［6］ Dafalias Y F.Bounding surface plasticity.Ⅰ：mathematical foundation and hypoplasticity［J］.Journal of Engineering Mechanics，1986，112（9）：966.

［7］ Dafalias Y F，Herrmann L R.Bounding surface plasticity.Ⅱ：application to isotropic cohesive soils ［J］.Journal of Engineering Mechanics，1986，112（12）：1263.

［8］ Simo J C，Taylor R L.Consistent tangent operators for rateindependent elastoplasticity［J］.Computer Methods in Applied Mechanics and Engineering，1985，48（3），101.

［9］ Manzari M T，Nour M A.On implicit integration of bounding surface plasticity models［J］.Computers ＆Structures，1997，3（3）：385－395.

［10］ 費(fèi)康，劉漢龍.邊界面模型在ABAQUS的開發(fā)應(yīng)用［J］.解放軍理工大學(xué)學(xué)報：自然科學(xué)版，2009，10（5）：447.FEI Kang，LIU Hanlong.Implementation and application of bounding surface model in ABAQUS［J］.Journal of PLA University of Science and Technology： Natural Science Edition，2009，10（5）：447.

［11］ Stipho A S A.Experimental and theoretical investigation of the behavior of anisotropically consolidated kaolin ［D］.Cardiff：Cardiff University，1978.

［12］ Ling H I，Yue D Y，Kaliakin V N，et al.Anisotropic elastoplastic bounding surface model for cohesive soils［J］.Journal of Engineering Mechanics，2002，128（7）：748.

［13］ Leoni M，Karstunen M，Vermeer P A.Anisotropic creep model for soft soils［J］.Geotechnique，2008，58（3）：215.

［14］ Yin Z Y，Chang C S，Karstuene，et al.An anisotropic elasticviscoplastic model for soft clays［J］.International Journal of Solids and Structures，2010，47（5）：665.