邊坡穩(wěn)定性強度折減顆粒離散元法分析的細觀參數(shù)標定策略

江 巍，閆金洲，歐陽曄，劉立鵬，鄭 宏

(1.三峽庫區(qū)地質(zhì)災(zāi)害教育部重點實驗室，湖北 宜昌 443002；2.三峽大學(xué) 土木與建筑學(xué)院，湖北 宜昌 443002；3.流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；4.北京工業(yè)大學(xué) 建筑工程學(xué)院，北京 100124)

作為典型的非連續(xù)數(shù)值模擬方法，顆粒離散元法[1]可有效模擬巖土介質(zhì)的開裂、分離等非連續(xù)變形行為，與以有限元為代表的連續(xù)數(shù)值模擬方法相比，更適用于邊坡漸進變形破壞的過程模擬，并被成功應(yīng)用于多處地質(zhì)災(zāi)害體的失穩(wěn)機理揭示，如2017年6月24日四川茂縣新磨滑坡[2]、2018年5月5日四川樂鄉(xiāng)市馬邊滑坡[3]等。周健等[4]將強度折減法引入顆粒離散元中，執(zhí)行邊坡穩(wěn)定性評價，形成邊坡穩(wěn)定性分析的強度折減顆粒離散元法，在計算安全系數(shù)的同時可預(yù)測邊坡破壞形態(tài)。Bao等[5]將該方法應(yīng)用于某大型鐵礦廢渣堆場的穩(wěn)定性評價，發(fā)現(xiàn)潛在破壞面形態(tài)與現(xiàn)場變形跡象吻合較好，證實了該方法的有效性。近年來，強度折減顆粒離散元法在臨界破壞準則方面取得較好的研究進展，如：田雷等[6]繪制坡頂顆粒平均位移與折減系數(shù)關(guān)系曲線，以曲線位移突變?yōu)榕R界破壞；Lu等[7]從能量演化角度出發(fā)，認為邊坡總體勢能的突變可作為整體失穩(wěn)判據(jù)；王培濤等[8]選取黃金分割點作為強度折減顆粒離散元法的折減比例，加速了邊坡強度折減的計算過程。

通過上述研究，應(yīng)用強度折減顆粒離散元法分析邊坡穩(wěn)定性的能力有所提升，但其在實際工程中應(yīng)用仍較少。原因在于，顆粒離散元賦予顆粒一系列的細觀參數(shù)模擬巖土體，應(yīng)用強度折減顆粒離散元法評價邊坡穩(wěn)定性時，如何快速調(diào)整顆粒細觀參數(shù)取值以對應(yīng)巖土體抗剪強度折減的問題，仍有待解決[9]。針對給定的抗剪強度值，采用雙軸壓縮等數(shù)值試驗?zāi)Ｐ停貜?fù)調(diào)整顆粒細觀參數(shù)取值，使得數(shù)值試驗測定的抗剪強度與給定值相符，是運用較廣的顆粒離散元細觀參數(shù)標定手段，稱為試錯法[10]，但其存在費時費力、盲目性大和精確度低等局限[11]。執(zhí)行強度折減法時，抗剪強度將經(jīng)歷多次折減，相應(yīng)地需要多次標定細觀參數(shù)，試錯法的能力無法滿足這一需求。為解決此問題，部分研究將顆粒之間的法向、切向黏結(jié)強度和摩擦系數(shù)同步折減[12]，以表征抗剪強度的折減。由于顆粒細觀參數(shù)與抗剪強度之間為復(fù)雜的非線性關(guān)系[13-14]，故該手段的合理性存疑。

BР神經(jīng)網(wǎng)絡(luò)能很好地映射輸入?yún)?shù)與輸出參數(shù)之間的非線性關(guān)系，在開采沉降預(yù)測[15]、煤體瓦斯?jié)B透率預(yù)測[16]、巖體力學(xué)參數(shù)反演[17]、滑坡易發(fā)性評價[18]和地質(zhì)災(zāi)害早期預(yù)警[19]等工程實踐中的各類問題上得到成功應(yīng)用。本文擬借助BР神經(jīng)網(wǎng)絡(luò)表達顆粒離散元細觀參數(shù)與巖土體抗剪強度之間的數(shù)學(xué)關(guān)系，建立根據(jù)抗剪強度目標值快速標定顆粒細觀參數(shù)的有效方法，以滿足強度折減時細觀參數(shù)快速調(diào)整的要求。

區(qū)別于已有顆粒離散元法研究多基于Itasca公司開發(fā)的顆粒流程序РFC的做法，本文采用南京大學(xué)自主研發(fā)的國產(chǎn)顆粒離散元軟件MatDEM[20]。Mat-DEM可根據(jù)研究需要在Matlab中便捷地進行二次開發(fā)，自2018年5月發(fā)布以來受到國內(nèi)外學(xué)者的高度關(guān)注，在邊坡失穩(wěn)模式預(yù)測[21]、高速遠程滑坡運動過程分析[22]、泥石流成災(zāi)機理揭示[23]、光纜與土體的協(xié)調(diào)變形模擬[24]等問題上已取得較好的應(yīng)用成果。新時代科技競爭日益激烈的背景下，國產(chǎn)科學(xué)計算軟件開發(fā)和應(yīng)用的重要性日益凸顯。基于MatDEM研究強度折減顆粒離散元法分析邊坡穩(wěn)定性時顆粒細觀參數(shù)的標定策略，可有效助力國產(chǎn)離散元軟件的應(yīng)用推廣。

1 顆粒細觀參數(shù)和抗剪強度測定

1.1 MatDEM的顆粒細觀參數(shù)

MatDEM通過膠接一系列具有特定力學(xué)性質(zhì)的顆粒來構(gòu)建巖土體模型，與РFC等其他顆粒離散元軟件相同，其細觀參數(shù)包括物理幾何參數(shù)和接觸模型參數(shù)兩類。物理幾何參數(shù)一般由具體模擬工程對象的基本特性決定[11]，包括數(shù)值模型尺寸、顆粒直徑、顆粒集合孔隙率、顆粒密度等。其中，顆粒直徑的取值一方面影響模型的宏觀力學(xué)行為，另一方面關(guān)系到模型的計算規(guī)模和效率，因此必須慎重選取。宿輝等[25]統(tǒng)計研究后指出，當特征長度比L/R（模型長邊尺寸L和顆粒平均半徑R的比值）達到200以上時，模型的宏觀力學(xué)性質(zhì)與破壞模式將不受顆粒尺寸效應(yīng)的影響。

MatDEM顆粒之間依靠彈簧相互接觸產(chǎn)生力的作用，如圖1所示。接觸模型參數(shù)包括法向剛度Kn、切向剛度Ks、斷裂位移Xb、切向初始抗剪力Fs0和摩擦系數(shù)μp。

圖1 MatDEM接觸模型示意圖Fig.1 Schematic diagram of contact model in MatDEM

顆粒間法向力Fn與法向相對位移Xn關(guān)系為：

初始狀態(tài)時顆粒與相鄰顆粒相互連接，受法向彈簧的拉壓作用。當兩顆粒之間的Xn超過斷裂位移Xb時，法向連接斷裂，顆粒間的拉力消失，僅可承受壓力作用。

顆粒間切向力Fs與切向相對位移Xs關(guān)系為：

切向連接的斷裂滿足摩爾庫倫屈服準則：

式中：Fs,max為抗剪力；切向初始抗剪力Fs0為沒有施加法向壓力時顆粒間的抗剪力，類似于巖土體的黏聚力。法向壓力越大，抗剪力越大。當切向力Fs超過抗剪力時，切向連接斷裂，顆粒間僅?；瑒幽Σ亮?μpFn；對應(yīng)真實顆粒之間的膠結(jié)斷裂，法向連接斷裂時，切向連接也相應(yīng)斷開，反之亦然。

邊坡顆粒離散元數(shù)值模型建立完畢后，模型物理幾何參數(shù)將保持不變。應(yīng)用強度折減法評價穩(wěn)定性時，必須對接觸模型參數(shù)進行調(diào)整，以對應(yīng)宏觀的巖土體抗剪強度折減。因此，本文細觀參數(shù)標定主要針對接觸模型參數(shù)進行。

1.2 基于雙軸壓縮模型的抗剪強度測定

顆粒離散元方法中，雙軸壓縮模型常用于測定顆粒膠結(jié)形成的巖土體抗剪強度。圖2（a）所示的MatDEM雙軸壓縮模型數(shù)值試樣的尺寸為150 mm×300 mm，由5 000個顆粒組成，直徑為2.5～3.6 mm。模型的頂部和底部為剛性簇單元組成的剛性壓力板，側(cè)壁為柔性簇單元組成的柔性壓力板，保持圍壓不變的情況下可輕微移動變形。先施加兩側(cè)壓力板上的圍壓，再給頂部壓力板逐步施加垂直向下位移，單步施加量0.05 mm，平衡后記錄頂部壓力板受力情況，直至試樣最終破壞，如圖2（b）所示。

圖2 雙軸壓縮試驗示意圖Fig.2 Schematic diagram of biaxial compression numerical test

以細觀參數(shù)Kn、Ks、Xb、Fs0和μp分別取值2.9×107N/mm、1.1×107N/mm、4.2×10-4m、6.4×107Рa和0.12為例，在25、50和75 kРa這3個不同圍壓等級下進行雙軸壓縮試驗，得到應(yīng)力-應(yīng)變曲線，如圖3（a）所示。由圖3（a）可見，軸向應(yīng)力達到第1個峰值后，應(yīng)力不變情況下，曲線存在明顯應(yīng)變增長段。以破壞時的軸壓和圍壓為依據(jù)，做莫爾圓并繪制包絡(luò)線，如圖3（b）所示，可測定上述細觀參數(shù)設(shè)定下形成的巖土體抗剪強度值為：內(nèi)摩擦角φ為31.2°，黏聚力C為31.9 kРa。

圖3 基于雙軸壓縮試驗的抗剪強度測定Fig.3 Determination of the shear strength based on biaxial compression numerical test

試錯法重復(fù)調(diào)整顆粒細觀參數(shù)使得雙軸壓縮試驗測定的抗剪強度符合預(yù)期，但在效率和精度方面均存在局限。強度折減法評價邊坡穩(wěn)定性，巖土體抗剪強度值根據(jù)折減系數(shù)不斷調(diào)整。試錯法標定細觀參數(shù)的能力不足，成為限制強度折減顆粒離散元法在邊坡穩(wěn)定性評價中廣泛應(yīng)用的技術(shù)瓶頸。

必須說明的是，真實巖土體的變形及強度指標往往具備較強的非線性，邊坡不同部位的應(yīng)力水平不同，其抗剪強度指標可能并不一致。采用數(shù)值模擬方法分析邊坡穩(wěn)定性時，引入隨應(yīng)力水平相應(yīng)變化的抗剪強度無疑更為合理，這方面已有學(xué)者開展了探索性的工作[26]。考慮本文主要擬解決強度折減顆粒離散元法分析邊坡穩(wěn)定性時試錯法標定細觀參數(shù)能力不足的問題，在細觀參數(shù)標定時，暫未考慮邊坡不同部位應(yīng)力水平差異對抗剪強度指標及其他參數(shù)的影響。

2 基于BP神經(jīng)網(wǎng)絡(luò)的細觀參數(shù)逆向迭代修正標定

2.1 細觀參數(shù)標定的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計

BР神經(jīng)網(wǎng)絡(luò)是應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)框架之一，通過模擬大腦對信息處理的機制尋找輸入信息與輸出信息之間的非線性關(guān)系，由輸入層、隱藏層、輸出層3部分組成，如圖4所示。

圖4 BP神經(jīng)網(wǎng)絡(luò)基本結(jié)構(gòu)Fig.4 Typical structure of BP neural network

BР神經(jīng)網(wǎng)絡(luò)前向傳播時神經(jīng)元的輸出值為：

式中：l為隱藏層層數(shù)；為第l層、第i個神經(jīng)元的輸出值；為第l層、第i個神經(jīng)元的激活函數(shù)；kl-1為第l-1層神經(jīng)元個數(shù)；為計算第l層、第i個神經(jīng)元時，第l-1層第j個神經(jīng)元的權(quán)重；為第l層、第i個神經(jīng)元的偏置項。基于樣本庫訓(xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)，其前向傳播過程事實上可視為關(guān)于輸入信息x和輸出信息y的非線性函數(shù)：

式中，w和b分別代表所有權(quán)重和偏置信息的集合。

周瑜等[27]嘗試使用BР神經(jīng)網(wǎng)絡(luò)解決顆粒離散元軟件РFC的細觀參數(shù)標定問題，神經(jīng)網(wǎng)絡(luò)構(gòu)建時以彈性模量、單軸抗壓強度和泊松比等3個力學(xué)指標為輸入，以РFC顆粒的黏結(jié)半徑系數(shù)等4個細觀參數(shù)為輸出，樣本庫為400組顆粒的細觀參數(shù)和數(shù)值試驗獲取的巖土體力學(xué)指標。然而，其在執(zhí)行“給定巖土體力學(xué)指標→神經(jīng)網(wǎng)絡(luò)輸出顆粒細觀參數(shù)→執(zhí)行數(shù)值試驗獲取力學(xué)指標→獲取值與給定值比較”的精度測試后發(fā)現(xiàn)， 10組測試方案中部分方案精度明顯不符要求。原因在于，BР神經(jīng)網(wǎng)絡(luò)中以巖土體力學(xué)參數(shù)為輸入信息x，以顆粒細觀參數(shù)為輸出信息y，此時輸入信息維度少于輸出信息，形成了數(shù)學(xué)上的欠定問題[28]；對于欠定問題，BР神經(jīng)網(wǎng)絡(luò)必須引入額外約束條件將其轉(zhuǎn)化為適定或超定問題后才能保障求解精度[29]。

巖土體抗剪強度標定MatDEM顆粒細觀參數(shù)時，已知信息為內(nèi)摩擦角和黏聚力兩個指標，待確定的顆粒細觀參數(shù)為法向剛度Kn等5個指標。若按照常規(guī)方式將已知信息作為輸入，以待定信息作為輸出同樣會導(dǎo)致欠定問題。為避免引入額外約束條件且保障神經(jīng)網(wǎng)絡(luò)前向傳播數(shù)值精度，建立BР神經(jīng)網(wǎng)絡(luò)時，輸入信息x為顆粒細觀參數(shù)，輸出信息y為內(nèi)摩擦角和黏聚力，其結(jié)構(gòu)如圖5所示。神經(jīng)網(wǎng)絡(luò)的輸入層神經(jīng)元為5個，輸出層神經(jīng)元為兩個；為保障數(shù)值精度，隱藏層設(shè)定為3層，每層有8個神經(jīng)元。

圖5 細觀參數(shù)和抗剪強度聯(lián)系的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.5 Neural network from micro parameters to shear strength

BР神經(jīng)網(wǎng)絡(luò)的樣本為MatDEM顆粒細觀參數(shù)和雙軸壓縮試驗獲取的相應(yīng)巖土體抗剪強度。參考細觀參數(shù)的取值范圍[20]，F(xiàn)s0設(shè)置5個水平，Xb設(shè)置3個水平，μp、Kn和Ks統(tǒng)一設(shè)置4個水平，各水平取值見表1。全面正交組合共計形成960組細觀參數(shù)進行雙軸壓縮數(shù)值試驗。顆粒密度ρ取值1 900 kg/m3，顆粒直徑取2.5～3.6 mm。

表1 細觀參數(shù)水平設(shè)定Tab.1 Level setting of micro parameters

將960組MatDEM顆粒細觀參數(shù)和雙軸壓縮數(shù)值試驗獲取的抗剪強度值作為神經(jīng)網(wǎng)絡(luò)樣本數(shù)據(jù)庫，執(zhí)行BР神經(jīng)網(wǎng)絡(luò)的學(xué)習訓(xùn)練。首先，對樣本數(shù)據(jù)進行歸一化處理，避免各項數(shù)據(jù)量綱/量級的不一致對結(jié)果產(chǎn)生影響。歸一化計算公式為：

式中，X、Xnorm分別為樣本歸一化之前和歸一化之后的數(shù)據(jù)，Xmin和Xmax分別為數(shù)據(jù)最小值和最大值。

MatDEM是在Matlab框架下編制的，其可便捷地在Matlab中進行二次開發(fā)；同時，Matlab神經(jīng)網(wǎng)絡(luò)工具箱又為BР神經(jīng)網(wǎng)絡(luò)的學(xué)習訓(xùn)練提供了大量的預(yù)置函數(shù)，因此，所有程序部分可在Matlab中實現(xiàn)有效集成。BР神經(jīng)網(wǎng)絡(luò)學(xué)習訓(xùn)練時，訓(xùn)練函數(shù)選用trainlm函數(shù)，輸出函數(shù)選用pureline函數(shù)，激活函數(shù)選擇tansig函數(shù)，誤差函數(shù)選用均方誤差RMSE。

2.2 細觀力學(xué)參數(shù)的逆向迭代修正標定

BР神經(jīng)網(wǎng)絡(luò)訓(xùn)練完畢后，式（5）中的所有權(quán)重和偏置信息均確定，其事實上也定義了黏聚力、內(nèi)摩擦角與細觀參數(shù)之間的復(fù)雜函數(shù)關(guān)系：

式中，f1和f2代表訓(xùn)練好的BР神經(jīng)網(wǎng)絡(luò)中由細觀參數(shù)逐步運算至黏聚力和內(nèi)摩擦角的過程，其直觀表達式由于權(quán)重、偏置項和激活函數(shù)等的引入很難直接給出，但其真實存在。

在BР神經(jīng)網(wǎng)絡(luò)框架下，當給定輸入信息時，可進行該輸入信息的梯度計算，因此，仿照非線性方程數(shù)學(xué)求解，可基于式（7）和（8）進行給定抗剪強度目標值下細觀參數(shù)取值的數(shù)學(xué)推演。首先，給予一組細觀參數(shù)初值，使用前向傳播獲取抗剪強度值，將該值與目標值的誤差定義為：

式中，C0和φ0分別為黏聚力和內(nèi)摩擦角的目標值。設(shè)定誤差限為RMSE0，若誤差不滿足要求，沿誤差梯度下降方向調(diào)整細觀參數(shù)，直至誤差滿足要求，如圖6所示。圖6中：n為循環(huán)次數(shù)；lr為學(xué)習率，其引進目的為執(zhí)行沿誤差梯度下降方向的步長搜索逼近，本文取0.01。區(qū)別于周瑜等[27]的做法，本文的BР神經(jīng)網(wǎng)絡(luò)以顆粒細觀參數(shù)為輸入，以巖土體抗剪強度為輸出，通過反推的方式獲取巖土體抗剪強度給定后的顆粒細觀參數(shù)，故稱之為逆向標定。

圖6 BP神經(jīng)網(wǎng)絡(luò)逆向標定細觀參數(shù)Fig.6 Micro parameters determined by BP neural network

對于逆向標定得到的細觀參數(shù)，執(zhí)行雙軸壓縮數(shù)值試驗，檢查抗剪強度與給定值的吻合程度，即精度檢查。精度檢查的結(jié)果可能并不理想，原因在于BР神經(jīng)網(wǎng)絡(luò)樣本庫質(zhì)量未能達到要求。訓(xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)前向傳播過程可視為雙軸壓縮數(shù)值試驗的一種替代，樣本庫的樣本數(shù)量越多，其替代的可靠程度越高。為保障效率，神經(jīng)網(wǎng)絡(luò)的樣本庫的樣本數(shù)量有限，故前向傳播過程替代雙軸壓縮數(shù)值試驗的可靠度未必足夠，導(dǎo)致逆向標定所得細觀參數(shù)精度不足。

對于非線性方程求解而言，已訓(xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)雖然精度不足，但提供了粗略的搜索方向，因此，逆向標定獲取的顆粒細觀參數(shù)將比給定的細觀參數(shù)初值更加趨近于最終解。鑒于此，當精度檢查結(jié)果不理想時，將逆向標定得到的細觀參數(shù)與其雙軸壓縮數(shù)值試驗獲取的抗剪強度作為一組新樣本添加到已有樣本庫中，重新進行BР神經(jīng)網(wǎng)絡(luò)的學(xué)習訓(xùn)練，此過程可稱為“樣本修正”。BР神經(jīng)網(wǎng)絡(luò)經(jīng)過重新學(xué)習訓(xùn)練后，式（7）和（8）中的函數(shù)f1和f2將相應(yīng)地調(diào)整變化，以前次標定獲得的細觀參數(shù)為初值，再次進行逆向標定和精度檢查。重復(fù)執(zhí)行逆向標定—精度檢查—樣本修正的過程，直至精度檢查順利通過，則顆粒細觀參數(shù)的標定過程完成，整個流程如圖7所示。

圖7 細觀參數(shù)標定的完整流程Fig.7 Whole calibration process of micro parameters

該流程中，每一輪的逆向標定均以上一輪標定獲得的細觀參數(shù)為初值，故顆粒細觀參數(shù)的最終取值獲取類似于數(shù)值求解中的迭代方式，因此，將這一總體流程稱為顆粒細觀參數(shù)的逆向迭代修正標定。

2.3 細觀參數(shù)標定策略的有效性測試

為測試逆向迭代修正標定策略的有效性，給定抗剪強度目標值C=3.0 kРa、φ=19.6°，執(zhí)行標定流程以獲取顆粒細觀參數(shù)。經(jīng)過4輪逆向標定后，細觀參數(shù)精度檢查通過，每輪逆向標定獲取的細觀參數(shù)和雙軸壓縮數(shù)值試驗獲取的抗剪強度值見表2，抗剪強度與目標值的誤差如圖8所示。

表2 每輪標定得到的細觀參數(shù)及測定的抗剪強度Tab.2 Micro parameters obtained in each round and its resulted shear strength

圖8 兩種神經(jīng)網(wǎng)絡(luò)誤差對比Fig.8 Comparison on results of the two neural networks

以抗剪強度給定值為標準，細觀參數(shù)逆向迭代修正標定執(zhí)行過程中，第1輪標定的黏聚力相對誤差為27.52%，摩擦系數(shù)相對誤差為59.70%；迭代1次后黏聚力相對誤差降至4.19%，摩擦系數(shù)相對誤差降至15.21%；4次迭代后，黏聚力相對誤差降至0.70%，摩擦系數(shù)相對誤差降至0.36%。因此，該方法在顆粒細觀參數(shù)的標定上具有可靠的數(shù)值精度。同時，必須注意的是，前面2輪標定得到的細觀參數(shù)精度不足。但是，逐輪的樣本修正后，標定結(jié)果的數(shù)值精度將迅速提升。其原因在于，如第2.2節(jié)所述，每輪逆向標定獲取的顆粒細觀參數(shù)將比上輪標定獲取的細觀參數(shù)更趨近于最終解，那么，樣本修正可理解為在最終解附近區(qū)域?qū)ι窠?jīng)網(wǎng)絡(luò)不斷加密。

為說明BР神經(jīng)網(wǎng)絡(luò)解決數(shù)學(xué)欠定問題時的局限性，圖8補充繪制了欠定BР神經(jīng)網(wǎng)絡(luò)輸出結(jié)果的誤差情況。該欠定BР神經(jīng)網(wǎng)絡(luò)以抗剪強度為輸入，以顆粒細觀參數(shù)為輸出，隱藏層層數(shù)和每層神經(jīng)元的布置仍如圖5所示。采用同樣的數(shù)據(jù)樣本庫，BР神經(jīng)網(wǎng)絡(luò)訓(xùn)練完畢后，直接輸入抗剪強度給定值得到顆粒細觀參數(shù)。用欠定BР神經(jīng)網(wǎng)絡(luò)輸出的顆粒細觀參數(shù)進行雙軸壓縮模型試驗，黏聚力的相對誤差為44.67%，摩擦系數(shù)相對誤差為45.85%。

3 雙軸壓縮數(shù)值模型并行試驗技術(shù)

BР神經(jīng)網(wǎng)絡(luò)標定細觀參數(shù)的全過程中，數(shù)據(jù)樣本庫建立是一項基礎(chǔ)而繁瑣的工作。為加快樣本庫建立過程，開發(fā)雙軸壓縮數(shù)值模型并行試驗技術(shù)，以充分發(fā)揮MatDEM中GРU計算模式的能力[30]。Mat-DEM中有GРU和CРU兩種計算模式，執(zhí)行不同功能時調(diào)用不同的計算模式。一般來說，數(shù)值模型建立和材料參數(shù)賦值為CРU計算模式。受力后的顆粒平衡迭代計算則會根據(jù)顆粒數(shù)量在GPU和CPU兩種模式中自動切換，單元數(shù)量較少時選擇CPU模式，單元數(shù)量增加到一定程度則選擇GPU模式。單個雙軸壓縮數(shù)值模型試驗盒中，顆粒數(shù)量一般在5 000左右，無法發(fā)揮GPU模式計算優(yōu)勢。

本文對MatDEM進行二次開發(fā)，實現(xiàn)雙軸壓縮數(shù)值模型的并行試驗；同時，構(gòu)建多個雙軸壓縮模型，每個模型的細觀參數(shù)單獨賦值，壓力板位移同步施加，顆粒受力后的平衡迭代同步進行，壓力板數(shù)據(jù)同步記錄。此時，參與同步計算的顆粒數(shù)量將達數(shù)十萬，利用專業(yè)工作站的高性能GPU，可有效提高計算效率。圖9為54個雙軸壓縮數(shù)值模型的排列分布情況和并行試驗獲取的試樣破壞時的位移云圖。

圖9 54個雙軸壓縮模型的排列分布和位移云圖Fig.9 Distribution and displacement diagrams of 54 biaxial compression models

試驗研究發(fā)現(xiàn)，雙軸壓縮數(shù)值模型并行試驗所需要的時間一方面取決于GPU的計算能力，另一方面取決于雙軸壓縮數(shù)值模型的排列分布。當雙軸壓縮數(shù)值模型的排列分布呈現(xiàn)圖9所示的正方形布局時，并行試驗的效率最高。其原因在于，以顆粒離散元為理論基礎(chǔ)的MatDEM，執(zhí)行逐時間步的增量計算，必然涉及動力系統(tǒng)的能量平衡和耗散問題。默認阻尼下，MatDEM數(shù)值模型完成一次標準平衡，即應(yīng)力波從模型一端傳至另一端所需的時間Ttotal為：

式中：Npack表示模型最長維度上堆積的顆粒數(shù)；ΔT為動力計算采用的時間增量步步長，取值需遠小于彈簧簡諧振動周期，以精確模擬數(shù)值模型彈性變形過程。多個雙軸壓縮數(shù)值模型并行試驗時，Mat-DEM會將同步試驗的數(shù)值模型視為整體。雙軸壓縮數(shù)值模型總個數(shù)一定時，采用正方形布局時，整體系統(tǒng)在兩個維度上的顆粒數(shù)基本接近，效率最高。

硬件配置Intel（R） Xeon（R） E-2224@3.40 GHz CPU，NIVIDIA Quadro RTX 4000顯卡的情況下，單個雙軸壓縮模型的試驗耗時約50 min。同樣配置條件下，54個雙軸壓縮模型并行試驗總耗時約480 min，節(jié)約82%的計算時間。因此，采用并行試驗技術(shù)進行雙軸壓縮數(shù)值模型試驗，可充分發(fā)揮GPU計算模式的能力。

4 邊坡穩(wěn)定性算例分析

以澳大利亞計算機應(yīng)用協(xié)會（簡稱ACADS）的兩個邊坡經(jīng)典考題為例，采用MatDEM執(zhí)行邊坡穩(wěn)定性的強度折減法分析，考察本文提出的顆粒細觀參數(shù)標定策略的實用性。穩(wěn)定性分析時，臨界破壞準則參照田雷等[6]的做法，記邊坡模型中所有顆粒的位移平均值為模型平均位移S，以其發(fā)生突變?yōu)檫吰逻_到極限狀態(tài)的判據(jù)。

4.1 算例1：均質(zhì)土質(zhì)邊坡

圖10為均質(zhì)土質(zhì)邊坡，土體力學(xué)參數(shù)為C=3.0 kPa，φ=19.6°，重度γ=20.0 kN/m3。ACADS給出的安全系數(shù)推薦值為1.00，對應(yīng)的滑動范圍如圖10橙色區(qū)域所示。

圖10 算例1模型及ACADS推薦解的滑動范圍Fig.10 Model of Example 1 and the failure region by ACADS

在MatDEM中建立該邊坡的顆粒離散元模型如圖11所示，模型含82 698個直徑67～96 mm的顆粒。此時，根據(jù)特征長度比L/R要求[25]將用于測定抗剪強度的雙軸壓縮數(shù)值模型尺寸調(diào)整為4 000 mm×8 000 mm，以避免顆粒尺寸效應(yīng)影響細觀參數(shù)標定。

圖11 邊坡的MatDEM模型和臨界失穩(wěn)狀態(tài)位移Fig.11 MatDEM model and its critical displacements

分析穩(wěn)定性時，初始折減系數(shù)設(shè)為0.70；若邊坡處于穩(wěn)定狀態(tài)，則將折減系數(shù)增加0.10繼續(xù)計算，折減系數(shù)達到1.10時邊坡失穩(wěn)破壞。在折減系數(shù)0.95到1.10范圍內(nèi)進一步細化分析，最終確定邊坡穩(wěn)定系數(shù)為1.02。邊坡模型平均位移S隨著強度折減系數(shù)Fsr增大的變化曲線如圖12所示。由圖12可知：Fsr達到1.02前，S基本呈水平直線；數(shù)據(jù)檢查發(fā)現(xiàn)，S隨著Fsr增加仍略有增加，表明抗剪強度下降會影響顆粒平衡后的位置，宏觀上體現(xiàn)為變形輕微增長；Fsr達到1.02后，S值隨Fsr增加而急劇變大，表明邊坡進入臨界失穩(wěn)狀態(tài)。邊坡達到臨界失穩(wěn)狀態(tài)時的位移，如圖11所示。此時邊坡自坡腳到坡頂已出現(xiàn)明顯貫通的滑動面，邊坡沿該滑動面將開始發(fā)生滑動。

圖12 算例1的S-Fsr曲線Fig.12 S-Fsr curves for Example 1

強度折減過程中，抗剪強度不同折減系數(shù)折減后對應(yīng)的細觀參數(shù)見表3。由表3可知，由于抗剪強度與細觀參數(shù)之間的非線性映射關(guān)系非常復(fù)雜，抗剪強度的折減很難通過直接對細觀參數(shù)進行折減來實現(xiàn)。

表3 不同折減系數(shù)下土體抗剪強度與細觀力學(xué)參數(shù)Tab.3 Micro parameters and soil shear strength corresponding to reduction coefficients

4.2 算例2：堆石壩模型

算例2為一由堆石體、過渡層、心墻和反濾層組成的堆石壩，壩高146 m，壩底寬648 m，其形態(tài)如圖13所示。

圖13 堆石壩模型及ACADS推薦解的滑動區(qū)域Fig.13 A rock fill dam and the failure region by ACADS

壩體模型各種材料的力學(xué)參數(shù)見表4，由于算例2是用于穩(wěn)定性分析，其參數(shù)與實際工程中各種材料的參數(shù)取值存在一定差異。ACADS對該堆石壩模型給出的安全系數(shù)推薦值為1.976，對應(yīng)的滑動面如圖13所示。

表4 堆石壩模型各材料力學(xué)參數(shù)Tab.4 Mechanical properties of materials in rockfill dam model

在MatDEM中建立該模型的顆粒離散元模型如圖14所示，模型含89 768個直徑為0.33～0.48 m的顆粒。顆粒細觀參數(shù)標定時，按照特征長度比L/R的要求[25]，將雙軸壓縮數(shù)值模型尺寸調(diào)整為20 m×40 m。

圖14 堆石壩的MatDEM模型Fig.14 MatDEM model for the rock fill dam

采用強度折減法分析該堆石壩模型的穩(wěn)定性，模型的平均位移S和強度折減系數(shù)Fsr關(guān)系曲線如圖15（a）所示。以S突變?yōu)槟Ｐ团R界失穩(wěn)判據(jù)，確定該堆石壩安全系數(shù)為1.975。模型達到臨界失穩(wěn)狀態(tài)時的位移如圖15（b）所示。圖15（a）顯示，即便模型處于穩(wěn)定狀態(tài)，模型平均位移S已超過2 m。其原因在于，顆粒離散元邊坡模型建立完畢之后，由于存在重力勢能，顆粒自平衡過程中會導(dǎo)致模型產(chǎn)生一定的變形，類似于重力作用下的沉降變形。事實上，算例1的邊坡模型處于穩(wěn)定狀態(tài)時，也存在0.002 m的模型平均位移S。算例2模型平均位移S偏大，其原因主要為，施加于模型上的邊界約束條件，堆石壩模型僅在底部邊界設(shè)置位移約束，而算例1邊坡模型在底部和左右兩側(cè)邊界均設(shè)置位移約束，因此堆石壩模型重力效應(yīng)可充分發(fā)揮，變形相對較大；此外，堆石壩模型高度達146 m，大尺寸規(guī)模下，顆粒自平衡過程中的模型變形也相對較大。

圖15 堆石壩的安全系數(shù)結(jié)果和臨界失穩(wěn)狀態(tài)位移Fig.15 Safety factors and the critical displacements for therock fill dam

強度折減過程中，材料顆粒細觀參數(shù)隨著折減系數(shù)變化而調(diào)整的情況見表5。

表5 不同折減系數(shù)下堆石壩模型中材料的抗剪強度與細觀力學(xué)參數(shù)Tab.5 Micro parameters and shear strength of material in rockfill dam model corresponding to reduction coefficients

兩個算例的執(zhí)行情況顯示，本文提出的細觀參數(shù)標定策略，可滿足MatDEM等顆粒離散元軟件采用強度折減法分析邊坡穩(wěn)定性時，迅速調(diào)整顆粒細觀參數(shù)取值，以對應(yīng)巖土體抗剪強度折減的需求。安全系數(shù)結(jié)果方面：算例1安全系數(shù)為1.02，與ACADS推薦的安全系數(shù)相差0.02；算例2安全系數(shù)為1.975，與ACADS推薦的系數(shù)相差0.001?；瑒用娼Y(jié)果方面：MatDEM獲取的臨界失穩(wěn)狀態(tài)下模型的滑動區(qū)域與ACADS推薦的滑動范圍基本吻合?？傮w上，MatDEM獲取的安全系數(shù)和滑動范圍與ACADS推薦解之間具有良好的可比性，表明采用本文策略標定的顆粒細觀參數(shù)，在對應(yīng)不同程度折減后的巖土體抗剪強度方面具有滿意的數(shù)值精度。因此，在應(yīng)用顆粒離散元強度折減法評價邊坡穩(wěn)定性時，本文提出的顆粒細觀參數(shù)逆向迭代修正標定策略可發(fā)揮積極作用。

5 結(jié) 論

針對以MatDEM為代表的顆粒離散元法軟件應(yīng)用強度折減技術(shù)評價邊坡穩(wěn)定性時的顆粒細觀參數(shù)標定問題，以顆粒細觀參數(shù)為輸入，抗剪強度為輸出構(gòu)造BР神經(jīng)網(wǎng)絡(luò)，建立顆粒細觀參數(shù)的逆向迭代修正標定方法，開發(fā)雙軸壓縮數(shù)值模型并行試驗技術(shù)以加快神經(jīng)網(wǎng)絡(luò)樣本庫樣本的獲取，并采用ACADS的兩個邊坡穩(wěn)定性分析經(jīng)典考題驗證了顆粒細觀參數(shù)逆向迭代修正標定方法的能力和精度?；谘芯窟^程中的發(fā)現(xiàn)，可得結(jié)論如下：

1）細觀參數(shù)逆向迭代修正標定方法可有效解決顆粒離散元強度折減法分析邊坡穩(wěn)定性時的顆粒細觀參數(shù)標定問題。標定得到的細觀參數(shù)可良好對應(yīng)巖土體的抗剪強度，數(shù)值算例計算結(jié)果與推薦解之間的可比性好。

2）MatDEM的二次開發(fā)在Matlab程序?qū)崿F(xiàn)便捷，使得MatDEM具備很好的研究潛力和應(yīng)用價值。本文開發(fā)的雙軸壓縮數(shù)值模型并行試驗技術(shù)有效加速了神經(jīng)網(wǎng)絡(luò)樣本庫的建立過程，可為其他小型模型的重復(fù)試驗工作提供參考。

3）使用BР神經(jīng)網(wǎng)絡(luò)執(zhí)行參數(shù)反演問題研究時，必須注意神經(jīng)網(wǎng)絡(luò)輸入信息和輸出信息維度的大小關(guān)系，以避免形成數(shù)學(xué)上的欠定問題。