夏海城,鄔愛清,盧 波,徐棟棟
(長江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室,武漢 430010)
巖石材料是由礦物顆粒、膠結(jié)材料和空隙等組成的具有一定缺陷的不均勻物質(zhì),非均質(zhì)性是巖石材料的本質(zhì)特征。這種細(xì)觀結(jié)構(gòu)的不均勻性導(dǎo)致巖石材料的彈性模量、強(qiáng)度、泊松比等宏觀力學(xué)參量具有一定的隨機(jī)性,同時(shí)也直接影響受載巖石的破裂演化過程。因此研究非均質(zhì)性對(duì)巖石材料力學(xué)性質(zhì)的影響機(jī)制具有重要意義。
國內(nèi)外學(xué)者對(duì)巖石非均質(zhì)問題展開了廣泛的研究,并取得諸多有價(jià)值的成果。巖石材料物理參數(shù)分布的不均勻性是其非均質(zhì)性的重要方面,胡小榮等[1]應(yīng)用地質(zhì)統(tǒng)計(jì)學(xué)方法對(duì)有限單元體力學(xué)參數(shù)的賦值進(jìn)行了研究;陳永強(qiáng)[2]采用不同的統(tǒng)計(jì)分布模擬材料的非均勻性,得到了不同特性的材料破壞過程應(yīng)力-應(yīng)變曲線;康健等[3]研究了非均質(zhì)隨機(jī)分布的熱膨脹系數(shù)和熱傳導(dǎo)系數(shù)對(duì)巖石熱破裂的影響;馮增朝等[4-5]研究了巖石細(xì)胞單元特性及其非均質(zhì)分布對(duì)巖石破壞全過程曲線的影響,巖石亞細(xì)觀組構(gòu)的非均質(zhì)性及其分布規(guī)律是決定巖石單軸壓縮全過程曲線的重要因素;張征等[6]分析了巖土參數(shù)不確定性的主要原因,并研究了巖土參數(shù)空間變異性分析的原理、方法和步驟。唐春安等[7-8]、Kaiser等[9]用Weibull分布函數(shù)模擬巖石材料參數(shù)的非均質(zhì)性,研究了二維巖石模型脆性破壞的特點(diǎn),提出“RFPA(Rock Failure Process Analysis)方法”。梁正召[10]開發(fā)了RFPA3D程序,模擬了巖石的破壞過程。尤名慶等[11]討論了巖石非均質(zhì)性和強(qiáng)度尺寸效應(yīng)之間的關(guān)系;羅榮等[12-13]研究了基于材料組分的非均質(zhì)巖土體數(shù)值建模方法及非均質(zhì)度對(duì)巖石力學(xué)特性的影響;柯長仁等[14]將非均質(zhì)性引入虛內(nèi)鍵模型,模擬了巖石式樣動(dòng)態(tài)破壞過程;Cil和Buscarnera[15]利用Weibull 分布的最弱環(huán)節(jié)理論,分析了不同粒徑尺寸對(duì)單軸最大抗壓強(qiáng)度的影響??偟膩碚f,在巖石非均質(zhì)的模擬方面的研究都采用基于某種統(tǒng)計(jì)分布對(duì)每個(gè)巖石單元賦予各不相同的物理參數(shù),使細(xì)觀單元的力學(xué)性質(zhì)各不相同,以模擬巖石宏觀上的非均質(zhì)性。
非均質(zhì)巖石材料力學(xué)特性的影響因素不僅在于巖石細(xì)觀材料參數(shù)之間的差異,這種差異形成的力學(xué)薄弱環(huán)節(jié)在空間的分布特征直接控制著受載巖石的破裂演化過程,從而影響巖石的強(qiáng)度和應(yīng)力-應(yīng)變?nèi)^程曲線的特征。上述研究主要關(guān)注因巖石成分不同導(dǎo)致的細(xì)觀材料參數(shù)的差異對(duì)巖石宏觀力學(xué)性狀的影響。事實(shí)上,因成巖機(jī)制和成巖過程的差異,巖石的成分、結(jié)構(gòu)和構(gòu)造在空間上呈現(xiàn)不均勻分布特征;巖石的彈性模型、強(qiáng)度等宏觀物理力學(xué)參量的離散性則是這種細(xì)觀不均勻性的綜合體現(xiàn)。既有研究成果未充分關(guān)注細(xì)觀不均勻性的空間分布特征對(duì)巖石宏觀力學(xué)特性的影響。本文采用具有伺服控制機(jī)制的非連續(xù)變形分析(Discontinuous Deformation Analysis,DDA)方法,根據(jù)巖石顆粒呈不規(guī)則多邊形的特點(diǎn),采用Voronoi結(jié)構(gòu)來模擬巖石礦物顆粒,構(gòu)建用于DDA分析的巖石數(shù)值試樣,通過對(duì)不同的Voronoi塊體單元賦予不同的物理力學(xué)參數(shù)的方式來實(shí)現(xiàn)對(duì)非均質(zhì)性的模擬。通過開展巖石單軸壓縮數(shù)值實(shí)驗(yàn),在驗(yàn)證非均質(zhì)度對(duì)巖石強(qiáng)度影響的基礎(chǔ)上,探究了巖石單元空間分布的影響作用,補(bǔ)充了主要結(jié)論。
DDA方法是石根華博士于20世紀(jì)80年代提出的平行于有限單元法的數(shù)值模擬方法,其單元或塊體形狀可以是任何凸?fàn)钚位虬紶钚蔚亩噙呅?。[16]相比有限單元法只能用標(biāo)準(zhǔn)形狀單元具有一定的優(yōu)越性,在DDA方法中,任意塊體的位移Di可用6個(gè)位移變量表示,即
Di=(u0v0r0εxεyγxy) ;
(1)
式中:(u0v0)是塊體內(nèi)特殊點(diǎn)(x0y0)的剛體位移;角r0是塊體繞轉(zhuǎn)動(dòng)中心(u0v0)的轉(zhuǎn)動(dòng)角,用弧度給出;(εxεyγxy)是該塊體的法向和切向應(yīng)變。
塊體中任意點(diǎn)(x,y)處的位移可由變形變量[Di]表示為
(2)
[Ti]=
(3)
式中[Ti]為位移轉(zhuǎn)換矩陣。
各連接塊體依靠接觸和位移約束形成塊體系統(tǒng)。n個(gè)塊體的聯(lián)立方程為
式中Kij是6×6子矩陣,與材料特性有關(guān);[Di]和[Fi]是6×1子矩陣,此處Di代表塊體i的變形變量,F(xiàn)i為塊體i上的6個(gè)變形變量的荷載。對(duì)于平衡方程(式(4)),由力和應(yīng)力產(chǎn)生的總勢(shì)能Π最小化來推導(dǎo),對(duì)Π求導(dǎo),得出的總勢(shì)能最小。子矩陣Kij為
(5)
子矩陣Fi是
(6)
式中d表示塊體的位移變量。
在對(duì)式(2)和式(4)求解后,可計(jì)算出塊體變形后的位移與應(yīng)變。由于DDA法中塊體無嵌入和張拉的條件,需在接觸位置加上或去掉剛性彈簧,修正總體平衡方程,并再次求解直至所有邊界滿足無嵌入和張拉,迭代完成。
楊懋偲[17]在DDA源代碼的基礎(chǔ)上進(jìn)行了二次開發(fā),初步實(shí)現(xiàn)了對(duì)電液伺服控制機(jī)制的模擬,得到了符合室內(nèi)常規(guī)實(shí)驗(yàn)結(jié)果的應(yīng)力-應(yīng)變?nèi)^程曲線。
本文將使用具有伺服控制機(jī)制的DDA程序開展非均質(zhì)度對(duì)巖石材料宏觀力學(xué)性質(zhì)的影響機(jī)制研究。
巖石材料物理參數(shù)的不均勻分布是引起巖石非線性變形的重要因素,在細(xì)觀數(shù)值實(shí)驗(yàn)中,將巖石試樣劃分為多個(gè)Voronoi單元,每個(gè)Voronoi單元具有不同的物理力學(xué)參數(shù),這種細(xì)觀單元物理力學(xué)參數(shù)的不均勻分布可以較好地模擬出巖石非均質(zhì)性的特點(diǎn)。本文采用正態(tài)分布函數(shù)模擬礦物顆粒物理參數(shù)分布的隨機(jī)性,并定義了量化巖石非均質(zhì)度大小的方法。通過調(diào)整巖石試樣非均質(zhì)度和期望的大小,比較分析了巖石力學(xué)性能的變化規(guī)律,并探究巖石力學(xué)性能的影響因素。
構(gòu)成巖石的礦物顆粒物理力學(xué)參數(shù)盡管呈現(xiàn)出強(qiáng)烈的隨機(jī)性特點(diǎn),但仍然滿足一定的統(tǒng)計(jì)學(xué)分布規(guī)律。根據(jù)以往的研究成果,采用不同的參數(shù)分布模型開展數(shù)值實(shí)驗(yàn)得到的破壞過程曲線也不一樣。目前,常用于模擬巖石非均質(zhì)性的隨機(jī)函數(shù)模型有Weibull分布模型和正態(tài)分布模型,其中正態(tài)分布模型得到的曲線更加符合巖石脆性破壞的特點(diǎn)[2]。本文采用正態(tài)分布函數(shù)對(duì)巖石非均質(zhì)進(jìn)行表征。正態(tài)分布函數(shù)的概率密度函數(shù)如式(7),其中μ為期望,反映了非均質(zhì)巖石材料的綜合性能;σ為方差,反映了巖石材料彈性模量的偏離值。
(7)
將服從正態(tài)分布的隨機(jī)數(shù)賦給巖石單元的彈性模量,以模擬巖石材料的非均質(zhì)性。正態(tài)分布函數(shù)中,標(biāo)準(zhǔn)差反映的巖石彈性模量偏離值是度量巖石非均質(zhì)性的重要參數(shù)。為準(zhǔn)確定義巖石的非均質(zhì)度,有效量化非均質(zhì)度的大小,本文將巖石樣本中彈性模量的標(biāo)準(zhǔn)差與期望之比Cv定義為變異系數(shù)來度量巖石非均質(zhì)度的大小,即
(8)
式中:變異系數(shù)Cv無量綱,反映了巖石單元隨機(jī)物理參數(shù)的偏離度,用來度量巖石的非均質(zhì)性;Ed為巖石單元彈性模量的標(biāo)準(zhǔn)差;Em為巖石單元彈性模量的平均值。
生成服從正態(tài)分布的隨機(jī)數(shù)組時(shí),通過調(diào)整標(biāo)準(zhǔn)差和期望的大小來控制變異系數(shù)。顯然在Cv=0時(shí),巖石為完全均質(zhì)材料,各細(xì)胞單元具有完全相同的物理參數(shù)。隨著Cv增大,巖石材料彈性模量的非均質(zhì)度增大,單元彈性模量離散度增強(qiáng),巖石均勻性變差。Cv越小,巖石單元的物理參數(shù)越趨向于期望,巖石整體的均勻性越好。
為聚焦巖石材料非均質(zhì)度對(duì)其力學(xué)性能的影響,排除非均勻度以外單元分布結(jié)構(gòu)的干擾,本文構(gòu)建2組不同彈性模量期望的非均質(zhì)巖石數(shù)值模型(如表1),使用正態(tài)分布函數(shù)產(chǎn)生的隨機(jī)數(shù)組對(duì)巖石單元彈性模量進(jìn)行賦值,每種變異系數(shù)生成6個(gè)不同的隨機(jī)數(shù)組(也即形成6個(gè)數(shù)值試樣),滿足相同的期望和標(biāo)準(zhǔn)差,并將6個(gè)試樣進(jìn)行相同條件的數(shù)值實(shí)驗(yàn)。不同變異系數(shù)的巖石樣本單元彈性模量統(tǒng)計(jì)數(shù)目直方圖如圖1所示。
表1 正態(tài)分布模型參數(shù)及非均質(zhì)度Table 1 Parameters of normal distribution modeland degree of heterogeneity
圖1 Voronoi單元統(tǒng)計(jì)數(shù)目直方圖Fig.1 Histogram of statistical number of Voronoielement of elastic modulus
圖2 非均質(zhì)巖石材料模型Fig.2 Heterogeneous rockmaterial model
圖2為本文所采用尺寸為50 cm×100 cm的平面數(shù)值模型,為真實(shí)模擬巖石細(xì)胞單元的形狀結(jié)構(gòu),采用Voronoi多邊形細(xì)觀單元形成的塊體系統(tǒng)進(jìn)行模擬,每個(gè)Voronoi單元的面積約為5 cm2,共計(jì)975個(gè)單元。
本文中數(shù)值實(shí)驗(yàn)為單軸壓縮實(shí)驗(yàn),對(duì)剛性加載板采用施加固定位移的方法對(duì)巖石模型進(jìn)行壓縮,每個(gè)時(shí)步的加載位移為5×10-8。在單軸壓縮實(shí)驗(yàn)中巖石的非均質(zhì)性主要是由巖石Voronoi單元之間物理參數(shù)不均勻引起的,為排除其他物理參數(shù)干擾的影響,本文所有巖石Voronoi單元泊松比為0.3,Voronoi單元之間的界面強(qiáng)度設(shè)置為統(tǒng)一值。
本文中數(shù)值加載模型主要可分為剛性加載板、非均質(zhì)巖石模型和剛性支架3部分,使用混合高階DDA數(shù)值實(shí)驗(yàn)平臺(tái),通過對(duì)剛性加載板施加豎向位移進(jìn)行伺服控制單軸壓縮數(shù)值實(shí)驗(yàn)。
圖3 全過程應(yīng)力-應(yīng)變曲線Fig.3 Complete stress-strain curve
圖3為巖石樣本應(yīng)力-應(yīng)變?nèi)^程關(guān)系曲線,巖石樣本軸向應(yīng)變通過加載板的位移計(jì)算,加載板在軸向的位移量在數(shù)值上等于巖石試樣的變形量。巖石塊體軸向應(yīng)力隨著加載板的不斷加載而增大,并在到達(dá)極限荷載后承載力迅速下降,下降到某一值后保留一定的殘余強(qiáng)度,符合真實(shí)的脆性巖石的承載規(guī)律。顯然巖石變形破壞曲線中應(yīng)力最大的點(diǎn)也是巖石極限承載力的點(diǎn),加載板施加的荷載力最大時(shí),巖石應(yīng)力達(dá)到最大。巖石的宏觀彈性模量為達(dá)到峰值應(yīng)力50%時(shí),應(yīng)力與對(duì)應(yīng)應(yīng)變的切線模量。
圖4 均質(zhì)巖石應(yīng)力-應(yīng)變?nèi)^程曲線Fig.4 Complete stress-strain curve of homo-geneous rock
巖石樣本為完全均質(zhì)時(shí),應(yīng)力-應(yīng)變?nèi)^程曲線如圖4所示,此時(shí)所有Voronoi單元的彈性模量均為30 GPa。以此為參照,研究材料非均質(zhì)度對(duì)巖石宏觀彈性模量和極限承載力的影響,對(duì)每組變異系數(shù)的巖石試樣均開展數(shù)值模擬單軸壓縮實(shí)驗(yàn),取巖石樣本全過程破壞曲線較為平均的一組對(duì)比。
圖5為彈性模量期望為30、50 GPa時(shí)不同變異系數(shù)下的應(yīng)力-應(yīng)變曲線。由圖5可知,隨著變異系數(shù)的增大,巖石的極限承載力逐漸下降,峰后曲線具有明顯的差異。對(duì)于巖石的宏觀彈性模量,初始加載階段不同變異系數(shù)的巖石樣本宏觀彈性模量幾乎相同,隨著加載的不斷進(jìn)行,變異系數(shù)較大的巖石宏觀彈性模量先出現(xiàn)下降,總體來看非均質(zhì)度對(duì)巖石樣本宏觀彈性模量的影響是較弱的。
圖5 不同非均質(zhì)度模型應(yīng)力-應(yīng)變曲線Fig.5 Stress-strain curves of rock models ofdifferent degrees of heterogeneity
為了有效研究非均質(zhì)度對(duì)巖石極限承載力的影響,消除相同變異系數(shù)下不同單元分布對(duì)變形曲線的干擾,聚焦巖石極限應(yīng)力隨非均質(zhì)度的變化趨勢(shì),對(duì)每個(gè)變異系數(shù)下6組彈性模量單元的數(shù)值實(shí)驗(yàn)結(jié)果取平均值,得到平均極限應(yīng)力,并做出變異系數(shù)-極限應(yīng)力的折線圖,結(jié)果如圖6所示。可以看出,服從正態(tài)分布的巖石單元彈性模量,期望值分別為30 GPa和50 GPa時(shí),巖石極限承載力均隨著材料變異系數(shù)的增大而減小,兩者具有較好的線性關(guān)系。
圖6 變異系數(shù)與巖石極限應(yīng)力的關(guān)系Fig.6 Relation between coefficient of variation andultimate stress of rock
值得注意的是,滿足相同變異系數(shù)和彈性模量期望的不同巖石樣本的應(yīng)力-應(yīng)變曲線并不完全相同(見圖7),在峰值和峰后階段具有較大差異。在給巖石單元的彈性模量賦值時(shí),服從某種正態(tài)分布的單元集合具有多種分布的可能性,這導(dǎo)致了在相同變異系數(shù)的巖石樣本中,每個(gè)樣本內(nèi)部單元的空間排列方式都是不同的,其巖石的破裂路徑也因此存在差異。所以,巖石材料的極限應(yīng)力不僅受到材料物理參數(shù)非均質(zhì)度的影響,樣本中巖石單元的排列方式即細(xì)觀組構(gòu)也對(duì)巖石力學(xué)特性存在重要作用。
圖7 彈性模量期望50 GPa和變異系數(shù)0.2時(shí)不同組構(gòu)的應(yīng)力-應(yīng)變曲線Fig.7 Stress-strain curves of rock material of differentfabrics with a modulus of elasticity of 50 GPa anda coefficient of variation of 0.2
“組構(gòu)”是描述巖石中幾何性與物理性要素的重要概念,巖石內(nèi)所有統(tǒng)計(jì)單元以某種排列方式遍布整個(gè)區(qū)域,其排列方式形成某種組構(gòu)。理論上講,滿足同一期望和標(biāo)準(zhǔn)差的彈性模量單元集合有無數(shù)個(gè),都能形成不同的組構(gòu)。巖石材料的非均質(zhì)性不僅在于其物理參數(shù)的不均勻,也在于材料組構(gòu)的隨機(jī)性排列。如圖7所示,不同組構(gòu)的巖石樣本在單軸壓縮荷載時(shí)有不同的應(yīng)力-應(yīng)變曲線。
巖石內(nèi)部不同礦物顆粒的彈性模量和承載能力存在差異。顯然,如果2個(gè)接觸單元的彈性模量和承載力接近,2單元的應(yīng)力和應(yīng)變應(yīng)基本保持一致,穩(wěn)定與失穩(wěn)也會(huì)同步;如果2相鄰單元彈性模量差異較大,2單元的應(yīng)變大小會(huì)存在較大差異,接觸面成為容易發(fā)生破壞的“力學(xué)薄弱環(huán)節(jié)”。在巖石樣本中,因不同礦物成分力學(xué)性狀差異形成的力學(xué)薄弱環(huán)節(jié),因礦物組構(gòu)不同而呈現(xiàn)隨機(jī)分布的特征,也即在相同的非均質(zhì)度下,不同的數(shù)值試樣“力學(xué)薄弱環(huán)節(jié)”分布也不盡相同。這些薄弱環(huán)節(jié)往往是巖石壓縮過程中率先開裂和破壞的面,其過程如圖8所示。
圖8 非均質(zhì)巖石單元接觸面剪切過程示意圖Fig.8 Schematic diagram of the shearing process ofthe contact surface of heterogeneous rock element
由圖8可知,率先失去強(qiáng)度的單元會(huì)轉(zhuǎn)移原來承擔(dān)的部分荷載,周圍單元承擔(dān)的荷載則相應(yīng)地增大,增幅取決于單元自身與失穩(wěn)單元彈性模量的差值。很顯然,相鄰單元之間的彈性模量差異越大,應(yīng)力梯度則越大,越容易發(fā)生開裂破壞。因此,具有相同非均質(zhì)度而細(xì)觀組構(gòu)不同的數(shù)值試樣,力學(xué)薄弱環(huán)節(jié)分布的不同會(huì)造成不同的應(yīng)力分異格局,從而導(dǎo)致巖石破裂路徑不同。
圖9 單元彈性模量分布柱狀圖Fig.9 Histogram of elastic modulus distribution
圖10 非均質(zhì)度為0.5時(shí)不同組構(gòu)的巖石樣本破裂過程Fig.10 Rupture process of rock specimens of differentfabrics with a degree of heterogeneity of 0.5
圖9展示了巖石非均質(zhì)度為0.5時(shí)Voronoi單元的彈性模量在巖石樣本中的空間分布,每個(gè)柱體的高度值代表了Voronoi單元彈性模量大小,顯然非均質(zhì)巖石單元的“高度值”是大小不一的。非均質(zhì)度為0.5、巖石彈性模量期望為30 GPa時(shí),2個(gè)不同組構(gòu)的樣本裂紋萌生與擴(kuò)展路徑不同(如圖10所示)。該現(xiàn)象反映出同等非均質(zhì)度不同組構(gòu)的樣本單元分布結(jié)構(gòu)會(huì)對(duì)巖石破裂路徑會(huì)產(chǎn)生干擾,使裂紋擴(kuò)展過程受到力學(xué)薄弱環(huán)節(jié)分布結(jié)構(gòu)的影響,進(jìn)而影響巖石的強(qiáng)度力學(xué)特性。彈性模量較小的單元接觸面易成為巖石中的薄弱環(huán)節(jié),筆者篩選出彈性模量<8 GPa的單元位置,并繪出其分布圖(如圖11),將其與巖石樣本破裂時(shí)裂紋分布圖(如圖12)對(duì)比,兩者具有較好的對(duì)應(yīng)關(guān)系。因此,巖石樣本的極限承載力與巖石彈性模量單元的空間分布的組構(gòu)特征存在密切聯(lián)系。
圖11 巖石樣本彈性模量<8 GPa的單元分布Fig.11 Element distribution of rock specimen withelastic modulus smaller than 8 GPa
圖12 巖石破壞時(shí)裂隙密集區(qū)域分布Fig.12 Fracture distribution when rock is destroyed
(1)本文使用DDA方法對(duì)巖石樣本進(jìn)行數(shù)值建模與計(jì)算,將巖石樣本劃分為Voronoi結(jié)構(gòu)單元,用服從正態(tài)分布的隨機(jī)數(shù)賦值給單元的彈性模量,將剛性加載板采用固定加載位移的方法進(jìn)行單軸壓縮實(shí)驗(yàn),最終得到非均值脆性巖石材料的單軸壓縮應(yīng)力-應(yīng)變?nèi)^程曲線。
(2)標(biāo)準(zhǔn)差是一個(gè)數(shù)據(jù)集合離散程度的反映,用巖石單元彈性模量的標(biāo)準(zhǔn)差與期望之比(變異系數(shù))描述巖石強(qiáng)度的非均質(zhì)性是合理的:比值越大,表明巖石強(qiáng)度的非均質(zhì)度越大,反之巖石強(qiáng)度的非均質(zhì)度越小。
(3)巖石材料的變異系數(shù)和巖石極限承載力之間存在較好的線性關(guān)系:隨著變異系數(shù)的增大,巖石極限承載力減小。巖石的宏觀彈性模量同變異系數(shù)之間沒有較為明顯的線性變化特征。
(4)物理參數(shù)的非均質(zhì)性并非巖石強(qiáng)度的唯一影響因素,巖石的組構(gòu)方式也對(duì)巖石強(qiáng)度存在重要影響作用。樣本中薄弱環(huán)節(jié)不同的分布結(jié)構(gòu)會(huì)改變巖石裂紋的萌生與擴(kuò)展路徑,從而導(dǎo)致相同非均質(zhì)度的巖石力學(xué)強(qiáng)度性質(zhì)有所不同。
(5)巖石試樣受荷載作用由變形到破壞,是連續(xù)小變形逐漸發(fā)展成宏觀斷裂的非連續(xù)大變形的過程。傳統(tǒng)基于連續(xù)介質(zhì)的數(shù)值方法(如FEM)受到自身理論框架的限制,在處理大量顆粒間復(fù)雜接觸行為時(shí),模擬效果與計(jì)算準(zhǔn)確度不如DDA方法。DDA方法由于具有完備的運(yùn)動(dòng)學(xué)理論和嚴(yán)格的接觸算法,可以更加正確地描述顆粒之間的接觸力學(xué)行為。但是DDA方法仍存在一些不足,如對(duì)某些物理參數(shù)取值較為敏感,大規(guī)模計(jì)算迭代收斂慢的問題。