彭霞，饒秋華，李卓，張杰

(中南大學(xué)土木工程學(xué)院，湖南長沙，410075)

在深部采礦、地?zé)衢_采、隧道開挖等巖體工程中，巖體內(nèi)部含有的多種缺陷(如孔隙、裂紋、節(jié)理等)在長期荷載作用下會發(fā)生損傷演變，造成工程安全隱患，因此，研究巖石破壞的宏?細(xì)觀機(jī)理對巖體工程的安全評定、防災(zāi)減災(zāi)等具有重大的理論指導(dǎo)意義。連續(xù)介質(zhì)力學(xué)難以解決多裂紋擴(kuò)展問題，且實驗研究存在一定的局限性。目前，數(shù)值方法已成為研究巖石細(xì)觀破壞機(jī)理的有效方法，其中，基于非連續(xù)性介質(zhì)理論的顆粒流軟件PFC(particle flow code)[1]已廣泛應(yīng)用于巖土工程領(lǐng)域。顆粒流軟件(PFC)是將材料離散成剛性顆粒后，通過模擬顆粒間的相互運動及其作用來研究真實材料的力學(xué)特性，其前提條件是選取合適的細(xì)觀參數(shù)。孟京京等[2?3]采用“試錯法”標(biāo)定PFC細(xì)觀參數(shù)，但由于細(xì)觀參數(shù)因素多且標(biāo)定具有一定的盲目性與不確定性，導(dǎo)致標(biāo)定周期長，模擬試驗次數(shù)多。YOON等[4]采用PB(Plackett-Burman)設(shè)計法測試細(xì)觀參數(shù)對宏觀參數(shù)(如單軸抗壓強度、彈性模量、泊松比等)的敏感性，選取影響最大的2個細(xì)觀參數(shù)，再通過CCD(central composite design)設(shè)計法建立宏?細(xì)觀參數(shù)之間的非線性關(guān)系，可以求得細(xì)觀參數(shù)，但每個宏觀參數(shù)僅由較少的細(xì)觀參數(shù)確定，沒有考慮其他細(xì)觀參數(shù)對宏觀參數(shù)的影響。周喻等[5]采用BP(back propagation)神經(jīng)網(wǎng)絡(luò)法方法，建立了宏?細(xì)觀參數(shù)的非線性模型，輸入宏觀力學(xué)參數(shù)，即可得到巖土體細(xì)觀力學(xué)參數(shù)，但BP神經(jīng)網(wǎng)絡(luò)模型的創(chuàng)建、學(xué)習(xí)、訓(xùn)練需要大量的隨機(jī)組合樣本。陳鵬宇等[6]對細(xì)觀參數(shù)進(jìn)行正交設(shè)計，以多因素分析、量綱分析研究宏?細(xì)觀參數(shù)之間的關(guān)系，確定的函數(shù)關(guān)系式只考慮單個細(xì)觀參數(shù)的影響，提出了試錯法標(biāo)定細(xì)觀參數(shù)的具體流程。目前，對于PFC 細(xì)觀參數(shù)的定量確定方法仍較少。為降低模擬試驗次數(shù)且提高精度，需引入一種新方法來定量建立PFC宏?細(xì)觀參數(shù)關(guān)系式。盛海林等[7]提出了球面對稱設(shè)計法，在球面上均勻、對稱地選取試驗點，通過試驗建立了自變量與因變量之間的關(guān)系，該方法具有試驗次數(shù)較少、精度較高等優(yōu)點，更適用于多因素、多水平試驗。本文采用球面對稱設(shè)計法，建立PFC2D宏?細(xì)觀參數(shù)定量關(guān)系式。依據(jù)材料宏觀力學(xué)參數(shù)定量確定其細(xì)觀參數(shù)，模擬計算單軸壓縮下完整巖石試件的應(yīng)力?應(yīng)變曲線和宏觀力學(xué)參數(shù)、含多裂紋巖石試件的斷裂軌跡，并與試驗結(jié)果進(jìn)行對比驗證，以便為PFC 細(xì)觀參數(shù)的定量分析提供一種有效方法。

1 球面對稱設(shè)計法

在球面對稱設(shè)計方法中，設(shè)因素個數(shù)為n，則球的半徑r為n的平方根，球的空間維數(shù)等于n，在該球面上均勻、對稱地選取試驗點?？紤]到任意因素對各因變量的影響均有一定范圍，將這些范圍以一定的代碼予以統(tǒng)一，規(guī)定每個因素有5個水平，分別用?r，?1，0，+1和+r共5個代碼表示，其中，各因素水平代碼?r，0和+r對應(yīng)的值分別為該因素范圍的最小值、平均值、最大值；根據(jù)任意2個水平代碼差的比值與該水平代碼對應(yīng)值差的比值相等的原則，可得到各因素水平代碼?1和+1 對應(yīng)的值。試驗次數(shù)包含3 部分：所有因素的水平代碼僅為?1或+1的試驗次數(shù)(2n次)；1個因素水平代碼為?r或+r、其余因素水平代碼為0的試驗次數(shù)(2n次)；各因素水平代碼為0的試驗次數(shù)(1次)?？傇囼灤螖?shù)為2n+2n+1次。

以5因素為例建立球面對稱設(shè)計表。設(shè)5個因素為X1，X2，X3，X4和X5，即n=5，球半徑r=5，則各因素的5個水平代碼為-5，?1，0，+1和+5。在設(shè)計表中，總列數(shù)為6(第1列為試驗次數(shù)號)，其余每一列元素為某一因素在每次試驗中的水平代碼；總行數(shù)與試驗總次數(shù)43(25+2×5+1)相同，每行元素為5個因素的某個水平代碼(5個水平代碼之一)。其中，第1～32行(2n次試驗)中的每一行是各因素的水平代碼?1 或+1的任意組合，且各行不重復(fù)；第33～42 行(2n次試驗)中的每一行是某一因素水平代碼為-5 或+5、其余因素水平代碼為0的任意組合，且各行不重復(fù)；第43 行為各因素水平代碼為0的組合。由此可得5因素球面對稱設(shè)計表，見表1。

2 細(xì)觀參數(shù)計算公式

2.1 PFC2D建模

在PFC2D 中，顆粒之間通過黏結(jié)產(chǎn)生相互作用，黏結(jié)模型包含接觸黏結(jié)模型、平行黏結(jié)模型、平直節(jié)理接觸模型3種[8?9]。接觸黏結(jié)模型中的圓形顆粒與顆粒之間為點接觸，因此，不能傳遞力矩，當(dāng)法向或切向力超過對應(yīng)的黏結(jié)強度時，黏結(jié)破壞。平行黏結(jié)模型中兩圓形顆粒間由平行鍵黏結(jié)，可以傳遞顆粒之間的力與力矩，該模型比接觸黏結(jié)模型更接近于巖石材料真實受力情況，應(yīng)用更廣。平直節(jié)理接觸模型是將多邊形顆粒代替圓形顆粒，能夠傳遞力與力矩，同時能抑制顆粒黏結(jié)破壞后的旋轉(zhuǎn)，但需要確定更多的細(xì)觀參數(shù)。

表1 5因素球面對稱設(shè)計表Table 1 Spherical symmetric design table of five factors

本文選取平行黏結(jié)模型，其細(xì)觀參數(shù)如表2所示。為簡化計算，進(jìn)行以下假設(shè)[9?11]：；。設(shè)平行黏結(jié)法向應(yīng)力強度與法向應(yīng)力強度之比σn/σs=m(m為比例系數(shù))，本文取m=1[12?13]，故待定的獨立細(xì)觀參數(shù)為6個，即Rmin，Ec，kn/ks，σn，μ和Rσ。

在PCF2D 建模中，先確定模型的最小顆粒半徑Rmin。研究結(jié)果表明[14?15]，當(dāng)Rmin小到一定程度時，宏觀力學(xué)參數(shù)計算結(jié)果變化不大，但Rmin過小會降低數(shù)值模擬計算效率。為得到合適的Rmin，以巖石單軸壓縮試驗為例說明。采用長為100 mm、寬為50 mm的長方形標(biāo)準(zhǔn)試樣，選取4 種不同的Rmin(0.1，0.2，0.3和0.4 mm)進(jìn)行建模，如圖1所示。設(shè)其他細(xì)觀參數(shù)為：，模擬計算不同Rmin下的巖石宏觀力學(xué)參數(shù)。由表3 可見計算得到的彈性模量E、泊松比v、單軸抗壓強度σc均在巖石的宏觀力學(xué)參數(shù)范圍內(nèi)[16]，E=10～70 GPa，v=0.15～0.30，σc=50～200 MPa，表明所選取的細(xì)觀參數(shù)是合理可行的；當(dāng)Rmin減少至0.1 mm時，E，v和σc的計算結(jié)果與Rmin=0.2時所得結(jié)果變化不大，但當(dāng)Rmin=0.2 mm 時，所需生成的顆粒數(shù)目比Rmin=0.1 mm 時大大減少，因而，選定Rmin=0.2 mm。待定的獨立細(xì)觀參數(shù)只有5個：Ec，kn/ks，σn，μ和Rσ。

表2 PFC2D平行黏結(jié)模型中的細(xì)觀參數(shù)Table 2 Microscopic parameters of PFC2D parallel bonded model

圖1 單軸壓縮下巖石試件PFC2D模型(Rmin=0.2 mm)Fig.1 PFC2D model of rock specimen under uniaxial compression(Rmin=0.2 mm)

表3 不同Rmin下巖石宏觀力學(xué)參數(shù)的模擬結(jié)果Table 3 Simulation results of macroscopic mechanical parameters of rock with different Rmin

2.2 細(xì)觀參數(shù)范圍確定

利用球面對稱設(shè)計法確定PFC2D細(xì)觀參數(shù)時，首先需根據(jù)宏觀參數(shù)的范圍確定細(xì)觀參數(shù)的取值范圍。對于巖石材料，其宏觀力學(xué)參數(shù)范圍[16]為：彈性模量E=10～70 GPa，泊松比v=0.15～0.30，單軸抗壓強度σc=50～200 MPa，起裂應(yīng)力與單軸抗壓強度σci/σc=0.2～0.6。在PFC2D 中，σci通過出現(xiàn)初始裂紋的應(yīng)力與峰值應(yīng)力的比值R設(shè)定，本文設(shè)定R為1%，即σci等于當(dāng)裂紋數(shù)達(dá)到峰值裂紋數(shù)的1%時所對應(yīng)的應(yīng)力。

圖2 細(xì)觀參數(shù)對v的影響Fig.2 Effect of mesoscopic parameters on v

圖3 細(xì)觀參數(shù)對E的影響Fig.3 Effect of mesoscopic parameters on E

圖4 細(xì)觀參數(shù)對σc的影響Fig.4 Effect of mesoscopic parameters on σc

圖5 細(xì)觀參數(shù)對σci/σc的影響Fig.5 Effect of mesoscopic parameters on σci/σc

2.3 宏?細(xì)觀參數(shù)關(guān)系式

考慮到待定的細(xì)觀參數(shù)有5個，采用5因素球面對稱設(shè)計表(表1)，5因素X1，X2，X3，X4和X5分別為Ec，kn/ks，σn，μ和Rσ。為確定表1中各因素水平代碼對應(yīng)的值，以因素Ec為例說明，Ec的范圍為7～51，則Ec的水平代碼，0和所對應(yīng)的值分別為7，29和51。令X1的水平代碼?1和1對應(yīng)的值分別為x和y(表4)，根據(jù)上述比值相等的原則，有

解得x=19.16，y=38.84。同理，可求得其他因素水平代碼對應(yīng)的值，如表5所示。將表5中各因素水平代碼對應(yīng)的值與表1中的值相對照，可得到具體的球面對稱設(shè)計試驗表，如表6所示。

根據(jù)球面對稱設(shè)計試驗表(表6)，選取每行中的5個細(xì)觀參數(shù)值，進(jìn)行PFC2D模擬計算。表7所示為模擬計算43 次試驗得到的巖石單軸壓縮試件宏觀力學(xué)參數(shù)(E，v，σc和σci/σc)，參數(shù)范圍為：E=10.16～72.77 GPa，v=0.1603～0.315，σc=49.70～198.98 GPa，σci/σc=0.244～0.674，基本符合巖石的力學(xué)參數(shù)范圍。

表4 因素Ec水平代碼對應(yīng)的值Table 4 Corresponding values for horizontal codes of factor Ec

表5 各因素水平代碼對應(yīng)的值Table 5 Corresponding values for horizontal codes of each factor

表6 球面對稱設(shè)計試驗表Table 6 Spherical symmetric design table

表7 球面對稱設(shè)計試驗?zāi)M結(jié)果Table 7 Simulation results of spherical symmetry design test

球面對稱設(shè)計的試驗結(jié)果可采用線性或非線性方法擬合[7]，本文采用多元線性回歸對表7 中模擬試驗結(jié)果進(jìn)行擬合，分別建立4個宏觀參數(shù)(E，v，σc和σci)與5個細(xì)觀參數(shù)(Ec，kn/ks，σn，μ和Rσ)之間的關(guān)系式：

相關(guān)系數(shù)R2均在90%以上，表明擬合精度較高，能夠準(zhǔn)確反映巖石宏?細(xì)觀參數(shù)的關(guān)聯(lián)性，其中摩擦因數(shù)主要影響峰值后的響應(yīng)，建議取μ=0.5[9]。式(1)適用于宏觀力學(xué)參數(shù)滿足E=10～70 GPa，v=0.15～0.30，σc=50～200 MPa，σci/σc=0.2～0.6的巖石材料。若已知巖石4個宏觀參數(shù)(E，v，σc和σci)，由式(1)可計算得到4個細(xì)觀參數(shù)。

3 實驗驗證

3.1 巖石單軸壓縮實驗

為驗證宏?細(xì)觀參數(shù)關(guān)系式(1)，選取7 種不同的巖石材料，根據(jù)其宏觀力學(xué)參數(shù)(見表8)，分別計算得到PFC2D 細(xì)觀參數(shù)值(Ec，kn/ks，σn和Rσ)，如表9所示。

基于不同巖石材料的PFC2D 細(xì)觀參數(shù)(表9)，采用標(biāo)準(zhǔn)巖石試樣(長×寬為50 mm×100 mm)進(jìn)行單軸壓縮模擬計算，得到應(yīng)力?應(yīng)變曲線(見圖6，以砂巖和BS 花崗巖為例)和宏觀力學(xué)參數(shù)(見表10)，并與試驗結(jié)果進(jìn)行對比。由圖6 可見：模擬應(yīng)力?應(yīng)變模擬曲線(OB′C′D′)可分為彈性(OB′)、損傷(B′C′)和破壞(C′D′)3個階段，但沒有出現(xiàn)類似試驗曲線(OABCD)的壓密過程(OA)，這是由于在壓縮試驗的初始階段，微孔、微裂縫等閉合引起非線性不可逆變形(OE)。計算得到的砂巖彈性模量E(OB′段斜率)和壓縮強度σc(峰值點C′應(yīng)力)分別為34.17 GPa和193.53 MPa，與實測值E=34.68 GPa，σc=190.8 MPa 相比相對誤差均較小，均約為2%；同理，計算得到的BS花崗巖E和σc與實測值相對誤差均小于2%。當(dāng)忽略壓密階段、試驗曲線ABCD向原點水平左移OE段后，所得曲線與模擬曲線基本重合，兩者吻合較好。

表8 不同巖石的宏觀力學(xué)參數(shù)(試驗結(jié)果)Table 8 Macroscopic mechanical parameters of different rocks(experimental results)

不同巖石試樣在單軸壓縮下宏觀力學(xué)參數(shù)的模擬結(jié)果與試驗結(jié)果對比見表10。從表10 可見：兩者相對誤差小于5%，精度較高，表明采用球面對稱設(shè)計建立的宏?細(xì)觀關(guān)系式(1)計算PFC2D 細(xì)觀參數(shù)的方法是合理的、有效的。

3.2 巖石多裂紋斷裂試驗

為進(jìn)一步驗證宏?細(xì)觀參數(shù)關(guān)系式(1)，模擬含3 條平行裂紋的砂巖試件(長160 mm、寬80 mm、厚30 mm)的單軸壓縮斷裂試驗[17]。如圖7所示，試件模型中從上至下依次為預(yù)制裂紋1、預(yù)制裂紋2、預(yù)制裂紋3，預(yù)制裂紋長為15 mm，寬為2.5 mm，預(yù)制裂紋1 與預(yù)制裂紋2和3 內(nèi)尖端的距離2b=20 mm，垂直于預(yù)制裂紋2和3的距離d=10 mm，預(yù)制裂紋1的水平夾角α=45°，預(yù)制裂紋1和3內(nèi)尖端連線與水平方向的夾角β=75°，實測的砂巖宏觀力學(xué)參數(shù)和基于關(guān)系式(1)計算得到的PFC2D 細(xì)觀參數(shù)見表8和表9。

表9 不同巖石的PFC2D細(xì)觀參數(shù)(計算結(jié)果)Table 9 PFC2D microscopic parameters of different rocks(calculated results)

圖6 單軸壓縮下完整巖石的應(yīng)力σ?應(yīng)變ε曲線Fig.6 Stress?strain curves of intact rocks under uniaxial compression

表10 不同巖石的宏觀力學(xué)參數(shù)(模擬與試驗結(jié)果對比)Table 10 Macroscopic mechanical parameters of different rocks(comparison of simulation and test results)

圖7 含3條平行裂紋砂巖試件模型Fig.7 Model of sandstone specimen with three parallel cracks

圖8 單軸壓縮下含3條平行裂紋的砂巖試件裂紋擴(kuò)展軌跡(模擬結(jié)果)Fig.8 Propagation trajectories of sandstone specimen with three parallel cracks under uniaxial compression(simulation results)

圖9 單軸壓縮下含3條平行裂紋的砂巖試件裂紋擴(kuò)展軌跡(試驗結(jié)果)[17]Fig.9 Propagation trajectories of sandstone specimen with three parallel cracks under uniaxial compression(test results)

在單軸壓縮下，含3條平行裂紋的砂巖試件裂紋擴(kuò)展軌跡其PFC2D 模擬結(jié)果和試驗結(jié)果分別如圖8和圖9所示。由模擬結(jié)果可知：3 條平行裂紋AB，CD和EF均有新裂紋分別從裂尖起裂(圖8(a))，沿著軸向加載方向擴(kuò)展(與含單一斜裂紋的試件擴(kuò)展形式相同)，其中，淺色和深色分別表示拉應(yīng)力和剪應(yīng)力超過平行鍵法向強度、剪切強度而形成的拉裂紋、剪裂紋，拉裂紋居多，起裂主要受拉應(yīng)力控制。從圖9(a)可見：只有2條平行裂紋AB和CD的尖端先出現(xiàn)新裂紋，這可能是實際巖石材料的非均質(zhì)性所致。從圖8(b)可見；隨著荷載增加，裂尖E的新裂紋因到裂紋AB和CD的距離較遠(yuǎn)，相互影響較小，仍沿原軸向加載方向擴(kuò)展，而裂尖F的新裂紋逐漸擴(kuò)展至裂尖C附近區(qū)域，最終與裂尖C貫通；裂尖B和裂紋C擴(kuò)展因受到裂尖D和A的應(yīng)力強度因子影響而改變方向，最終與裂尖D和A貫通尖端；裂尖A和D新裂紋沿軸向加載方向雙向擴(kuò)展，裂尖A和D逐漸貫通，與圖9(b)所示實驗結(jié)果相符，此時，裂尖A，C和D均有2 條裂紋相互貫通，應(yīng)力集中得到釋放，裂尖B只有1條新裂紋，故產(chǎn)生次生裂紋(圖8(c)，9(c))。從圖8(d)和圖9(d)可見：隨著進(jìn)一步加載，裂尖E的新裂紋逐步向試件上端擴(kuò)展時受到壓密阻礙，又分叉衍生出向下擴(kuò)展的次生裂紋。但模擬中次生裂紋趨向于與試件右側(cè)邊緣貫通，試驗中次生裂紋軸向加載方向延伸，其中差異可能是實際巖石材料的非均質(zhì)性所致。由此可見：在單軸壓縮下，含3條平行裂紋的砂巖試件裂紋擴(kuò)展軌跡其PFC2D 模擬結(jié)果與試驗結(jié)果較吻合，從而驗證了采用球面對稱設(shè)計法標(biāo)定PFC2D 細(xì)觀參數(shù)的合理性和可靠性。

4 結(jié)論

1)采用球面對稱設(shè)計法，通過宏觀參數(shù)范圍確定細(xì)觀參數(shù)取值范圍，得到多因素球面對稱設(shè)計試驗表，利用試驗?zāi)M計算結(jié)果和多元回歸分析，建立了PFC2D宏?細(xì)觀參數(shù)關(guān)系式。

2)采用基于宏?細(xì)觀參數(shù)關(guān)系式計算得到的PFC2D 細(xì)觀參數(shù)，模擬計算單軸壓縮下完整巖石的應(yīng)力?應(yīng)變曲線和宏觀力學(xué)參數(shù)、含多裂紋巖石試件斷裂軌跡，計算結(jié)果與試驗結(jié)果較吻合，從而驗證了該宏?細(xì)觀參數(shù)關(guān)系式的合理性和可行性。

3)采用球面對稱設(shè)計建立PFC2D宏?細(xì)觀參數(shù)關(guān)系式的方法可推廣應(yīng)用于其他脆性材料，為PFC定量分析復(fù)雜荷載條件下的材料細(xì)觀破壞機(jī)理提供參考。