基于中心組合設計的顆粒流平直節(jié)理模型宏-細觀參數(shù)相關性研究

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

在礦產(chǎn)資源開采、邊坡、隧道、石油開采等巖體工程中，天然巖石材料存在的細觀缺陷(如裂紋或孔隙)極大地影響著巖石的宏觀力學性能。在外載作用下，裂紋或孔隙的演變發(fā)展形成宏觀缺陷會引起巖體工程災害隱患，因此，研究巖石細觀缺陷發(fā)展至宏觀破壞過程及其機理，對巖體工程的強度計算、安全評估和防災減災等具有重要理論和實際意義。由于現(xiàn)有試驗測試手段的局限性和巖石內(nèi)部缺陷的復雜性，人們難以在巖體工程中實時監(jiān)測到巖石細觀損傷積累至宏觀破壞的全過程，故多借助于數(shù)值模擬方法。目前，研究巖石裂紋擴展的數(shù)值模擬方法主要有離散元法、有限元法、邊界元法等，其中，基于非連續(xù)介質(zhì)力學理論的離散元法在模擬巖石裂隙擴展時具有獨特的優(yōu)勢，因為該方法既能避免選取復雜的本構(gòu)關系和裂尖奇異性，又能夠從細觀層面上模擬巖石宏觀力學特性并揭示其破壞機理。離散元顆粒流軟件PFC(particle flow code)在巖石類材料的數(shù)值模擬中應用較廣泛[1]。PFC 軟件是將巖石離散為剛性顆粒，通過定義顆粒接觸本構(gòu)關系來描述顆粒之間的相互作用，以此模擬巖石的宏觀力學特性。

目前，PFC 顆粒黏結(jié)模型大多為接觸黏結(jié)(linear contact bond)、平行黏結(jié)(parallel bond)和平直節(jié)理模型(flat jointed model)。在接觸黏結(jié)模型中，黏結(jié)斷裂時顆粒接觸剛度仍然存在，這與巖石材料不相符[2]；在平行黏結(jié)模型中，模擬得到的巖石抗拉強度與抗壓強度的比值偏大，與巖石材料的力學性能也不相符[3-4]，通常采用簇平行和新的膠接模型來解決該問題[5-6]。在平直節(jié)理模型中[7]，圓形顆粒被構(gòu)造成多邊形從而抑制了顆粒的旋轉(zhuǎn)，使得該模型計算得到的單軸抗拉強度和抗壓強度之比減小，更符合巖石力學特性。采用顆粒流PFC2D軟件模擬巖石的宏觀力學特性時，首先需要確定顆粒黏結(jié)模型的細觀參數(shù)。目前，采用較多的是試錯法，如文獻[8-9]基于試錯法標定了平行黏結(jié)的細觀參數(shù)；文獻[10]基于試錯法標定了平直節(jié)理接觸模型的細觀參數(shù)，并與平行黏結(jié)模型的細觀參數(shù)進行對比，分析了平直節(jié)理模型的優(yōu)勢，但未建立宏-細觀參數(shù)之間的關系式。由于試錯法具有一定的盲目性，細觀參數(shù)較多且標定過程較繁瑣，從而導致試驗次數(shù)很多，因此，研究者尋求更有效的細觀參數(shù)定量確定方法，如：文獻[11]采用單因素和多因素分析方法，分析了平行黏結(jié)模型的細觀參數(shù)對宏觀參數(shù)的影響，建立了平行黏結(jié)模型的宏-細觀參數(shù)關系式；文獻[12]采用BP神經(jīng)網(wǎng)絡法，建立了簇平行黏結(jié)模型的宏-細觀參數(shù)的非線性模型，但BP 神經(jīng)網(wǎng)絡需要大量的隨機組合樣本才能具有較高的精度；文獻[13]采用球面對稱設計法建立了平行黏結(jié)宏-細觀參數(shù)的關系式，但考慮的宏觀參數(shù)有限；文獻[14-15]通過對平直節(jié)理模型細觀參數(shù)進行正交分析，得到了細觀參數(shù)對宏觀參數(shù)的影響規(guī)律，并利用回歸分析建立了宏-細觀參數(shù)關系式，但因試驗次數(shù)較少導致精度較低?？傊壳瓣P于PFC2D平直節(jié)理接觸模型的細觀參數(shù)定量確定方法偏少。為提高精度、減少試驗次數(shù)，需要尋求一種更合適的方法來建立PFC2D平直節(jié)理接觸模型的宏-細觀參數(shù)關系式。中心組合設計(又稱星點設計，central composite design)[16-17]是一種多因素五水平的試驗設計方法，相對于正交試驗和均勻化試驗，它在提高精度的同時不需要大量增加實驗次數(shù)，特別適用于多因素、多水平試驗情形。

本文作者基于PFC2D平直節(jié)理模型，通過控制變量法得到細觀參數(shù)的取值范圍和對宏觀參數(shù)的影響規(guī)律，采用中心組合設計和回歸分析，建立巖石PFC2D宏-細觀參數(shù)關系式以定量確定細觀參數(shù)，模擬計算巖石在單軸和雙軸壓縮下完整巖石試件的應力-應變曲線和含多裂紋巖石試件的斷裂軌跡，并與試驗結(jié)果進行對比，以驗證所建立的巖石PFC2D宏-細觀參數(shù)關系式的有效性。

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

平直節(jié)理模型是針對硬巖提出的一種模型，如圖1所示。該模型將圓形顆粒構(gòu)造成多邊形，從而限制了顆粒的旋轉(zhuǎn)，致使模型單軸抗壓強度和抗拉強度的比值增大，更加符合巖石類材料的力學性能。平直節(jié)理模型采用扁平面將2 個顆粒黏結(jié)，黏結(jié)面由片段的單元組成。每個單元具有黏結(jié)和未黏結(jié)2種模式，當黏結(jié)界面上的荷載單元達到黏接強度時，黏結(jié)單元便會破壞成非黏結(jié)單元，其本構(gòu)關系如圖2所示。表1所示為平直節(jié)理模型的細觀參數(shù)，包括顆粒細觀參數(shù)(6個)和平直節(jié)理參數(shù)(7個)。

圖1 PFC2D平直節(jié)理模型[7]Fig.1 Flat joint model of PFC2D[7]

圖2 PFC2D平直節(jié)理模型的本構(gòu)模型[7]Fig.2 Constitutive relation of PFC2D flat joint model[7]

2 PFC細觀參數(shù)確定方法

2.1 細觀參數(shù)取值范圍的確定

為了模擬巖石類材料宏觀力學特性，首先需要確定PFC2D平直節(jié)理模型的細觀參數(shù)。由于平直節(jié)理模型的細觀參數(shù)較多(共13個)，為了減少標定次數(shù)，根據(jù)前人的研究成果[13]，進行以下假設：

表1 PFC2D平直節(jié)理模型的細觀參數(shù)[7]Table 1 Microscopic parameters of PFC2D flat joint model[7]

1)顆粒接觸模量Eb等于平直節(jié)理接觸模量Ec；

2)顆粒剛度比(kn/ks)b等于平直節(jié)理接觸剛度比kn/ks；

3)最大、最小顆粒半徑之比Rmax/Rmin=1.66；

4)顆粒密度ρ=2 300 kg/m3；

5)平直節(jié)理半徑系數(shù)-λ=1；

6)交界面段數(shù)N=4；

7)平直節(jié)理抗拉強度小于抗剪強度。

為了確定PFC2D細觀參數(shù)，需要分析細觀參數(shù)對巖石材料宏觀力學參數(shù)(含彈性模量E、泊松比ν、抗壓強度σf、抗拉強度σt、黏聚力C2、摩擦角φ)的影響，以此確定細觀參數(shù)的范圍。針對剩余的7個細觀參數(shù)采用控制變量法分析，選取適當?shù)某跏技氂^參數(shù)，即Ec=40 GPa，kn/ks=1，C1=30 MPa，σc=3 MPa，μb=0.1，μ=0.1，R=0.2 mm[18]，在改變單一細觀參數(shù)的情況下開展巖石材料的單軸壓縮(E，ν，σf)、單軸拉伸(σt)和雙軸壓縮(C2，φ)數(shù)值模擬試驗，計算得到的單一細觀參數(shù)對宏觀力學參數(shù)的影響結(jié)果如圖3～9 所示，其中橫坐標為PFC2D細觀參數(shù)，縱坐標為歸一化處理后的巖石宏觀參數(shù)，即細觀參數(shù)改變后通過模擬計算得到的宏觀參數(shù)與由初始細觀參數(shù)模擬計算得到的宏觀參數(shù)之比。

從圖3可知：接觸模量Ec對宏觀彈性模量E影響最大，對其他宏觀參數(shù)幾乎沒有影響(在圖中表現(xiàn)為其他參數(shù)曲線基本重合)。由圖4可知：剛度比kn/ks主要影響泊松比ν，隨著kn/ks的增加，宏觀彈性模量E、抗拉強度σt略減少(在圖中表現(xiàn)為E，σt和C2曲線基本重合)，其他宏觀參數(shù)基本保持不變。由圖5可知：平直節(jié)理黏聚力C1主要影響σf和黏聚力C2，對泊松比ν和摩擦角φ的影響較小。由圖6可知：平直節(jié)理抗拉強度σc主要影響σt，隨著σc增加，σf略增加。由圖7和圖8可知：顆粒摩擦因素μb和平直節(jié)理內(nèi)摩擦因數(shù)μ主要影響摩擦角φ，且μb對φ的影響更顯著；隨著μ增加，σf略微增加，C2略微減少，對其他宏觀參數(shù)影響較小(在圖中表現(xiàn)為E、σt曲線基本重合)。由于μb和μ對宏觀參數(shù)的影響具有相同的規(guī)律，可取μ=μb[7]。由圖9可知，顆粒半徑R對宏觀參數(shù)影響的幅度較小，這與文獻[19]中的結(jié)論一致?？紤]到計算效率和精度，本文選取最小半徑Rmin=0.3 mm。因此，待定的獨立細觀參數(shù)只有5個：Ec，kn/ks，C1，σc和μb。

圖3 細觀參數(shù)Ec對宏觀參數(shù)的影響Fig.3 Effect of microscopic parameters Ec on macroscopic parameters

圖4 細觀參數(shù)kn/ks對宏觀參數(shù)的影響Fig.4 Effect of microscopic parameters kn/ks on macroscopic parameters

圖5 細觀參數(shù)C1對宏觀參數(shù)的影響Fig.5 Effect of microscopic parameters C1 on macroscopic parameters

圖6 細觀參數(shù)σc對宏觀參數(shù)的影響Fig.6 Effect of microscopic parameters σc on macroscopic parameters

圖7 細觀參數(shù)μb對宏觀參數(shù)的影響Fig.7 Effect of microscopic parameters μb on macroscopic parameters

圖8 細觀參數(shù)μ對宏觀參數(shù)的影響Fig.8 Effect of microscopic parameters μ on macroscopic parameters

圖9 細觀參數(shù)R對宏觀參數(shù)的影響Fig.9 Effect of microscopic parameters R on macroscopic parameters

由上述分析可知，細觀參數(shù)Ec，kn/ks，σc和μb與宏觀參數(shù)關聯(lián)度最大的分別是E，v，σt和φ，可由對應的宏觀力學參數(shù)范圍來確定該細觀參數(shù)范圍。細觀參數(shù)C1與宏觀參數(shù)σf和C2關聯(lián)度均較大，考慮到巖石類材料的黏聚力C2差異較大，變化范圍難以確定，故C1的取值范圍由σf的取值范圍來確定。對于大多數(shù)巖石類材料，宏觀力學參數(shù)的大致范圍為[20-21]：E=30～120 GPa，v=0.10～0.35，σf=50～300 MPa，σt=3～25 MPa，φ=30°～65°。故由圖3～9 可大致確定細觀參數(shù)的取值范圍為：Ec=35～110 GPa，kn/ks=1.6～5.3，C1=25～150 MPa，σc=5.3～30 MPa，μb=0.1～0.7。

2.2 基于中心組合設計樣本的細觀參數(shù)計算公式

中心組合設計又稱星點設計，是一種多因素多水平的試驗設計方法，它是在二水平析因設計方法的基礎之上加上極值點和中心點組成的。可見，中心組合設計相對于正交、均勻設計具有較高的實驗精度，并且也不過多地增加實驗次數(shù)，更適合于多水平多因素的試驗設計。為了定量確定PFC2D細觀參數(shù)，采用中心組合設計法建立細觀參數(shù)與宏觀參數(shù)之間的關系式。考慮到待定的細觀參數(shù)有5個(Ec，kn/ks，C1，σc和μb)，利用Design-Expert 軟件生成5 因素5 水平的中心組合設計代碼表(表2)，X1，X2，X3，X4和X5分別代表Ec，kn/ks，C1，σc和μb。每個元素均有5 個水平代碼：-2.38，-1，0，+1 和+2.38，其中，代碼-2.38，0 和+2.38對應的值分別為該因素范圍的最小值、平均值、最大值，代碼-1和+1對應的值可通過“任意2個水平代碼之差的比值與該水平代碼對應值之差的比值相等”原則計算求得。以Ec為例，其范圍為35～110 GPa，則Ec的水平代碼-2.38，0和+2.38所對應的值分別為35，72.5和110。設其水平代碼-1和+1對應的值分別為x和y，根據(jù)上述原則，有

經(jīng)求解得到x=56.74，y=88.26。同理，可求得其他因素水平代碼對應的值，見表3。將表3中各因素水平代碼對應的值代入表2，可得到中心組合設計試驗表，如表4所示。

根據(jù)上述控制變量法分析(圖3～9)、中心組合設計試驗表(表4)可知：E受Ec和kn/ks的影響，v主要受kn/ks的影響，σf受C1，σc和μb的影響，σt受σc和kn/ks的影響，C2主要受C1和μb的影響，φ主要受C1和μb的影響。對6個數(shù)值模型宏觀力學性能(E，v，σf，σt，C2，φ)和5 個細觀參數(shù)(Ec，kn/ks，C1，σc，μb)，利用SPSS軟件進行非線性回歸，得到模型宏觀力學參數(shù)與PFC2D細觀參數(shù)之間的關系式：

表2 中心組合設計代碼表Table 2 Central composite design code table

表3 細觀參數(shù)的中心組合設計代碼值Table 3 CCD code value of microscopic parameters

其中，R2為相關系數(shù)。可見，相關系數(shù)R2均較高，表明擬合精度較高，能夠準確地反映巖石宏-細觀參數(shù)的關聯(lián)性。若已知巖石的宏觀力學參數(shù)，則由式(2)可計算出PFC2D細觀參數(shù)，其中Ec，kn/ks和σc由E，v和σt的關系式聯(lián)立求解，C1和μb通過相關系數(shù)R2較高的σf和φ關系式求得。

3 試驗驗證

3.1 巖石單軸和雙軸壓縮試驗

為了驗證宏-細觀參數(shù)關系式(2)，選取產(chǎn)自于云南的紅砂巖開展試驗研究。表5所示為本試驗測得的紅砂巖以及其他4 種巖石材料[22-23]的宏觀力學參數(shù)。在TRW-30001 軸伺服試驗機上開展巖石的單軸(直徑×高度為50 mm×100 mm的圓柱體試件)、雙軸(長×寬×高為100 mm×100 mm×40 mm 的立方體試件)壓縮試驗，其中雙軸壓縮中側(cè)向壓力為6 MPa，加載速率為500 N/s。在試驗過程中自動記錄應力-應變曲線。

根據(jù)不同巖石的宏觀力學參數(shù)(表5)和模型宏觀力學參數(shù)與PFC2D細觀參數(shù)之間的關系式，可分別計算得到相應的PFC2D細觀參數(shù)(Ec，kn/ks，C1，σc和μb)，如表6所示。采用該細觀參數(shù)，分別開展不同巖石的單軸壓縮(E，ν，σf)、單軸拉伸(σt)和雙軸壓縮(C2，φ)模擬計算，得到的巖石宏觀力學參數(shù)(E，v，σf，σt，C2和φ)見表7，并與試驗結(jié)果進行對比。從表7可見：各參數(shù)的實驗結(jié)果與模擬結(jié)果相對誤差均小于10%，模擬精度較高，從而驗證了采用中心組合設計建立的模型宏觀力學參數(shù)與PFC2D細觀參數(shù)之間的關系式是合理、有效的。

表4 中心組合設計試驗表Table 4 Central composite design table

表5 不同巖石的宏觀力學參數(shù)Table 5 Macroscopic mechanical parameters of different rocks

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

圖10所示為紅砂巖在單軸、雙軸壓縮下數(shù)值模擬計算得到的應力-應變曲線，并與試驗結(jié)果進行對比。從圖10可見：在單軸、雙軸壓縮下，實測的紅砂巖應力-應變曲線(OABCD)可分為壓密階段(OA)、彈性變形階段(AB)、損傷階段(BC)和破壞階段(CD)。模擬結(jié)果(OB′C′D曲線)中沒有出現(xiàn)試驗曲線的壓密階段，這是因為PFC2D建模時顆粒為連續(xù)體，而實際巖石存在有天然的孔隙、裂紋等缺陷，在壓縮時會有微孔、微裂縫等閉合壓密階段(OA)。計算得到的紅砂巖單軸壓縮彈性模量、抗壓強度分別為11.38 GPa 和59.45 MPa，實測的單軸彈性模量、抗壓強度分別為10.69 GPa 和57.20 MPa，相對誤差分別為6.45%和3.93%；計算得到的雙軸壓縮彈性模量和抗壓強度分別為10.83 GPa和78.20 MPa，實測的雙軸彈性模量和抗壓強度分別為11.82 GPa 和82.40 MPa，相對誤差分別為8.38%和5.10%。可見，彈性模量和抗壓強度模擬結(jié)果(未考慮初始缺陷)均略高于實驗結(jié)果，雙軸壓縮比單軸壓縮的抗壓強度更高。忽略初始壓密階段，單、雙軸壓縮試驗曲線與模擬曲線基本重合，彈性模量、抗壓強度均較吻合，從而驗證了采用宏-細觀參數(shù)關系式(2)計算PFC2D 細觀參數(shù)方法的合理性和有效性。

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

圖10 紅砂巖單、雙軸壓縮應力-應變曲線Fig.10 Stress-strain curves of sandstone

3.2 巖石多裂紋斷裂試驗

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

為了進一步驗證宏-細觀參數(shù)關系式(2)，對含3條平行裂紋(裂紋長度為20 mm，寬度為1 mm)的紅砂巖試件進行斷裂試驗。采用與上述壓縮試驗相同的紅砂巖，宏觀力學參數(shù)、細觀參數(shù)與表5和表6所示的相同。試件是長×寬×高為100 mm×100 mm×20 mm的立方體，預制3條等長平行裂紋AB，CD和EF，裂紋長度為20 mm，寬度為1 mm，傾角α=45°，間距為15 mm，如圖11所示。在DNS100 伺服加載系統(tǒng)上進行單軸壓縮斷裂試驗，最大荷載為1 000 kN，以位移加載速率為0.1 mm/s進行加載，采用攝像機記錄試件斷裂破壞過程。

圖12和圖13所示分別為單軸壓縮下紅砂巖試件細觀裂紋演變和宏觀裂紋擴展軌跡的數(shù)值模擬結(jié)果，其中，黑色、紅色、白色分別代表顆粒、細觀裂紋、宏觀裂紋。從圖12和圖13可見：在初始加載時，3 條裂紋的6 個裂紋尖端都會產(chǎn)生初始損傷而出現(xiàn)微裂紋(圖12(a))，此時沒有宏觀裂紋(圖13(a))；隨著進一步加載，各裂紋尖端的損傷不斷積累會產(chǎn)生越來越多的微裂紋并集聚(圖12(b))，由于3條裂紋之間的損傷區(qū)相互作用，裂紋之間的微裂紋首先貫通而形成了最初的2 條宏觀裂紋BC和DE(圖13(b))。當外載不斷增加到破壞荷載時，上部裂紋AB的尖端A和下部裂紋EF的尖端F附近微裂紋急劇增多(圖12(c))，最后形成了2條沿軸向加載方向的宏觀主裂紋(尖端A向上延伸、尖端F向下延伸)以及衍生的次生裂紋，破壞軌跡如圖13(c)所示。

圖12 單軸壓縮下含3條平行裂紋的紅砂巖試件細觀裂紋演變過程(模擬結(jié)果)Fig.12 Micro-crack evolutions of sandstone specimen with three parallel cracks under uniaxial compression(simulation results)

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

圖14所示為攝像機記錄的單軸壓縮下紅砂巖試件裂紋擴展軌跡試驗結(jié)果。從圖14可見：在加載初期，各裂紋尖端會有局部損傷產(chǎn)生微裂紋，但沒有出現(xiàn)宏觀裂紋(圖14(a))；隨著外載增加，3條裂紋之間的微裂紋聚集、演變、貫通而產(chǎn)生了2條宏觀裂紋；當外載不斷增加到破壞荷載時，裂紋尖端A和尖端F的微裂紋不斷發(fā)育，最后形成了沿軸向方向延伸的2 條主裂紋(尖端A向上延伸、尖端F向下延伸)和衍生的次生裂紋，破壞軌跡如圖14(c)所示，其中，試件的上部碎片是試驗機卸載時能量釋放沖壓所致。

比較圖13和圖14可見：除了最后破壞階段的次生裂紋有細微差別外(這是實際巖石材料的非均質(zhì)性所致)，單軸壓縮下含3 條平行裂紋的紅砂巖試件裂紋擴展軌跡試驗結(jié)果與PFC2D模擬結(jié)果較吻合，從而進一步驗證了采用宏-細觀參數(shù)關系式(2)計算PFC2D細觀參數(shù)方法的合理性和可行性。

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

4 結(jié)論

1)采用中心組合設計法和回歸分析，建立了巖石的PFC2D平直節(jié)理模型宏-細觀參數(shù)關系式，可定量計算巖石的細觀參數(shù)。與現(xiàn)有的試錯法和正交設計方法相比，該方法試驗次數(shù)較少，精度較高。

2)PFC2D細觀參數(shù)對巖石的宏觀力學參數(shù)的影響規(guī)律為：Ec對宏觀彈性模量E影響最大，kn/ks對泊松比ν影響最大，C1對σf和C2影響最大，μb對φ影響最大，它們各自呈正相關關系。

3)通過開展單、雙軸壓縮試驗和斷裂試驗，測得了完整巖石試件的應力-應變曲線、宏觀力學參數(shù)和含多裂紋巖石試件的斷裂軌跡，該試驗結(jié)果與PFC2D模擬計算結(jié)果較吻合，從而驗證了所建立的PFC2D平直節(jié)理模型宏-細觀參數(shù)關系式的有效性。

4)采用中心組合設計來確定巖石PFC2D平直節(jié)理模型宏-細觀參數(shù)關系的方法，可推廣到其他脆性材料的細觀參數(shù)確定，為PFC 定量分析多裂紋多場耦合斷裂的細觀機理提供可靠依據(jù)。