基于球面對稱設計的離散元細觀參數(shù)定量確定方法

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

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

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

1 球面對稱設計法

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

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

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

2.1 PFC2D建模

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

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

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

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

表2 PFC2D平行黏結模型中的細觀參數(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下巖石宏觀力學參數(shù)的模擬結果Table 3 Simulation results of macroscopic mechanical parameters of rock with different Rmin

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

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

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

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

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

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

2.3 宏?細觀參數(shù)關系式

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

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

根據球面對稱設計試驗表(表6)，選取每行中的5個細觀參數(shù)值，進行PFC2D模擬計算。表7所示為模擬計算43 次試驗得到的巖石單軸壓縮試件宏觀力學參數(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，基本符合巖石的力學參數(shù)范圍。

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

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

表6 球面對稱設計試驗表Table 6 Spherical symmetric design table

表7 球面對稱設計試驗模擬結果Table 7 Simulation results of spherical symmetry design test

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

相關系數(shù)R2均在90%以上，表明擬合精度較高，能夠準確反映巖石宏?細觀參數(shù)的關聯(lián)性，其中摩擦因數(shù)主要影響峰值后的響應，建議取μ=0.5[9]。式(1)適用于宏觀力學參數(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個細觀參數(shù)。

3 實驗驗證

3.1 巖石單軸壓縮實驗

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

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

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

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

3.2 巖石多裂紋斷裂試驗

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

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

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

表10 不同巖石的宏觀力學參數(shù)(模擬與試驗結果對比)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條平行裂紋的砂巖試件裂紋擴展軌跡(模擬結果)Fig.8 Propagation trajectories of sandstone specimen with three parallel cracks under uniaxial compression(simulation results)

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

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

4 結論

1)采用球面對稱設計法，通過宏觀參數(shù)范圍確定細觀參數(shù)取值范圍，得到多因素球面對稱設計試驗表，利用試驗模擬計算結果和多元回歸分析，建立了PFC2D宏?細觀參數(shù)關系式。

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

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