含應(yīng)變軟化本構(gòu)關(guān)系的巖-土接觸元直剪試驗(yàn)數(shù)值模擬

鄭立寧，康景文，謝 強(qiáng)，陳 云，李可一，賈 鵬

（1.中建地下空間有限公司，成都 610081；2.中國(guó)建筑西南勘察設(shè)計(jì)研究院有限公司，成都 610052；3.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，成都 610031）

1 引言

直剪試驗(yàn)是獲取巖-土接觸界面剪切強(qiáng)度參數(shù)的主要途徑。大量的實(shí)踐表明[1]，在直剪過程中接觸面內(nèi)部受力并非呈理想特征，但由于目前儀器水平限制，無法對(duì)剪切過程中接觸面內(nèi)部的應(yīng)力-應(yīng)變特征進(jìn)行準(zhǔn)確測(cè)試。數(shù)值技術(shù)的發(fā)展為該問題提供了解決途徑。接觸元數(shù)值解法中的本構(gòu)關(guān)系可采用雙曲線模型、剛塑性模型、非線性彈性模型及理想彈塑性模型[2–4]。實(shí)踐表明，一般巖-土接觸界面均具有明顯的應(yīng)變軟化特征，王恭先等[5]試驗(yàn)得到某剪切帶剪切位移達(dá)到5 mm 左右，強(qiáng)度衰減至穩(wěn)定，另一種剪切位移則需達(dá)到10～20 mm；任光明等[6]研究認(rèn)為，某剪切帶位移僅3 mm 強(qiáng)度即進(jìn)入殘余狀態(tài)。但由于應(yīng)變軟化接觸本構(gòu)關(guān)系數(shù)值方法實(shí)現(xiàn)較為困難，在直剪過程中考慮接觸元應(yīng)變軟化特征的研究一直較為少見。

FLAC 軟件提供了基于FISH 語言修改的具有鏈表結(jié)構(gòu)的各種對(duì)象單元類（如接觸單元鏈表int.fin），為實(shí)現(xiàn)各種自定義本構(gòu)模型提供了空間[7]。本文通過對(duì)接觸單元鏈表數(shù)據(jù)的訪問，完善具有應(yīng)變軟化接觸本構(gòu)數(shù)值模型，并對(duì)含應(yīng)變軟化接觸本構(gòu)數(shù)值模型的直剪試驗(yàn)進(jìn)行模擬，深入研究接觸面直剪漸進(jìn)性破壞過程以及剪切接觸面內(nèi)部應(yīng)力變化特征。另外，系統(tǒng)總結(jié)了接觸單元剪切強(qiáng)度參數(shù)演化對(duì)剪切過程力學(xué)行為特征的影響，為深入揭示接觸面直剪內(nèi)部機(jī)制提供了理論解釋。

2 應(yīng)變軟化接觸本構(gòu)數(shù)值實(shí)現(xiàn)

2.1 FLAC 接觸單元數(shù)據(jù)的存儲(chǔ)及訪問

在FLAC 中，所有重要信息（網(wǎng)格坐標(biāo)、材料參數(shù)等）均采用內(nèi)存堆棧技術(shù)進(jìn)行存貯，即調(diào)用FISH 函數(shù)的參數(shù)地址存貯函數(shù)imem（i）和參數(shù)存貯函數(shù)fmem（i+j）（其中i為地址指針，j為地址指針偏移量）。一個(gè)地址指向?qū)?yīng)的位置，用于存放數(shù)據(jù)；指針變量為整數(shù)變量，而數(shù)據(jù)變量為浮點(diǎn)數(shù)變量[7–8]。每個(gè)接觸單元包括若干接觸節(jié)點(diǎn)，單元或節(jié)點(diǎn)的其他參數(shù)依靠地址指針偏移量的變化進(jìn)行存儲(chǔ)或訪問。

2.2 FLAC 中含應(yīng)變軟化本構(gòu)關(guān)系接觸單元的實(shí)現(xiàn)

本文建立的應(yīng)變軟化接觸本構(gòu)模型以三線型應(yīng)變軟化曲線為基礎(chǔ)[9]。當(dāng)接觸單元所受剪應(yīng)力較小時(shí)，單元僅產(chǎn)生線彈性變形，彈性變形階段材料強(qiáng)度參數(shù)維持峰值（φp、cp）不變。當(dāng)剪應(yīng)力達(dá)到最大剪應(yīng)力值τp（τp按摩爾-庫侖準(zhǔn)則計(jì)算）后，模型進(jìn)入線性軟化階段，其強(qiáng)度參數(shù)（φ、c）隨塑性剪切位移L 的增加線性衰減；當(dāng)塑性剪切位移達(dá)到Lp時(shí)，剪應(yīng)力達(dá)到最小值τr，接觸單元強(qiáng)度參數(shù)降至殘余強(qiáng)度（φr、cr）且保持不變，變化過程見圖1、2。

在FLAC 程序各計(jì)算時(shí)步中，調(diào)用FISH 函數(shù)的參數(shù)地址存貯函數(shù)imem（i）和參數(shù)存貯函數(shù)fmem（i+j），提取各接觸單元的最大剪應(yīng)力及接觸節(jié)點(diǎn)對(duì)的相對(duì)剪切應(yīng)變，對(duì)比預(yù)先確定的應(yīng)變軟化曲線及剪切強(qiáng)度參數(shù)衰減折線（見圖1、2），如接觸單元接觸節(jié)點(diǎn)進(jìn)入塑性剪切階段，則進(jìn)行相應(yīng)的強(qiáng)度參數(shù)衰減，由此實(shí)現(xiàn)接觸單元的應(yīng)變軟化接觸本構(gòu)模型[7]。

圖1 3 線型應(yīng)變軟化曲線Fig.1 Simplified polyline type of strain softening curve

圖2 剪切強(qiáng)度參數(shù)衰減折線Fig.2 Simplified polyline type for shear strength parameters attenuation

3 直剪試驗(yàn)數(shù)值模擬

3.1 模型的建立

參照直剪試驗(yàn)，模型長(zhǎng)為6 cm，高為2 cm，中間為接觸單元，共20個(gè)（編號(hào)從左至右為1～20），每個(gè)接觸單元長(zhǎng)為0.3 cm；接觸單元上、下兩側(cè)為實(shí)體單元，各200個(gè)。模型上表面受均勻豎向載荷（100 kPa），在左側(cè)上部表面施加速度為v=0.002 mm/時(shí)步的水平荷載，模型如圖3 所示。

圖3 直剪試驗(yàn)數(shù)值模型Fig.3 Element model of direct shear test

結(jié)合常見巖土體特征，模型上下塊體采用摩爾-庫侖模型，材料參數(shù)見表1。

表1 塊體模型參數(shù)Table 1 Parameters of solid models

表2 接觸單元模型參數(shù)Table 2 Parameters of contact elements models

模型接觸單元考慮兩大類6 種本構(gòu)模型分別計(jì)算，其一為強(qiáng)度不變型（傳統(tǒng)摩爾-庫侖模型），不同剪切強(qiáng)度值的3 種（B1、B2、B3）；其二為應(yīng)變軟化型3 種，包括摩擦角衰減型（F1）、黏聚力衰減型（C1）、摩擦角與黏聚力綜合衰減型（FC1），具體參數(shù)見表2。

3.2 剪應(yīng)力-位移關(guān)系曲線模擬結(jié)果

計(jì)算得到的不同接觸單元模型如圖4 所示。

圖4 不同模型剪應(yīng)力-位移關(guān)系曲線Fig.4 Shear stress-displacement curves of different models

從圖可知，當(dāng)模型接觸單元為強(qiáng)度不變型時(shí)，計(jì)算得到的剪切位移-剪應(yīng)力曲線為峰值后穩(wěn)定型；當(dāng)模型接觸單元為應(yīng)變軟化型時(shí)，剪切位移-剪應(yīng)力曲線為峰值后線性衰減型，不同的應(yīng)變軟化參數(shù)，衰減的變化斜率及殘余值均不一致。

3.3 直剪破壞過程及內(nèi)部剪應(yīng)力變化特征

以摩擦角與黏聚力綜合軟化模型（FC1）為例，分析直剪漸進(jìn)性破壞過程及接觸面內(nèi)部剪應(yīng)力的變化特征，不同計(jì)算時(shí)步模型的剪應(yīng)力分布云圖見圖5，不同計(jì)算時(shí)步各接觸單元剪應(yīng)力-變化曲線特征見圖6。

由圖5、6 可知，接觸單元的剪應(yīng)力值在剪切作用初期，受單元切向剛度控制產(chǎn)生彈性增加達(dá)到極值，極值過后進(jìn)入應(yīng)變軟化階段，剪應(yīng)力產(chǎn)生衰減直至穩(wěn)定。整個(gè)剪切過程具有漸進(jìn)性破壞特征，剪切面的剪應(yīng)力集中區(qū)從推力作用區(qū)開始逐步向模型中后部轉(zhuǎn)移，最終在模型后部形成剪應(yīng)力最大區(qū)域，且由于法向應(yīng)力的推后集中及側(cè)推力的疊加作用，最大剪應(yīng)力值約為140 kPa。隨著剪切作用的持續(xù)發(fā)展，所有接觸單元均進(jìn)入應(yīng)變軟化階段，內(nèi)部剪應(yīng)力不斷降低且向整體均布過渡。

不同計(jì)算時(shí)步各接觸單元剪切強(qiáng)度參數(shù)變化特征曲線見圖7。

圖7 剪切強(qiáng)度參數(shù)變化曲線Fig 7 The variation curves of shear strength parameters

由圖可見，各接觸單元自身的強(qiáng)度參數(shù)值受控于塑性剪切變形的增加而逐步衰減，直至殘余值。在整個(gè)剪切過程中，靠近推力作用區(qū)的接觸單元先期產(chǎn)生強(qiáng)度衰減，后逐步過渡至其他單元。

4 接觸單元剪切強(qiáng)度參數(shù)演化與剪切應(yīng)力-位移關(guān)系研究

參照文獻(xiàn)[10]中剪切強(qiáng)度參數(shù)應(yīng)變硬/軟化變化特征，分別建立接觸單元摩擦角與黏聚力均軟化的綜合軟化型及摩擦角增長(zhǎng)、黏聚力衰減的摩擦角增長(zhǎng)型兩類模型，每種模型考慮5 種剪切塑性變形特征值。接觸單元模型參數(shù)見表3，計(jì)算剪應(yīng)力-位移關(guān)系曲線見圖8、9。

表3 接觸單元模型參數(shù)Table 3 The parameters of contact elements models

4.1 接觸單元剪切強(qiáng)度參數(shù)演化

由圖8、9 可見，接觸單元剪切強(qiáng)度參數(shù)的演化直接導(dǎo)致在剪切過程中剪應(yīng)力-位移關(guān)系曲線呈線性變化關(guān)系，不同的應(yīng)變硬/軟化參數(shù)對(duì)剪應(yīng)力-位移關(guān)系曲線形態(tài)特征影響較大。結(jié)合本文計(jì)算結(jié)果并綜合室內(nèi)試驗(yàn)經(jīng)驗(yàn)，考慮巖土材料在剪切過程中黏聚力值以衰減為主，剪切過程中摩擦角的變化特征主控直剪試驗(yàn)剪應(yīng)力-位移曲線的不同形式特征。

圖8 綜合軟化型剪應(yīng)力-位移關(guān)系曲線Fig.8 Shear stress-displacement curves of cohesion and friction angle attenuation type

圖9 摩擦角增長(zhǎng)型剪應(yīng)力-位移關(guān)系曲線Fig.9 Shear stress-displacement curves of friction angle increasing type

4.1 接觸單元剪應(yīng)力-位移曲線分類

由此將剪應(yīng)力-位移曲線類型分為4 種：綜合軟化型、整體穩(wěn)定型、前期硬化型、后期硬化型，見圖10。

圖10 剪應(yīng)力-位移曲線分類Fig.10 Classification of shear stress-displacement curves

綜合軟化型一般表示巖土材料在剪切過程中，摩擦角及黏聚力均快速衰減，破壞過程以脆性破壞為主，如超固結(jié)土、結(jié)構(gòu)性土等；整體穩(wěn)定型一般表示巖土材料在剪切過程中，摩擦角及黏聚力均變化不大，或剪切強(qiáng)度參數(shù)較低；前期硬化型一般表示巖土材料在剪切過程中，剪切增加了接觸面的粗糙程度，摩擦角具有短期增加現(xiàn)象；后期硬化型一般表示巖土材料在剪切過程中，前期產(chǎn)生較大的黏聚力衰減，后期摩擦角逐步增大并主導(dǎo)了后續(xù)剪切過程。

綜上可見，在直剪過程中如產(chǎn)生剪應(yīng)力-位移曲線上升，主要是由內(nèi)摩擦角增加所致；如產(chǎn)生曲線下行，則致因較多，需綜合分析判斷。

5 結(jié)論

（1）基于FLAC 軟件中FISH 程序語言，可實(shí)現(xiàn)接觸單元的應(yīng)變軟化接觸本構(gòu)模型。

（2）通過含應(yīng)變軟化接觸單元直剪模型數(shù)值模擬，表明直剪試驗(yàn)具有明顯的漸進(jìn)性破壞特征，模型的剪應(yīng)力集中區(qū)域及破壞點(diǎn)從推力側(cè)啟動(dòng)，后逐步貫通整個(gè)接觸面。

（3）受應(yīng)變軟化接觸單元影響，各接觸點(diǎn)的應(yīng)力狀態(tài)均經(jīng)歷線彈性增長(zhǎng)、峰值后衰減和逐步穩(wěn)定均布的應(yīng)變軟化過程。

（4）在剪切過程中一般接觸面黏聚力以衰減為主，摩擦角的增減變化特征對(duì)直剪試驗(yàn)剪切力-位移曲線的形式影響較大，據(jù)此可將曲線類型分為綜合軟化型、整體穩(wěn)定型、前期硬化型和后期硬化型4種。

（5）由于巖土材料的復(fù)雜性，不同類型剪切應(yīng)力-位移曲線均存在多解性，在直剪過程中如產(chǎn)生剪應(yīng)力位移曲線上升現(xiàn)象，可判斷其主要是由內(nèi)摩擦角增加所致，如產(chǎn)生曲線下行，則影響原因較多，需綜合分析判斷。

[1]蔣明鏡,沈珠江.結(jié)構(gòu)性黏土剪切帶的徽觀分析[J].巖土工程學(xué)報(bào),1998,20(2):102－108.JIANG Ming-jing,SHEN Zhu-jiang.Microscopic analysis of shear band in structured clay[J].Chinese Journal of Geotechnical Engineering,1998,20(2):102－108.

[2]CLOUGH G W,DUNEAN J M.Finite element analysis of retaining wall behavior[J].Journal of Soil Mechanics&Foundation Division,ASCE,1971,97(12):1657－1674.

[3]BOULON M,NOVA R.Modelling of soil structure interface behavior:A comparison between elastoplastic and rate-type laws[J].Computers and Geotechnics,1990,9:21－46.

[4]胡黎明,馬杰,張丙印.直剪試驗(yàn)中接觸面漸進(jìn)破壞的數(shù)值模擬[J].清華大學(xué)學(xué)報(bào)(自然科學(xué)版).2008,48(6):943－946.HU Li-ming,MA Jie,ZHANG Bing-yin.Numerical simulation of interface failure during direct shear tests[J]Journal of Tsinghua University(Science and Technology),2008,48(6):943－946.

[5]王恭先,劉進(jìn)舉.順層滑坡薄層滑帶的抗剪強(qiáng)度特征與指標(biāo)的試驗(yàn)和選擇[J].水文地質(zhì)工程地質(zhì),1980,(6):6－9.

[6]任光明.聶德新,左三勝.滑帶土結(jié)構(gòu)強(qiáng)度再生研究[J].地質(zhì)災(zāi)害與環(huán)境保護(hù),1996,7(3):7－12.REN Guang-ming,NIE De-xin,ZUO San-sheng.Study on cohesion restoration of slide zone soil[J].Chinese Journal of Geological Hazards and Environment Preservation,1996,7(3):7－13.

[7]Itasca.Fish in Flac,User Manual of FLAC Version 5.0[M].[s.l.]:[s.n.],2003.

[8]漆泰岳,陸士良,高波.FLAC 錨桿單元模型的修正及其應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2004,23(13):2197－2200.QI Tai-yue,LU Shi-liang,GAO Bo.A mended model of bolt element in flac and its application[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(13):2197－2200.

[9]川本眺萬.ひずみ軟化を考慮した巖盤掘削の解析[C]//土木學(xué)會(huì)論文集.[s.l.]:[s.n.],1981,107－117.

[10]王水林,王威,吳振君.巖土材料峰值后區(qū)強(qiáng)度參數(shù)演化與應(yīng)力-應(yīng)變曲線關(guān)系研究[J].巖石力學(xué)與工程學(xué)報(bào),2010,29(8):1524－1528.WANG Shui-lin,WANG Wei,WU Zhen-jun.Study of relationship between evolution of post-peak strength parameters and stress-strain curves of geomaterials[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(8):1524－1528.