夏海城，鄔愛清，盧 波，徐棟棟

(長江科學(xué)院 水利部巖土力學(xué)與工程重點實驗室，武漢 430010)

1 研究背景

巖石材料是由礦物顆粒、膠結(jié)材料和空隙等組成的具有一定缺陷的不均勻物質(zhì)，非均質(zhì)性是巖石材料的本質(zhì)特征。這種細觀結(jié)構(gòu)的不均勻性導(dǎo)致巖石材料的彈性模量、強度、泊松比等宏觀力學(xué)參量具有一定的隨機性，同時也直接影響受載巖石的破裂演化過程。因此研究非均質(zhì)性對巖石材料力學(xué)性質(zhì)的影響機制具有重要意義。

國內(nèi)外學(xué)者對巖石非均質(zhì)問題展開了廣泛的研究，并取得諸多有價值的成果。巖石材料物理參數(shù)分布的不均勻性是其非均質(zhì)性的重要方面，胡小榮等[1]應(yīng)用地質(zhì)統(tǒng)計學(xué)方法對有限單元體力學(xué)參數(shù)的賦值進行了研究；陳永強[2]采用不同的統(tǒng)計分布模擬材料的非均勻性，得到了不同特性的材料破壞過程應(yīng)力-應(yīng)變曲線；康健等[3]研究了非均質(zhì)隨機分布的熱膨脹系數(shù)和熱傳導(dǎo)系數(shù)對巖石熱破裂的影響；馮增朝等[4-5]研究了巖石細胞單元特性及其非均質(zhì)分布對巖石破壞全過程曲線的影響，巖石亞細觀組構(gòu)的非均質(zhì)性及其分布規(guī)律是決定巖石單軸壓縮全過程曲線的重要因素；張征等[6]分析了巖土參數(shù)不確定性的主要原因，并研究了巖土參數(shù)空間變異性分析的原理、方法和步驟。唐春安等[7-8]、Kaiser等[9]用Weibull分布函數(shù)模擬巖石材料參數(shù)的非均質(zhì)性，研究了二維巖石模型脆性破壞的特點，提出“RFPA(Rock Failure Process Analysis)方法”。梁正召[10]開發(fā)了RFPA3D程序，模擬了巖石的破壞過程。尤名慶等[11]討論了巖石非均質(zhì)性和強度尺寸效應(yīng)之間的關(guān)系；羅榮等[12-13]研究了基于材料組分的非均質(zhì)巖土體數(shù)值建模方法及非均質(zhì)度對巖石力學(xué)特性的影響；柯長仁等[14]將非均質(zhì)性引入虛內(nèi)鍵模型，模擬了巖石式樣動態(tài)破壞過程；Cil和Buscarnera[15]利用Weibull 分布的最弱環(huán)節(jié)理論，分析了不同粒徑尺寸對單軸最大抗壓強度的影響?？偟膩碚f，在巖石非均質(zhì)的模擬方面的研究都采用基于某種統(tǒng)計分布對每個巖石單元賦予各不相同的物理參數(shù)，使細觀單元的力學(xué)性質(zhì)各不相同，以模擬巖石宏觀上的非均質(zhì)性。

非均質(zhì)巖石材料力學(xué)特性的影響因素不僅在于巖石細觀材料參數(shù)之間的差異，這種差異形成的力學(xué)薄弱環(huán)節(jié)在空間的分布特征直接控制著受載巖石的破裂演化過程，從而影響巖石的強度和應(yīng)力-應(yīng)變?nèi)^程曲線的特征。上述研究主要關(guān)注因巖石成分不同導(dǎo)致的細觀材料參數(shù)的差異對巖石宏觀力學(xué)性狀的影響。事實上，因成巖機制和成巖過程的差異，巖石的成分、結(jié)構(gòu)和構(gòu)造在空間上呈現(xiàn)不均勻分布特征；巖石的彈性模型、強度等宏觀物理力學(xué)參量的離散性則是這種細觀不均勻性的綜合體現(xiàn)。既有研究成果未充分關(guān)注細觀不均勻性的空間分布特征對巖石宏觀力學(xué)特性的影響。本文采用具有伺服控制機制的非連續(xù)變形分析(Discontinuous Deformation Analysis,DDA)方法，根據(jù)巖石顆粒呈不規(guī)則多邊形的特點，采用Voronoi結(jié)構(gòu)來模擬巖石礦物顆粒，構(gòu)建用于DDA分析的巖石數(shù)值試樣，通過對不同的Voronoi塊體單元賦予不同的物理力學(xué)參數(shù)的方式來實現(xiàn)對非均質(zhì)性的模擬。通過開展巖石單軸壓縮數(shù)值實驗，在驗證非均質(zhì)度對巖石強度影響的基礎(chǔ)上，探究了巖石單元空間分布的影響作用，補充了主要結(jié)論。

2 基于DDA的細觀數(shù)值實驗方法及數(shù)值模型構(gòu)建

2.1 DDA基本理論

DDA方法是石根華博士于20世紀80年代提出的平行于有限單元法的數(shù)值模擬方法，其單元或塊體形狀可以是任何凸狀形或凹狀形的多邊形。[16]相比有限單元法只能用標準形狀單元具有一定的優(yōu)越性，在DDA方法中，任意塊體的位移Di可用6個位移變量表示，即

Di=(u0v0r0εxεyγxy) ；

(1)

式中：(u0v0)是塊體內(nèi)特殊點(x0y0)的剛體位移；角r0是塊體繞轉(zhuǎn)動中心(u0v0)的轉(zhuǎn)動角，用弧度給出；(εxεyγxy)是該塊體的法向和切向應(yīng)變。

塊體中任意點(x，y)處的位移可由變形變量[Di]表示為

(2)

[Ti]=

(3)

式中[Ti]為位移轉(zhuǎn)換矩陣。

各連接塊體依靠接觸和位移約束形成塊體系統(tǒng)。n個塊體的聯(lián)立方程為

式中Kij是6×6子矩陣，與材料特性有關(guān)；[Di]和[Fi]是6×1子矩陣，此處Di代表塊體i的變形變量，F(xiàn)i為塊體i上的6個變形變量的荷載。對于平衡方程(式(4))，由力和應(yīng)力產(chǎn)生的總勢能Π最小化來推導(dǎo)，對Π求導(dǎo)，得出的總勢能最小。子矩陣Kij為

(5)

子矩陣Fi是

(6)

式中d表示塊體的位移變量。

在對式(2)和式(4)求解后，可計算出塊體變形后的位移與應(yīng)變。由于DDA法中塊體無嵌入和張拉的條件，需在接觸位置加上或去掉剛性彈簧，修正總體平衡方程，并再次求解直至所有邊界滿足無嵌入和張拉，迭代完成。

楊懋偲[17]在DDA源代碼的基礎(chǔ)上進行了二次開發(fā)，初步實現(xiàn)了對電液伺服控制機制的模擬，得到了符合室內(nèi)常規(guī)實驗結(jié)果的應(yīng)力-應(yīng)變?nèi)^程曲線。

本文將使用具有伺服控制機制的DDA程序開展非均質(zhì)度對巖石材料宏觀力學(xué)性質(zhì)的影響機制研究。

2.2 正態(tài)分布函數(shù)及非均質(zhì)度

巖石材料物理參數(shù)的不均勻分布是引起巖石非線性變形的重要因素，在細觀數(shù)值實驗中，將巖石試樣劃分為多個Voronoi單元，每個Voronoi單元具有不同的物理力學(xué)參數(shù)，這種細觀單元物理力學(xué)參數(shù)的不均勻分布可以較好地模擬出巖石非均質(zhì)性的特點。本文采用正態(tài)分布函數(shù)模擬礦物顆粒物理參數(shù)分布的隨機性，并定義了量化巖石非均質(zhì)度大小的方法。通過調(diào)整巖石試樣非均質(zhì)度和期望的大小，比較分析了巖石力學(xué)性能的變化規(guī)律，并探究巖石力學(xué)性能的影響因素。

構(gòu)成巖石的礦物顆粒物理力學(xué)參數(shù)盡管呈現(xiàn)出強烈的隨機性特點，但仍然滿足一定的統(tǒng)計學(xué)分布規(guī)律。根據(jù)以往的研究成果，采用不同的參數(shù)分布模型開展數(shù)值實驗得到的破壞過程曲線也不一樣。目前，常用于模擬巖石非均質(zhì)性的隨機函數(shù)模型有Weibull分布模型和正態(tài)分布模型，其中正態(tài)分布模型得到的曲線更加符合巖石脆性破壞的特點[2]。本文采用正態(tài)分布函數(shù)對巖石非均質(zhì)進行表征。正態(tài)分布函數(shù)的概率密度函數(shù)如式(7)，其中μ為期望，反映了非均質(zhì)巖石材料的綜合性能；σ為方差，反映了巖石材料彈性模量的偏離值。

(7)

將服從正態(tài)分布的隨機數(shù)賦給巖石單元的彈性模量，以模擬巖石材料的非均質(zhì)性。正態(tài)分布函數(shù)中，標準差反映的巖石彈性模量偏離值是度量巖石非均質(zhì)性的重要參數(shù)。為準確定義巖石的非均質(zhì)度，有效量化非均質(zhì)度的大小，本文將巖石樣本中彈性模量的標準差與期望之比Cv定義為變異系數(shù)來度量巖石非均質(zhì)度的大小，即

(8)

式中：變異系數(shù)Cv無量綱，反映了巖石單元隨機物理參數(shù)的偏離度，用來度量巖石的非均質(zhì)性；Ed為巖石單元彈性模量的標準差；Em為巖石單元彈性模量的平均值。

生成服從正態(tài)分布的隨機數(shù)組時，通過調(diào)整標準差和期望的大小來控制變異系數(shù)。顯然在Cv=0時，巖石為完全均質(zhì)材料，各細胞單元具有完全相同的物理參數(shù)。隨著Cv增大，巖石材料彈性模量的非均質(zhì)度增大，單元彈性模量離散度增強，巖石均勻性變差。Cv越小，巖石單元的物理參數(shù)越趨向于期望，巖石整體的均勻性越好。

2.3 非均質(zhì)材料數(shù)值模型與參數(shù)賦值

為聚焦巖石材料非均質(zhì)度對其力學(xué)性能的影響，排除非均勻度以外單元分布結(jié)構(gòu)的干擾，本文構(gòu)建2組不同彈性模量期望的非均質(zhì)巖石數(shù)值模型(如表1)，使用正態(tài)分布函數(shù)產(chǎn)生的隨機數(shù)組對巖石單元彈性模量進行賦值，每種變異系數(shù)生成6個不同的隨機數(shù)組(也即形成6個數(shù)值試樣)，滿足相同的期望和標準差，并將6個試樣進行相同條件的數(shù)值實驗。不同變異系數(shù)的巖石樣本單元彈性模量統(tǒng)計數(shù)目直方圖如圖1所示。

表1 正態(tài)分布模型參數(shù)及非均質(zhì)度Table 1 Parameters of normal distribution modeland degree of heterogeneity

圖1 Voronoi單元統(tǒng)計數(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ù)值模型，為真實模擬巖石細胞單元的形狀結(jié)構(gòu)，采用Voronoi多邊形細觀單元形成的塊體系統(tǒng)進行模擬，每個Voronoi單元的面積約為5 cm2，共計975個單元。

本文中數(shù)值實驗為單軸壓縮實驗，對剛性加載板采用施加固定位移的方法對巖石模型進行壓縮，每個時步的加載位移為5×10-8。在單軸壓縮實驗中巖石的非均質(zhì)性主要是由巖石Voronoi單元之間物理參數(shù)不均勻引起的，為排除其他物理參數(shù)干擾的影響，本文所有巖石Voronoi單元泊松比為0.3，Voronoi單元之間的界面強度設(shè)置為統(tǒng)一值。

2.4 非均質(zhì)巖石材料破壞過程曲線

本文中數(shù)值加載模型主要可分為剛性加載板、非均質(zhì)巖石模型和剛性支架3部分，使用混合高階DDA數(shù)值實驗平臺，通過對剛性加載板施加豎向位移進行伺服控制單軸壓縮數(shù)值實驗。

圖3 全過程應(yīng)力-應(yīng)變曲線Fig.3 Complete stress-strain curve

圖3為巖石樣本應(yīng)力-應(yīng)變?nèi)^程關(guān)系曲線，巖石樣本軸向應(yīng)變通過加載板的位移計算，加載板在軸向的位移量在數(shù)值上等于巖石試樣的變形量。巖石塊體軸向應(yīng)力隨著加載板的不斷加載而增大，并在到達極限荷載后承載力迅速下降，下降到某一值后保留一定的殘余強度，符合真實的脆性巖石的承載規(guī)律。顯然巖石變形破壞曲線中應(yīng)力最大的點也是巖石極限承載力的點，加載板施加的荷載力最大時，巖石應(yīng)力達到最大。巖石的宏觀彈性模量為達到峰值應(yīng)力50%時，應(yīng)力與對應(yīng)應(yīng)變的切線模量。

圖4 均質(zhì)巖石應(yīng)力-應(yīng)變?nèi)^程曲線Fig.4 Complete stress-strain curve of homo-geneous rock

3 巖石材料強度非均質(zhì)度對其力學(xué)性能的影響

巖石樣本為完全均質(zhì)時，應(yīng)力-應(yīng)變?nèi)^程曲線如圖4所示，此時所有Voronoi單元的彈性模量均為30 GPa。以此為參照，研究材料非均質(zhì)度對巖石宏觀彈性模量和極限承載力的影響，對每組變異系數(shù)的巖石試樣均開展數(shù)值模擬單軸壓縮實驗，取巖石樣本全過程破壞曲線較為平均的一組對比。

圖5為彈性模量期望為30、50 GPa時不同變異系數(shù)下的應(yīng)力-應(yīng)變曲線。由圖5可知，隨著變異系數(shù)的增大，巖石的極限承載力逐漸下降，峰后曲線具有明顯的差異。對于巖石的宏觀彈性模量，初始加載階段不同變異系數(shù)的巖石樣本宏觀彈性模量幾乎相同，隨著加載的不斷進行，變異系數(shù)較大的巖石宏觀彈性模量先出現(xiàn)下降，總體來看非均質(zhì)度對巖石樣本宏觀彈性模量的影響是較弱的。

圖5 不同非均質(zhì)度模型應(yīng)力-應(yīng)變曲線Fig.5 Stress-strain curves of rock models ofdifferent degrees of heterogeneity

為了有效研究非均質(zhì)度對巖石極限承載力的影響，消除相同變異系數(shù)下不同單元分布對變形曲線的干擾，聚焦巖石極限應(yīng)力隨非均質(zhì)度的變化趨勢，對每個變異系數(shù)下6組彈性模量單元的數(shù)值實驗結(jié)果取平均值，得到平均極限應(yīng)力，并做出變異系數(shù)-極限應(yīng)力的折線圖，結(jié)果如圖6所示。可以看出，服從正態(tài)分布的巖石單元彈性模量，期望值分別為30 GPa和50 GPa時，巖石極限承載力均隨著材料變異系數(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)，在峰值和峰后階段具有較大差異。在給巖石單元的彈性模量賦值時，服從某種正態(tài)分布的單元集合具有多種分布的可能性，這導(dǎo)致了在相同變異系數(shù)的巖石樣本中，每個樣本內(nèi)部單元的空間排列方式都是不同的，其巖石的破裂路徑也因此存在差異。所以，巖石材料的極限應(yīng)力不僅受到材料物理參數(shù)非均質(zhì)度的影響，樣本中巖石單元的排列方式即細觀組構(gòu)也對巖石力學(xué)特性存在重要作用。

圖7 彈性模量期望50 GPa和變異系數(shù)0.2時不同組構(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

4 材料空間組構(gòu)對巖石力學(xué)性質(zhì)的影響

“組構(gòu)”是描述巖石中幾何性與物理性要素的重要概念，巖石內(nèi)所有統(tǒng)計單元以某種排列方式遍布整個區(qū)域，其排列方式形成某種組構(gòu)。理論上講，滿足同一期望和標準差的彈性模量單元集合有無數(shù)個，都能形成不同的組構(gòu)。巖石材料的非均質(zhì)性不僅在于其物理參數(shù)的不均勻，也在于材料組構(gòu)的隨機性排列。如圖7所示，不同組構(gòu)的巖石樣本在單軸壓縮荷載時有不同的應(yīng)力-應(yīng)變曲線。

巖石內(nèi)部不同礦物顆粒的彈性模量和承載能力存在差異。顯然，如果2個接觸單元的彈性模量和承載力接近，2單元的應(yīng)力和應(yīng)變應(yīng)基本保持一致，穩(wěn)定與失穩(wěn)也會同步；如果2相鄰單元彈性模量差異較大，2單元的應(yīng)變大小會存在較大差異，接觸面成為容易發(fā)生破壞的“力學(xué)薄弱環(huán)節(jié)”。在巖石樣本中，因不同礦物成分力學(xué)性狀差異形成的力學(xué)薄弱環(huán)節(jié)，因礦物組構(gòu)不同而呈現(xiàn)隨機分布的特征，也即在相同的非均質(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可知，率先失去強度的單元會轉(zhuǎn)移原來承擔(dān)的部分荷載，周圍單元承擔(dān)的荷載則相應(yīng)地增大，增幅取決于單元自身與失穩(wěn)單元彈性模量的差值。很顯然，相鄰單元之間的彈性模量差異越大，應(yīng)力梯度則越大，越容易發(fā)生開裂破壞。因此，具有相同非均質(zhì)度而細觀組構(gòu)不同的數(shù)值試樣，力學(xué)薄弱環(huán)節(jié)分布的不同會造成不同的應(yīng)力分異格局，從而導(dǎo)致巖石破裂路徑不同。

圖9 單元彈性模量分布柱狀圖Fig.9 Histogram of elastic modulus distribution

圖10 非均質(zhì)度為0.5時不同組構(gòu)的巖石樣本破裂過程Fig.10 Rupture process of rock specimens of differentfabrics with a degree of heterogeneity of 0.5

圖9展示了巖石非均質(zhì)度為0.5時Voronoi單元的彈性模量在巖石樣本中的空間分布，每個柱體的高度值代表了Voronoi單元彈性模量大小，顯然非均質(zhì)巖石單元的“高度值”是大小不一的。非均質(zhì)度為0.5、巖石彈性模量期望為30 GPa時，2個不同組構(gòu)的樣本裂紋萌生與擴展路徑不同(如圖10所示)。該現(xiàn)象反映出同等非均質(zhì)度不同組構(gòu)的樣本單元分布結(jié)構(gòu)會對巖石破裂路徑會產(chǎn)生干擾，使裂紋擴展過程受到力學(xué)薄弱環(huán)節(jié)分布結(jié)構(gòu)的影響，進而影響巖石的強度力學(xué)特性。彈性模量較小的單元接觸面易成為巖石中的薄弱環(huán)節(jié)，筆者篩選出彈性模量<8 GPa的單元位置，并繪出其分布圖(如圖11)，將其與巖石樣本破裂時裂紋分布圖(如圖12)對比，兩者具有較好的對應(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 巖石破壞時裂隙密集區(qū)域分布Fig.12 Fracture distribution when rock is destroyed

5 結(jié) 論

(1)本文使用DDA方法對巖石樣本進行數(shù)值建模與計算，將巖石樣本劃分為Voronoi結(jié)構(gòu)單元，用服從正態(tài)分布的隨機數(shù)賦值給單元的彈性模量，將剛性加載板采用固定加載位移的方法進行單軸壓縮實驗，最終得到非均值脆性巖石材料的單軸壓縮應(yīng)力-應(yīng)變?nèi)^程曲線。

(2)標準差是一個數(shù)據(jù)集合離散程度的反映，用巖石單元彈性模量的標準差與期望之比(變異系數(shù))描述巖石強度的非均質(zhì)性是合理的：比值越大，表明巖石強度的非均質(zhì)度越大，反之巖石強度的非均質(zhì)度越小。

(3)巖石材料的變異系數(shù)和巖石極限承載力之間存在較好的線性關(guān)系：隨著變異系數(shù)的增大，巖石極限承載力減小。巖石的宏觀彈性模量同變異系數(shù)之間沒有較為明顯的線性變化特征。

(4)物理參數(shù)的非均質(zhì)性并非巖石強度的唯一影響因素，巖石的組構(gòu)方式也對巖石強度存在重要影響作用。樣本中薄弱環(huán)節(jié)不同的分布結(jié)構(gòu)會改變巖石裂紋的萌生與擴展路徑，從而導(dǎo)致相同非均質(zhì)度的巖石力學(xué)強度性質(zhì)有所不同。

(5)巖石試樣受荷載作用由變形到破壞，是連續(xù)小變形逐漸發(fā)展成宏觀斷裂的非連續(xù)大變形的過程。傳統(tǒng)基于連續(xù)介質(zhì)的數(shù)值方法(如FEM)受到自身理論框架的限制，在處理大量顆粒間復(fù)雜接觸行為時，模擬效果與計算準確度不如DDA方法。DDA方法由于具有完備的運動學(xué)理論和嚴格的接觸算法，可以更加正確地描述顆粒之間的接觸力學(xué)行為。但是DDA方法仍存在一些不足，如對某些物理參數(shù)取值較為敏感，大規(guī)模計算迭代收斂慢的問題。