考慮滾動阻抗的線性接觸模型離散元宏細觀參數(shù)敏感性研究*

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

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

0 引言

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

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

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

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

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

1.1 計算參數(shù)選取

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

表1 計算參數(shù)

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

圖1 原始級配與模擬級配曲線

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

1.2 模擬結(jié)果

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

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

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

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

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

(1)

(2)

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

表3 模擬試驗方案

2.1 有效模量的影響

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

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

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

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

2.2 剛度比的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 模型建立過程

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

3.2 模型調(diào)參

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

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

3.3 權(quán)重分析

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

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

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

4 結(jié)論

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

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

3)隨機森林可以較好地表征宏細觀參數(shù)中的影響關(guān)系。對于試樣初始彈性模量而言，主要因素的影響權(quán)重排序為有效模量(0.695 14)>剛度比(0.201 57)；對于試樣泊松比而言，主要因素的影響權(quán)重排序為剛度比(0.751)>有效模量(0.193 64)；對于試樣峰值而言，主要因素的影響權(quán)重排序為抗轉(zhuǎn)動系數(shù)(0.564 72)>摩擦系數(shù)(0.427 5)。