考慮微元強度韋伯分布與裂紋長度冪律排布的砂巖強度預測模型

石浩,張后全,吳疆宇,宋雷,李明,榮傳新,陸鵬舉

(1.安徽理工大學 礦山建設工程安徽省高校重點實驗,安徽 淮南 232001;2.中國礦業大學 深部巖土力學與地下工程國家重點實驗室,江蘇 徐州 221116;3.中國礦業大學 力學與土木工程學院,江蘇 徐州 221116; 4.安徽理工大學 土木建筑學院,安徽 淮南 232001;5.山東科技大學 省部共建礦山巖層智能控制與綠色開采國家重點實驗室培育基地,山東 青島 266590;6.中煤礦山建設集團有限責任公司,安徽 合肥 230000)

隨著煤炭開采工藝的進步、機械化水平的提高以及信息技術的發展,以數字化、信息化為前提和基礎的智慧礦山建設已經取得了階段性進展;相信在不遠的未來,我國必然會實現深部煤炭的無人化開采[1]。智慧礦山的建設對巖石力學特別是強度的準確判識要求極高[2],而巖石內部的各類損傷對其強度特性的影響占重要地位,是引起巖石強度離散的關鍵[3-6]。因此,為了保障智慧礦山的建設以及井下人員、設備的安全,有必要掌握巖石內部的各類損傷對其強度特性的影響規律。

目前,關于損傷缺陷影響巖石強度特性的研究已取得了十分豐碩的成果[7-27],并且根據研究對象的不同可分為2類。第一類研究主要針對巖石微元體強度劣化等微細觀損傷缺陷[7-12]。此類研究基于自然巖石為非均質材料,認為其宏觀破壞對應微細觀破壞由量變到質變的過程[7,9]。專家學者多從唯象學角度出發,假定微元體損傷參量服從韋伯統計分布規律,基于D-P準則或M-C準則建立與溫度、蠕變等特定工況相關的巖石損傷本構模型。相關成果對于豐富巖石損傷描述方法及揭示巖石損傷內在機制具有重要意義,但眾多研究主要是基于理論方法開展的,重在分析韋伯分布參數對應力-應變關系(曲線)的影響[7-10,12]。在本質上所研究巖石試件仍以黑箱的形式存在,無法設定微元體強度或損傷參量的分布規律,導致理論分布假設的合理性得不到充分驗證。第二類研究主要針對裂紋、裂隙等宏觀損傷缺陷[13-27]。此類研究預制缺陷數量差異較大,多數研究預制缺陷數量小于4,也有少量研究預制缺陷數量成百上千。對于預制缺陷數量小于4的研究,研究者主要分析單個裂紋、裂隙的傾角、數量、位置、張開程度以及多個裂紋、裂隙的組合形式等對巖石強度的影響[13-19]。相關成果對于了解缺陷巖石的力學性能具有重要意義,但對缺陷信息(如傾角等)的設置過于理想,也導致研究結果的普適性較弱[20-22]。對于預制缺陷數量成百上千的研究,研究者通常注重對缺陷生成方法(包括概率分布或利用地層掃描和計算機信息提取相結合)的提出和改善[23-26],而對缺陷如何影響巖石強度特性的分析相對不足[25-26]。另外,裂紋長度服從冪律排布的規律也未被充分考慮。趙洪寶等[22]回顧了關于裂隙(紋)巖石的研究,指出應當考慮分布裂紋數量及裂紋間的相互作用并建立與裂紋分布相關的強度理論,才能真正提高研究成果的工程應用價值。以趙洪寶等人的觀點作為指引,本課題組[27]根據Mori-Tanaka理論進行推導,首先建立了與預制裂紋分布信息相關的損傷砂巖強度模型。

鑒于自然界的巖石試件其內部宏細微觀缺陷并存,僅對宏觀或細微觀的缺陷進行單獨研究難以全面地掌握巖石承載失效特性,而目前統籌考慮巖石宏細微觀缺陷對其強度特性影響的研究較少。因此,本文綜合考慮巖石微元強度的韋伯分布與裂紋長度的冪律排布規律,采用數值模擬與理論推導方法分別構建損傷巖石試件的離散元計算模型與強度預測模型,并利用數值計算結果對理論模型的合理性進行驗證。

1 巖石微元強度韋伯分布與裂紋長度冪律排布理論

1.1 巖石微元強度韋伯分布理論

巖石可視為由眾多基元介質組成的有機整體,巖石破壞過程對應基元體損傷由量變到質變,基元介質的力學性質描述適用于統計的方法[7-9]。和眾多專家學者的研究一致,本文假設巖石微元破壞的概率隨微元強度F=f(σa)的分布密度P為[7]

P=[f(σa)]

(1)

式中,σa為試件的許用強度。

根據(1)式可得損傷度D與微元破壞的概率密度的關系為

(2)

根據(2)式可知,確定巖石的微元強度與破壞概率是建立巖石損傷模型抑或巖石強度預測模型的前提。目前,韋伯函數已被證明適合用于描述巖石微元強度或損傷的統計分布規律[7-12],對應的巖石微元破壞概率密度方程為

(3)

式中:n,F0均為韋伯分布參數;n為均質系數,用以反映微元體強度等力學參數分布的密集程度;F0代表尺度參數。

1.2 巖石裂紋長度冪律排布理論

已有研究表明[25,27-28],巖石中的裂紋長度服從冪律排布,對應的分布方程為

(4)

式中:l表示裂紋長度;x表示有關裂紋長度分布的自變量;另外,a為大于0的常數且a=zeta(1/b-1),而b為大于1的冪律指數。在實際應用中,需要預先給出預制裂紋的最大長度lmax、最小長度lmin和冪律指數b,從而確定特定的冪律排布形式。然后,根據(4)式的增減性可以建立(5)式所示關系

(5)

式中,xmax,xmin分別為x的上界和下界。另外,考慮到(5)式為超靜定方程,a,xmin及xmax均為未知參數,需要提前給一個參數賦值,故假定(6)式成立

a=1或xmin=1或xmax=1

(6)

根據(5)～(6)式可以確定裂紋長度的冪律排布形式。在此基礎上,通過給定預制裂紋數量n,具體的裂紋長度可以利用線性插值(理論)或概率分配(數值模擬)的方式獲得。

2 考慮微元強度韋伯分布與裂紋長度冪律排布的砂巖強度模擬研究

無論是巖石微元強度的韋伯分布抑或裂紋長度的冪律排布,對其研究都難以通過室內實驗開展,但理論模型的建立需有數據支撐?？紤]到韋伯分布及冪律排布在本質上都是數學運算[7,27],完全可以通過數值方法編程實現,故在此采用數值模擬方法分析微元強度韋伯分布與裂紋長度冪律排布對巖石強度的影響。

2.1 顆粒黏結模型

PFC2D(particle flow code 2D)作為離散元顆粒流軟件,被廣泛應用于巖土力學方面的研究[27]。該軟件優勢契合本文研究主題,能夠較好地模擬微元強度及裂紋長度等分布對介質模型強度的影響,故被本文研究選用。PFC2D軟件對巖土介質的表征主要采用顆粒及黏結兩部分,其中顆粒為剛體,通常用圓盤表示;而對于顆粒間的黏結,適用于巖石材料的模型主要有2種,分別為接觸黏結模型和平行黏結模型,兩模型對應的顆粒及黏結的力學行為如圖1所示。接觸黏結模型(見圖1a))相對簡單,僅能表征由顆粒接觸、碰撞產生的法向力和切向力;而平行黏結模型(見圖1b))既能反映顆粒間的接觸作用,亦能反映顆粒間的彎矩作用。通常認為,此2種模型在巖石內部均存在[27,29-30],故本文研究采用包含這2種模型的顆粒黏結模型(PBM)。

圖1 黏結模型及其微觀力學行為示意圖[27,29-30]

2.2 無初始損傷砂巖數值模型及細觀參數的確定

研究巖石微元強度韋伯分布及裂紋長度冪律排布對其強度的影響,理想狀態下需獲取無初始損傷試件的強度特性[14,27]。但在現實條件下,所有可獲得的巖石試件都不可避免地存在各類損傷且難以量化,專家學者一般選擇無明顯損傷的試件作為無初始損傷試件開展研究[14]。按照此種思路,筆者基于標準尺寸青砂巖的單軸壓縮室內實驗,在先前的研究中已經建立了相應的PFC2D數值模型并標定了模型細觀參數[27]。模型參數已在表1中列出,本文研究均基于該數值模型及標定的細觀參數。需要特別指出的是,所提青砂巖試件其單軸壓縮實驗強度約為52.22 MPa[27]。

表1 PFC2D介質細觀參數[27]

2.3 微元強度韋伯分布與裂紋長度冪律排布的實現及數值計算結果

2.3.1 微元強度韋伯分布的實現

對于表1中列出的眾多模型細觀參數,用以控制微元強度的主要為黏結黏聚力及黏結拉伸強度[27]。所以對于數值模型微元強度韋伯分布的實現,僅需利用PFC2D軟件自帶的FISH語言編程生成特定形式的韋伯分布隨機數,隨機數數量與黏結數量相同,然后將韋伯分布隨機數分別與黏結的黏聚力和拉伸強度同時結合(取乘積),生成新的平行黏結黏聚力及平行黏結拉伸強度數值。本文綜合考慮均質系數n對Weibull分布覆蓋區間范圍及區間內分布梯度的影響,設置其值分別為1,1.5,3,5,7。通過對模型顆粒間的黏結信息進行統計,得到的顆粒黏結強度σb模擬與理論分布之間的對比情況如圖2所示?？梢钥闯?理論分布與模擬分布高度一致,顆粒黏結強度均呈減速降低趨勢。并且,均質系數越小,顆粒黏結強度參數極大值越大;當均質系數為1,1.5,3,5,7時,對應的黏結強度參數極大值分別為401.79,172.87,74.38,53.08,45.93 MPa。另外,各均質系數下的曲線交匯點對應的黏結強度值為32 MPa,與表1中黏結強度值一致。

圖2 顆粒黏結強度模擬分布與理論分布對比

2.3.2 裂紋長度冪律排布的實現

在PFC2D軟件中,裂紋以構件的形式存在。長度服從冪律排布的裂紋是利用FISH語言編譯軟件中的DFN(discrete fracture network)模塊生成的,具體步驟為:①設置裂紋的長度范圍、長度分布形式、傾角范圍、位置分布形式;②給定裂紋的數量以及在試件中的分布范圍,由此生成裂紋;③賦值裂紋力學屬性參數。本文設置的裂紋參數與文獻[27]一致,裂紋長度范圍為0.5～10 mm,服從冪律排布;傾角范圍為0～180°;位置分布隨機,分布范圍與試件面積范圍一致[31]。預制裂紋力學參數賦值基于平滑節理模型(smooth-joint model),模型參數已在表2中列出[14,27]?？梢钥闯?裂紋存在的主要作用是破壞裂紋面兩側黏結的黏聚力和拉伸強度。

表2 平滑節理模型的力學參數[14,27]

溫度應力、風化、開采擾動、海水侵蝕以及地質作用程度的不同都會對巖石損傷產生較大影響,導致其內部裂紋信息存在較大差異?；诖朔N考量,本文建立的數值模型根據預制裂紋數量及裂紋分布冪律指數可以分為4類。具體而言,預制裂紋數量分別設置為50,100,200,400,而冪律指數b分別設置為1.1,2,3,6,該參數同樣是在考慮冪律排布的覆蓋區間范圍及區間內分布梯度變化的情況下確定的。當預制裂紋數量為50,預制裂紋長度lp在不同冪律指數下的理論分布與模擬分布對比如圖3所示?？梢钥闯?模擬與理論分布形式基本一致,均呈減速降低趨勢。并且,冪律指數越小,裂紋長度降低速率越平緩(見圖3中冪律指數為1.1對應的分布曲線),對應長裂紋占比越高;冪律指數越大,裂紋長度降低更接近于兩段式,先急速降低再平穩變化(見圖3中冪律指數為6對應的分布曲線),對應長裂紋占比越低。

2.3.3 數值計算結果

根據2.3.1節及2.3.2節可知,本文共計建立了5均質系數×4裂紋數量×4冪律指數=80個參數組合。為了盡可能降低試件強度離散性的影響,每個參數組合包含5個試件。需要指出的是,每個參數組合5個試件的所有參數信息都是一致的,但5個試件對應的微元強度韋伯分布模塊及DFN模塊調用次數分別為1～5。當微元強度韋伯分布模塊及DFN模塊調用次數一定時,試件內部的微元強度和裂紋分布信息是完全一致的,即試件的模擬計算是可重復的;而當調用次數變化時,微元強度和單個裂紋的實際存在位置及形式將再次發生隨機分布[27]。由上述介紹可知,本文針對損傷巖石共計建立并計算了400個數值模型試件,基本能夠實現對微元強度、預制裂紋信息及巖石強度間關系的統計研究[27]。

以預制裂紋數量np為x軸(第一維),均值系數n為y軸(第二維),冪律指數b為z軸(第三維),單軸抗壓強度σ1c為顏色軸(第四維),根據模擬結果構建得到的試件單軸強度四維散點σ1c(np,n,b),如圖4所示?？梢钥闯?各參數組合試件1～5的強度分布規律基本一致,抗壓強度最大值超過26 MPa,最小值接近1 MPa;另外,抗壓強度具有明晰的定向流動規律,其最小值均對應預制裂紋數量較多、均質系數較小且冪律指數較小處,抗壓強度最大值均對應預制裂紋數量較少、均質系數較大且冪律指數較大處。

圖4 400個試件抗壓強度模擬結果四維空間散點

3 考慮微元強度韋伯分布與裂紋長度冪律排布的砂巖強度預測模型

3.1 考慮微元強度韋伯分布的砂巖強度預測模型及驗證

(2)～(3)式分別給出了損傷變量與微元破壞概率密度的關系以及巖石微元強度與破壞概率的影響因素?？紤]到巖石最終破壞是損傷不斷增加的連續過程所引發的結果,故將(3)式代入(2)式進行積分,得到了損傷度與微元強度概率間的對應關系為

(7)

試件的有效應力σe與損傷程度的關系為

σe=σa(1-D)

(8)

式中,σa為無初始損傷試件的最大承載應力,數值上取為無初始損傷試件的單軸壓縮強度[27],根據2.2節知該值為52.22 MPa。

根據F的定義,其數學期望值E(F)、標準差SD均可表示為n,F0的函數,為

式中,符號Γ為伽馬函數,其在實數域上的定義為

(11)

(12)

考慮試件強度的離散性及(12)式的增減特性,將實際應力矩陣以區間范圍的形式表示為

(13)

實際應用中需要給出應力矩陣與均質系數之間的確切關系,如此便需確定方差系數?？紤]到試件(非)均質離散性隨均質系數的變化而變化,假定方差系數服從表達形式覆蓋廣泛的Logistic函數[3],由此可得實際應力矩陣與均質系數之間的關系

(14)

式中,ξ1,ξ2,ξ3和ξ4均為方差系數表達式的待定常數。

利用(14)式擬合模擬結果,如圖5所示,可以看出將方差系數假定為均質系數的Logistic函數能夠較好地描述有效應力隨均質系數的變化規律,相關系數R2超0.999,說明該假定是合理的。并且計算可得,在n的變化區間內,ξ1,ξ2,ξ3和ξ4的組合數值介于-1～1之間,滿足要求。

圖5 理論與模擬得到的巖石抗壓強度隨均質系數的變化規律

3.2 微元強度韋伯分布與裂紋長度冪律排布耦合影響下的砂巖強度預測模型及驗證

當僅考慮裂紋對試件強度的影響時,筆者已根據Mori-Tanaka方法建立了有效應力與裂紋密度之間的關系[27]

(15)

同樣利用Logistic函數對文獻[27]的數據進行回歸,可得η和d與冪律指數b的關系分別為[32-34]

相關系數R2均超0.999。根據(15)式中單位面積裂紋數量m與預制裂紋數量n的關系以及裂紋平均半長d與冪律指數b之間的關系,可得有效應力σe與預制裂紋數量n及冪律指數b之間的關系為

(16)

裂紋作為結構單元,其存在直接導致所在位置黏結的拉伸強度及黏結強度降低為0(見表2)。雖然裂紋對試件強度的作用在本質上也是通過影響(“殺死”)微元強度實現的,但其對黏結的破壞程度決定了受裂紋影響的黏結將可能不再受微元強度韋伯分布的影響,即裂紋與微元強度韋伯分布對黏結的影響互不干擾。據此,結合(14)及(16)式,疊加考慮微元強度韋伯分布及裂紋長度冪律排布對試件強度的影響,可建立12參數的砂巖強度預測模型為

(17)

結合圖5的擬合結果以及文獻[27]中的數據可得模型涉及的12個參數數值,在表3中列出。

表3 砂巖強度預測模型參數數值

理論模型建立后,需要對其合理性進行驗證。本文首先利用(17)式計算砂巖試件在不同預制裂紋數量、裂紋分布冪律指數及微元強度均值系數影響下的單軸抗壓強度值,然后構建得到了砂巖單軸抗壓強度理論值的四維空間散點,其與模擬結果平均值的對比如圖6所示?？梢灾庇^看出,抗壓強度理論分布(見圖6b))與模擬分布(見圖6a))基本一致,即抗壓強度最小值均對應預制裂紋數量較多、均質系數較小且冪律指數較小處;抗壓強度最大值均對應預制裂紋數量較少、均質系數較大且冪律指數較大處。

圖6 抗壓強度模擬結果平均值與理論計算值四維空間散點

上述是對理論結果與數值結果分布相似性的主觀描述,為了定量描述理論結果與模擬結果的吻合程度,依照數值模擬設置的參數組合,利用(17)式計算獲得了砂巖試件在相應預制裂紋數量、裂紋分布冪律指數及微元強度均值系數影響下的單軸抗壓強度,其與模擬結果平均值對比如圖7所示。圖中,橫坐標表示80個參數組合按照所得砂巖強度模擬均值自小至大依次排列時對應的序號?？梢钥闯?理論結果與模擬結果高度一致,并且計算得相關系數R2=0.991,說明利用所建強度預測模型能夠較好地描述微元強度韋伯分布與裂紋長度冪律排布對試件強度的影響。

圖7 理論計算值與模擬數據的對比

4 結 論

本文綜合考慮巖石微元強度的韋伯分布與裂紋長度的冪律排布規律,構建了損傷砂巖試件的PFC2D計算模型及強度預測理論模型,并利用數值計算結果對理論模型的合理性進行了驗證。主要結論為:

1) 通過編程同時實現了PFC2D計算模型中微元強度的韋伯分布與裂紋長度的冪律排布,分析了巖石中的宏、細微觀損傷與相應分布參數間的定量對應關系。

2) 數值建立了400個同時考慮宏、細微觀損傷的砂巖試件并對其進行了模擬加載,實現了對微元強度、預制裂紋信息及巖石強度間關系的統計分析。根據模擬結果構建了砂巖試件單軸抗壓強度四維空間散點σ1c(np,n,b),得到了多損傷參量影響下的試件抗壓強度流動規律。

3) 聯合Mori-Tanaka方法及巖石損傷概率分布理論,推導建立了12參數的巖石強度預測模型,經與數值模擬結果對比分析,該模型被證明能夠有效地描述細微觀微元強度及宏觀裂紋信息對巖石強度的影響。