馬少坤,黃 驍,韋榕寬,劉 瑩,段智博

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004;2.廣西大學(xué) 工程防災(zāi)與結(jié)構(gòu)安全重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004)

0 引言

巖土體作為1種天然的工程材料，其復(fù)雜的工程特性在生產(chǎn)實(shí)踐中常常引發(fā)不可預(yù)知的安全事故，例如邊坡失穩(wěn)[1]、隧道塌方[2]等。而Cundall于20世紀(jì)70年代提出的離散元法為研究巖土的相關(guān)特性提供1種有效手段[3]，對(duì)工程安全具有重要意義。離散元法的核心接觸模型區(qū)別于傳統(tǒng)土力學(xué)理論中的宏觀本構(gòu)模型，主要由顆粒與顆粒間的細(xì)觀接觸特性及參數(shù)所表征[3-4]。然而在實(shí)際工程和試驗(yàn)中所測(cè)得巖土體的變形或強(qiáng)度特性均以宏觀參數(shù)來(lái)體現(xiàn)，這使得在利用離散元法對(duì)土體進(jìn)行仿真模擬時(shí)，選取與其匹配的接觸模型及細(xì)觀參數(shù)較困難。

大量學(xué)者針對(duì)接觸模型細(xì)觀參數(shù)與土體宏觀力學(xué)參數(shù)的關(guān)系進(jìn)行定量研究。曾遠(yuǎn)等[5]利用PFC2D離散元分析軟件，進(jìn)行多組無(wú)黏性土雙軸模擬試驗(yàn)，詳細(xì)分析細(xì)觀參數(shù)對(duì)土體宏觀特性的影響;陳亞?wèn)|等[6]建立多組不同細(xì)觀參數(shù)的模擬真三軸數(shù)值試驗(yàn)，分析無(wú)黏性土細(xì)觀參數(shù)與宏觀參數(shù)的關(guān)系，并提出針對(duì)無(wú)黏性土的細(xì)觀參數(shù)確定方法；徐小敏等[7]開(kāi)展一系列基于線性接觸模型的模擬試驗(yàn)，分析顆粒細(xì)觀參數(shù)與試樣宏觀參數(shù)的關(guān)系，并建立經(jīng)驗(yàn)公式。上述基于無(wú)黏性土的細(xì)觀參數(shù)研究均以“ball”為基礎(chǔ)單元，鮮有考慮真實(shí)土顆粒的形狀效應(yīng)，即忽略土顆粒之間的抗轉(zhuǎn)動(dòng)作用；而Iwashita等[8]通過(guò)試驗(yàn)發(fā)現(xiàn)顆粒間的抗轉(zhuǎn)動(dòng)作用對(duì)顆粒的力學(xué)性質(zhì)有顯著的影響。因此，部分學(xué)者也對(duì)抗轉(zhuǎn)動(dòng)模型進(jìn)行研究[9]。

此外，為探討宏細(xì)觀參數(shù)間的非線性關(guān)系，利用計(jì)算機(jī)算法對(duì)模擬數(shù)據(jù)進(jìn)行分析。周喻等[10]采用BP神經(jīng)網(wǎng)絡(luò)方法對(duì)400 組細(xì)觀力學(xué)參數(shù)隨機(jī)組合的三軸模擬試驗(yàn)數(shù)據(jù)建立宏觀力學(xué)參數(shù)與細(xì)觀力學(xué)參數(shù)的非線性模型，結(jié)果表明此模型可快速、準(zhǔn)確地反演巖土體的非線性關(guān)系；李澄清等[11]利用25組雙軸壓縮離散元數(shù)值試驗(yàn)?zāi)M結(jié)果建立BP人工神經(jīng)網(wǎng)絡(luò)反演系統(tǒng)，結(jié)果表明模型反演性能較好；崔洋洋[12]開(kāi)展多組基于線性接觸模型的雙軸壓縮試驗(yàn)，利用BP人工神經(jīng)網(wǎng)絡(luò)分析方法尋找宏-細(xì)觀參數(shù)之間的映射關(guān)系。

綜上所述，本文基于抗轉(zhuǎn)動(dòng)線性接觸模型，采用控制變量法開(kāi)展一系列三軸試驗(yàn)數(shù)值分析，通過(guò)對(duì)比不同細(xì)觀參數(shù)下的模擬試驗(yàn)結(jié)果，研究模型的不同細(xì)觀參數(shù)對(duì)宏觀參數(shù)的影響；針對(duì)對(duì)比試驗(yàn)無(wú)法體現(xiàn)出參數(shù)交互作用這一難題，引入隨機(jī)森林機(jī)器學(xué)習(xí)法，詳細(xì)探究宏細(xì)觀參數(shù)的權(quán)重關(guān)系。

1 模型驗(yàn)證及結(jié)果

1.1 計(jì)算參數(shù)選取

為研究模型細(xì)觀參數(shù)的敏感性，基于已有試驗(yàn)數(shù)據(jù)開(kāi)展相應(yīng)的數(shù)值試驗(yàn)以獲得1組合理的土體基本細(xì)觀參數(shù)?；谖墨I(xiàn)[13]，選取豐浦砂為模擬對(duì)象，選取相關(guān)計(jì)算參數(shù)，見(jiàn)表1。

表1 計(jì)算參數(shù)

此外，采用擴(kuò)大顆粒粒徑的方法保證計(jì)算效率。通過(guò)多次試算，最終選取由3 018個(gè)顆粒組成的試樣，該數(shù)量足以反應(yīng)材料特性[7],級(jí)配曲線對(duì)比如圖1所示。最終標(biāo)定結(jié)果見(jiàn)表2。

圖1 原始級(jí)配與模擬級(jí)配曲線

表2 豐浦砂標(biāo)定參數(shù)

1.2 模擬結(jié)果

豐浦砂三軸固結(jié)排水壓縮室內(nèi)試驗(yàn)與數(shù)值試驗(yàn)的結(jié)果對(duì)比如圖2～3所示。由圖2～3可知，在試驗(yàn)初期，式樣的總體積減小，而偏應(yīng)力-軸應(yīng)變曲線近似為直線。隨著試驗(yàn)不斷進(jìn)行，試樣出現(xiàn)剪脹的趨勢(shì)，而偏應(yīng)力的增長(zhǎng)不斷放緩，直至到達(dá)峰值后，偏應(yīng)力逐漸下降?？傮w而言，模擬結(jié)果與室內(nèi)試驗(yàn)的結(jié)果契合度均較好，整體誤差較小，證明利用該組細(xì)觀參數(shù)可較好地表征豐浦砂。

圖2 豐浦砂室內(nèi)試驗(yàn)與數(shù)值試驗(yàn)偏應(yīng)力-軸應(yīng)變曲線對(duì)比

圖3 豐浦砂室內(nèi)試驗(yàn)與數(shù)值試驗(yàn)的體應(yīng)變-軸應(yīng)變曲線對(duì)比

2 參數(shù)敏感性分析

在土體參數(shù)標(biāo)定中，將土體宏觀參數(shù)(包括初始彈性模量、泊松比、峰值強(qiáng)度等)與模型細(xì)觀參數(shù)相對(duì)應(yīng)。根據(jù)徐小敏等[7]的研究，在離散元三軸模擬試驗(yàn)中，當(dāng)試樣的軸向應(yīng)變處于峰值應(yīng)力所對(duì)應(yīng)軸向應(yīng)變的1%范圍內(nèi)時(shí)，可視為試樣處于彈性形變階段。可根據(jù)此范圍內(nèi)的偏應(yīng)力-軸應(yīng)變曲線與體應(yīng)變-軸應(yīng)變曲線計(jì)算出試樣所對(duì)應(yīng)的初始彈性模量與泊松比，如式(1)～(2)所示：

(1)

(2)

為研究和表征出各個(gè)細(xì)觀參數(shù)對(duì)宏觀參數(shù)的影響，基于本文所得豐浦砂對(duì)應(yīng)的模型參數(shù)，為各個(gè)細(xì)觀參數(shù)設(shè)計(jì)進(jìn)行對(duì)比試驗(yàn)，模擬試驗(yàn)方案見(jiàn)表3。

表3 模擬試驗(yàn)方案

2.1 有效模量的影響

有效模量會(huì)顯著影響試樣的初始彈性模量。不同圍壓下有效模量與試樣初始彈性模量的變化規(guī)律如圖4所示。由圖4可知，二者呈正相關(guān)關(guān)系。即試樣的初始彈性模量隨著有效模量的增大而增大。這是因?yàn)楫?dāng)顆粒間有效模量增大時(shí)，顆粒間法向剛度與切向剛度隨之增大，而三軸試驗(yàn)中上下加載板的加載速率保持恒定不變，因此在試樣加載過(guò)程中，顆粒間的變形等速率增加，而顆粒間法向剛度與切向剛度增大導(dǎo)致其接觸力變大，故在試樣處于較小的且相同的軸應(yīng)變時(shí)，有效模量越大的試樣的軸應(yīng)力越大，從而在宏觀上體現(xiàn)為顆粒間有效模量的增加會(huì)顯著增大試樣的初始彈性模量。

圖4 不同圍壓下初始彈性模量隨有效模量變化規(guī)律曲線

此外，有效模量的變化還會(huì)對(duì)試樣的泊松比產(chǎn)生影響。不同圍壓下泊松比隨有效模量變化規(guī)律曲線如圖5所示。由圖5可知，泊松比會(huì)隨著有效模量的增大而增大。這是因?yàn)楫?dāng)顆粒間有效模量增大時(shí)，顆粒間法向剛度和切向剛度隨之增大，而試樣處于軸向加載速率不變的條件下，顆粒間的徑向合力會(huì)增大，此時(shí)伺服側(cè)壁需要更大的徑向位移以保證圍壓的恒定，在宏觀參數(shù)上表現(xiàn)為泊松比隨之增大。

圖5 不同圍壓下泊松比隨有效模量變化規(guī)律曲線

2.2 剛度比的影響

剛度比的變化會(huì)顯著影響土體初始抗變形能力。不同圍壓下泊松比隨剛度比變化關(guān)系曲線如圖6所示。由圖6可知，試樣的泊松比隨剛度比的增大而增大。這是因?yàn)楫?dāng)顆粒的有效模量不變時(shí)，剛度比的增大使得顆粒間的切向剛度不斷減小，剪切帶內(nèi)部顆?？商峁┑慕佑|力徑向合力不足以與圍壓保持平衡。因此試樣內(nèi)部需調(diào)動(dòng)更大范圍的顆粒來(lái)平衡外部圍壓，造成更大范圍內(nèi)顆粒的相對(duì)滑移，從而使試樣產(chǎn)生更大的徑向變形。因此剛度比的增大會(huì)導(dǎo)致試樣整體的切向抗變形能力下降，在宏觀參數(shù)上反映為泊松比隨剛度比的增大而增大。

圖6 不同圍壓下泊松比隨剛度比變化關(guān)系曲線

此外，剛度比還對(duì)試樣的初始彈性模量產(chǎn)生影響。不同圍壓下初始彈性模量隨剛度比變化關(guān)系曲線如圖7所示。由圖7可知，初始彈性模量會(huì)隨著剛度比的增大而不斷減小，并最終趨于穩(wěn)定。這是因?yàn)楫?dāng)剛度比較小時(shí)，試樣具有較高的切向抗變形能力，剪切帶內(nèi)顆?？杀ＷC在較小的相對(duì)位移下提供與側(cè)向圍壓保持平衡的力，此時(shí)試樣總體積減小，內(nèi)部顆粒重疊量增大，軸向應(yīng)力較高，從而具有較高的初始彈性模量；而隨著剛度比的增加，其切向抗變形能力急劇下降，使得試樣內(nèi)部調(diào)動(dòng)更多的顆粒相對(duì)移動(dòng)以保持與側(cè)向圍壓的平衡，導(dǎo)致在試驗(yàn)初期產(chǎn)生較大側(cè)向變形。試樣總體積增加，內(nèi)部顆粒重疊量變小，軸向應(yīng)力下降，從而造成初始彈性模量下降。而當(dāng)顆粒的切向剛度減小到一定程度時(shí)，試樣已完全調(diào)動(dòng)內(nèi)部顆粒以抵擋外部荷載，此時(shí)繼續(xù)增加剛度比對(duì)初始彈性模量影響不大，因此在曲線上體現(xiàn)為初始彈性模量隨著剛度比的增大而趨于定值。

圖7 不同圍壓下初始彈性模量隨剛度比變化關(guān)系曲線

2.3 摩擦系數(shù)的影響

摩擦系數(shù)主要對(duì)土體的強(qiáng)度產(chǎn)生影響。因3種不同圍壓條件下，不同摩擦系數(shù)下的偏應(yīng)力-軸應(yīng)變關(guān)系曲線的規(guī)律相似，此處選取100 kPa為例進(jìn)行分析，規(guī)律曲線如圖8所示，其中μ為摩擦系數(shù)。由圖8可知，摩擦系數(shù)與峰值強(qiáng)度呈正相關(guān)，而殘余強(qiáng)度會(huì)隨著摩擦系數(shù)的增大而趨于定值。

圖8 100 kPa圍壓時(shí)不同摩擦系數(shù)下的偏應(yīng)力-軸應(yīng)變關(guān)系曲線

隨著摩擦系數(shù)的增加，偏應(yīng)力-軸應(yīng)變曲線的變化可大致分為2個(gè)階段。當(dāng)摩擦系數(shù)較小時(shí)，土體偏應(yīng)力-軸應(yīng)變曲線整體呈現(xiàn)出“松砂型”的應(yīng)變硬化趨勢(shì)。而當(dāng)摩擦系數(shù)較大時(shí)，偏應(yīng)力-軸應(yīng)變曲線呈現(xiàn)出“密砂型”的應(yīng)變軟化趨勢(shì)。顯然，試樣的細(xì)觀摩擦系數(shù)和所模擬土的宏觀密實(shí)度之間存在相關(guān)性，即在同一粒徑分布下，細(xì)觀摩擦系數(shù)與土體的宏觀壓實(shí)度存在一定關(guān)系，這與文獻(xiàn)[14]類似。

為考察摩擦系數(shù)對(duì)試樣黏聚力的影響，繪制不同摩擦系數(shù)下強(qiáng)度包線及內(nèi)摩擦角關(guān)系圖，如圖9所示。由圖9可知，隨著摩擦系數(shù)的增大，試樣的內(nèi)摩擦角顯著增大，而增長(zhǎng)速率隨之減?。欢鴱膹?qiáng)度包線明顯看出，摩擦系數(shù)的增加不會(huì)使試樣產(chǎn)生黏聚力。

圖9 摩擦系數(shù)對(duì)強(qiáng)度包線及內(nèi)摩擦角的影響

2.4 抗轉(zhuǎn)動(dòng)系數(shù)的影響

抗轉(zhuǎn)動(dòng)系數(shù)會(huì)顯著影響土體的抗剪強(qiáng)度。因3種不同圍壓條件下，不同抗轉(zhuǎn)動(dòng)系數(shù)下的偏應(yīng)力-軸應(yīng)變關(guān)系曲線的規(guī)律相似，此處選取100 Pka為例進(jìn)行分析，100 kPa圍壓時(shí)不同抗轉(zhuǎn)動(dòng)系數(shù)下的偏應(yīng)力-軸應(yīng)變關(guān)系曲線如圖10所示，其中，μr為抗轉(zhuǎn)動(dòng)系數(shù)。通過(guò)對(duì)模擬結(jié)果進(jìn)行對(duì)比與分析，發(fā)現(xiàn)抗轉(zhuǎn)動(dòng)系數(shù)的取值存在作用區(qū)間，在此區(qū)間內(nèi)，試樣的峰值強(qiáng)度和殘余強(qiáng)度均會(huì)隨著系數(shù)的增大而增大，呈現(xiàn)出整體上升的趨勢(shì)。當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)接近閾值時(shí)，土體偏應(yīng)力-軸應(yīng)變曲線基本不發(fā)生變化。其主要原因是當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)較大時(shí)，顆粒間臨界轉(zhuǎn)動(dòng)阻力矩較大，此時(shí)試樣中的顆粒間的相對(duì)轉(zhuǎn)動(dòng)矩未達(dá)到臨界轉(zhuǎn)動(dòng)阻力矩，導(dǎo)致顆粒間發(fā)生的相對(duì)轉(zhuǎn)動(dòng)較小，宏觀上表現(xiàn)為試樣強(qiáng)度和變形參數(shù)幾乎不隨之變化。

圖10 100 kPa圍壓時(shí)不同抗轉(zhuǎn)動(dòng)系數(shù)下偏應(yīng)力-軸應(yīng)變關(guān)系曲線

為考察抗轉(zhuǎn)動(dòng)系數(shù)對(duì)試樣黏聚力的影響，繪制不同抗轉(zhuǎn)動(dòng)系數(shù)下強(qiáng)度包線與內(nèi)摩擦角關(guān)系圖，如圖11所示。由圖11可知，抗轉(zhuǎn)動(dòng)系數(shù)不會(huì)使試樣產(chǎn)生黏聚力，進(jìn)一步驗(yàn)證抗轉(zhuǎn)動(dòng)模型對(duì)無(wú)黏性土的適用性。

圖11 抗轉(zhuǎn)動(dòng)系數(shù)對(duì)強(qiáng)度包線及內(nèi)摩擦角的影響

此外，抗轉(zhuǎn)動(dòng)系數(shù)也對(duì)試樣的剪脹性有顯著影響。不同抗轉(zhuǎn)動(dòng)系數(shù)下試樣的體應(yīng)變-軸應(yīng)變變化曲線如圖12所示。由圖12可知，隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大，試樣的體應(yīng)變?cè)诩虞d初期的彈性階段并未產(chǎn)生較大的變化，即泊松比不被抗轉(zhuǎn)動(dòng)系數(shù)所影響。而隨著加載不斷進(jìn)行，曲線發(fā)生分化，抗轉(zhuǎn)動(dòng)系數(shù)越大的試樣剪脹趨勢(shì)越強(qiáng)，且當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)接近閾值時(shí)趨于穩(wěn)定，也進(jìn)一步證明抗轉(zhuǎn)動(dòng)系數(shù)的取值存在主要作用區(qū)間這一結(jié)論。

圖12 100 kPa圍壓時(shí)不同抗轉(zhuǎn)動(dòng)系數(shù)下的體應(yīng)變-軸應(yīng)變曲線

3 基于隨機(jī)森林回歸模型的參數(shù)權(quán)重分析

由于模型的參數(shù)間存在一定的交互作用，即多個(gè)細(xì)觀參數(shù)會(huì)對(duì)單個(gè)宏觀參數(shù)共同產(chǎn)生影響。而以上控制變量試驗(yàn)無(wú)法精確體現(xiàn)出影響的權(quán)重，因此引入隨機(jī)森林機(jī)器學(xué)習(xí)法作為工具，對(duì)參數(shù)進(jìn)行進(jìn)一步的精細(xì)化分析。

隨機(jī)森林算法是廣泛使用于計(jì)算機(jī)領(lǐng)域的1種機(jī)器學(xué)習(xí)算法，其基于決策樹(shù)與bagging思想，通過(guò)建立多棵決策樹(shù)以提高預(yù)測(cè)性能，其預(yù)測(cè)值由各棵樹(shù)輸出的平均數(shù)所決定[15]。該算法具有權(quán)重分析功能，可以考察各個(gè)自變量特征的重要性。

本文基于抗轉(zhuǎn)動(dòng)線性接觸模型開(kāi)展300組隨機(jī)細(xì)觀參數(shù)的三軸試驗(yàn)，并利用python語(yǔ)言sklearn庫(kù)基于隨機(jī)森林算法建立預(yù)測(cè)模型，根據(jù)文獻(xiàn)[16]與本文模擬結(jié)果，確定主要參數(shù)范圍，見(jiàn)表4。

表4 隨機(jī)參數(shù)范圍

3.1 模型建立過(guò)程

以模型的4個(gè)自有參數(shù)作為輸入層參數(shù)，以初始彈性模量、泊松比與最大偏應(yīng)力3個(gè)參數(shù)作為輸出層參數(shù)，并隨機(jī)選取300組數(shù)據(jù)中的240組作為訓(xùn)練集，60組作為測(cè)試集，以擬合度R2作為模型好壞評(píng)估的指標(biāo)。

3.2 模型調(diào)參

隨機(jī)森林模型具有較多的模型參數(shù)，正確的參數(shù)選擇能夠減少模型過(guò)擬合并保證其具有較高的準(zhǔn)確度和泛化性能。其中回歸樹(shù)數(shù)量n_estmators與最大樹(shù)深max_depth2個(gè)參數(shù)對(duì)模型優(yōu)劣起著重要作用。采用網(wǎng)格搜索法找到最佳模型參數(shù)，所建立模型的優(yōu)劣以擬合優(yōu)度R2為考察指標(biāo)，所得參數(shù)和模擬結(jié)果見(jiàn)表5，結(jié)果表明，訓(xùn)練的模型具有較高的準(zhǔn)確性。

表5 隨機(jī)森林模型參數(shù)與擬合度

3.3 權(quán)重分析

由于隨機(jī)森林模型存在隨機(jī)性，每次訓(xùn)練均可能得到規(guī)律相同但權(quán)重系數(shù)有一定差別的結(jié)果，因此為消除其影響，用相同的模型參數(shù)訓(xùn)練10 000次模型再對(duì)權(quán)重參數(shù)取平均值作為最終結(jié)果。

隨機(jī)森林輸出的宏細(xì)觀參數(shù)權(quán)重柱狀圖如圖13所示。由圖13可知，權(quán)重值越大代表影響程度越大，其總體規(guī)律與本文對(duì)比試驗(yàn)結(jié)果相似。其中，抗轉(zhuǎn)動(dòng)系數(shù)(0.564 72)與摩擦系數(shù)(0.427 5)顯著影響試樣的最大偏應(yīng)力且前者的影響偏大；剛度比主要對(duì)泊松比(0.751)產(chǎn)生決定性影響，對(duì)初始彈性模量(0.201 57)起次要作用；有效模量對(duì)初始彈性模量(0.695 14)起主要作用，而對(duì)泊松比(0.193 64)起次要作用。

圖13 隨機(jī)森林輸出的宏細(xì)觀參數(shù)權(quán)重柱狀圖

4 結(jié)論

1)試樣的變形參數(shù)主要受顆粒有效模量與剛度比的影響。有效模量的增加會(huì)顯著增加試樣的初始彈性模量且增加試樣的泊松比。剛度比的增加會(huì)顯著降低試樣的初始彈性模量且增加泊松比。

2)試樣的強(qiáng)度參數(shù)主要受顆粒間摩擦系數(shù)與抗轉(zhuǎn)動(dòng)系數(shù)的影響。摩擦系數(shù)的增加會(huì)顯著增加試樣的峰值強(qiáng)度而幾乎不影響試樣的殘余強(qiáng)度；而抗轉(zhuǎn)動(dòng)系數(shù)的取值存在作用區(qū)間，當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)處于區(qū)間內(nèi)時(shí)，試樣的峰值強(qiáng)度與殘余強(qiáng)度會(huì)隨著抗轉(zhuǎn)動(dòng)系數(shù)的增大而增大，而當(dāng)抗轉(zhuǎn)動(dòng)系數(shù)取值大于區(qū)間時(shí)，試樣的峰值強(qiáng)度與殘余強(qiáng)度將不再增加。此外，抗轉(zhuǎn)動(dòng)系數(shù)與試樣的剪脹性顯著相關(guān)。

3)隨機(jī)森林可以較好地表征宏細(xì)觀參數(shù)中的影響關(guān)系。對(duì)于試樣初始彈性模量而言，主要因素的影響權(quán)重排序?yàn)橛行Ａ?0.695 14)>剛度比(0.201 57)；對(duì)于試樣泊松比而言，主要因素的影響權(quán)重排序?yàn)閯偠缺?0.751)>有效模量(0.193 64)；對(duì)于試樣峰值而言，主要因素的影響權(quán)重排序?yàn)榭罐D(zhuǎn)動(dòng)系數(shù)(0.564 72)>摩擦系數(shù)(0.427 5)。