佟 安 張軍徽 武 娜

(北方工業(yè)大學(xué)土木工程學(xué)院，北京100144)

近年來，隨著我國城市化水平的提高，城市地下空間開發(fā)利用得到大幅增長，同時也不可避免地發(fā)生著各類地質(zhì)災(zāi)害。研究巖石的破壞機制對于地下工程施工和地質(zhì)災(zāi)害治理等具有極其重要的意義[1]。經(jīng)典的連續(xù)介質(zhì)方法建立的巖石強度理論，難以解釋巖石自身表現(xiàn)出來的非連續(xù)性(存在節(jié)理裂隙)、大變形以及各向異性等特征，以“均勻連續(xù)、小變形”為基本假設(shè)，建立的有限元軟件不再適用[2]。1971 年，Cundall 最早提出適用于巖石特性的離散元法，1986年，王泳嘉將離散元方法引入國內(nèi)，并不斷發(fā)展至今。顆粒流分析程序(particle flow code，PFC)是一種基于離散元方法開發(fā)的用于模擬圓形顆粒介質(zhì)的運動及其相互作用的數(shù)值分析軟件，可以很好地模擬巖石的宏觀力學(xué)行為[3]。

目前，PFC存在的主要問題是細觀參數(shù)難以標定，原因在于極難獲得顆粒細觀參數(shù)與室內(nèi)試驗獲取的宏觀力學(xué)參數(shù)之間的關(guān)系，兩者之間存在諸多難以控制的因素，并且這些因素之間的非線性關(guān)系極為明顯[4]。多位學(xué)者使用“試湊法”解決此問題，但是盲目性較大且耗費大量時間。

反向傳播算法(back propagation，BP)神經(jīng)網(wǎng)絡(luò)具有優(yōu)秀的非線性動態(tài)處理能力，可實現(xiàn)高度非線性映射，具有較強的學(xué)習(xí)、存儲、計算能力和容錯特性，能實現(xiàn)輸入樣本和輸出樣本的非線性映射。本文提出以PFC2D 為基礎(chǔ)獲取宏細觀參數(shù)樣本，利用BP 神經(jīng)網(wǎng)絡(luò)建立宏細觀參數(shù)映射關(guān)系的方式代替“試湊法”。

1 PFC2D砂巖單軸壓縮試驗?zāi)M

PFC2D 是一種常被用于模擬巖石力學(xué)行為的二維離散元軟件，使用剛性圓盤模擬巖樣中的顆粒[5]。在選擇接觸本構(gòu)模型時，平行粘結(jié)模型是最常用的接觸本構(gòu)模型。但此本構(gòu)模型存在著明顯的問題，其建立的巖石模型不能同時匹配抗壓強度(uniaxia compressive strength，UCS)和抗拉強度(tensile strength，TS)，即UCS/TS值(一般為5左右)小于巖石的強度比值(10～20)[6]，原因在于此模型構(gòu)造圓形顆粒進行模擬，不能提供足夠的抗扭轉(zhuǎn)力，尤其是在模型產(chǎn)生裂縫后。平直節(jié)理模型可以解決此類問題，由圖1 可以看出，該模型將圓形顆粒構(gòu)造成多邊形顆粒，以此抑制顆粒破壞后的旋轉(zhuǎn)，使得UCS/TS值顯著增大。

平直節(jié)理模型是一種描述有限尺寸、線彈性且考慮局部損傷的粘結(jié)模型，將交界面進行離散化，每個交界面都有粘結(jié)和未粘結(jié)2種狀態(tài)，見圖1。本文選用平直節(jié)理模型用于模擬紅砂巖，試件大小選取50 mm×100 mm 的標準試樣。為減少計算時間，力學(xué)模型盡量簡化，由于巖石結(jié)構(gòu)的復(fù)雜性，無法做到模擬顆粒與真實顆粒大小一致，當顆粒數(shù)超過3000后，模型力學(xué)性能基本保持穩(wěn)定[7]，本文的PFC 模型顆?？倲?shù)為3602 個，滿足要求?；A(chǔ)參數(shù)選取參考文獻[8]，具體取值見表1。

圖1 平直節(jié)理模型

表1 PFC 數(shù)值試樣的基礎(chǔ)細觀參數(shù)

隨著時間步長的增加，上下加載板勻速移動，擠壓試件，當加載至應(yīng)力下降到峰值應(yīng)力的70%時，停止加載，此時砂巖試件已明顯破壞，試驗結(jié)束。

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

本文針對紅砂巖抗壓強度偏小的特點，選取合適的細觀參數(shù)范圍進行單軸壓縮試驗的敏感性分析，來確定各個細觀參數(shù)對試驗的影響規(guī)律及程度，從而分類篩選出對砂巖單軸壓縮應(yīng)力–應(yīng)變曲線影響程度較大的主要細觀參數(shù)，為下一步選取BP 神經(jīng)網(wǎng)絡(luò)參數(shù)做準備工作。經(jīng)過大量的敏感性分析，篩選出以下幾個主要的參數(shù)：顆粒有效模量(E?)、顆粒剛度比(κ?)、平直節(jié)理有效模量平直節(jié)理剛度比抗拉強度(σ)和黏聚力(C)、平直節(jié)理摩擦系數(shù)(μ)，根據(jù)文獻[9]，將與值保持一致，取值列于表2中。

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

2.1 彈性模量、剛度比敏感性分析

本文通過控制單一變量實現(xiàn)敏感性分析，保持表2 中各個初始值不變，從左至右依次同時由初始值改變的數(shù)值。同理，保持其他數(shù)值不變，調(diào)整分別分析接觸彈性模量、剛度比對應(yīng)力–應(yīng)變曲線的影響規(guī)律及程度。

由圖2可知，隨著彈性模量的增加，單位應(yīng)變條件下應(yīng)力增長幅度呈減小趨勢，表現(xiàn)出應(yīng)力–應(yīng)變曲線初始端越緩、斜率越小即宏觀的彈性模量減小，但彈性模量的增加對峰值應(yīng)力影響不大。同理，由圖3可知，隨著剛度比的增加，應(yīng)力–應(yīng)變曲線的彈性模量和峰值應(yīng)力均增大。彈性模量與剛度比的變化對試樣的應(yīng)力–應(yīng)變曲線作用顯著。

圖2 彈性模量敏感性分析

圖3 剛度比敏感性分析

2.2 抗拉強度、黏聚力敏感性分析

從初始值開始依次向右改變表2 中的C取值，其他細觀參數(shù)取值不變，同理改變σ的取值。分析黏聚力、抗拉強度對應(yīng)力–應(yīng)變曲線的影響規(guī)律及程度。

由圖4 和圖5 可知，黏聚力和抗拉強度的敏感性分析結(jié)果類似，隨著數(shù)值的增加，峰值強度顯著增大，而彈性模量基本保持不變。黏聚力、抗拉強度對試樣的應(yīng)力–應(yīng)變曲線作用顯著。

圖4 黏聚力敏感性分析

圖5 抗拉強度敏感性分析

2.3 摩擦系數(shù)敏感性分析

從初始值開始依次向右改變表2 中μ的取值，其他細觀參數(shù)取值不變。分析μ對模型應(yīng)力–應(yīng)變曲線的影響規(guī)律及程度。

從圖6 中可以看出，平直節(jié)理摩擦系數(shù)對應(yīng)力–應(yīng)變曲線的峰值應(yīng)力和彈性模量影響并不明顯，因此，將此變量取固定值，不作為BP 神經(jīng)網(wǎng)絡(luò)的輸出參數(shù)。

圖6 摩擦系數(shù)敏感性分析

3 建立宏細觀參數(shù)映射關(guān)系

3.1 BP神經(jīng)網(wǎng)絡(luò)基本原理

本文采用的BP 神經(jīng)網(wǎng)絡(luò)是一種按誤差逆向傳播算法訓(xùn)練得到的多層前饋網(wǎng)絡(luò)，在巖土工程領(lǐng)域應(yīng)用較為廣泛[10]。此方法的優(yōu)勢在于能夠存儲與學(xué)習(xí)數(shù)量龐大的輸入輸出映射關(guān)系，且無需映射關(guān)系之間有著明確的數(shù)學(xué)表達式。BP 神經(jīng)網(wǎng)絡(luò)的應(yīng)用已十分廣泛，其原理此處不多贅述。

3.2 BP神經(jīng)網(wǎng)絡(luò)建立的總體流程

從PFC2D 模擬紅砂巖單壓縮試驗得到的應(yīng)力–應(yīng)變曲線上獲取的宏觀力學(xué)參數(shù)與細觀參數(shù)建立BP 神經(jīng)網(wǎng)絡(luò)關(guān)系，主要分為6 個工作流程：(1)獲取與細觀參數(shù)相對應(yīng)的宏觀力學(xué)參數(shù)；(2)宏觀力學(xué)參數(shù)作為輸入樣本，細觀參數(shù)作為輸出樣本，用于訓(xùn)練BP神經(jīng)網(wǎng)絡(luò)；(3) 調(diào)整BP神經(jīng)網(wǎng)絡(luò)程序，使宏細觀參數(shù)建立精確的映射關(guān)系；(4)建立BP 神經(jīng)網(wǎng)絡(luò)映射關(guān)系后，輸入專門用于測試的宏觀力學(xué)參數(shù)，反演相應(yīng)的細觀參數(shù)；(5)將反演出的細觀參數(shù)輸入到PFC2D中，進行數(shù)值模擬試驗，并將數(shù)值模擬計算所得到的宏觀力學(xué)參數(shù)與第4 步輸入的宏觀力學(xué)參數(shù)進行對比，驗證BP 神經(jīng)網(wǎng)絡(luò)的精準性；(6)不斷重復(fù)上述工作流程，直到建立出較為精準高效的BP神經(jīng)網(wǎng)絡(luò)。

3.3 BP神經(jīng)網(wǎng)絡(luò)的創(chuàng)建

首先，選取建立網(wǎng)絡(luò)的宏細觀參數(shù)。采用“試湊法”基于模擬相似應(yīng)力–應(yīng)變曲線的原理，試用4 個宏觀參數(shù)描述應(yīng)力–應(yīng)變曲線作為網(wǎng)絡(luò)的輸入樣本，分別為：彈性模量、峰值處的應(yīng)力、應(yīng)變值以及破壞斜率，如圖7。此4 個宏觀參數(shù)是針對PFC2D 模擬紅砂巖單軸壓縮試驗應(yīng)力–應(yīng)變曲線的特點總結(jié)，其彈性段和最終破壞段，都具有相當明顯的線性增長(下降)的趨勢，便于精確獲取且利于網(wǎng)絡(luò)尋找宏細觀參數(shù)映射關(guān)系。使用6 個細觀參數(shù)作為神經(jīng)網(wǎng)絡(luò)的輸出樣本，分別為：

圖7 基于應(yīng)力–應(yīng)變曲線選取宏觀力學(xué)參數(shù)

利用PFC2D內(nèi)置的FISH以及Python語言，二次開發(fā)紅砂巖單軸壓縮程序，使其可以自動在6 組細觀參數(shù)的取值范圍內(nèi)產(chǎn)生隨機數(shù)，見表3。并記錄相應(yīng)的宏細觀參數(shù)，程序自動運行40 次作為網(wǎng)絡(luò)輸入輸出的樣本庫。

表3 PFC 細觀參數(shù)取值范圍

由于用于訓(xùn)練的輸入樣本和輸出樣本不同，需要進行歸一化處理，便于快速尋找映射關(guān)系，歸一化公式如式(1)所示。利用此式將訓(xùn)練樣本數(shù)據(jù)，全部歸一化至[0,1]之間。

隨后創(chuàng)建BP 神經(jīng)網(wǎng)絡(luò)需要考慮網(wǎng)絡(luò)層數(shù)及中間層神經(jīng)元個數(shù)，一般情況下具有一個中間層的BP神經(jīng)網(wǎng)絡(luò)可以很好地模擬任何有理函數(shù)，故本文選取具有一個中間層的三層BP神經(jīng)網(wǎng)絡(luò)。4個宏觀力學(xué)參數(shù)作為輸入層樣本參數(shù)，6 個細觀參數(shù)作為輸出層樣本參數(shù)，其中取相同值，因此，輸出層樣本參數(shù)為4 個。最終，輸入節(jié)點為4，輸出層節(jié)點數(shù)為4。

中間層神經(jīng)元個數(shù)的選取對整體性能影響較大，數(shù)目過少，網(wǎng)絡(luò)將不能建立復(fù)雜的判斷界，難以訓(xùn)練出合理的網(wǎng)絡(luò)，兼容性較差；數(shù)目過多會導(dǎo)致訓(xùn)練時間較長，網(wǎng)絡(luò)的泛化能力下降，計算誤差可能變大。因此，選取合適的神經(jīng)元個數(shù)至關(guān)重要，但目前并沒有確定的科學(xué)辦法選取神經(jīng)元個數(shù)，一般根據(jù)如下經(jīng)驗公式[11]來確定

其中，p為中間層神經(jīng)元個數(shù)，n為輸入層節(jié)點數(shù)。本文根據(jù)式(2)選神經(jīng)元個數(shù)為5個左右，在此范圍上下浮動5 個神經(jīng)元，選出期望值與反演輸出值均方差最小的神經(jīng)元個數(shù)，最終確定神經(jīng)元個數(shù)為9，如圖8所示。

圖8 神經(jīng)元個數(shù)選取對比

最后選取建立網(wǎng)絡(luò)的必要函數(shù)，網(wǎng)絡(luò)輸入層與隱藏層及隱藏層與輸出層之間的傳遞函數(shù)采用logsig 函數(shù)和purelin 函數(shù)；考慮到網(wǎng)絡(luò)的規(guī)模和學(xué)習(xí)時間，選用trainlm 函數(shù)，trainlm 函數(shù)使用Levenberg–marquardt算法，訓(xùn)練網(wǎng)絡(luò)。

3.4 BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練與反演

利用之前獲取的40 組宏細觀參數(shù)作為訓(xùn)練樣本集，訓(xùn)練次數(shù)為2000，訓(xùn)練的目標誤差為0.008，學(xué)習(xí)率設(shè)為0.05。將訓(xùn)練好的網(wǎng)絡(luò)用sim 函數(shù)進行反演模擬，并不斷調(diào)試。

3.5 BP神經(jīng)網(wǎng)絡(luò)反演能力驗證

再次使用3.3節(jié)中的方法，獲取3組不同于訓(xùn)練樣本集的宏細觀參數(shù)來驗證BP 神經(jīng)網(wǎng)絡(luò)的反演能力，如表4所示。

表4 BP 神經(jīng)網(wǎng)絡(luò)反演值與期望值對比

表4 中反演輸出值與期望值之間總體誤差保持在目標誤差0.01 以內(nèi)。并將第1 組反演輸出結(jié)果輸入至PFC2D中進行模擬，對比結(jié)果如圖9所示，可以看出兩曲線高度相似，表明創(chuàng)建的BP 神經(jīng)網(wǎng)絡(luò)具有較強的反演能力，且精度較高。

圖9 模擬曲線與原模擬曲線對比

4 BP神經(jīng)網(wǎng)絡(luò)的室內(nèi)試驗驗證

紅砂巖試樣取自武漢一采石場，嚴格按照規(guī)范[12]要求加工成50 mm×100 mm的長方體，并確保試樣沒有較大裂縫。本試驗選用吉林省金力試驗技術(shù)有限公司研制的YAW-2000 微機控制電液伺服壓力試驗機，加載方式為位移控制，加載速率為0.005 m/s,以上條件均滿足PFC2D的參數(shù)選取。試驗獲取多組應(yīng)力–應(yīng)變曲線，來驗證BP 神經(jīng)網(wǎng)絡(luò)的可行性。使用本文方法獲取試驗曲線的4 個宏觀參數(shù)作為輸入樣本，輸入到上述創(chuàng)建成功的BP 神經(jīng)網(wǎng)絡(luò)中，網(wǎng)絡(luò)的反演輸出得到的細觀參數(shù)如表5 所示。再將得到的細觀參數(shù)輸入至PFC2D 單軸壓縮程序中，模擬得到的應(yīng)力–應(yīng)變曲線與試驗曲線進行對比。從圖10中可以看出，PFC模擬曲線與試驗曲線十分接近(由于砂巖試件內(nèi)部結(jié)構(gòu)復(fù)雜，分布并不均勻，試驗曲線彈性段未能體現(xiàn)出線性增長)，進一步從數(shù)據(jù)上看，如表6 所示，輸出的宏觀參數(shù)值與室內(nèi)單軸壓縮試驗所獲取的宏觀參數(shù)值極為接近，同時如圖11所示兩破壞形式也很相近。由此證明，本文所創(chuàng)建的BP神經(jīng)網(wǎng)絡(luò)系統(tǒng)是合理可行的。

表5 BP 神經(jīng)網(wǎng)絡(luò)反演的PFC 細觀參數(shù)值

圖10 PFC 模擬曲線與試驗曲線對比

表6 PFC 模擬獲取的宏觀參數(shù)與單軸壓縮試驗參數(shù)對比

圖11 PFC 試樣與試件破壞形式對比

5 UCS/TS值驗證

圖12 為上述神經(jīng)網(wǎng)絡(luò)模擬室內(nèi)試驗最終識別出的紅砂巖單軸壓縮曲線，其UCS值為34.45 MPa。同樣將網(wǎng)絡(luò)輸出的細觀參數(shù)輸入至PFC2D 紅砂巖拉伸模擬試驗中，獲得的單軸拉伸應(yīng)力–應(yīng)變曲線如圖13 所示，其TS 值為3.18 MPa。UCS/TS 值約為11，在合理范圍內(nèi)，表明平直節(jié)理模型可以正常模擬出巖石UCS/TS 的真實值，彌補了平行粘結(jié)模型在此方面的缺陷。

圖12 PFC 單軸壓縮曲線

圖13 PFC 拉伸曲線

6 結(jié)論與展望

本文基于BP 神經(jīng)網(wǎng)絡(luò)并二次開發(fā)PFC2D 巖石單軸壓縮程序來解決參數(shù)難以標定的問題，在BP 神經(jīng)網(wǎng)絡(luò)中輸入應(yīng)力–應(yīng)變曲線上的宏觀力學(xué)參數(shù)，反演出細觀力學(xué)參數(shù)，再輸入至PFC2D 中進行模擬，得到與試驗曲線高度吻合的模擬曲線，得出以下結(jié)論：

(1)通過上述對比可以看出，BP 神經(jīng)網(wǎng)絡(luò)具有非線性映射能力、泛化能力等優(yōu)勢，可以精確地建立宏細觀參數(shù)映射關(guān)系，反演得到的應(yīng)力–應(yīng)變曲線與室內(nèi)試驗曲線高度吻合，且本文創(chuàng)建的網(wǎng)絡(luò)僅需要40組樣本數(shù)據(jù)，較“試湊法”節(jié)省大量時間，表明BP神經(jīng)網(wǎng)絡(luò)可以很好地解決細觀參數(shù)標定難題。

(2)本文選取的平直節(jié)理模型有效地解決平行粘結(jié)模型關(guān)于UCS/TS 過小問題，展現(xiàn)出平直節(jié)理模型模擬巖石的優(yōu)勢。

(3) PFC2D 與BP 神經(jīng)網(wǎng)絡(luò)的聯(lián)合使用效果顯著，但是過程中依舊存在一些不確定因素，比如神經(jīng)元個數(shù)的選取依舊難以尋找最優(yōu)解，目前只有經(jīng)驗公式可供選擇，有待解決。

(4) PFC2D 模擬紅砂巖的單軸壓縮試驗所得到的應(yīng)力–應(yīng)變曲線與室內(nèi)試驗獲得的曲線仍存在差別，彈性段為直線，過于理想，程序有待再度開發(fā)。