羅 榮,曾亞武
(武漢大學土木建筑工程學院,武漢 430072)
巖石作為一種地質材料,非均質性是其基本特性之一[1]。造成巖石的非均質特性的主要原因是在漫長的地質年代里,巖石經歷結晶過程的演化、高溫高壓的影響、以及各種地質營力的作用,使巖石具有不同的細觀結構(礦物顆粒組成和各類缺陷)。巖石的非均質性對巖石的力學特性及受荷載后的力學行為有很大的影響,研究表明:巖石內部微觀介質參數的不均勻性是造成巖石力學性質非均勻和非線性的主要原因[2]。
國內外學者就巖石的非均質特征對巖石力學特性的影響開展了一系列工作,并取得了大量研究成果。Song 和 Kim[3],Napier和 Dede[4],Li[5],楊強[6]等通過隨機指定格構單元的強度和形狀描述巖石材料的非均質性模擬了原巖的破壞過程。張征[7]、譚文輝[8]、趙紅亮[9]等基于地質統計學方法的區域化變量理論,考慮巖石力學參數的空間變異二重性特征研究了巖石的非均質特性。唐春安和Kaiser[10-11]考慮巖石材料參數的非均勻性,假定單元參數服從某種隨機分布,對二維脆性破裂進行了詳細研究,唐春安[12-13]在此基礎上進一步開發了RFPA模型模擬巖石破壞過程。馬志濤等[14]利用Weibull分布對巖石材料的非均質性進行描述,從能量角度建立模擬細觀層次上巖石破壞演化的Mh-PCA模型。馮增朝等[15]利用Weibull分布描述巖石的非均質參數,并研究了不同的非均質統計分布對巖石全曲線性態的影響。岳中琦等[16]提出了一種基于數字圖像的非均質巖土工程結構的二維數值分析方法,利用數字圖像技術獲取非均質巖石的細觀結構,并結合有限差分方法分析巖石的破壞過程。總體來說,都是通過對組成巖石模型的各細觀單元賦予不同的力學參數模擬巖石的非均質特性。
巖石是由不同礦物顆粒集合體和膠結材料組成的非均質體,巖石中的礦物成分及其含量是決定巖石物理力學性質的主要因素,在某些條件下甚至會產生決定性影響[17],巖石的組成非均質性主要為巖石的礦物顆粒的非均質構成。一般來說,在巖石局部范圍內,相同的礦物顆粒集合體具有相近的力學性質,不同的礦物顆粒集合體具有不同的力學性質。本文利用基于巖石礦物組成及含量的巖石礦物細胞元隨機性參數賦值方法描述巖石的非均質性,通過數值試驗研究巖石非均質構成對巖石宏觀力學特性的影響。
在數值模擬中,將巖石劃分為若干單元,單元尺寸為0.2~0.5 mm,并假設每一個單元覆蓋一個巖石顆粒,可將這些細觀的巖石組構稱之為巖石的細胞元[15]。巖石礦物細胞元隨機性參數賦值方法:假設每一個巖石細胞元僅由一種礦物構成,并將其稱為礦物細胞元;假定同屬性細胞元在局部范圍內具有相同的物理力學性質,不同屬性的細胞元物理力學性質不同;利用巖石的礦物含量定義細胞元類別的判定區間,根據Monte Carlo方法對各個細胞元進行類別判定和參數賦值,從而得出非均質巖石的數值計算模型。由于礦物細胞元隨機性參數賦值方法考慮了巖石的組成礦物種類及含量,因此可很好地用于研究巖石的非均質構成對巖石宏觀力學特性的影響。
當i=m時,ui=1,ui為一組單調遞增的數列。定義區間 Ai=[ui-1,ui](其中 i=1,2,…m - 1;當i=m,取右閉區間)為非均質巖石礦物細胞元的類別判定區間,其中子區間Ai即為第i類別礦物細胞元的判定區間。由于區間Ai在[0,1]區間所占比例為 ui-ui-1=ni,ni又為礦物 i的實際含量,所以定義的礦物細胞元類別判定區間Ai與礦物i的含量具有一一對應關系。
巖石礦物細胞元隨機性參數賦值方法利用Monte Carlo方法對非均質巖石單元進行隨機判定。根據非均質巖石剖分的有限元單元數,產生一組服從[0,1]均勻分布的隨機數序列 ξk,選用[0,1]的均勻分布作為隨機序列目標分布,概率分布函數為F(x)=x,概率密度函數為f(x)=1。根據 Monte Carlo方法,對于初始隨機數 ξk,令F(xk)= ξk,即得到用以進行礦物細胞元類別判定的數序列xk=ξk,利用所得類別判定隨機數序列依次對各單元進行礦物細胞元類別判定。
礦物細胞元類別判定過程如圖1所示,對于第k號單元,判定類別隨機數xk所屬區間,若xk∈Ai,則判定該單元為第i類別礦物細胞元,并對該礦物細胞元賦予相應的i類別物理力學參數。依次對所有的單元體進行類別判定和賦值,即可完成有限元模型的非均質賦值。
礦物細胞元類別判定數xk是服從[0,1]均勻分布的隨機數序列,故單元判定為第i類別礦物細胞元的概率 Pi=uiui-1=ni,其中ni等于巖石第i類別礦物的含量,因此利用該礦物細胞元判定區間進行單元類別判定可使得模型中各類別礦物細胞元的含量與巖石實際礦物含量相等;同時,由于類別判定隨機數是服從[0,1]的均勻分布的隨機數序列,所以賦值完成后的模型各礦物細胞元是隨機均勻混合的。

圖1 礦物細胞元類別判定示意圖Fig.1 Schematic diagram of the sort judgment for mineral cell unit
對于地表巖石結構,形成材料宏觀非線性的主要原因并不是微觀介質的非線性,而是微觀介質參數的不均勻性[2]。因此本文選用理想彈塑性模型描述礦物細胞元的變形特征,并假定在同一模型內,礦物細胞元的強度隨彈性模量增大而增大,且礦物細胞元強度、彈性模量的變化是成比例的[15],即強度越高,彈性模量越大;反之,都小。即

式中:E',R'分別為礦物細胞元的彈性模量和屈服強度;E0,R0分別為對應數值計算模型的所有細胞元的彈性模量、屈服強度的平均值(期望值)。
數值試驗采用平面應變模型,幾何尺寸為20 cm×10 cm,有限元單元(即巖石礦物細胞元)尺寸為0.5 mm ×0.5 mm,共 劃 分400×200個單元,數值試驗采用位移加載模式,計算簡圖如圖2所示,其中不同灰度代表不同的礦物細胞元類別。

圖2 模型的幾何特征及邊界條件Fig.2 Geometry and boundary conditions of the model
本文所有數值模擬計算中,模型的各種物理力學參數取相同的期望值,彈性模量E0=50 GPa,屈服強度R0=50 MPa,由于各類巖石材料泊松比差別不大,數值試驗泊松比均取ν=0.3。
本文取礦物細胞元含量相同、力學參數不同的非均質巖石構成建立非均質巖石模型,分別對2種礦物細胞元類型(n1=n2=0.5)和3種礦物細胞元類型(n1=n2=n3=0.333)的混合模型進行單軸壓縮試驗,其各類別礦物細胞元參數分別見表1、表2。

表1 2類礦物細胞元混合模型計算參數Table 1 Calculation parameters of hybrid model with two mineral cell units

表2 3類礦物細胞元混合模型計算參數Table 2 Calculation parameters of hybrid model with three mineral cell units
數值試驗采用位移加載方式,巖石宏觀軸向應變按下式計算:

式中:Δh為試件壓縮變形;h為試件壓縮變形前的高度。
巖石宏觀軸向應力通過頂部邊界的約束反力來表示,即

式中:n為頂部邊界節點數;Fi為頂部邊界上i節點的約束反力;l為模型的寬度;d為模型的厚度(取單位厚度)。
根據表1、表2所述參數分別對2類礦物細胞元和3類礦物細胞元混合的非均質巖石模型進行加載試驗,各試件宏觀應力-應變關系曲線結果整理如圖3、圖4所示。

圖3 2類礦物細胞元混合模型的應力-應變曲線Fig.3 Stress-strain curves of hybrid model with two cell units

圖4 3類礦物細胞元混合模型的應力-應變曲線Fig.4 Stress-strain curves of hybrid model with three cell units
根據圖3、圖4所示結果,可得出以下結論:
(1)試件1的各類礦物細胞元參數相同,為均質巖石模型,巖石宏觀應力-應變關系曲線表現為礦物細胞元的變形特征——理想彈塑性,從線彈性階段到理想塑性階段具有明顯的屈服轉折點;對于其它非均質巖石模型,宏觀應力-應變關系曲線除具有線彈性階段和理想塑性階段外,還具有明顯的非線性屈服段,宏觀應力-應變關系曲線由礦物細胞元的變形特征和巖石的非均質構成共同作用決定。
(2)在低應力狀態下,由于各礦物細胞元均處于彈性變形階段,各試件也表現為宏觀線彈性變形特征。由于非均質巖石由不同類別的礦物細胞元混合組成,隨著荷載的增加,當強度較低的礦物細胞元達到屈服,應力將發生轉移和重分布,隨荷載的繼續增加,進入屈服的礦物細胞元將逐步增多,直到所有礦物細胞元達到屈服,荷載不再增加,表現出理想塑性的變形特征。礦物細胞元逐步發生屈服的過程在宏觀應力-應變關系曲線上表現出明顯的非線性特征。
(3)各組非均質巖石模型的宏觀彈性模量、極限強度均不相同,且都比均質巖石宏觀彈性模量和極限強度低,說明巖石的非均質構成對巖石的宏觀力學性能有較大的影響。
為研究非均質巖石構成差異性對巖石宏觀受力性能的影響,本文利用巖石的礦物構成定義巖石的非均質度。
因巖石礦物細胞元隨機性參數賦值方法具有隨機性特征,可利用無量綱量標準差率Cv表示模型中所有礦物細胞元參數的隨機離散程度,即巖石的非均質度。根據礦物細胞元隨機性參數賦值方法,任一單元被判定為第i類別礦物細胞元的概率等于組成巖石的第i類礦物的含量ni,以彈性模量E為例,即礦物細胞元參數為Ei的概率為ni,利用各礦物細胞元彈性模量計算巖石的非均質度,

式中σ(Ei),E(Ei)分別為模型各礦物細胞元彈性模量值的標準差、期望值。
由式(2)可知,礦物細胞元強度、彈性模量是成比例的,因此對于本文由礦物細胞元強度計算的非均質度與彈性模量計算得到的結果是一致的。
巖石礦物細胞元隨機性參數賦值方法基于非均質巖石的真實礦物構成,且非均質度的計算僅與非均質巖石的各礦物種類及含量有關,因此該非均質度的計算公式對于真實巖石依然成立。根據非均質度的定義可知,當Cv=0時,巖石模型的各礦物細胞元參數值的標準差為零,即各礦物細胞元參數均相等,為均質巖石;當Cv>0時,為非均質巖石;在巖石模型所有細胞元的參數均值不變時,Cv越大,各礦物細胞元參數標準差越大,巖石均質性越差。
非均質巖石的宏觀彈性模量是構成非均質巖石的各類礦物細胞元均未進入屈服狀態所表現的彈性變形力學特征。根據第3節數值試驗參數及計算結果,分別計算各試件的非均質度和宏觀彈性模量,結果見圖5。
通過圖5可知,對于2類礦物細胞元模型和3類礦物細胞元模型,當非均質度等于0時,均質巖石宏觀彈性模量即為礦物細胞元的彈性模量E0;當非均質度較小時(Cv<0.3),宏觀彈性模量變化較小,基本與非均質巖石模型的所有礦物細胞元彈模均值E0相等;當非均質度較大(Cv>0.3)時,宏觀彈模逐漸變小,由礦物細胞元參數和巖石非均質度共同決定。對于2類礦物細胞元模型和3類礦物細胞元模型,其宏觀彈性模量隨巖石非均質度的變化具有近似相同的曲線規律。

圖5 宏觀彈模隨非均質度的變化曲線Fig.5 Variation of elastic modulus of model against inhomogeneity degree
本文采用理想彈塑性模型描述巖石礦物細胞元的變形特征,當構成非均質巖石的所有礦物細胞元均達到屈服后,表現出理想塑性的變形特征,此時應力即為非均質巖石的極限強度。根據第3節數值試驗參數及計算結果,分別計算各試件的非均質度和極限強度,結果見圖6。

圖6 極限強度隨非均質度的變化曲線Fig.6 Variation of ultimate strength of model against inhomogeneity degree
通過圖6可知,當非均質度等于0時,均質巖石的極限強度即為礦物細胞元的極限強度;隨著非均質增大,礦物細胞元力學特性差異性增大,極限強度降低,這是因為在外荷載作用下強度較低的細胞元發生屈服后,將作為承載缺陷引起巖石宏觀承載能力的下降。對于非均質巖石,極限強度將由礦物細胞元強度和巖石非均質度共同決定。對于2類礦物細胞元模型和3類礦物細胞元模型,極限強度隨巖石非均質度的變化具有相同的曲線規律。
巖石的宏觀彈模和極限強度分別反映巖石在外荷載作用下的彈性變形特征和整體抵抗屈服破壞的能力。根據以上試驗結果,巖石非均質度是影響巖石宏觀力學參數的重要因素,巖石非均質度越小,巖石各礦物細胞元之間力學參數差異越小,巖石越均勻,宏觀總體力學參數也越接近礦物細胞元參數期望值;反之,非均質度越大,各礦物細胞元之間力學參數差異越大,巖石越不均勻,宏觀總體力學參數越小,即巖石非均質性對巖石力學特性具有弱化效應。比較分析2類礦物混合模型和3類礦物混合模型的試驗結果,可以看出,非均質度是描述巖石宏觀力學參數受巖石非均質性弱化影響的重要指標參數。總體上說,當巖石為均質體時,其宏觀力學參數由礦物細胞元參數決定,當巖石為非均質體時,其宏觀力學參數由礦物細胞元參數和非均質度共同決定。
巖石是由不同礦物顆粒組成的非均質體,本文利用巖石礦物細胞元隨機性參數賦值方法研究了非均質巖石構成對其宏觀應力-應變關系曲線的影響;并利用巖石礦物構成定義了巖石的非均質度,研究了巖石非均質度對巖石宏觀力學參數的影響,得出以下結論:
(1)均質巖石的宏觀應力-應變關系曲線與構成巖石的礦物細胞元一致,表現為理想彈塑性變形特征;非均質巖石宏觀應力-應變關系曲線由細胞元變形特征和巖石非均質構成共同決定,表現為線彈性、逐步屈服和理想塑性3個階段;逐步屈服階段表現出非線性特征,非均質度越大時,巖石均質性越差,非線性特征越明顯。
(2)非均質性對巖石宏觀力學參數具有弱化影響。根據非均質巖石構成定義的巖石非均質度反映了非均質巖石構成的各礦物顆粒力學性能的差異性大小,非均質度越小,巖石越均勻。當非均質度不大時(Cv<0.3),非均質性對宏觀彈模的弱化作用較小,當非均質度較大時(Cv>0.3),宏觀彈模受非均質性影響較大;非均質巖石的極限強度受非均質特性的弱化影響與彈性模量一致。
(3)通過對2類礦物細胞元混合模型和3類礦物細胞元混合模型的試驗結果比較分析,非均質巖石模型宏觀力學參數隨非均質度的變化規律一致,非均質度是描述巖石宏觀力學參數受巖石非均質性弱化影響的重要指標參數。
綜上所述,非均質巖石與均質巖石的力學特征具有明顯區別:對于均質巖石,其力學特性由非均質巖石構成礦物細胞元力學特性決定;對于非均質巖石,其力學特性由細胞元的力學特性和巖石非均質特性組成共同決定。
[1]唐春安.巖石破裂過程數值試驗[M].北京:科學出版社,2003.(TANG Chun-an.Numerical Test of Rock Failure Process[M].Beijing:Science Press,2003.(in Chinese))
[2]唐春安.巖石聲發射規律數值模擬初探[J].巖石力學與工程學報,1997,16(4):368 -374.(TANG Chunan.Numerical Simulation of AE in Rock Failure[J].Chinese Journal of Rock Mechanics and Engineering,1997,16(4):368 -374.(in Chinese))
[3]SONG J,KIM K.Micromechanical Modeling of the Dynamic Fracture Process During Rock Blasting[J].International Journal of Rock Mechanics and Mining Sciences& Geomechanics Abstracts,1996,33(4):387 -394.
[4]NAPIER J A L,DDEDE T.A Comparison Between Random Mesh Schemes and Explicit Growth Rules for Rock Fracture Simulation[J].International Journal of Rock Mechanics and Mining Sciences&Geomechanics Abstracts,1997,34(3 -4):356.
[5]LI L H,BAI Y L,XIA M F,et al.Damage Localizations as a Possible Precursor of Earthquake Rupture[J].Pure and Applied Geophysics,2000,157:1929 -1943.
[6]楊 強,張 浩,吳榮宗.二維格構模型在巖石類材料開裂模擬中的應用[J].巖石力學與工程學報,2000,19(增1):941 -945.(YANG Qiang,ZHANG Hao,WU Rong-zong.Lattice Model Simulation of Fracture in Disordered Brittle Materials[J].Chinese Journal of Rock Mechanics and Engineering,2000,19(Sup.1):941 - 945.(in Chinese))
[7]張 征,劉淑春,鞠碩華.巖土參數空間變異性分析原理與最優估計模型[J].巖土工程學報,1996,18(4):40 - 47.(ZHANG Zheng,LIU Shu-chun,JU Shuo-hua.The Optimum Estimation Model and the Principle of Spatial Variability Analysis of Rock and Soil Parameters[J].Chinese Journal of Geotechnical Engineering,1996,18(4):40 -47.(in Chinese))
[8]譚文輝,王家臣,周汝弟.巖體強度參數空間變異性分析[J].巖石力學與工程學報,1999,18(5):546-549.(TAN Wen-hui,WANG Jia-chen,ZHOU Ru-di.Analysis on Spatial Variability of Strength Parameter of Rock Mass[J].Chinese Journal of Rock Mechanics and Engineering,1999,18(5):497 -502.(in Chinese))
[9]趙紅亮,馮夏庭,張東曉,等.巖土力學參數空間變異性的集合卡爾曼濾波估值[J].巖土力學,2007,28(10):2219 -2223.(ZHAO Hong-liang,FENG Xia-ting,ZHANG Dong-xiao,et al.Spatial Variability of Geomechanical Parameter Estimation via Ensemble Kalman Filter[J].Rock and Soil Mechanics,2007,28(10):2219 -2223.(in Chinese))
[10]TANG Chun-an,KAISER P K.Numerical Simulation of Cumulative Damage and Seismic Energy Release During Brittle Rock Failure-Part I:Fundamentals[J].International Journal of Rock Mechanics and Mining Sciences,1998,35(2):113 -121.
[11]KAISER P K,TANG Chun-an.Numerical Simulation of Damage Accumulation and Seismic Energy Release During Brittle Rock Failure-PartⅡ:Rib Pillar Collapse[J].International Journal of Rock Mechanics and Mining Sciences,1998,35(2):123 -134.
[12]傅宇方,梁正召,唐春安.巖石介質細觀非均勻性對宏觀破裂過程的影響[J].巖土工程學報,2000,22(6):705 - 710.(FU Yu-fang,LIANG Zheng-zhao,TANG Chun-an.Numerical Simulation on Influence of Mesoscopic Heterogeneity on Macroscopic Behavior of Rock Failure[J].Chinese Journal of Geotechnical Engineering,2000,22(6):705 -710.(in Chinese))
[13]梁正召,楊天鴻,唐春安,等.非均勻性巖石破壞過程的三維損傷軟化模型與數值模擬[J].巖土工程學報,2005,27(12):1447 -1452.(LIANG Zheng-zhao,YANG Tian-hong,TANG Chun-an,et al. Three-dimensional Damage Soften Model for Failure Process of Heterogeneous Rocks and Associated Numerical Simulation[J].Chinese Journal of Geotechnical Engineering,2005,27(12):1447 -1451.(in Chinese))
[14]MA Zhi-tao,TAN Yun-liang,ZHANG Ting.Modeling of Rock Failure Based on Physical Vellular Automata[J].Journal of Southeast University(English Edition),2005,21(3):348-352.
[15]馮增朝,趙陽升,段康廉.巖石的細胞元特性及其非均質分布對巖石全曲線性態的影響[J].巖石力學與工程學報,2004,23(11):1819 -1823.(FENG Zeng-chao,ZHAO Yang-sheng,DUAN Kang-lian.Influence of Rock Cell Characteristics and Rock Inhomogeneity Parameter on Complete Curve of Stress-Strain[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(11):1819 -1823.(in Chinese))
[16]陳 沙,岳中琦,譚國煥.基于數字圖像的非均質巖土工程材料的數值分析方法[J].巖土工程學報,2005,27(8):956 - 964.(CHEN S,YUE Z Q,THAM L G.Digital Image Based Numerical Modeling Method for Heterogeneous Geomaterials[J].Chinese Journal of Geotechnical Engineering,2005,27(8):956 - 964.(in Chinese))
[17]陶振宇,潘別桐.巖石力學原理與方法[M].武漢:中國地質大學出版社,1991.(TAO Zhen-yu,PAN Bietong.Principle and Method of Rock Mechanics[M].Wuhan:China University of Geosciences Press,1991.(in Chinese ))