紅砂巖單軸壓縮宏細(xì)觀參數(shù)映射關(guān)系研究1)

佟 安 張軍徽 武 娜

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 彈性模量敏感性分析

圖3 剛度比敏感性分析

2.2 抗拉強(qiáng)度、黏聚力敏感性分析

從初始值開(kāi)始依次向右改變表2 中的C取值，其他細(xì)觀參數(shù)取值不變，同理改變?chǔ)业娜≈?。分析黏聚力、抗拉?qiáng)度對(duì)應(yīng)力–應(yīng)變曲線的影響規(guī)律及程度。

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

圖4 黏聚力敏感性分析

圖5 抗拉強(qiáng)度敏感性分析

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

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

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

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

3 建立宏細(xì)觀參數(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)勢(shì)在于能夠存儲(chǔ)與學(xué)習(xí)數(shù)量龐大的輸入輸出映射關(guān)系，且無(wú)需映射關(guān)系之間有著明確的數(shù)學(xué)表達(dá)式。BP 神經(jīng)網(wǎng)絡(luò)的應(yīng)用已十分廣泛，其原理此處不多贅述。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 模擬曲線與原模擬曲線對(duì)比

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

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

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

圖10 PFC 模擬曲線與試驗(yàn)曲線對(duì)比

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

圖11 PFC 試樣與試件破壞形式對(duì)比

5 UCS/TS值驗(yàn)證

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

圖12 PFC 單軸壓縮曲線

圖13 PFC 拉伸曲線

6 結(jié)論與展望

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

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

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

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

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