李新平 黃明智 王 剛 徐 坤 劉婷婷

(武漢理工大學(xué)土木工程與建筑學(xué)院，武漢430070)

基于離散單元法發(fā)展起來的顆粒流程序(particle flow code，PFC)，將巖、土等材料介質(zhì)離散為抽象的顆粒單元組成的集合體，通過賦予顆粒及顆粒間接觸細(xì)觀參數(shù)表征類似礦物晶體及晶體結(jié)構(gòu)間的黏結(jié)作用，結(jié)合顆粒間牛頓運(yùn)動定律和接觸間力與位移定律來模擬真實材料宏觀力學(xué)特性，可以從細(xì)觀角度展現(xiàn)裂紋孕育、擴(kuò)展直至貫通的時空演化規(guī)律，在巖石力學(xué)特性數(shù)值模擬仿真中已得到大量應(yīng)用[1-2]。

PFC中需定義顆粒接觸的本構(gòu)關(guān)系來模擬巖石宏觀力學(xué)性質(zhì)，在程序發(fā)展初期定義了兩種接觸模型：接觸黏結(jié)模型和平行黏結(jié)模型，但用這兩種顆粒黏結(jié)模型進(jìn)行數(shù)值模擬過程中發(fā)現(xiàn)兩者都存在著一定的缺陷，即模擬巖石材料室內(nèi)試驗所得出的單軸抗壓強(qiáng)度與單軸抗拉強(qiáng)度之比偏低，與大多數(shù)巖石材料實際值不符[3]。為解決這一問題，許多學(xué)者從不同的角度進(jìn)行了大量的改進(jìn)研究，如：Cho等[4]通過改進(jìn)平行黏結(jié)顆粒間接觸方式而建立簇平行黏結(jié)模型；Potyondy[5]則從另一角度出發(fā)，通過改變顆粒形狀研究出適用于硬質(zhì)巖石的平直節(jié)理接觸模型，在此接觸模型下巖石材料的單軸抗壓強(qiáng)度與單軸抗拉強(qiáng)度之比值得以明顯增大。

在顆粒流數(shù)值模擬中，顆粒和黏結(jié)模型代表了模型材料的本構(gòu)關(guān)系，其細(xì)觀參數(shù)決定了數(shù)值模型的宏觀力學(xué)參數(shù)，對于細(xì)觀參數(shù)的標(biāo)定，目前已有學(xué)者在不同的接觸模型下進(jìn)行了初步的研究分析。文獻(xiàn)[6-9]以平行黏結(jié)接觸模型為基礎(chǔ)，通過數(shù)理統(tǒng)計方法與常用數(shù)值模擬試驗探討了顆粒材料中細(xì)觀參數(shù)對宏觀參數(shù)的影響關(guān)系，并嘗試建立宏細(xì)觀參數(shù)間的線性表達(dá)式；文獻(xiàn)[10-11]則研究了平直節(jié)理接觸模型顆粒材料中細(xì)觀參數(shù)對宏觀參數(shù)的影響規(guī)律及其排序結(jié)果，并根據(jù)線性回歸分析得出了宏細(xì)觀參數(shù)之間近似的擬合關(guān)系。

上述研究初步揭示了巖石材料宏細(xì)觀參數(shù)之間復(fù)雜的映射關(guān)系，為細(xì)觀參數(shù)的準(zhǔn)確標(biāo)定提供了一定參考，但研究成果中所討論的主要力學(xué)參數(shù)偏少、未對宏觀力學(xué)參數(shù)內(nèi)聚力c和內(nèi)摩擦角φ值進(jìn)行標(biāo)定,且進(jìn)行驗證的數(shù)值模擬與室內(nèi)試驗較少，擬合出的宏細(xì)觀參數(shù)表達(dá)式不夠準(zhǔn)確且需多次迭代調(diào)整?；谝陨涎芯恐械牟蛔?，為了能更真實地反映出巖石的宏觀力學(xué)特性，本文以平直節(jié)理接觸模型為基礎(chǔ)，通過單軸壓縮、雙軸壓縮和直接拉伸數(shù)值模擬試驗探究巖石材料宏細(xì)觀參數(shù)之間的關(guān)系，進(jìn)而采用善于處理復(fù)雜非線性映射關(guān)系的神經(jīng)網(wǎng)絡(luò)算法對細(xì)觀參數(shù)進(jìn)行標(biāo)定研究，結(jié)合花崗巖室內(nèi)常規(guī)試驗對標(biāo)定結(jié)果的準(zhǔn)確性進(jìn)行驗證分析，為巖石力學(xué)數(shù)值仿真研究提供一定的參考依據(jù)。

1 PFC 2D基本理論

1.1 平直節(jié)理接觸模型

平直節(jié)理接觸模型與平行黏結(jié)模型、接觸黏結(jié)模型最大的區(qū)別在于能夠抑制顆粒黏結(jié)破壞后的旋轉(zhuǎn)，如圖1所示顆粒形狀構(gòu)造成多邊形并相互“咬合”在一起，當(dāng)顆粒黏結(jié)破壞后，由于顆粒之間的相互“咬合”作用，顆粒單元無法發(fā)生自由旋轉(zhuǎn)運(yùn)動，只能整體間滑動或脫落，與實際巖石材料內(nèi)部微觀構(gòu)造更為接近。

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

平直節(jié)理接觸模型有未黏結(jié)和黏結(jié)兩種模式，兩者的本構(gòu)關(guān)系不相同。對于未黏結(jié)部分，剪切強(qiáng)度為

式中，μ是摩擦系數(shù)。對于黏結(jié)部分，剪切強(qiáng)度為

式中，c和φ分別為顆粒間黏聚力和內(nèi)摩擦角。

1.2 巖石數(shù)值模擬試驗

在進(jìn)行PFC建模時要使模型能夠較真實地反映出巖石材料的力學(xué)特性，首要考慮的就是細(xì)觀參數(shù)標(biāo)定，以此實現(xiàn)數(shù)值模型與實際巖石材料特性相匹配，細(xì)觀參數(shù)標(biāo)定思路就是通過對比分析巖石的數(shù)值模擬試驗結(jié)果與相對應(yīng)的室內(nèi)試驗結(jié)果。與室內(nèi)試驗獲取試樣宏觀力學(xué)參數(shù)類似，在顆粒流中單軸壓縮數(shù)值模擬試驗(圖2)也是獲取顆粒材料宏觀參數(shù)的常用方法，根據(jù)模擬試驗結(jié)果可得到試樣在平面應(yīng)力狀態(tài)下的單軸抗壓強(qiáng)度σf、彈性模量E和泊松比ν等力學(xué)參數(shù)；直接拉伸數(shù)值模擬試驗(圖3)可直接獲取模型試樣的抗拉強(qiáng)度σt，圖中黑色部分代表試樣內(nèi)部微裂紋分布。

圖2 單軸壓縮

圖3 直接拉伸

除以上四個基本力學(xué)參數(shù)外，為確定巖石材料另外兩個重要力學(xué)參數(shù)內(nèi)聚力c和內(nèi)摩擦角φ，可通過如圖4所示的雙軸壓縮數(shù)值模擬試驗(對應(yīng)巖石常規(guī)三軸壓縮試驗)，其試驗原理為：首先生成四道邊界墻及離散元顆粒模型，待顆粒達(dá)到平衡后固定所有邊界墻并通過伺服系統(tǒng)程序?qū)δＰ褪┘訃鷫?，達(dá)到預(yù)定值后通過伺服程序不斷調(diào)整左右邊墻位移速度來實現(xiàn)施加恒定圍壓，最后在豎直方向上通過上下邊界墻恒定的加載速率對數(shù)值試樣模型施加軸向載荷，直至試樣破壞為止。在不同梯度的圍壓下進(jìn)行上述數(shù)值模擬試驗，結(jié)合莫爾庫侖屈服準(zhǔn)則分析即可計算出內(nèi)聚力c和內(nèi)摩擦角φ的值。

圖4 雙軸壓縮

2 顆粒流細(xì)觀力學(xué)參數(shù)敏感性

為了研究顆粒流模型中細(xì)觀力學(xué)參數(shù)對宏觀參數(shù)的影響效果，以平直節(jié)理接觸模型為研究對象，主要細(xì)觀參數(shù)涉及：N，λ，kn/ks，μ，σc，cc，φ，Ec，其中：N為顆粒徑向單元數(shù)；λ為顆粒最大半徑與最小半徑比；kn/ks為平直節(jié)理接觸模型法切向剛度比；μ為接觸模型摩擦系數(shù)；σc為接觸模型抗拉強(qiáng)度；cc為接觸模型黏聚力；φ為接觸模型內(nèi)摩擦角；Ec為接觸模型有效模量[12]。由于平直節(jié)理接觸模型包含的細(xì)觀參數(shù)較多，如果盲目地逐個進(jìn)行調(diào)試會導(dǎo)致大量的數(shù)值試驗和繁瑣的計算，可能無法建立與實際情況相匹配的數(shù)值模型，因此可先探究接觸模型各細(xì)觀參數(shù)對顆粒材料宏觀力學(xué)參數(shù)的顯著性影響，以此作為指導(dǎo)依據(jù)，有利于減少細(xì)觀參數(shù)標(biāo)定過程中的盲目性，快捷準(zhǔn)確地建立數(shù)值模型。

2.1 正交試驗設(shè)計

當(dāng)遇到試驗結(jié)果受多方面因素共同影響時，若采用全面試驗的方法，將導(dǎo)致試驗次數(shù)大大增加，而正交試驗恰好可以解決試驗次數(shù)過多的問題。正交試驗設(shè)計就是通過特定模式來優(yōu)化試驗組合的方式來開展試驗研究，利用正交表可以科學(xué)地安排各因素組合，能夠在明顯減少試驗次數(shù)的前提下得到可靠性較強(qiáng)的試驗結(jié)果[13-14]。

考慮到平直節(jié)理接觸模型細(xì)觀參數(shù)中一些非特征性參數(shù)對材料宏觀特性影響較小，可根據(jù)相關(guān)研究進(jìn)行簡化。由文獻(xiàn)[5,10,15-16]可知，在一定范圍內(nèi)，顆粒半徑大小對顆粒流數(shù)值模型模擬結(jié)果的影響微乎其微，可作以下假設(shè)：平直節(jié)理接觸模型中抗拉強(qiáng)度要小于抗剪強(qiáng)度，以更符合巖石實際力學(xué)特性，同時可獲得與實際巖石相近的拉壓強(qiáng)度比；按程序默認(rèn)數(shù)值取λ=1，N=4；顆粒密度ρ取常用一般值2700 kg/m3；顆粒半徑比Rmax/Rmin取定值1.66，且最小粒徑Rmin=0.3 mm；顆粒接觸模量、剛度比和摩擦系數(shù)與平直節(jié)理接觸模型的相同。

除去已經(jīng)確定的細(xì)觀參數(shù)，剩余待標(biāo)定的細(xì)觀參數(shù)有6個，而進(jìn)行單軸壓縮、直接拉伸和雙軸壓縮(圍壓分別為5 MPa，10 MPa，15 MPa和20 MPa)數(shù)值模擬試驗(模擬試樣高100 mm，寬50 mm)可得到的宏觀參數(shù)同樣有6個，分別是：彈性模量E，泊松比ν，單軸抗壓強(qiáng)度σf，單軸抗拉強(qiáng)度σt，內(nèi)聚力c和內(nèi)摩擦角φ。由上述細(xì)觀參數(shù)在不同因素水平下建立的正交試驗設(shè)計表如表1所示，其中x為因素水平，y為對應(yīng)的細(xì)觀參數(shù)取值，與之相對應(yīng)的平直節(jié)理接觸模型正交設(shè)計矩陣序列如表2所示。

表2 平直節(jié)理接觸模型正交設(shè)計矩陣序列

2.2 細(xì)觀參數(shù)敏感性

當(dāng)涉及到多個因素對試驗結(jié)果產(chǎn)生影響時，可選用多因素方差法進(jìn)行分析計算，基于假設(shè)檢驗的原理，多因素方差分析法除可判別多個因素自身是否對因變量具有顯著性影響外，還可研究各因素之間共同作用對因變量產(chǎn)生的影響[17]，以此方法可探究多個細(xì)觀參數(shù)因素及它們之間的相互作用對各宏觀參數(shù)變量的顯著性影響關(guān)系，進(jìn)而分析細(xì)觀參數(shù)的敏感性。由表2中細(xì)觀參數(shù)進(jìn)行數(shù)值模擬試驗后，所獲得的材料宏觀力學(xué)參數(shù)結(jié)果如表3所示。多因素方差分析中考慮到各細(xì)觀因素的主效應(yīng)可取假設(shè)檢驗的顯著性水平為α=0.05，如果相伴概率值Sig≤α，則可認(rèn)為對應(yīng)因素對因變量產(chǎn)生顯著性影響，如果相伴概率值Sig>α，則無顯著性影響。

表3 正交設(shè)計矩陣序列對應(yīng)的宏觀參數(shù)計算結(jié)果

正交設(shè)計序列中σc，cc/σc和cc之間具有明顯的線性相關(guān)性，根據(jù)宏細(xì)觀參數(shù)之間的關(guān)聯(lián)性程度，在多因素方差分析中選擇σc，cc/σc和cc其中的一個進(jìn)行分析。針對不同的宏觀變量分析時選取的細(xì)觀參數(shù)組合略有差異，大量數(shù)值模擬試驗分析表明，宏觀參數(shù)中的ν和σt/σf與cc/σc的關(guān)聯(lián)性較強(qiáng)，所以當(dāng)以這兩個宏觀參數(shù)為因變量進(jìn)行多因素方差分析時，可選用細(xì)觀參數(shù)cc/σc，而對于其他宏觀參數(shù)皆采用cc進(jìn)行多因素方差分析。將正交設(shè)計矩陣表2和表3中數(shù)據(jù)導(dǎo)入SPSS數(shù)據(jù)處理軟件中，采用多因素方差法分析計算，根據(jù)多因素方差分析的F統(tǒng)計量和相伴概率值Sig的計算結(jié)果可知，宏觀參數(shù)的顯著性影響因素及各細(xì)觀參數(shù)的敏感性排序結(jié)果為：E：Ec>kn/ks>μ；ν：kn/ks>cc/σc>μ>φ>σc>Ec；σf：σc>cc>φ>kn/ks>Ec；σt：σc；σt/σf：cc/σc>φ>μ>Ec；c：σc>cc>kn/ks>φ>μ；φ：μ>φ>Ec>kn/ks>cc>σc。

圖5為SPSS多因素方差分析的F統(tǒng)計量結(jié)果及相伴概率Sig的值，根據(jù)圖中結(jié)果可得以下結(jié)論：

圖5 多因素方差分析的F統(tǒng)計量及相伴概率Sig(續(xù))

圖5 多因素方差分析的F統(tǒng)計量及相伴概率Sig

(1)彈性模量E主要受平直節(jié)理接觸模型有效模量Ec的影響，而接觸模型剛度比kn/ks和摩擦系數(shù)μ對其影響很小，幾乎可以忽略不計，其余細(xì)觀參數(shù)對其無影響。

(2)泊松比ν基本上受所研究的細(xì)觀因素共同影響，其中影響最大的因素是接觸模型法切向剛度比kn/ks，接觸模型的強(qiáng)度比cc/σc影響次之，接觸模型摩擦系數(shù)μ和內(nèi)摩擦角φ對其影響較小，其他細(xì)觀參數(shù)對其影響可忽略不計。

(3)單軸抗壓強(qiáng)度σf受兩大主要細(xì)觀因素影響，分別是接觸模型抗拉強(qiáng)度σc和接觸模型黏聚力cc，其中以σc影響更大，此外接觸模型內(nèi)摩擦角φ及其余細(xì)觀參數(shù)對其影響很小。

(4)影響單軸抗拉強(qiáng)度σt的細(xì)觀因素基本上只有接觸模型抗拉強(qiáng)度σc，使其標(biāo)定起來最為簡單。

(5)拉壓強(qiáng)度比σt/σf主要受接觸模型的強(qiáng)度比cc/σc的影響，接觸模型內(nèi)摩擦角φ對其影響程度相對較小，其他細(xì)觀參數(shù)對其影響極小。

(6)影響內(nèi)聚力c的因素比較多，基本上所有細(xì)觀參數(shù)都對其有影響，其中以接觸模型抗拉強(qiáng)度σc和黏聚力cc對其影響最大，接觸模型法切向剛度比kn/ks和內(nèi)摩擦角φ對其影響次之，接觸模型摩擦系數(shù)μ對其影響較小，其余細(xì)觀參數(shù)對其影響極小。

(7)內(nèi)摩擦角φ受所有細(xì)觀因素共同影響，其中對其影響最大的因素是接觸模型摩擦系數(shù)μ，而接觸模型內(nèi)摩擦角φ次之，接觸模型有效模量Ec對其影響僅次于內(nèi)摩擦角φ，其余3個細(xì)觀參數(shù)都對其有一定的影響，但影響程度相對來說很小。

由以上分析可知，顆粒流數(shù)值模擬中各細(xì)觀參數(shù)對宏觀參數(shù)的影響顯著性存在較大區(qū)別，宏細(xì)觀參數(shù)之間具有高度非線性特征，目前大多數(shù)研究者普遍采用“試錯法”進(jìn)行標(biāo)定[18-19],但調(diào)整過程盲目性大，費(fèi)時費(fèi)力。趙國彥等[7]、叢宇等[9]、陳鵬宇等[10]將數(shù)學(xué)分析與數(shù)值模擬試驗相結(jié)合，根據(jù)線性回歸分析擬合出宏細(xì)觀參數(shù)之間的近似對應(yīng)關(guān)系表達(dá)式，并以此對細(xì)觀參數(shù)進(jìn)行標(biāo)定，但基于線性擬合法的標(biāo)定過程誤差較大且后期調(diào)整過程工作量繁雜，因此，亟需尋求一種更加高效、快捷的標(biāo)定方法來確定細(xì)觀參數(shù)。而神經(jīng)網(wǎng)絡(luò)在處理非線性映射關(guān)系方面能力突出且無需事前揭示各變量之間明確的函數(shù)關(guān)系，它具有高度非線性映射能力、極強(qiáng)的自學(xué)習(xí)能力和良好的自適應(yīng)性，能夠通過機(jī)器的自我學(xué)習(xí)和訓(xùn)練建立輸入樣本向量和輸出樣本向量之間的非線性預(yù)測關(guān)系，是研究已知量與未知量之間較為復(fù)雜非線性映射關(guān)系的有力工具，應(yīng)用于細(xì)觀參數(shù)的標(biāo)定研究極具優(yōu)勢[20-21]。

3 細(xì)觀力學(xué)參數(shù)反演及結(jié)果分析

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

基于誤差逆向傳播而建立起的BP(back propagation)神經(jīng)網(wǎng)絡(luò)模型采用前反饋多層結(jié)構(gòu)體系，其最大的優(yōu)勢在于不需人為的干擾，只需通過大量數(shù)據(jù)訓(xùn)練來自我學(xué)習(xí)和儲存，便可建立樣本數(shù)據(jù)對應(yīng)的輸入輸出參數(shù)映射關(guān)系模型，且事前無需明確的數(shù)學(xué)表達(dá)式。其中自我學(xué)習(xí)機(jī)制采用梯度下降法，根據(jù)反向傳播的誤差值反饋來不斷調(diào)整網(wǎng)絡(luò)結(jié)構(gòu)并優(yōu)化模型內(nèi)部權(quán)值和閾值，以達(dá)到網(wǎng)絡(luò)模型的誤差平方和最小，從而保證模型預(yù)測的準(zhǔn)確性。原則上只要樣本訓(xùn)練數(shù)據(jù)足夠且網(wǎng)絡(luò)模型在自我學(xué)習(xí)訓(xùn)練中能達(dá)到收斂條件，就可建立起輸入變量與輸出變量之間的映射關(guān)系模型，并以此應(yīng)用于數(shù)據(jù)的預(yù)測分析，常用的BP神經(jīng)網(wǎng)絡(luò)模型拓?fù)浣Y(jié)構(gòu)如圖6所示，主要包含以下三層結(jié)構(gòu)。

圖6 前饋BP神經(jīng)網(wǎng)絡(luò)壓縮結(jié)構(gòu)模型

3.2 BP神經(jīng)網(wǎng)絡(luò)模型的建立與訓(xùn)練

基于BP神經(jīng)網(wǎng)絡(luò)分析方法，針對平直節(jié)理接觸模型，建立人工神經(jīng)網(wǎng)絡(luò)細(xì)觀參數(shù)反演系統(tǒng)。選取細(xì)觀力學(xué)參數(shù)特征值：kn/ks，μ，σc，cc，φ，Ec；選取宏觀力學(xué)參數(shù)特征值：E，ν，σf，σt，c，φ。由于要對細(xì)觀參數(shù)進(jìn)行標(biāo)定，即根據(jù)巖石實際宏觀力學(xué)參數(shù)反演出數(shù)值模型中相應(yīng)的細(xì)觀參數(shù)，可選取宏觀力學(xué)參數(shù)作為輸入層參數(shù)，而細(xì)觀力學(xué)參數(shù)作為輸出層參數(shù)，以此建立神經(jīng)網(wǎng)絡(luò)模型。參照正交設(shè)計表(25組)中細(xì)觀參數(shù)的取值范圍，隨機(jī)生成75組細(xì)觀參數(shù)并進(jìn)行相應(yīng)的數(shù)值模擬試驗，測試出與之對應(yīng)的宏觀參數(shù)。選取這100組不同宏細(xì)觀參數(shù)組合作為細(xì)觀參數(shù)反演系統(tǒng)的訓(xùn)練樣本，表4為上述宏細(xì)觀參數(shù)的取值范圍。由于選取的樣本數(shù)據(jù)中宏細(xì)觀參數(shù)都為6個，所以選取輸入層和輸出層神經(jīng)元數(shù)目也均為6個，調(diào)試神經(jīng)網(wǎng)絡(luò)程序內(nèi)部參數(shù)后導(dǎo)入樣本數(shù)據(jù)進(jìn)行訓(xùn)練學(xué)習(xí)，建立神經(jīng)網(wǎng)絡(luò)模型。

3.3 BP神經(jīng)網(wǎng)絡(luò)標(biāo)定法可靠性驗證

基于以上學(xué)習(xí)樣本數(shù)據(jù)建立起B(yǎng)P神經(jīng)網(wǎng)絡(luò)模型，為進(jìn)一步檢驗其可靠性，參照表4中宏觀力學(xué)參數(shù)的取值范圍隨機(jī)生成5組宏觀力學(xué)參數(shù)的測試組合樣本(表5)，用來測試神經(jīng)網(wǎng)絡(luò)的預(yù)測反演能力，將其導(dǎo)入已經(jīng)建立的網(wǎng)絡(luò)模型預(yù)測對應(yīng)的細(xì)觀參數(shù)。為校核細(xì)觀參數(shù)反演結(jié)果的準(zhǔn)確性，將反演結(jié)果輸入PFC中進(jìn)行數(shù)值模擬，根據(jù)模擬結(jié)果計算出所需的宏觀力學(xué)參數(shù)(表6)，并與導(dǎo)入網(wǎng)絡(luò)模型中的宏觀力學(xué)參數(shù)測試樣本值進(jìn)行對比分析。

表4 宏觀參數(shù)和細(xì)觀參數(shù)取值范圍

表5 BP神經(jīng)網(wǎng)絡(luò)測試樣本

表6 細(xì)觀參數(shù)反演后進(jìn)行數(shù)值模擬測試出的模擬值

數(shù)值模擬測試出的宏觀力學(xué)參數(shù)與隨機(jī)生成的宏觀力學(xué)參數(shù)測試樣本進(jìn)行對比分析，可驗算出反演結(jié)果的精度，以此衡量反演網(wǎng)絡(luò)系統(tǒng)的準(zhǔn)確性。精度T的定義為

式中，a為將反演后的細(xì)觀參數(shù)輸入PFC中數(shù)值模擬測試出的宏觀參數(shù)模擬值，A為隨機(jī)生成的宏觀參數(shù)測試樣本值。各宏觀參數(shù)的精度計算值如圖7所示。

圖7 測試樣本各宏觀參數(shù)計算精度值

根據(jù)測試樣本預(yù)測結(jié)果來評價BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測反演性能，可設(shè)實際隨機(jī)測試樣本數(shù)據(jù)序列為μm(i)，數(shù)值模擬計算數(shù)值序列為μn(i)，其殘差為

平均殘差為

則殘差均方差為

計算結(jié)果如表7所示。

由圖7可知，測試樣本反演結(jié)果精度普遍在90%以上，滿足建模精度要求。由表7可知，測試結(jié)果平均殘差值較小，說明預(yù)測誤差總體比較??；殘差均方差極小，表明預(yù)測結(jié)果較為穩(wěn)定，建立的神經(jīng)網(wǎng)絡(luò)模型預(yù)測精度值較高且模型結(jié)構(gòu)較為可靠。

表7 各宏觀參數(shù)殘差及均方差

3.4 BP神經(jīng)網(wǎng)絡(luò)參數(shù)反演標(biāo)定細(xì)觀參數(shù)

采用四川雙江口水電站地下廠房爆破開挖后的斑狀花崗巖，礦物成分以石英、長石居多并夾雜著黑云母和白云母等，質(zhì)地較為堅硬，水工切割打磨成直徑50 mm，高100 mm的標(biāo)準(zhǔn)試樣，每組試樣3個，在全自動化MTS(mechanical test system)試驗機(jī)上進(jìn)行單軸壓縮、巴西劈裂和常規(guī)三軸壓縮(圍壓分別為20 MPa，40 MPa，60 MPa和80 MPa)試驗，取測試結(jié)果平均值。正交試驗中數(shù)值模擬是通過直接拉伸模擬試驗獲得抗拉強(qiáng)度，而室內(nèi)試驗中則采用巴西劈裂試驗換算出抗拉強(qiáng)度，兩者雖略有不同，但分析正交試驗?zāi)M結(jié)果可知，直接拉伸數(shù)值模擬與巴西劈裂數(shù)值模擬測試出的抗拉強(qiáng)度相差較小，對細(xì)觀參數(shù)標(biāo)定的影響微乎其微。經(jīng)數(shù)據(jù)處理后可得室內(nèi)常規(guī)試驗巖石力學(xué)參數(shù)的實際值，以此作為預(yù)測輸入樣本，導(dǎo)入上述建立的神經(jīng)網(wǎng)絡(luò)模型可反演出相應(yīng)的細(xì)觀參數(shù)，再輸入PFC中進(jìn)行數(shù)值模擬可測試出巖石材料宏觀參數(shù)模擬值，兩者對比如表8所示。

表8 巖石試樣宏觀力學(xué)參數(shù)的試驗值和模擬值

由表8數(shù)據(jù)可以看出，數(shù)值試驗?zāi)M值與室內(nèi)試驗實際值非常接近，各宏觀參數(shù)的模擬精度都在90%以上，說明通過實際宏觀力學(xué)參數(shù)反演出的細(xì)觀參數(shù)(見表9)比較準(zhǔn)確，建立的顆粒流模型與實際巖石力學(xué)特性比較吻合，表明BP神經(jīng)網(wǎng)絡(luò)法標(biāo)定平直節(jié)理接觸模型細(xì)觀參數(shù)是可行的。

表9 巖石顆粒試樣細(xì)觀參數(shù)標(biāo)定結(jié)果

4 實例分析

根據(jù)表9中標(biāo)定的平直節(jié)理接觸模型細(xì)觀參數(shù)建立起花崗巖的顆粒流數(shù)值試驗?zāi)Ｐ?，分別進(jìn)行單軸壓縮、巴西劈裂和雙軸壓縮等數(shù)值模擬試驗，并與花崗巖對應(yīng)的室內(nèi)試驗結(jié)果進(jìn)行對比分析來檢驗所建立數(shù)值模型的準(zhǔn)確性。

4.1 應(yīng)力-應(yīng)變曲線的對比

室內(nèi)試驗與數(shù)值模擬獲得的單軸壓縮和常規(guī)三軸壓縮應(yīng)力?應(yīng)變曲線對比如圖8所示，從圖8(a)中可以看出，單軸壓縮室內(nèi)試驗和數(shù)值模擬試驗的應(yīng)力?應(yīng)變曲線相似度較高，彈性變形階段直線近似平行，峰后階段曲線近似重合且都反映出花崗巖的脆性破壞特征。初始階段兩者在曲線形態(tài)上雖略有差異，這是因為顆粒流程序建模過程中顆粒間未設(shè)置初始缺陷，應(yīng)力?應(yīng)變曲線在試樣破壞前是一條經(jīng)過原點(diǎn)的傾斜直線；而室內(nèi)試驗曲線由于天然花崗巖內(nèi)初始微裂隙的存在，導(dǎo)致初期應(yīng)力?應(yīng)變呈非線性發(fā)展即壓密階段；但這并未影響兩者在力學(xué)參數(shù)與破壞形式上的吻合度。由圖8(b)可知，不同圍壓(20 MPa，40 MPa，60 MPa)下的數(shù)值模擬試驗應(yīng)力?應(yīng)變曲線與室內(nèi)常規(guī)三軸壓縮應(yīng)力?應(yīng)變曲線契合度較高，且隨著圍壓的增大，巖石的峰值強(qiáng)度增大、塑性特征增強(qiáng)。初始階段數(shù)值模擬各曲線近似平行，說明不同圍壓對脆硬性花崗巖彈性模量影響不大，室內(nèi)試驗各曲線在加載初期近似重合，反映的花崗巖彈性變形特征與數(shù)值模擬結(jié)果相吻合。由此可見，基于神經(jīng)網(wǎng)絡(luò)反演出的細(xì)觀參數(shù)所建立的顆粒流模型與室內(nèi)試驗匹配度較高，兩者所得巖石的力學(xué)參數(shù)和應(yīng)力?應(yīng)變曲線極為接近，由于對實際復(fù)雜問題需要進(jìn)行簡化以方便建模，與實際情況雖存在略微差異，但大體上可以表征巖石的真實力學(xué)特性。

圖8 應(yīng)力?應(yīng)變曲線對比

4.2 試樣破壞形式對比

單軸壓縮和巴西劈裂試樣與相對應(yīng)數(shù)值模型試樣破壞模式對比如圖9所示，可知室內(nèi)試驗與數(shù)值模擬巖樣的破壞模式基本相同，單軸壓縮下試樣裂隙發(fā)展較多且較為分散，破碎較為完全，但主要表現(xiàn)為沿軸向的脆性劈裂破壞，局部發(fā)生輕微剪切破壞并伴隨有表面層狀剝落，宏觀裂隙方向與最大主應(yīng)力方向近似平行，破裂角大多在80°～90°之間。巴西劈裂試驗下，試樣首先從兩端起裂，而后沿徑向方向擴(kuò)展直至貫通整個圓盤，最后形成一條沿徑向方向的主裂縫，裂縫方向與軸向載荷加載方向基本平行；在主裂隙擴(kuò)展貫通過程中，試樣兩端還衍生出其他次生裂隙，并沿軸向加載方向擴(kuò)展，但擴(kuò)展長度較小，未貫通試樣。

圖9 單軸壓縮和巴西劈裂巖石破壞形態(tài)

常規(guī)三軸壓縮試驗與雙軸壓縮數(shù)值試驗結(jié)果對比如圖10和圖11所示。由圖10可知，20～80 MPa圍壓下試件以整體剪切破壞為主，形成單一宏觀破裂面，在圍壓為20 MPa和40 MPa時，破裂面起于試樣頂部，破裂角分別為72°和70°；圍壓為60 MPa時，破裂面從試樣端部延伸，破裂角為64°；圍壓為80 MPa時，破裂面下移，破裂角為58°。由此可知，隨著圍壓的增大，試樣整體剪切破壞現(xiàn)象愈發(fā)明顯，且滑移剪切帶的位置也在不斷發(fā)生變化，試件破裂角逐漸減小。圖11中數(shù)值模擬試驗試樣破碎分區(qū)圖(不同黑色部分代表不同破碎塊體)也反映了類似規(guī)律，模型試樣破壞形式基本上為宏觀剪切破壞，且隨圍壓增大，試樣整體剪切破壞越明顯，同時破裂角不斷減小。兩者對比結(jié)果說明標(biāo)定的細(xì)觀參數(shù)比較可靠，建立的顆粒流模型能夠較真實地反映花崗巖力學(xué)特性。

圖10 室內(nèi)常規(guī)三軸壓縮試驗不同圍壓下巖樣的破壞模式

圖11 雙軸數(shù)值模擬試驗不同圍壓下巖樣的破壞模式(續(xù))

5 結(jié)論

本文基于花崗巖室內(nèi)單軸、三軸壓縮及巴西劈裂試驗，結(jié)合離散元軟件PFC2D，探究了平直節(jié)理接觸模型中宏細(xì)觀參數(shù)之間的影響關(guān)系，并通過多因素方差分析對細(xì)觀參數(shù)敏感性進(jìn)行排序，進(jìn)一步運(yùn)用BP神經(jīng)網(wǎng)絡(luò)法對細(xì)觀參數(shù)標(biāo)定進(jìn)行應(yīng)用研究，得出的主要結(jié)論如下：

(1)運(yùn)用正交試驗設(shè)計科學(xué)安排數(shù)值試驗并采用多因素方差分析法探究各因素對因變量的顯著性影響關(guān)系，以此研究了平直節(jié)理接觸模型中宏細(xì)觀參數(shù)之間的關(guān)系并對細(xì)觀參數(shù)敏感性排序，結(jié)果表明各細(xì)觀參數(shù)對宏觀參數(shù)的影響顯著性存在較大區(qū)別，宏細(xì)觀參數(shù)之間的關(guān)系具有高度非線性特征。

(2)建立神經(jīng)網(wǎng)絡(luò)預(yù)測反演模型，并對平直節(jié)理接觸模型中細(xì)觀參數(shù)進(jìn)行標(biāo)定，標(biāo)定結(jié)果表明巖石材料各宏觀參數(shù)模擬值精度都在90%以上，且模擬值總體誤差較小，神經(jīng)網(wǎng)絡(luò)預(yù)測模型反演出的細(xì)觀參數(shù)較為可靠；神經(jīng)網(wǎng)絡(luò)標(biāo)定法在平直節(jié)理接觸模型中可作為一種高效、便捷的建模方法。

(3)以花崗巖室內(nèi)單軸壓縮、巴西劈裂和三軸壓縮試驗為基礎(chǔ)，采用神經(jīng)網(wǎng)絡(luò)標(biāo)定法建立花崗巖顆粒流模型，進(jìn)行單軸壓縮、巴西劈裂和雙軸壓縮數(shù)值模擬試驗，模擬結(jié)果與室內(nèi)試驗結(jié)果對比表明兩者獲得的巖石材料宏觀力學(xué)參數(shù)和應(yīng)力?應(yīng)變曲線基本一致，且兩者試樣的破壞模式極為接近，驗證了所建立花崗巖顆粒流模型的準(zhǔn)確性，即可表征巖石材料的真實力學(xué)特性。