含原生隱微裂隙巖石顆粒流模型構建及細觀參數標定方法研究*

李志強 劉國鋒 晏長根 董 凱 崔 潔

(①長安大學,公路學院,西安 710064,中國)(②海南大學,土木建筑工程學院,?？?570228,中國)

0 引 言

原生隱微裂隙是火成巖、沉積巖成巖過程的伴生物,在玄武巖、大理巖等硬、脆性巖石中十分常見。與工程尺度的斷層或斷裂、巖體統計尺度的節理相比,隱微裂隙尺度小,且分布隱蔽,無規律。以往建立巖石體力學模型或進行數值仿真時,較少考慮隱微裂隙的影響效應,實際上其對巖石的力學特性和破壞模式有著顯著影響,具體表現為:工程擾動下,裂隙尖端局部應力急劇增大,加劇次生裂隙的擴展演化過程,導致巖石的宏觀特征強度大幅降低(魏元龍等,2015; 唐禮忠等,2019; 張傳慶等,2019); 且巖石的破壞模式與裂隙的分布特征如位置、角度、長度等密切相關(Xie et al.,2011; 袁廣祥等,2019; 趙建軍等,2019; 張杰等,2021)。在巖體工程領域,含裂隙巖石力學特性、破裂機理及裂隙作用機制的理論和應用研究一直是熱點問題之一,目前常用的研究手段主要為室內試驗與數值仿真,室內力學試驗多以預制含裂隙試樣為研究對象(Liu et al.,2018; Jiang et al.,2019; Song et al.,2019; 陳偉等,2021),然而預制裂隙工藝復雜、制作困難且難以準確反映真實裂隙的物理力學特性。以顆粒流分析程序PFC為代表的離散元仿真作為一種簡單、便捷、行之有效的研究手段愈發受到研究者的青睞。

PFC能夠從細觀結構角度進行介質材料的力學特性的描述和受力變形的分析,可直觀顯現荷載作用下巖石內次生裂隙萌生、擴展直至貫通的全過程,在研究含裂隙巖石的力學特性和裂紋擴展特征方面具有明顯優勢(周喻等,2013; Itasca Consulting Group,Inc.2014)。陳鵬宇(2018)從接觸模型選取、細觀參數標定等方面詳盡介紹了PFC2D模擬裂隙巖石的主流方法。陳秀云(2015)基于PFC構建了單一中心閉合裂隙巖石顆粒流模型,開展多工況數值實驗,探究閉合裂隙特征參數,如寬度、長度、傾角等對裂隙巖石破壞特征的影響。黃丹等(2020),于輝等(2020)分別針對性地研究了閉合裂隙傾角與含裂隙大理巖、砂巖強度和變形特性的響應關系,其中于輝等(2020)同時考慮了圍壓對裂隙砂巖力學特性和破壞特征的影響。黃達等(2020)開展了巖石拉-壓數值模擬的試驗,實現了巖石高壓拉強度比的離散元模擬,并分析了其破壞機制。鄧清海等(2017)以直切槽式圓盤試樣為研究對象,基于PFC進行多組數值巴西試驗,探究裂隙傾角、長度對直切槽式圓盤試樣裂紋擴展規律的影響。唐世斌等(2021)重復開展多次數值巴西劈裂試驗,解析裂隙類型(張開、閉合、填充)、裂隙傾角等參數對次生裂隙演化過程的影響。李玉成等(2019),張晗等(2021)利用PFC開展含平行雙裂隙巖樣的數值壓縮試驗,揭示單軸加載條件下含平行雙裂隙巖樣的破壞規律及裂隙分布對該規律的影響。黃彥華等(2014)考慮巖橋傾角的差異構建多種含斷續節理巖樣顆粒流模型,并在不同圍壓下進行離散元仿真,分析圍壓和巖橋傾角對紅砂巖力學特性和破壞模式的影響。何棟梁等(2019)建立了多個具有不同錨固角度的錨固節理顆粒流模型,并在不同法向荷載下進行直剪,探究錨固節理的力學特性及節理面的破壞特征。方前程等(2014)通過類巖石材料模型試驗將研究對象尺度擴展為含斷續節理巖體,基于PFC仿真,研究不同圍壓下含不同傾角斷續節理巖體的力學響應。周喻等(2016)通過PFC構建斷續節理巖質邊坡顆粒流模型,分析工程尺度巖質邊坡的細觀破壞機制。上述有關裂隙巖石體的研究中,無論對象為小尺度的巖塊或是大尺度的巖體,PFC仿真均可作為一項重要的研究手段。

與其他能夠指定模型介質的宏觀本構關系并賦予物理力學參數的數值軟件不同,PFC采用顆粒局部接觸和演化特征來反映宏觀問題(李東東,2018),顆粒流模型的準確性及細觀參數的合理性直接決定著模擬成功與否,然而多數研究者利用其進行裂隙巖石數值分析時,主流思路仍是以預設規則裂隙顆粒流巖石模型為研究對象,考慮原生隱微裂隙真實分布形態構建真實裂隙巖石顆粒流模型的研究不多見,此外含原生裂隙顆粒流模型所涉及的多個接觸本構,其細觀參數的標定流程同樣尚不夠明確。鑒于此,本文以含原生隱微裂隙巖石為研究對象,對其真實顆粒流模型的構建方法及仿真所需細觀參數的標定流程展開系統研究,研究成果不僅極大增強了含原生隱微裂隙巖石顆粒流模型構建的準確性及細觀參數取值的透明性和合理性,同時可為內置多接觸本構顆粒流模型細觀參數的標定提供重要參考。

1 隱微裂隙巖石建模及參數標定方法

根據適用維度的不同,離散元仿真軟件PFC又細分為2D及3D版(Itasca Consulting Group,Inc.,2014)。相較于3D版,2D版計算效率高,硬件要求低,結果清晰、直觀便于分析,在模擬裂隙巖石過程中更為常用,本文建模及參數確定方法主要針對PFC2D,但并不局限于2D版,結果同樣可為3D版提供參考。

1.1 隱微裂隙巖石顆粒流模型構建

隱微裂隙巖石顆粒流模型的構建遵循“先構建完整顆粒模型,后導入數字化裂隙”的原則,首先建立與巖樣同尺寸的閉合“墻體”,作為模型邊界; 其后在“墻體”內部填充顆粒單元,運算平衡,完整顆粒流模型即構建完成。在此,需要明確顆粒最小粒徑、最大最小粒徑比、孔隙率、密度等顆粒幾何參數。眾多學者(Jensen et al.,1999; 周喻等,2011; 石崇等,2015; Chen,2016; 黃彥華等,2016; 陳鵬宇,2018)針對如何確定顆粒單元幾何參數開展了大量工作,參照其研究結果同時結合室內物理力學測試即可確定上述參數。

完整模型建立成功并驗證無誤后,后續工序便是將隱微裂隙準確地導入其中,如何更精確對巖樣內隱微裂隙的分布形態進行刻畫顯得尤為關鍵。朱珍德等(2007)綜合利用電鏡掃描、數字圖像處理及數理統計等分析手段對大理巖內原生微裂隙的分布規律進行研究。鐘志彬等(2017),汪權明等(2020)通過普通低倍放大鏡分別對流紋巖、白云巖側表面的微裂隙進行觀察和描繪,借助CAD實現裂隙幾何分布信息的量化。袁廣祥等(2019)憑借偏光顯微鏡觀測、拍照及數字化處理獲取花崗巖內原生裂隙的分布規律。Jing et al.(2021)基于CT掃描和數字圖像處理技術實現煤巖三維裂隙網絡模型的重構及裂隙結構的定量化提取。以上諸多方法中,CT掃描能夠在無損條件下獲得巖樣內隱微裂隙更為準確、全面的分布特征。利用該特性,提出圖1所示含隱微裂隙巖石顆粒流模型構建方法。該方法包括:試樣CT掃描,裂隙識別、切片導入CAD,坐標系(與完整模型相同)建立,裂隙臨摹,幾何坐標獲取,fish裂隙生成、導入完整模型等步序。相對于拍照、臨摹等方式所建立的含裂隙巖石離散元模型,上述方法所得模型更為準確。

圖1 隱微裂隙巖石顆粒流模型構建方法

1.2 巖石及隱微裂隙細觀參數標定

PFC軟件內置多個接觸模型用以描述介質顆粒間接觸的力學行為,多位學者(Lee et al.,2011; Mas et al.,2011; 胡波等,2019)研究表明:平行黏結模型、平直節理模型可以很好地模擬完整巖石的力學特性,節理裂隙用光滑節理模型模擬較為合適,影響上述3種接觸力學行為的主要細觀參數分別如表1、表2和表3所示。

表1 平行黏結模型主要細觀參數

表2 平直節理模型主要細觀參數

表3 光滑節理模型主要細觀參數

表中各細觀參數并沒有實際的物理意義(Mehranpour et al.,2017),標定其數值時常采用大量試算建立模擬結果(宏觀力學參數)與所用細觀參數的關系,依據試驗所得宏觀參數反推細觀參數的方式。通過歸納、總結前人研究,提出了圖2所示完整顆粒流模型細觀參數的標定方法,平行黏結模型及平直節理模型均可用于描述完整模型顆粒單元接觸間的力學行為,對于不同類型裂隙巖石,“平行黏結模型+光滑節理模型”和“平直節理模型+光滑節理模型”兩種組合在描述裂隙與巖石的交互影響時各有千秋,可通過手動試算確定最佳接觸模型,也可參考他人研究進行選擇。整體上,完整模型細觀參數的標定流程可概括為:模型選擇、參數組合、批量模擬、敏感性分析、參數關系構建、參數反算和結果微調。

圖2 完整顆粒流模型參數確定流程圖

光滑節理模型模擬節理裂隙時,節理經過的顆粒,其原本的接觸本構會被光滑節理接觸替換,采用該方式模擬節理裂隙遠比移除節理位置的接觸或降低黏結強度的方式更為合理。結合前人研究,確定裂隙與顆粒單元間接觸的細觀參數時,首先需在考慮節理裂隙影響的條件下,分析光滑節理模型各細觀參數與模擬所得宏觀參數的響應關系。為此應先預建含不同傾角裂隙巖石的顆粒流模型,完整模型及參數上節已定保持不變,指定裂隙接觸為光滑節理模型,設計多種參數組合方案,批量模擬構建結果樣本,定性分析模擬結果與細觀參數間的映射關系,隨后根據CT掃描結果,建立能夠反映裂隙真實分布形態的隱微裂隙巖石顆粒流模型,結合試驗所得裂隙巖石力學參數及宏觀結果與細觀參數的關系,反復調算,即可確定裂隙顆粒流模型所需細觀參數,詳細流程見圖3。

圖3 隱微裂隙顆粒流模型參數確定流程圖

2 方法應用與驗證

隱晶質玄武巖是我國西南山區白鶴灘水電站地下廠房施工期間最廣泛揭露的圍巖之一,新鮮、完整、質地堅硬、巖性極好,見圖4。然而廠房施工期間,隱晶質玄武巖洞壁卻出現了大量片幫剝落、卸荷開裂等破壞現象,如圖5所示。室內實驗結果顯示,隱晶質玄武巖的單軸強度呈現出極強的離散性,最高可達351MPa,最低僅為107MPa(江權等,2017)。劉國鋒等(2016)發現其內部隨機發育有大量分布隱蔽、無規律,細小閉合的隱裂隙。多位學者(江權等,2017; 張傳慶等,2019)均認為正是由于原生隱裂隙的存在才導致部分隱晶質玄武巖試樣強度極低,進而導致宏觀強度表現出較高的離散性。

圖4 隱裂隙發育的隱晶質玄武巖

圖5 洞室施工期間隱晶質玄武巖典型破壞(片幫剝落)

以原生隱微裂隙發育的隱晶質玄武巖為例進行裂隙巖石顆粒流模型構建及參數標定方法的應用和驗證。實驗所用試樣為國際巖石力學學會ISRM建議的標準圓柱樣(直徑50mm,高100mm),利用CT掃描測定內部隱裂隙分布形態,判定試樣為完整試樣或含隱裂隙試樣,其后在RMT150-C剛性試驗機上進行力學特性實驗,采用位移控制模式,以0.002mm·s-1的速度進行加載,直至試樣破壞。對于完整試樣而言,當加載強度即將達到峰值強度之前,試樣會斷續發出“噼啪”聲響,一旦加載達到試樣破壞所需的峰值強度,試樣會隨即產生脆性爆裂破壞,伴隨著巨大聲響,炸裂成塊狀或片狀并彈射出去。而對于隱含裂隙的試樣而言,其破壞過程與完整試樣則明顯不同,加載過程較為平靜,峰值破壞后,試樣仍可保持相對較好的完整性。典型試樣宏觀力學參數及內部隱裂隙數量如表4所示。

表4 隱晶質玄武巖典型試樣宏觀力學參數及隱裂隙數量

2.1 完整模型構建及參數標定

通過指定“墻體”的端點和法線方向生成長100mm,寬50mm的矩形閉合“墻體”,采用“半徑擴大法”在“墻體”內填充顆粒單元,調整內應力,刪除懸浮顆粒,完整模型即構建完成。顆粒單元幾何參數取值見表5。

表5 顆粒單元幾何參數

顆粒分布滿足預定要求后,即可按圖2所示流程進行參數標定。首先指定顆粒間的接觸本構,這里以平行黏結模型為例。關于平行黏結模型,眾多學者(阿比爾的等,2018; 胡訓健等,2021; Wu et al.,2021)已經開展了大量研究,其細觀參數與宏觀計算結果間的關系較為明晰,主要結論有:宏觀彈性模量(E)對平行黏結有效模量(pb＿emod)最為敏感; 法向剛度(pb＿kn)和切向剛度(pb＿ks)的比值(pb＿kratio)對泊松比(ν)的影響最大; 而法向和切向黏結強度的比值(pb-ten/pb＿shear)則主要影響試樣的破壞模式。保證結果可靠性的前提下,為簡化平行黏結模型細觀參數的標定過程,預先做出兩個基本假設,即顆粒接觸模量等于平行黏結有效模量,顆粒法向和切向的剛度比等于平行黏結剛度比,同時固定對計算結果影響不大的參數值,顆粒間的摩擦系數(fric)取0.5,黏結半徑系數(pb＿rad)為1.0。

2.1.1 平行黏結有效模量(pb＿emod)的確定

分別取黏結有效模量為5GPa、10GPa、15GPa、20GPa、25GPa、30GPa和35GPa,固定其他參數的取值不變,重復進行7次單軸壓縮模擬,模擬所得彈性模量與黏結有效模量的關系如圖6所示。

圖6 pb＿emod與宏觀彈性模量的關系

彈性模量隨黏結有效模量取值的增加而增大,根據表4,完整巖石的彈性模量為53.45GPa,初步確定平行黏結有效模量等于21.87GPa。

2.1.2 平行黏結剛度比(pb＿kratio)的確定

固定顆粒接觸模量和黏結有效模量為21.87GPa,設置黏結剛度比為1.0、1.5、2.0、2.5、3.0、3.5和4.0,其他參數保持不變,7次重復試驗所得泊松比隨黏結剛度比的變化見圖7。

圖7 pb＿kratio與泊松比的關系

隨著黏結剛度比的增加,泊松比也在不斷增大,但是增大的速率在逐漸變緩。將表4所示完整試樣的泊松比(0.23)帶入圖中對數函數擬合公式,所得黏結剛度比為2.06。

2.1.3 法向黏結強度(pb＿ten)和切向黏結強度(pb＿shear)的確定

黏結強度比(pb＿ten/pb＿shear)影響模型最終的破壞模式,確定黏結強度前,應先根據試樣的破壞模式,確定黏結強度比。人為設定切向黏結強度為10MPa,取黏結強度比為0.2、0.6、1.0、1.4、1.8和2.2,其他參數不變,重復進行6次單軸壓縮數值模擬,典型破壞形式見圖8。

圖8 不同黏結強度比pb＿ten/pb＿shear對應的破壞形式(藍色為張拉裂隙,紅色為剪切裂隙)

圖8a為黏結強度比取0.2的模擬結果,典型的脆性劈裂破壞,破壞產生的裂隙主要為張拉裂隙,未出現剪切裂隙。隨著黏結強度比的增大,剪切破壞產生的裂隙逐漸增多,試樣的破壞形式也在由劈裂破壞向剪切破壞轉換。玄武巖屬于硬脆性巖石,完整試樣破壞形式以脆性劈裂破壞為主,反復試算后,取黏結強度比為0.29,則對應的法向黏結強度為2.9MPa,設置黏結強度變化系數k為0.5、1、2、4、6、8和10,得到7組比值均為0.29的法向黏結強度和切向黏結強度,展開模擬,峰值強度與黏結強度變化系數的關系如圖9所示。

圖9 k與峰值強度的關系

峰值強度σ隨著黏結強度變換系數k近似線性變化,依據表4所示完整試樣的峰值強度(331.81MPa),黏結強度變化系數計算結果為13.69,則法向黏結強度和切向黏結強度分別為39.7MPa和136.9MPa。

以上述參數進行數值單軸壓縮試驗,所得力學參數與力學試驗結果存在較大誤差,見表6。

表6 初始參數模擬結果

上述問題是由于模型各細觀參數對宏觀計算結果存在交叉影響(石崇等,2018),即同一個細觀參數影響多個宏觀參數的計算結果,同一個宏觀參數的計算結果受到多個細觀參數的影響,需要結合已建立的宏、細觀參數間的關系按照式(1)對已確定的初始細觀參數進行微調。

(1)

式中:plt表示室內實驗所得宏觀參數;psm為數值模擬的結果;a代表宏、細觀參數函數關系的比例系數;pbi指第一次確定的初始細觀參數;pbc表示調整后的細觀參數。

杜朗的徒弟陳洋接替了他在羅恬身邊的位置——不光是要替代杜朗在富人圈子里表演魔術，還要替代杜朗做羅恬的情人。

調整后的參數如表7所示,該參數模擬結果見表8,此時模擬與實驗結果的誤差皆已控制在5%以內。數值模擬和室內試驗所得應力-應變曲線及破壞形態的對比見圖10,兩者吻合度高,模擬較為成功。

表7 PFC數值仿真所用最終參數

表8 最終參數模擬結果

圖10 數值模擬和室內試驗結果的對比

2.2 原生隱微裂隙巖石模型構建及參數確定

裂隙試樣中的裂隙“導入”過程如圖11所示,該裂隙試樣的力學參數及裂隙特征見表4。細觀參數按圖3所示流程進行標定,預設規則裂隙巖石顆粒流模型,其中預設裂隙長為1.5cm,傾角φf分別為0°、15°、30°、45°、60°、75°和90°,且預設裂隙位于模型中心位置。部分模型如圖11所示,指定裂隙與顆粒單元間接觸的數值本構為光滑節理模型,細觀參數組合通過折減平行黏結模型對應細觀參數及人為給定的方式生成(具體方案見表9),控制變量重復進行119組數值單軸壓縮實驗,模擬所得力學參數、設定的節理參數、裂隙傾角三者間的關系見圖12。

表9 光滑節理模型參數組合方案

圖11 裂隙導入流程

圖12 含預設規則裂隙巖石顆粒流模型

對于裂隙傾角為90°的模型試樣而言,改變法向、切向剛度和摩擦系數,模擬所得力學參數均未發生明顯改變,表明其對圖示3個細觀參數并不敏感。其他裂隙試樣的力學參數則在細觀參數改變時均發生相應變化。

圖13a、圖13b、圖13c所示分別為節理法向剛度對裂隙顆粒流模型峰值強度、彈性模量和泊松比的影響。當裂隙傾角為0°～30°時,隨著節理法向剛度的減小,裂隙模型峰值強度明顯逐漸降低。傾角為45°時,其峰值強度變化趨勢同上,但變化值大幅減小。傾角為60°時,隨法向剛度減小,峰值強度先增大后減小往復變化,傾角為75°時,模型試樣的峰值強度恢復隨法向剛度折減而減小的變化規律,但減小量較小。對于0°和15°傾角的裂隙模型,法向剛度依次折減,彈性模量則不斷降低。傾角為30°時,彈性模量則先增加后降低。傾角為45°的試樣,彈性模量隨法向剛度減小而增加,且只有法向剛度減小到一定程度后,彈性模量才開始有明顯增加。對于60°傾角裂隙模型,其彈性模量隨法向剛度的變化規律恰好與45°模型相反,隨著法向剛度的減小,彈性模量先減小后變大。75°傾角模型,彈性模量并未隨法向剛度改變產生明顯變化。對于0°～30°傾角的裂隙模型,節理法向剛度減小,其泊松比對應增加,且只有法向剛度減小到一定程度后,泊松比開始有明顯增加。裂隙傾角為45°和60°時,隨著法向剛度的減小,泊松比分別呈現先減小后增大與先增大后減小的相反規律。75°傾角試樣,其泊松比不隨法向剛度的變化而變化,這一點與彈性模量相似。

圖13 模擬所得力學參數、設定的節理參數、預設裂隙傾角三者間的關系

節理切向剛度與裂隙模型峰值強度、彈性模量和泊松比的關系分別見圖13d、圖13e、圖13f。當裂隙傾角為0°、60°和75°時,改變節理切向剛度,裂隙巖石的力學參數并未產生明顯變化。對于15°和30°傾角的裂隙模型,降低節理切向剛度,峰值強度顯著降低。含45°傾角裂隙模型,峰值強度同樣隨切向剛度降低,但降幅較小。對于傾角為15°的裂隙模型,改變節理切向剛度,模型彈性模量發生小范圍變化。對于45°裂隙模型,其彈性模量隨切向剛度減小呈“臺階形”增長,即切向剛度折減到一定程度后,彈性模量才會發生明顯變化,泊松比變化規律類似。含30°傾角裂隙模型,當節理切向剛度折減1/50后,其泊松比顯著增加。

圖11所示為典型的含裂隙試樣,其裂隙傾角分別為30°、45°和60°。對于30°傾角裂隙試樣,降低節理法向、切向剛度,減小摩擦系數,可以降低其峰值強度。需要降低彈性模量時,則可通過大幅降低節理法向剛度和摩擦系數來實現。欲使泊松比增大時,則可采用減小節理法向、切向剛度及摩擦系數的方式。含45°傾角裂隙試樣,增大節理法向剛度和摩擦系數,可以使試樣的峰值強度增大; 減小節理法向、切向剛度可以增大試樣彈性模量; 增大節理切向剛度和摩擦系數,可以減小裂隙試樣的泊松比。對于含60°傾角裂隙試樣,通過折減節理法向剛度的方式可以減小試樣的峰值強度; 改變節理法向剛度和摩擦系數可以使試樣的宏觀彈性模量發生變化,但對應規律需要具體分析; 增加節理法向剛度可以增大試樣的泊松比。

結合上述定性分析,經反復試算,適用于表4所示典型含裂隙試樣的光滑節理接觸細觀參數如表10所示?？紤]到隱裂隙為閉合裂隙,特增加了法向黏結強度和切向黏結強度兩個參數,表征閉合裂隙的黏結作用?；诒?0中參數及圖11所示真實的裂隙巖石顆粒流模型進行數值單軸壓縮實驗,結果見表11,數值模擬與室內試驗獲得的力學參數十分接近。圖14所示為模擬和試驗所得應力-應變曲線的對比,兩者十分吻合,可見,該方法所建模型及所確定的細觀參數較為合理。

表10 仿真所用細觀參數

表11 裂隙巖石顆粒流模型模擬結果

圖14 模擬與試驗結果的對比

3 主要結論

本文結合前人研究及CT掃描、圖像數字化技術建立了原生隱微裂隙巖石顆粒流模型構建和細觀參數標定的方法,并以隱晶質玄武巖為例,對方法進行應用和驗證。主要結論如下:

(1)CT掃描配合圖像數字化技術可實現試樣內裂隙形態的精確刻畫,基于此提出了能夠反映裂隙真實分布形態的含隱微裂隙巖石顆粒流模型的精細化構建方法。

(2)平行黏結模型細觀參數與宏觀模擬結果(力學參數)的關系較為簡明,但考慮各參數的交互影響,初始參數仍需微調后方可使用。光滑節理模型細觀參數與宏觀模擬結果的關系相對較為復雜,確定其參數時,應先根據已建模型定性分析不同傾角下宏觀結果與細觀參數間的響應。

(3)“平行黏結+光滑節理”的數值本構組合能夠很好地描述含隱微裂隙隱晶質玄武巖顆粒流模型接觸間的力學性質。

(4)通過典型案例,驗證了本文所提出的含隱微裂隙巖石模型構建和參數標定方法,模擬與試驗結果的對比分析證實了本文所提方法的合理性和可行性。