劉富有，陳鵬宇，余宏明

(1.中國地質(zhì)大學(武漢) 工程學院，武漢　430074；2.河南省地質(zhì)礦產(chǎn)勘查開發(fā)局 第二地質(zhì)礦產(chǎn)調(diào)查院，鄭州　450000；3.內(nèi)江師范學院　地理與資源科學學院，四川 內(nèi)江　641100)

?

基于Flatjoint接觸模型的巖石單軸壓縮和巴西劈裂顆粒流模擬研究

劉富有1,2，陳鵬宇3，余宏明1

(1.中國地質(zhì)大學(武漢) 工程學院，武漢430074；2.河南省地質(zhì)礦產(chǎn)勘查開發(fā)局 第二地質(zhì)礦產(chǎn)調(diào)查院，鄭州450000；3.內(nèi)江師范學院地理與資源科學學院，四川 內(nèi)江641100)

基于顆粒流平臺，選擇Flatjoint接觸模型作為巖石模擬的基本模型，用單軸壓縮和巴西劈裂數(shù)值試驗測試巖石宏觀參數(shù)，對數(shù)值試驗進行正交設(shè)計，采用多因素方差分析和回歸分析研究宏細觀參數(shù)之間的關(guān)系。在此基礎(chǔ)上，提出了基于巖石單軸壓縮和巴西劈裂數(shù)值試驗的細觀參數(shù)匹配方法。以灰?guī)r的室內(nèi)試驗為基礎(chǔ)，對灰?guī)r細觀參數(shù)進行了標定，實現(xiàn)了灰?guī)r單軸壓縮和巴西劈裂的顆粒流模擬，模擬結(jié)果與試驗結(jié)果相接近，驗證了本文方法的有效性。

巖石；顆粒流；Flatjoint接觸模型；單軸壓縮；巴西劈裂

1　研究背景

顆粒流程序(particleflowcode，PFC)是一種離散單元法，其通過圓形顆粒之間的相互作用來模擬材料的宏觀力學性質(zhì)。最初，顆粒流程序中只有接觸黏結(jié)模型和平行黏結(jié)模型2種顆粒接觸黏結(jié)模型[1]。在實際應用中發(fā)現(xiàn)這2種模型計算得出的單軸抗拉強度都較高，導致單軸抗壓強度和單軸抗拉強度比值UCS/TS值為3～4。而對于很多巖石，UCS/TS常常超過10。為此，N.Cho等[2]提出了簇平行黏結(jié)模型，它將多個顆粒固結(jié)成簇，形成一個不規(guī)則的大顆粒，簇中單個顆粒的旋轉(zhuǎn)被抑制，這使UCS/TS值顯著增大。為了更好地解決這個問題，D.O.Potyondy[3]又提出了一種適用于硬質(zhì)巖石的平直節(jié)理顆粒黏結(jié)模型(Flat-JointedBonded-ParticleMaterial,FJBPM)，平直節(jié)理顆粒黏結(jié)模型將圓形顆粒構(gòu)造成多邊形顆粒，顆粒破壞后的旋轉(zhuǎn)被抑制。隨著顆粒流數(shù)值模擬方法的發(fā)展，其應用范圍也越來越廣泛，特別是巖石破壞的細觀模擬分析[4-6]。但是上述方法大都采用顆粒流程序中的平行黏結(jié)模型。為此，本文選用FJBPM作為基本模型，研究宏細觀參數(shù)之間的關(guān)系，提出了細觀參數(shù)的匹配方法，在此基礎(chǔ)上進行了灰?guī)r單軸壓縮與巴西劈裂破壞的顆粒流模擬。

2　顆粒流基本理論

PFC2D中通過設(shè)置每個接觸的本構(gòu)模型，從而實現(xiàn)對巖石材料宏觀本構(gòu)性質(zhì)的模擬。Flatjoint接觸模型能夠抑制黏結(jié)破壞后顆粒的旋轉(zhuǎn)。典型的Flatjoint接觸模型如圖1所示。

(a)假設(shè)的多邊形顆粒形態(tài)

(b)平直節(jié)理交界面形態(tài)圖1　平直節(jié)理接觸模型[3]Fig.1　Flatjoint contact model[3]

Flatjoint接觸模型有未黏結(jié)和黏結(jié)2種模式，兩者的本構(gòu)關(guān)系不同，具體可參見文獻[3]。

選定Flatjoint接觸模型，可通過巖石數(shù)值試驗與室內(nèi)物理試驗的對比進行細觀參數(shù)的標定。單軸壓縮數(shù)值試驗是經(jīng)宏觀參數(shù)得到細觀參數(shù)的最重要的途徑之一，該數(shù)值試驗主要可以得到單軸抗壓強度、變形模量和泊松比等參數(shù)；另一個確定細觀參數(shù)的重要數(shù)值試驗就是巴西劈裂數(shù)值試驗。根據(jù)數(shù)值試樣破裂時的峰值作用力Ff，可得數(shù)值試樣的抗拉強度。本文選擇這2個數(shù)值試驗作為細觀參數(shù)標定的基礎(chǔ)。

3　細觀參數(shù)的標定

3.1正交試驗設(shè)計

在進行正交試驗設(shè)計之前，考慮到接觸本構(gòu)模型細觀參數(shù)較多，增加了數(shù)值試驗的數(shù)量和細觀參數(shù)標定的難度，因此，有必要進行一定的假設(shè)，以減少細觀參數(shù)的數(shù)量。Flatjoint接觸模型的細觀參數(shù)包括：N，Ec，kn/ks，μb，λ，σb，cb，φb。其中：N為交界面段數(shù)；Ec為平直節(jié)理模量；kn/ks為平直節(jié)理剛度比；μb為平直節(jié)理摩擦系數(shù)；λ為平直節(jié)理兩端較小顆粒的半徑比；σb為平直節(jié)理抗拉強度；cb為平直節(jié)理黏聚力；φb為平直節(jié)理內(nèi)摩擦角。

參考D.O.Potyondy[3]，B.A.Poulsenn等[7]的研究，作以下假設(shè)：①平直節(jié)理抗拉強度σb小于抗剪強度τb，即σb<τb=cb+σtanφb，σ為平直節(jié)理所受正應力，可獲得符合實際的巖石拉壓強度比；②采用簡化的平直節(jié)理破壞準則，即τb=cb，φb=0，可以有效控制平直節(jié)理抗剪強度與抗拉強度比；③平直節(jié)理抗拉強度和抗剪強度兩者的標準差與平均值的比值相等，取SD/Mean=0.25；④平直節(jié)理兩端較小顆粒的半徑比λ=1；⑤取最小半徑Rmin=0.4mm，顆粒半徑比固定為1.66；⑥平直節(jié)理摩擦系數(shù)μ=0.5。

由此，根據(jù)所需確定的細觀參數(shù)建立正交試驗設(shè)計表1，設(shè)計的正交矩陣序列如表2所示。分別進行單軸壓縮數(shù)值試驗和直接拉伸數(shù)值試驗(試樣寬50mm，高100mm)，由此得到宏觀參數(shù)：單軸抗壓強度σf、彈性模量E、泊松比ν和抗拉強度σt，所得結(jié)果如表2所示，可以看出該結(jié)果符合大部分巖石的宏觀參數(shù)取值范圍，并且平直節(jié)理模型的壓拉比σf/σt大部分滿足10～20。

3.2多因素方差分析

采用SPSS軟件對表2的數(shù)據(jù)進行多因素方差分析。圖2給出了多因素方差分析的F統(tǒng)計量的結(jié)果，其中X1表示Ec，X2表示kn/ks，X3表示σb，X4表示τb/σb，X5表示N。

表1　Flatjoint接觸模型細觀參數(shù)正交試驗設(shè)計值Table 1　Orthogonal experimental design of mesoscopicparameters of Flatjoint model

表2　Flatjoint接觸模型細觀參數(shù)正交設(shè)計的矩陣序列和 宏觀參數(shù)計算結(jié)果Table 2　Orthogonal matrix sequence of mesoscopic parametersof Flat-jointed model and the corresponding macroscopicparameters

根據(jù)F統(tǒng)計量和相伴概率值Sig.的結(jié)果(取顯著性水平為α=0.05)可得到宏觀參數(shù)的顯著性影響因素及其排序結(jié)果如下：①彈性模量E,Ec>kn/ks>N；②泊松比υ,kn/ks>τb/σb;③單軸抗壓強度σf,σb>τb/σb>N;④抗拉強度σt,σb>kn/ks;⑤壓拉強度比σf/σt,τb/σb>kn/ks。

根據(jù)圖2的結(jié)果，可以得到以下結(jié)論：

(1)對于平直節(jié)理模型，彈性模量E主要受到節(jié)理模量Ec、節(jié)理剛度比kn/ks和交界面段數(shù)N的影響，其中以節(jié)理模量Ec最為顯著，節(jié)理剛度比kn/ks次之，交界面段數(shù)N的影響相對較弱，可以忽略。

(a)彈性模量E

(b)泊松比υ

(c)單軸抗壓強度σf

(d)抗拉強度σt

(e)壓拉強度比σf/σt圖2　平直節(jié)理模型宏細觀參數(shù)多因素方差分析的F統(tǒng)計量Fig.2　F statistics of multi-factor analysis of variance formacro-and-meso-scopic parameters of Flatjoint model

(2)泊松比υ主要受到節(jié)理剛度比kn/ks和節(jié)理強度比τb/σb的影響，節(jié)理剛度比kn/ks的影響遠大于節(jié)理強度比τb/σb。

(3)單軸抗壓強度主要受到節(jié)理抗拉強度σb、強度比τb/σb和交界面段數(shù)N的影響，以節(jié)理抗拉強度σb最為顯著，節(jié)理強度比τb/σb次之，交界面段數(shù)N影響最弱。

(4)抗拉強度σt主要受到節(jié)理抗拉強度σb和節(jié)理剛度比kn/ks的影響，但是節(jié)理剛度比kn/ks相比于節(jié)理抗拉強度影響較弱，可以忽略。

(5)壓拉強度比σf/σt主要受到節(jié)理強度比τb/σb和節(jié)理剛度比kn/ks的影響，節(jié)理強度比τb/σb的影響遠大于節(jié)理剛度比kn/ks。

3.3回歸分析

根據(jù)上述正交試驗計算結(jié)果可以建立宏觀參數(shù)與其主要顯著影響因素之間的關(guān)系式如表3所示，其中由于交界面段數(shù)N對宏觀參數(shù)的影響相對于其他參數(shù)較小，并未考慮該因素。可以看出，宏觀參數(shù)擬合公式的R2值在0.917～0.992之間，說明擬合效果良好，能夠比較準確反映宏觀參數(shù)與細觀參數(shù)之間的關(guān)系。其中，彈性模量E與節(jié)理模量Ec呈正相關(guān)關(guān)系，與節(jié)理剛度比kn/ks的自然對數(shù)呈負相關(guān)關(guān)系；泊松比v與節(jié)理剛度比kn/ks和強度比τb/σb呈正相關(guān)關(guān)系；單軸抗壓強度σf與節(jié)理抗拉強度σb和強度比τb/σb呈正相關(guān)關(guān)系；抗拉強度σt與節(jié)理抗拉強度σb呈正相關(guān)關(guān)系；壓拉強度比σf/σt與節(jié)理剛度比kn/ks和強度比τb/σb呈正相關(guān)關(guān)系。

表3　宏觀參數(shù)與細觀參數(shù)之間的關(guān)系式及擬合度Table 3　Regression equations of macro-and-meso-scopicparameters and corresponding degrees of fitting

根據(jù)表3中的擬合公式，在已知宏觀參數(shù)的情況下，即可反算對應的細觀參數(shù)，由于σf/σt，σf，σt三者之間具有相關(guān)性，選擇其中2個公式即可，為了計算方便，選擇σf，σt兩者的擬合公式，得到反算公式如表4所示。至于交界面段數(shù)N，可選擇為3～4[3,7]，根據(jù)反算公式初步計算細觀參數(shù)，然后進行數(shù)值試驗計算宏觀參數(shù)，對比計算宏觀參數(shù)與實際宏觀參數(shù)之間的差別，再根據(jù)擬合公式中所反映的宏細觀參數(shù)之間的趨勢性關(guān)系，可對細觀參數(shù)進行適當微調(diào)，直到達到合理的精度范圍。

表4　細觀參數(shù)的反算公式Table 4　Back calculation formulas for mesoscopicparameters

4　實例分析

4.1室內(nèi)試驗

采用焦作市龍寺廢棄礦山奧陶系中統(tǒng)上馬家溝組(O2s)厚層狀灰?guī)r，室內(nèi)試驗是在INSTORN-1346電液伺服巖石力學測試系統(tǒng)上進行，可自動計算宏觀力學參數(shù)。

4.1.1單軸壓縮試驗

圖3為試驗過程中測試系統(tǒng)自動記錄的部分應力-應變曲線，從中可以看出后區(qū)曲線較陡，說明灰?guī)r脆性明顯。圖4為巖石單軸壓縮試驗巖樣破壞形態(tài)，巖樣主要發(fā)生垂向的脆性劈裂破壞以及表層剝離，破壞面方向與軸向荷載加載的方向近乎平行。表5中給出了巖石單軸壓縮試驗成果。

(a)Y-1編號

(b)Y-2編號圖3　巖石單軸壓縮應力-應變曲線Fig.3　Stress-strain curves of rock uniaxial compressive test

圖4　巖石單軸壓縮試驗巖樣破壞形態(tài)Fig.4　Failure patterns of rock samples under uniaxialcompression test表5　巖石單軸壓縮試驗成果Table 5　Result of rock uniaxial compression test

編號抗壓強度/MPa試驗結(jié)果均值彈性模量/GPa試驗結(jié)果均值泊松比試驗結(jié)果均值Y-1111.98Y-2133.55Y-3122.17Y-4119.96121.9172.9668.4051.8751.9061.280.2020.2050.2190.2070.208

4.1.2巖石巴西劈裂試驗

巴西破裂試驗同樣是在INSTORN-1346電液伺服巖石力學測試系統(tǒng)上進行，可自動計算巖樣抗拉強度。圖5為試驗過程中測試系統(tǒng)自動記錄的部分應力-應變曲線，從中可以看出曲線呈倒V形，反映了灰?guī)r的硬脆性特征。圖6為巖石巴西劈裂試驗巖樣破壞形態(tài)，巖樣主要發(fā)生徑向的脆性劈裂破壞，破壞面方向與軸向荷載加載的方向平行。表6給出了巖石巴西劈裂試驗成果。

(a) L-1

(b)L-2圖5　巖石巴西劈裂應力-應變曲線Fig.5　Stress-strain curves of rock brazilian splitting test

圖6　巖石巴西劈裂試驗巖樣破壞形態(tài)Fig.6　Failure patterns of rock samples under Braziliansplitting test表6　巖石巴西劈裂試驗成果Table 6　Result of rock Brazilian splitting test

編號狀態(tài)抗拉強度/MPa試驗值均值L-1L-2L-3L-4天然4.5005.5035.2315.8145.262

4.2巖石破壞的顆粒流數(shù)值模擬

4.2.1巖石細觀參數(shù)的標定

根據(jù)上述試驗結(jié)果，即可進行巖石細觀參數(shù)的標定，宏觀參數(shù)選用巖石試驗結(jié)果的平均值。生成顆粒流試樣時，采用的空隙率為0.1，由此可以根據(jù)巖石密度確定顆粒密度，采用如下公式：

ρs=ρ/(1-r)。

(1)

式中：ρs為顆粒密度；ρ為巖石密度；r為數(shù)值試樣空隙率，取0.1。

在正交試驗時，本文采用的是直接拉伸試驗確定抗拉強度，而巖石試驗采用的是巴西劈裂試驗。因此，在細觀參數(shù)匹配時，采用的是巴西劈裂數(shù)值試驗。巴西劈裂試驗與直接拉伸試驗確定的抗拉強度存在一定差異，但是差異并不是特別大，所以仍可以用表4中公式反算細觀參數(shù)。然后，根據(jù)巖石單軸壓縮數(shù)值試驗對細觀參數(shù)進行微調(diào)，以達到合理的精度范圍。根據(jù)表5和表6的試驗結(jié)果標定各項細觀參數(shù)如表7所示，對應數(shù)值試樣宏觀參數(shù)如表8所示。

表7　巖石細觀參數(shù)標定結(jié)果Table 7　Calibration result of rock mesoscopic parameters

表8　巖石數(shù)值試樣宏觀參數(shù)Table 8　Macroscopic parameters of rock numerical sample

對比表5，表6和表8 的宏觀參數(shù)，從中可以看出，巖石數(shù)值試樣的宏觀參數(shù)和室內(nèi)試驗結(jié)果的平均值是非常接近的，誤差僅在0.48%～2.51%之間。

4.2.2數(shù)值模擬結(jié)果分析

4.2.2.1單軸壓縮數(shù)值試驗

圖7是PFC2D模擬的灰?guī)r單軸壓縮破壞過程的數(shù)值模擬結(jié)果。從圖中可以看出，數(shù)值試樣的破壞主要表現(xiàn)為軸向劈裂以及表層剝離，這與巖石試驗結(jié)果是非常一致的。從應力-應變曲線中可以看出后區(qū)曲線較陡，這與室內(nèi)試驗結(jié)果一致，說明數(shù)值試樣能夠很好地體現(xiàn)灰?guī)r的脆性特征。

(a)破壞試樣

(b)應力-應變曲線圖7　灰?guī)r單軸壓縮數(shù)值模擬結(jié)果Fig.7　Numerical simulation result of uniaxial compressiontest of limestone

4.2.2.2巴西劈裂數(shù)值試驗

圖8是PFC2D模擬的灰?guī)r巴西劈裂破壞過程的數(shù)值模擬結(jié)果。從圖中可以看出，數(shù)值試樣的破壞主要表現(xiàn)為徑向的劈裂破壞，這與巖石試驗結(jié)果是非常一致的。從應力-應變曲線中可以看出后區(qū)曲線較陡，這與室內(nèi)試驗結(jié)果一致，說明數(shù)值試樣能夠很好地體現(xiàn)灰?guī)r的脆性特征。

(a)破壞試樣

(b)荷載-應變曲線圖8　灰?guī)r巴西劈裂數(shù)值模擬結(jié)果Fig.8　Numerical simulation result of Brazilian splittingtest of limestone

5　結(jié)　論

(1)采用正交試驗設(shè)計和多因素方差分析研究巖石細觀參數(shù)和宏觀參數(shù)之間的關(guān)系，得到如下結(jié)論：彈性模量E主要受節(jié)理模量Ec、節(jié)理剛度比kn/ks和交界面段數(shù)N的影響，顯著性排序為Ec>kn/ks>N；泊松比v主要受節(jié)理剛度比kn/ks和節(jié)理強度比τb/σb的影響，顯著性排序為kn/ks>τb/σb；單軸抗壓強度主要受節(jié)理抗拉強度σb、強度比τb/σb和交界面段數(shù)N的影響，顯著性排序為σb>τb/σb>N；抗拉強度σt主要受節(jié)理抗拉強度和節(jié)理剛度比的影響，顯著性排序為σb>kn/ks；壓拉強度比σf/σt主要受到節(jié)理強度比τb/σb和節(jié)理剛度比kn/ks的影響，顯著性排序為τb/σb>kn/ks。

(2)根據(jù)正交試驗計算結(jié)果可以建立宏觀參數(shù)與其主要顯著影響因素之間的關(guān)系式，發(fā)現(xiàn)彈性模量E與節(jié)理模量Ec呈正相關(guān)關(guān)系，與節(jié)理剛度比kn/ks的自然對數(shù)呈負相關(guān)關(guān)系；泊松比v與節(jié)理剛度比kn/ks和節(jié)理強度比τb/σb呈正相關(guān)關(guān)系；單軸抗壓強度σf與節(jié)理抗拉強度σb和節(jié)理強度比τb/σb呈正相關(guān)關(guān)系；抗拉強度σt與節(jié)理抗拉強度σb呈正相關(guān)關(guān)系；壓拉強度比σf/σt與節(jié)理剛度比kn/ks和節(jié)理強度比τb/σb呈正相關(guān)關(guān)系。根據(jù)宏細觀參數(shù)之間的擬合公式得到了細觀參數(shù)的反算公式，由此建立了巖石細觀參數(shù)的標定方法。

(3)進行了灰?guī)r的單軸壓縮與巴西劈裂試驗，試驗結(jié)果顯示灰?guī)r的脆性特征明顯。單軸壓縮主要發(fā)生垂向的脆性劈裂破壞以及表層剝離，破壞面方向與軸向荷載加載的方向近乎平行。巴西劈裂主要發(fā)生徑向的脆性劈裂破壞，破壞面方向與軸向荷載加載的方向平行。建立了灰?guī)r的單軸壓縮與巴西劈裂數(shù)值模型，模擬結(jié)果顯示FJBPM可以很好地體現(xiàn)灰?guī)r的脆性破壞特征，數(shù)值模擬結(jié)果與室內(nèi)試驗結(jié)果比較接近，說明了本文方法的有效性。

[1]POTYONDYDO,CUNDALLPA.ABonded-particleModelforRock[J].InternationalJournalofRockMechanicsandMiningSciences,2004,41(8):1329-1364.

[2]CHON,MARTINCD,SEGODC.AClumpedParticleModelforRock[J].InternationalJournalofRockMechanicsandMiningSciences,2007,44(7):997-1010.

[3]POTYONDYDO.AFlat-jointedBonded-particleMaterialforHardRock[C]//AmericanRockMechanicsAssociation,Proceedingsofthe46thUSRockMechanics/Geome-chanicsSymposium,Chicago,IL,USA,June24-27,2012:1510-1519.

[4]張學朋，王剛，蔣宇靜，等.基于顆粒離散元模型的花崗巖壓縮試驗模擬研究[J].巖土力學，2014，35(增1)：99-105.

[5]劉寧，張春生，褚衛(wèi)江，等.深埋大理巖脆性破裂細觀特征分析[J].巖石力學與工程學報，2012，31(增2)：3557-3565.

[6]唐謙，李云安.圍壓對巖石裂紋擴展影響的顆粒流模擬研究[J].長江科學院院報，2015，32(4)：81-85.

[7]POULSENNBA,ADHIKARYDP.ANumericalStudyoftheScaleEffectinCoalStrength[J].InternationalJournalofRockMechanics&MiningSciences,2013,63:62-71.

(編輯：劉運飛)

PFC Simulation of Uniaxial Compression and Brazilian Splitting Test ofRock Based on Flat-jointed Bonded-particle Material

LIU Fu-you1,2,CHEN Peng-yu3,YU Hong-ming1

(1.FacultyofEngineering,ChinaUniversityofGeosciences,Wuhan430074,China;2.No.2InstituteofGeological&MineralResourcesofHenanProvince,Zhengzhou450000,China;3.SchoolofGeographyandResourcesScience,NeijiangNormalUniversity,Neijiang641100,China)

TheFlat-jointbonded-particlematerialwasselectedasthebasicmodelforrocksimulation.Themacro-parametersofrockwerecalculatedthroughnumericaluniaxialcompressiontestandBraziliansplittingtestinPFC.Orthogonaldesignwasadoptedtoconductnumericaltestsandmulti-factoranalysisofvarianceandregressionanalysiswereutilizedtoanalyzetherelationshipbetweenmacro-parametersandmeso-parameters.Onthisbasis,methodofmatchingthemeso-parametersfornumericaluniaxialcompressiontestandBraziliansplittingtestwasproposed.Accordingtolaboratorytestsoflimestone,themeso-parametersoflimestonewereobtainedandthePFCsimulationofuniaxialcompressionandBraziliansplittingtestsoflimestonewereconducted.Thesimulationresultsareclosetothelaboratorytestresults,whichvalidatedtheeffectivenessoftheproposedmethod.

rock;PFC;Flatjointcontactmodel;uniaxialcompression;Braziliansplitting

2015-07-07;

2015-10-12

國家自然科學基金項目(41272377，41302278)

劉富有(1968-),男,河南鄧州人,教授級高級工程師,博士，主要從事工程地質(zhì)、水文地質(zhì)和環(huán)境地質(zhì)方面的工作，(電話)0371-55156812(電子信箱)79012983@qq.com。

10.11988/ckyyb.20150566

2016,33(09):60-65

TU457

A

1001-5485(2016)09-0060-06