成晨，林杭

(中南大學(xué)資源與安全工程學(xué)院，湖南長(zhǎng)沙 410083)

在巖體工程的研究中，研究對(duì)象常常是一些大型的巖石構(gòu)造體，如大型邊坡穩(wěn)定性分析，深部隧道開(kāi)挖，地下開(kāi)采等，這些工程很難通過(guò)室內(nèi)試驗(yàn)進(jìn)行研究。隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，利用數(shù)值模擬可實(shí)現(xiàn)對(duì)巖體數(shù)值模型施加室內(nèi)試驗(yàn)無(wú)法實(shí)現(xiàn)的復(fù)雜應(yīng)力條件，因此該方法在巖體工程領(lǐng)域得到了廣泛的應(yīng)用。離散元法(DEM)是解決粒狀和不連續(xù)材料工程問(wèn)題的一種有效數(shù)值方法,常用于模擬巖石材料從微觀破裂到宏觀破裂的演化過(guò)程，可直觀地展現(xiàn)材料的破壞行為以及裂隙的擴(kuò)展[1]。利用顆粒離散元數(shù)值模擬研究巖體的力學(xué)性質(zhì)時(shí)，需要進(jìn)行全面的標(biāo)定工作，標(biāo)定過(guò)程包括變形和強(qiáng)度的參數(shù)修正。黏結(jié)的顆粒模型需要按照給定的材料確定的具體微觀參數(shù)進(jìn)行制備，以便再現(xiàn)特定的宏觀材料性能，這一過(guò)程直接影響最終模擬結(jié)果的可靠性。另一方面，接觸模型的選擇也會(huì)影響微觀參數(shù)匹配的效果和顆粒模型的宏觀力學(xué)表現(xiàn)，所以數(shù)值模型的標(biāo)定研究對(duì)巖石數(shù)值模擬試驗(yàn)有著重要意義。模型微觀參數(shù)一般采用試錯(cuò)法或者反演計(jì)算法標(biāo)定，該方法通常耗時(shí)費(fèi)力并且盲目性大。因此，國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了簡(jiǎn)化。POTYONDY 等[2]分析了顆粒與接觸的微觀變形強(qiáng)度參數(shù)對(duì)顆粒黏結(jié)模型雙軸、三軸、巴西劈裂試驗(yàn)?zāi)M結(jié)果的敏感性。YOON等[3-4]在使用PFC 模擬巖石單軸壓縮試驗(yàn)時(shí)分別測(cè)試了線性平行黏結(jié)模型的微觀參數(shù)對(duì)單軸抗壓強(qiáng)度，楊氏模量和泊松比的敏感性，通過(guò)統(tǒng)計(jì)中央復(fù)合設(shè)計(jì)方法與Plackett-Burman(P-B)設(shè)計(jì)方法相結(jié)合應(yīng)用計(jì)算機(jī)程序生成了最佳微觀參數(shù)集。李坤蒙等[5]提出了基于物理試驗(yàn)響應(yīng)的參數(shù)確定方法，利用物理試驗(yàn)響應(yīng)反演計(jì)算顆粒間的接觸本結(jié)構(gòu)參數(shù)。鄧樹新等[6]利用P-B 設(shè)計(jì)和響應(yīng)曲面法結(jié)合數(shù)學(xué)規(guī)劃建立了平行黏結(jié)模型的接觸模量，黏結(jié)強(qiáng)度等模型微觀參數(shù)的標(biāo)定方法。鄧樹新等[6]利用P-B設(shè)計(jì)和響應(yīng)曲面法結(jié)合數(shù)學(xué)規(guī)劃建立了平行黏結(jié)模型的接觸模量，黏結(jié)強(qiáng)度等模型微觀參數(shù)的標(biāo)定方法。WANG 等[7]通過(guò)敏感性分析和全局優(yōu)化過(guò)程來(lái)標(biāo)定黏結(jié)離散元的宏觀力學(xué)性能。POTYONDY[8]提出考慮顆粒的多邊形結(jié)構(gòu)的平板模型，用于取代平行黏結(jié)模型模擬典型硬巖，模擬結(jié)果表明平板顆粒黏結(jié)模型更加符合巖石破壞的微觀和宏觀機(jī)制。WU 等[9]基于黏結(jié)顆粒模型(BPM)的3 個(gè)內(nèi)在問(wèn)題，提出利用平板模型解決黏結(jié)顆粒模型的抗壓拉強(qiáng)度比、內(nèi)摩擦角的校正問(wèn)題。上述研究在黏結(jié)顆粒模型微觀參數(shù)對(duì)宏觀力學(xué)參數(shù)的影響，以及接觸微觀參數(shù)標(biāo)定方法上有較為豐富的成果，但存在以下問(wèn)題：首先，多數(shù)巖石模型參數(shù)標(biāo)定研究沒(méi)有考慮巖石的抗拉強(qiáng)度，強(qiáng)度包絡(luò)線等巖石力學(xué)宏觀特性的校正；其次，目前大多參數(shù)標(biāo)定是基于平行黏結(jié)模型，而平行黏結(jié)顆粒模型已被多數(shù)學(xué)者證實(shí)有抗拉強(qiáng)度失真的問(wèn)題，能解決該問(wèn)題的平板模型目前研究較少，其部分接觸微觀參數(shù)與宏觀相應(yīng)之間的的關(guān)系尚不明確，需要進(jìn)一步研究。針對(duì)上述問(wèn)題，本文采用Box-Behnken設(shè)計(jì)試驗(yàn)法，以多元線性模型擬合模擬結(jié)果，得到能夠反映平板接觸模型微觀參數(shù)對(duì)宏觀力學(xué)參數(shù)敏感性的關(guān)系方程，通過(guò)反解關(guān)系方程組，提出平板模型4個(gè)主要微觀參數(shù)的算法方程。

1 數(shù)值模擬

1.1 模型選擇

在PFC 軟件內(nèi)置的接觸模型中，接觸黏結(jié)模型(CBM)和平行黏結(jié)模型(PBM)常被用于模擬完整巖石，2 個(gè)模型的差異主要體現(xiàn)在力的相互作用和傳遞上。接觸黏結(jié)模型只在顆粒接觸點(diǎn)之間起作用而線性平行黏結(jié)模型是在接觸面上通過(guò)以接觸點(diǎn)為中心的有限矩形或圓盤橫截面上活動(dòng)，并且接觸黏結(jié)模型只能傳遞力的作用，無(wú)法抵抗彎矩或者滾動(dòng)，而線性平行黏結(jié)模型能夠同時(shí)傳遞顆粒之間的力與力矩。當(dāng)接觸破壞時(shí)，線性平行黏結(jié)模型接觸剛度能立即下降，而接觸黏結(jié)模型的接觸剛度則會(huì)一直保持作用。因?yàn)檫@些優(yōu)點(diǎn)，目前許多模擬研究選擇平行黏結(jié)模型來(lái)模擬完整巖體。而SHARRKOCK[10]對(duì)人造節(jié)理巖石的研究中明確指出，無(wú)論是對(duì)大尺寸試塊還是小尺寸試塊，平行黏結(jié)模型都無(wú)法實(shí)現(xiàn)對(duì)其抗拉強(qiáng)度的校準(zhǔn)(模擬抗拉強(qiáng)度約是實(shí)際抗拉強(qiáng)度的5 倍)。為更好地解決拉壓比的失真問(wèn)題，POTYONDY[8]提出了考慮顆粒的多邊形結(jié)構(gòu)平板模型(FJM)，當(dāng)顆粒相接觸時(shí)，會(huì)產(chǎn)生一個(gè)抽象的接觸表面，該接觸面被離散化為元素并且可以產(chǎn)生相對(duì)旋轉(zhuǎn)。接觸面的每個(gè)元素都具有獨(dú)立的狀態(tài)。3 種模型黏結(jié)示意圖見(jiàn)圖1。吳順川對(duì)平板模型材料的研究里證明了平板模型能夠解決顆粒黏結(jié)模型的拉壓比失真的問(wèn)題，所以本文選擇平板黏結(jié)顆粒模型作為主要研究對(duì)象。

圖1 顆粒接觸模型示意圖Fig.1 Schematic diagram of particle contact model

1.2 模型建立

為獲得模型宏觀參數(shù),首先需要建立一個(gè)與標(biāo)準(zhǔn)試驗(yàn)巖體樣品尺寸相同的顆粒黏結(jié)模型(BPM)，參考POTYONGDY 等[2]的模型建立方法，首先生成顆粒體系、之后按順序進(jìn)行伺服壓密、處理游離顆粒以及平衡系統(tǒng)等步驟(見(jiàn)圖2)，最后建立如圖3 所示的尺寸為φ50×100 mm3，共包含10 563 個(gè)圓形顆粒的顆粒黏結(jié)模型試樣。

1.3 單軸壓縮試驗(yàn)

試驗(yàn)加載前在模型中心設(shè)定一個(gè)測(cè)量圓，用于記錄每個(gè)計(jì)算步數(shù)的橫向應(yīng)變速率與縱向應(yīng)變速率以及材料中心的平均應(yīng)力。接著給上下墻設(shè)置一個(gè)0.01 m/s的加載速率，通過(guò)記錄上墻體收到的壓力來(lái)記錄模型受到的軸向應(yīng)力，在瞬時(shí)軸向應(yīng)力降為其峰值的70%時(shí)終止循環(huán)，輸出計(jì)算過(guò)程中的軸向應(yīng)力-應(yīng)變曲線。模型的泊松比可由測(cè)量圓測(cè)得的橫向應(yīng)變速率比縱向應(yīng)變速率計(jì)算得到。

1.4 直接拉伸模擬試驗(yàn)

直接拉伸試驗(yàn)?zāi)M采用與單軸壓縮實(shí)試驗(yàn)相同的材料模型，刪除上部與下部墻壁，并將材料頂部和底部的一排顆粒設(shè)置為“抓手”(見(jiàn)圖3)。頂部手柄向上以0.01 m/s的速度移動(dòng)，同時(shí)底部手柄向下以相同的速度移動(dòng)，以施加拉伸載荷。利用中心設(shè)置的測(cè)量圓的平均圓內(nèi)接觸力來(lái)計(jì)算模型中心拉應(yīng)力。

圖3 單軸壓縮模擬試驗(yàn)與直拉模擬試驗(yàn)Fig.3 Uniaxial compression simulation experiment and direct tension simulation experiment

2 微觀參數(shù)匹配

單元試件模型建立后，以試驗(yàn)獲得的巖石宏觀力學(xué)參數(shù)為依據(jù)，標(biāo)定模型的微觀本構(gòu)參數(shù)。為解決顆粒黏結(jié)模型模擬巖石抗拉強(qiáng)度失真的問(wèn)題，本文選取抗拉強(qiáng)度σt，單軸抗壓強(qiáng)度σc，彈性模量E和泊松比λ作為宏觀標(biāo)定參數(shù)。

2.1 常規(guī)匹配流程

參考Gutiérrez-Ch J.G(2018)[5]對(duì)完整巖石的標(biāo)定過(guò)程，分別標(biāo)定平板模型的法向與切向接觸剛度比kˉ*，有效模量Eˉ*，黏結(jié)抗拉強(qiáng)度σˉc與黏結(jié)內(nèi)聚力cˉ。遵循控制變量的原則，每組模擬試驗(yàn)改變一個(gè)微觀參數(shù)的值，保持其他參數(shù)不變，取取值范圍中間值。進(jìn)行一系列的單軸壓縮和拉伸模擬試驗(yàn)后得到了如圖5所示的微觀參數(shù)與宏觀參數(shù)之間的關(guān)系。以靈敏度分析的結(jié)果為指導(dǎo)，在GUTIéRREZ-Ch J.G 所提出的標(biāo)定過(guò)程中加入對(duì)抗拉強(qiáng)度的校準(zhǔn)，新的校準(zhǔn)順序遵循圖4 所示的流程圖。

圖4 微觀參數(shù)的分步迭代過(guò)程Fig.4 Iterative process of microscopic parameter

對(duì)于平行黏結(jié)模型參數(shù)的匹配，通常令模型的黏結(jié)抗拉強(qiáng)度與黏結(jié)內(nèi)聚力相等，BAHAADDINI 等[4]用平行黏結(jié)顆粒模型標(biāo)定Hawkesbury 砂巖的宏觀力學(xué)參數(shù)[10]時(shí)，得到的顆粒模型抗拉強(qiáng)度遠(yuǎn)大于真實(shí)巖體的抗拉強(qiáng)度，本文通過(guò)考慮抗拉強(qiáng)度校正的參數(shù)匹配方法(圖4)再次對(duì)平行黏結(jié)顆粒材料和平板黏結(jié)顆粒材料進(jìn)行標(biāo)定，其標(biāo)定得到的顆粒黏結(jié)材料宏觀屬性見(jiàn)表1。

由表1中本文的標(biāo)定結(jié)果可知，平行黏結(jié)模型依然不能同時(shí)滿足抗拉強(qiáng)度與抗壓強(qiáng)度的校正，這是因?yàn)榭估瓘?qiáng)度校正完成后，平行黏結(jié)模型材料能達(dá)到的單軸抗壓強(qiáng)度會(huì)產(chǎn)生一個(gè)最大限制，這個(gè)限制值遠(yuǎn)小于他所對(duì)應(yīng)真實(shí)巖體的抗壓強(qiáng)度，而平板模型材料不會(huì)存在這樣的限制(見(jiàn)圖6)。

圖6 接觸模型內(nèi)聚力對(duì)宏觀單軸壓縮強(qiáng)度的影響Fig.6 Effect of cohesion of contact model on uniaxial compressive strength of the macroscopic material

表1 平行黏結(jié)模型材料和平板模型材料的宏觀屬性Table 1 Macroscopic properties of parallel bond model materials and flat joint model materials

2.2 優(yōu)化試驗(yàn)設(shè)計(jì)及數(shù)值試驗(yàn)結(jié)果

雖然上述的迭代匹配方法能夠得到較為準(zhǔn)確的平板模型參數(shù)匹配結(jié)果，但是匹配過(guò)程需要嚴(yán)格遵守迭代的順序，只有在前一個(gè)參數(shù)校準(zhǔn)成功后，才能進(jìn)行下一個(gè)參數(shù)的校準(zhǔn)。并且若要對(duì)其他不同類型的巖石進(jìn)行標(biāo)定，整個(gè)標(biāo)定過(guò)程需要重新開(kāi)始，這將是個(gè)十分耗時(shí)并且繁瑣的工作。

目前試驗(yàn)設(shè)計(jì)方法(DOE)已被用于開(kāi)發(fā)更有效的DEM 標(biāo)定方法[12]。DOE 被定義為一種結(jié)構(gòu)化的、有組織的方法，用于確定影響過(guò)程的因素與該過(guò)程的輸出之間的關(guān)系[3]。故本文通過(guò)數(shù)理統(tǒng)計(jì)學(xué)中的試驗(yàn)設(shè)計(jì)，基于前文得到的微觀參數(shù)與宏觀參數(shù)之間的關(guān)系，針對(duì)多種巖石數(shù)值模擬研究中的微觀參數(shù)標(biāo)定問(wèn)題進(jìn)行進(jìn)一步的探討?，F(xiàn)已有許多不同的DOE 方法，方法選擇主要取決于分析的目的和要調(diào)查的因素的數(shù)量。Box-Behnken 設(shè)計(jì)試驗(yàn)法(BBD)可以對(duì)3～7 個(gè)因子進(jìn)行檢驗(yàn)，檢驗(yàn)次數(shù)為15～62 次。當(dāng)因子數(shù)目相同時(shí)所需的試驗(yàn)次數(shù)少于中心組合設(shè)計(jì)的試驗(yàn)次數(shù)，另外BBD法不存在軸向點(diǎn)，所以在設(shè)置因子水平范圍時(shí)可避免出現(xiàn)負(fù)值[13]。故選擇BBD 法簡(jiǎn)化參數(shù)標(biāo)定過(guò)程，設(shè)計(jì)過(guò)程如下：首先選取單軸抗壓強(qiáng)度USC，彈性模量E，泊松比PR和抗拉強(qiáng)度DTS作為試驗(yàn)指標(biāo)。若不考慮顆粒自身屬性和尺寸效應(yīng)的影響，只考慮接觸模型的變形與強(qiáng)度參數(shù)，選取黏結(jié)有效模量Eˉ(EC)，黏結(jié)剛度kˉn/kˉs(KK)，黏結(jié)抗拉強(qiáng)度σc(TS)，黏結(jié)內(nèi)聚力cˉ(CS)4 個(gè)平板模型的主要微觀參數(shù)作為試驗(yàn)因素。

最后計(jì)算試驗(yàn)因素范圍，根據(jù)3.1 數(shù)值試驗(yàn)結(jié)果進(jìn)行線性擬合：得到宏觀彈性模量與EC的關(guān)系式：y=0.648x-0.006。根據(jù)目標(biāo)標(biāo)定對(duì)象宏觀參數(shù)值，設(shè)計(jì)的巖石宏觀抗拉強(qiáng)度范圍為5～11 MPa。將y=2 和y=10 代入上式，得到x=7.79 和x=17.18，因此黏結(jié)抗拉強(qiáng)度TS 編碼范圍為7.5～17.5 MPa，同理，設(shè)計(jì)彈性模量范圍取50～80 GPa，泊松比范圍取0.15～0.3，單軸壓縮強(qiáng)度范圍取110～220 MPa，分別代入到圖5 中宏觀參數(shù)對(duì)應(yīng)線性擬合公式，得到黏結(jié)剛度KK，黏結(jié)有效模量EC以及黏結(jié)內(nèi)聚力CS的編碼范圍，具體取值范圍以及編碼公式見(jiàn)表2。

表2 微觀參數(shù)及其標(biāo)定范圍Table 2 Microscopic parameters and calibration ranges

圖5 平板模型的微觀參數(shù)與模型宏觀力學(xué)參數(shù)之間的敏感性分析Fig.5 Sensitivity analysis between the microscopic parameters of the flat joint model and the macro-mechanical parameters of the model

根據(jù)設(shè)計(jì)方案進(jìn)行模擬試驗(yàn)，基于單軸壓縮數(shù)值試驗(yàn)與直接拉伸數(shù)值試驗(yàn)的結(jié)果，計(jì)算得到每組數(shù)值試驗(yàn)的單軸抗壓強(qiáng)度UCS，彈性模量E，泊松比PR和抗拉強(qiáng)度DTS，方案和計(jì)算結(jié)果見(jiàn)表3。

表3 BBD設(shè)計(jì)矩陣以及數(shù)值模擬結(jié)果Table 3 BBD design matrice and numerical simulation results

3 結(jié)果與討論

3.1 統(tǒng)計(jì)分析

為簡(jiǎn)化求解，本文均采用如式(1)所示的線型模型，使用Design Expert 軟件別對(duì)泊松比、楊氏模量、單軸抗壓強(qiáng)度和抗拉強(qiáng)度4個(gè)宏觀力學(xué)響應(yīng)進(jìn)行擬合分析，擬合結(jié)果見(jiàn)表4，R2與精密度等統(tǒng)計(jì)特征通常作為模型選擇的依據(jù)。R2越接近1，說(shuō)明模型擬合效果越好，并且要求精密度比率大于4，模型才是合理有效的。由擬合結(jié)果可知R2與精度均滿足擬合條件，4 個(gè)宏觀參數(shù)均可以采用線性模型來(lái)擬合。最后分別計(jì)算每個(gè)宏觀參數(shù)對(duì)應(yīng)微觀參數(shù)的影響系數(shù)見(jiàn)表4。

表4 泊松比、彈性模量、抗拉強(qiáng)度和單軸壓縮強(qiáng)度的BBD設(shè)計(jì)多元線性模型分析結(jié)果Table 4 Analysis results of multiple linear model of BBD design of poisson’s ratio,elastic modulus,tensile strength and uniaxial compressive strength

Y=I0+A×KK+B×EC+C×TS+D×CS(1)

式中：I0表示常量，A，B，C和D為方程系數(shù)。

3.2 求解方程

由分析結(jié)果得到4 個(gè)宏觀因素與4 個(gè)微觀因素之間的矩陣關(guān)系方程：

式中：

故由方程(2)的解X=A-1(Y-b)，化簡(jiǎn)得到黏結(jié)剛度比、有效模量、黏結(jié)抗拉強(qiáng)度、黏結(jié)內(nèi)聚力的表達(dá)方程為：

3.3 方法驗(yàn)證

根據(jù)4種經(jīng)典硬巖的物理試驗(yàn)結(jié)果進(jìn)行微觀參數(shù)的標(biāo)定，代入前文提出的方程解析式(3)得到4種巖石的微觀參數(shù)，如表5所示。

將表5 中的微觀參數(shù)賦予到PFC2D 顆粒黏結(jié)模型中，得到單軸壓縮和拉伸條件下數(shù)值模擬結(jié)果，與物理試驗(yàn)宏觀參數(shù)的結(jié)果進(jìn)行對(duì)比，見(jiàn)圖7[8,14-16]，4 個(gè)指標(biāo)的相對(duì)誤差大部分都在5%左右。證明式(3)對(duì)同一尺寸模型的參數(shù)匹配具有通用性。

表5 平板模型黏結(jié)顆粒微觀參數(shù)Table 5 Microscopic parameters of flat joint bonded model

圖7 巖石宏觀力學(xué)參數(shù)模擬值與試驗(yàn)值對(duì)比Fig.7 Comparison of numerical and experimental values of rock macro-mechanical parameters

4 結(jié)論

1) 顆粒黏結(jié)模型的宏觀力學(xué)參數(shù)與接觸模型微觀參數(shù)存在一些對(duì)應(yīng)影響關(guān)系，并且宏觀力學(xué)參數(shù)隨著主要影響因素的增大而增大。顆粒黏結(jié)材料的泊松比主要受黏結(jié)接剛度比影響，有效模量主要與黏結(jié)有效模量有關(guān)，抗拉強(qiáng)度的變化主要由黏結(jié)抗拉強(qiáng)度控制，單軸壓縮強(qiáng)度受黏結(jié)內(nèi)聚力影響的同時(shí)也受黏結(jié)抗拉強(qiáng)度的控制。

2) 通過(guò)分步迭代調(diào)整接觸模型的黏結(jié)抗拉強(qiáng)度和黏結(jié)內(nèi)聚力時(shí)，平行黏結(jié)模型材料的單軸抗壓強(qiáng)度與抗拉強(qiáng)度之比只能達(dá)到一個(gè)較小的值，而平板模型材料能夠與較高抗壓拉強(qiáng)度比的巖體宏觀響應(yīng)相匹配。

3) 標(biāo)定同一尺寸多種巖石的顆粒黏結(jié)模型參數(shù)時(shí)，可利用BBD試驗(yàn)設(shè)計(jì)方法建立巖石泊松比、彈性模量、抗拉強(qiáng)度和單軸壓縮強(qiáng)度與接觸模型微觀參數(shù)黏結(jié)剛度比、有效模量、黏結(jié)抗拉強(qiáng)度和黏結(jié)內(nèi)聚力之間的線性關(guān)系方程簡(jiǎn)化模擬參數(shù)匹配的步驟。