巖石礦物細胞元隨機性參數賦值方法研究

羅 榮，曾亞武

（武漢大學 土木建筑工程學院，武漢 430072）

1 引 言

巖石是經過多種地質作用形成的天然礦物集合體，在漫長的地質年代經歷了結晶過程的演化、高溫、高壓的影響以及各種地質營力的作用，形成了極其復雜的結構和構造特性，其中最基本的特性之一就是巖石的非均質性[1]，即巖石是由不同的礦物顆粒集合體和膠結材料組成的非均質體。大多數巖石是由幾種礦物組成的，巖石中的礦物成分及其含量，是決定巖石物理力學性質的主要因素，在某些條件下甚至會產生決定性影響[2]。

對于巖石類非均質材料，在外力作用或環境因素（如溫度）改變等條件下，材料內部的應力和變形分布是相當復雜的，主要表現就是其物理力學性質的高度非均勻性和非線性，而微觀介質（可以理解為礦物顆粒）參數的不均勻性是造成巖石力學性質非均勻和非線性的主要原因[3－4]。大量的巖石力學試驗結果也表明巖石試樣的物理力學性質參數存在較大的離散性，這種試驗結果的離散性并不完全是由試驗條件等外部因素引起的，巖石自身的非均質性是一個重要的影響因素[5]。巖石的非均質性主要是由構成巖石的礦物顆粒的物理力學性質的差異性決定的，礦物顆粒的力學性質差異越小巖石均質性越好；反之則巖石均質性越差。巖石力學的分析方法要成功地解決實際巖土工程中的問題，必須考慮巖石的非均質特性[6]。

事實上，國內外學者已經在如何考慮巖石的非均質特性方面開展了一系列工作，并取得了大量的研究成果。Song和 Kim[7]，Napier和 Dede[8]，Li等[9]應用格構模型分析了原巖的破壞過程，巖石材料的非均質性通過隨機指定格構單元的強度和形狀實現；唐春安[3,10－11]和 Kaiser[10－11]利用 Weibull分布函數描述巖石材料的非均質參數，對二維問題做了較詳細地研究；陳永強等[12]對三維非均質脆性材料的破壞進行了研究和數值模擬；王學濱[13]對含初始隨機材料缺陷的巖樣破壞受峰后脆性的影響進行了數值模擬；康健等[14－15]利用隨機分布參數研究了巖石熱破裂現象；許湘華等[16]在進行邊坡可靠度分析時利用節理面的隨機分布考慮了節理幾何參數的不確定性。總體來說，上述研究中對巖石非均質特性的處理方法主要為對網格模型中的各單元體賦予不同的參數。因為巖石力學參數是巖石工程設計分析的關鍵性參數，其取值決定了計算結果的準確性、客觀性和實用性，因此，準確地反映巖石的力學參數以及不均勻的分布規律，對于巖土工程設計、施工和分析評價都具有重要的意義[17]。因此，研究巖石非均質性從微觀參數的不均勻分布著手，研究方向無疑是正確的。但上述研究中，對材料非均質性的處理采用了統計分析加隨機賦值的方法，即先按一定的網格尺寸將材料作網格劃分，假定網格單元的材料力學參數服從某一統計分布規律（如正態分布或Weibull分布等），然后再由確定的分布規律產生一個離散的材料力學參數值序列，并按 Monte Carlo隨機方法投放到各網格單元中去。這類方法對網格單元力學參數的賦值是隨機的，沒有考慮巖石的結構和構造特征。

從力學角度來看，巖石的非均質性還表現為力學參數的空間變異性，這種變異并不是完全隨機的，而是具有隨機性和結構性的二重性。張征[18]、譚文輝[19]、胡小榮[20]等基于地質統計學方法的區域化變量理論，考慮巖體力學參數的空間變異二重性特征研究了巖石的非均質賦值方法，該種分析處理并不采用傳統的概率統計方法，而是通過對樣品數據進行結構分析并借助變異函數這一工具來反映參數變量所具有的二重性特征，但這種方法由于本身存在的整體相關性和平滑化效應問題，單元賦值的精度不高，還有待進一步研究[20－21]。

由此可知，在考慮巖石的非均質性時，除了考慮微觀力學參數的隨機性外，更重要的是必須考慮巖石組成種類及含量所產生的結構性。本文提出了一種新的描述巖石非均質特性的有限元參數賦值方法——巖石礦物細胞元隨機性參數賦值方法，即將每一個有限元網格視為一個組成巖石礦物細胞元，基于組成巖石的礦物種類和含量，定義單元類別判定區間，對各個礦物細胞元進行礦物類別判定，并進行參數賦值，完成賦值過程。該方法通過組成巖石的礦物種類和含量，描述巖石非均質參數的結構性；不同種類的礦物顆粒集合體在巖石中的分布是隨機的，也不失隨機性。

2 巖石礦物細胞元隨機性參數賦值方法

假定組成巖石的礦物顆粒在局部范圍內具有相同的物理力學性質，這樣根據有限元計算要求剖分所得的網格單元，都可以代表一個局部范圍的礦物顆粒，并將其稱之為礦物細胞元[4]，具有相應的物理力學性質。巖石礦物細胞元隨機性參數賦值方法就是要判定每一個礦物細胞元的礦物屬性，同一類別的礦物細胞元賦予相同的物理力學參數，其參數賦值流程如圖1所示，步驟如下：

圖1 礦物細胞元隨機性參數賦值流程圖Fig.1 Flow chart of random parameter assignment of rock mineral cell unit

（1）根據研究對象尺度確定合適的有限元網格剖分尺寸，并進行有限元網格剖分。如對于巖石材料類問題，可以考慮細觀尺度的有限元網格剖分；

（2）根據巖石的礦物含量或工程問題的非均質特征定義礦物細胞元類別判定區間；

（3）利用Monte Carlo方法確定的一組離散的隨機數，逐個進行礦物細胞元類別判定；

（4）根據各礦物細胞元類別屬性，賦予相應的物理力學參數。

通過上述過程，即可建立起描述巖石非均質性的有限元模型。

2.1 礦物細胞元類別區間定義

圖2 巖石礦物含量示意圖Fig.2 Contents of the rock minerals

令

定義區間 Ai=[ui-1,ui)（其中i=1，2，…m-1；當i=m，取右閉區間）為巖石礦物細胞元類別判定區間，即將[0，1]區間劃分為m個區間，區間 Ai即為第i類別礦物的細胞元判定區間。由于區間 Ai在[0，1]區間所占比例為 ui-ui-1=ni，ni又為礦物 i的實際含量，所以判定區間 Ai與礦物i的含量有一一對應關系。

2.2 礦物細胞元類別判定

利用Monte Carlo方法，產生一組隨機數序列，該系列隨機數的個數與礦物細胞元數相同。依次利用該組隨機數序列，根據定義的礦物細胞元類別判定區間對細胞元進行類別判定。

2.2.1 Monte Carlo方法

Monte Carlo方法的理論依據是概率論中的大數定理，是一種通過隨機變量的數字模擬和統計分析求取數學物理工程技術問題近似解的數值方法。其基本思想是：首先為所要處理的問題建立一個概率模型，使所求問題的解剛好是該模型的參數、特征量或相關量；然后，通過統計試驗產生該問題的統計抽樣樣本；最后分析這些樣本的特性，并以此作為原問題的近似解。

Monte Carlo方法的主要手段是用隨機數進行統計試驗，產生符合某種統計規律的隨機數是應用Monte Carlo方法的基礎。首先產生[0，1]區間均勻分布的隨機數，然后利用這些隨機數產生服從其他分布的隨機數或統計量。設隨機變量x服從概率密度函數為f(x)的分布，其概率分布函數為F(x)。

圖3、4分別為隨機變量x的概率分布函數曲線和概率密度函數曲線示意圖。首先產生一組服從[0，1]區間均勻分布的隨機數序列ξi，根據累計分布函數式，對于 F(xi)=ξi利用其反函數求出相對應的xi，對于xi，根據概率密度函數可計算 f(xi)。這樣，由均勻隨機數序列ξi映射為一組特殊序列xi，由大數定理可知，當n取得足夠大時，xi服從概率密度為f(x)的分布，即由[0，1]區間的均勻隨機數映射得到一組服從給定分布的隨機數。

圖3 概率分布函數曲線示意圖Fig.3 Curve of probability distribution function

圖4 概率密度函數曲線示意圖Fig.4 Curve of probability density distribution function

2.2.2 礦物細胞元類別判定

選用[0，1]的均勻分布作為隨機序列目標分布，概率分布函數為F(x)=x，概率密度函數為f(x)=1。根據Monte Carlo方法，對于初始隨機數ξk，令F(xk)=ξk，即得到用以進行礦物類別判定數序列xk=ξk。

根據定義的細胞元類別判定區間 Ai，依次判別xk的所屬區間，如圖5所示，若 xk∈Ai，則判定該單元為第i類別礦物細胞元，并對該細胞元賦予相應的i類別物理力學參數。依次對所有的礦物細胞元進行類別判定和賦值，即完成整個模型的賦值。

圖5 礦物細胞元類別判定示意圖Fig.5 Schematic diagram of sort judgment for rock mineral cells

由于礦物細胞元類別判定數xk是服從[0，1]均勻分布的隨機數序列，細胞元類別判定為第i類別礦物細胞元的概率 Pi=ui-ui-1=ni，又ni為第 i類別礦物的含量，故細胞元的類別判定概率與礦物的真實類別含量相等。

礦物細胞元混合模型通過真實地反映巖石的礦物含量來描述巖石的非均質性，且模型使各礦物細胞元處于隨機均勻混合狀態，因此，該模型進行巖石非均質參數賦值時既考慮了巖石參數的結構性特征，又考慮了礦物細胞元的隨機分布特征。

3 礦物細胞元分布的隨機性特征

本文提出的巖石礦物細胞元隨機性參數賦值方法的主要思想是利用組成巖石的礦物種類及其含量定義礦物細胞元類別判定區間。利用Monte Carlo方法，產生服從[0，1]均勻分布的隨機數序列作為類別判定數列，利用定義的判定區間依次對細胞元進行類別判定并相應賦值，使各類別礦物細胞元均勻隨機混合。

礦物細胞元的類別判定中采用的服從[0，1]均勻分布的隨機數序列決定了礦物細胞元分布具有隨機性。雖然細胞元的類別判定區間使得各類別細胞元含量在概率上與巖石的礦物含量相等，但對于確定的細胞元類別判定區間，隨機數取[0，1]區間任何一個數的概率都相等，因此，服從[0，1]均勻分布的隨機數序列不是惟一確定的，判定數序列的隨機性引起礦物細胞元的類別判定具有隨機性，從而導致礦物細胞元的分布具有隨機性特征。

對同一類別礦物細胞元組成的巖石有限元模型，僅有[0，1]一個細胞元類別判定區間，對任意服從[0，1]均勻分布的隨機數序列，所有的礦物細胞元都能被惟一地確定為同一類別，即等同為均質模型，所有的礦物細胞元具有相同的物理力學參數，此時礦物細胞元的分布不具有隨機性特征。

對由多種礦物細胞元組成的巖石有限元模型，各單元依次進行礦物細胞元類別判定并對不同類別的礦物細胞元賦予不同的參數值，此時，類別判定數序列作為礦物細胞元類別判定指標具有的隨機性特征將直接影響類別判定結果，使得描述巖石非均質性的有限元模型不具有惟一確定性，依賴于隨機數序列具有隨機性特征。如圖6、7分別為兩礦物細胞元賦值混合模型和三礦物細胞元賦值混合模型示意圖，相同的巖石礦物含量定義相同的判定區間，但不同的隨機數序列進行礦物細胞元類別判定后獲得的巖石非均質模型具有明顯的隨機性特征。

圖6 兩礦物細胞元賦值Fig.6 Assignment of two mineral cell units

圖7 三礦物細胞元賦值Fig.7 Assignment of three mineral cell units

4 巖石礦物細胞元隨機性參數賦值數值試驗

巖石的宏觀力學特征是由組成巖石的各類礦物細胞元共同作用的結果[2－3]。為研究巖石礦物細胞元分布的隨機性特征對巖石宏觀力學特征的影響，本文在彈性范圍內，對相同的巖石礦物種類和含量，利用不同的隨機數序列進行礦物細胞元類別判定并賦值后建立的有限元模型進行了數值試驗研究。

4.1 數值試驗模型

為簡化計算，在進行有限元數值試驗時，采用線彈性本構關系。考慮巖石由兩礦物組成或三礦物組成，各選擇6組不同物理參數或礦物含量的試件，對每組試件進行10次礦物細胞元類別判定，建立礦物細胞元混合模型，用以研究礦物細胞元分布的隨機性特征對巖石宏觀彈性模量的影響。

數值試驗采用平面應變模型，試件幾何尺寸為20 cm×10 cm，有限元單元（即巖石礦物細胞元）尺寸為0.5 mm×0.5 mm[4]，共劃分400×200個單元，數值試驗采用位移加載模式。由于各類巖石材料的泊松比一般為0.3左右，差異性不大，故不考慮泊松比的差異性，計算泊松比取v=0.3；巖石材料彈性模量為1～100 GPa不等，故數值試驗中，力學參數僅考慮彈性模量的差異性。對于兩礦物和三礦物混合模型，各組試件礦物細胞元彈性模量賦值參數E及礦物含量n取值分別見表1、2。

表1 兩礦物混合模型數值計算參數Table1 Parameters of two mineral cell units

4.2 數值試驗結果分析

對每一組試件，僅考慮礦物細胞元類別判定的隨機性影響，進行10次隨機類別判定及賦值，分別計算得到的等效宏觀彈性模量，如圖 8、9及表 3所示。

圖9 三礦物混合模型隨機性試驗結果圖Fig.9 Experimental results of random characteristic(three mineral cell units)

表3 兩礦物、三礦物混合模型等效彈性模量計算結果Table3 Statistical results of equivalent elastic modulus of two mineral cell units and three mineral cell units

結果表明：

（1）無論是兩礦物或三礦物組成的巖石試件，其等效宏觀彈性模量計算結果幾乎不受礦物細胞元隨機分布的影響。

（2）只要礦物種類和含量確定，建立的礦物細胞元混合模型的等效宏觀力學參數都非常穩定，各組試件數值試驗結果的標準差、變異系數和極差都非常小。

（3）不同類別礦物細胞元之間的彈性模量相差越小，礦物細胞元混合模型的等效宏觀彈性模量的標準差、變異系數和極差也越小。

（4）相同礦物含量的巖石礦物細胞元混合模型的等效宏觀彈性模量均值隨礦物細胞元彈性模量的提高而提高，說明巖石材料宏觀力學參數受礦物細胞元力學參數的影響。

4.3 數值試驗結果比較

本文采用文獻[4]所述的巖石非均質參數賦值方法，對巖石單元的彈性模量進行賦值，進行相關數值試驗，計算該方法等效宏觀彈性模量的穩定性，單元彈性模量的Weibull分布密度函數表達式為

式中：β為非均質參數；E0為所有細胞元的彈性模量平均值；φ為彈性模量的密度函數值。

為便于比較，E0取巖石礦物細胞元隨機性參數方法中兩礦物混合模型第1組試件的所有細胞元彈性模量平均值，即取E0=44 GPa；非均質參數β分別取3、6、9、12、15、18。計算結果見表4。

表4 Weibull分布非均質模型等效彈性模量計算結果Table4 Statistical results of equivalent elastic modulus of inhomogeneity model(Weibull distribution)

通過比較可知，采用本文方法計算所得非均質巖石的等效宏觀彈性模量值與利用 Weibull分布函數賦值方法（給定非均質參數β）的計算結果在精度和穩定性方面基本一致，說明采用本文方法描述巖石力學參數的非均勻分布是可以接受的。

但與 Weibull分布函數賦值方法相比，本文方法具有以下兩個明顯特征：

（1）考慮了巖石組成的結構特征，即礦物組成及其含量。只要巖石的礦物組成及其含量確定，就不需要其他的非均質參數即可獲得非均質巖石模型，從而避免了非均質參數的選擇對巖石力學性能造成的影響。

（2）本文方法可完全退化為均質模型，即該方法可同時用于非均質巖石和均質巖石的數值模擬。

5 結 論

（1）巖石礦物細胞元隨機性參數賦值方法考慮了巖石組成的結構特征，在賦值過程中引入巖石自身的礦物種類及含量建立判定函數，在不失隨機性的基礎上考慮了巖石參數賦值的結構性。

（2）采用簡單的[0，1]均勻分布隨機數來描述巖石礦物細胞元的隨機分布特點，不依賴于其他的隨機性參數，避免了采用隨機分布（如正態分布，Weibull分布等）時分布參數的選擇帶來的不確定性影響。

（3）利用Monte Carlo方法進行礦物細胞元類別判定，對單一種類的礦物細胞元構成的模型等同于均質模型，無隨機性特征；而對于多種礦物細胞元構成的混合模型，礦物細胞元的隨機分布對巖石宏觀力學參數的影響非常小，這與利用 Weibull分布函數進行參數賦值得出的結論一致。并且本文方法統一了均質巖石和非均質巖石的數值分析賦值過程，比傳統的采用隨機分布函數對巖石有限元單元參數進行賦值的方法具有優越性。

（4）組成巖石的各類礦物細胞元力學參數差別越小，巖石越均勻，宏觀力學參數的變異性也越小，與利用 Weibull分布函數進行參數賦值得出的結論一致，但本文方法具有自身鮮明的特點。

（5）巖石等效宏觀彈性模量隨礦物細胞元彈性模量的提高而提高，反映了礦物類型和含量對巖石力學性質的結構性影響。

本文通過數值試驗主要研究了巖石礦物細胞元隨機分布對巖石宏觀彈性模量的影響，至于巖石宏觀應力-應變關系、強度及變形特征、裂紋擴展等與礦物細胞元參數及礦物細胞元分布之間的關系還有待進一步研究。

[1]唐春安. 巖石破裂過程數值試驗[M]. 北京: 科學出版社,2003.

[2]陶振宇,潘別桐. 巖石力學原理與方法[M]. 武漢: 中國地質大學出版社,1991.

[3]唐春安. 巖石聲發射規律數值模擬初探[J]. 巖石力學與工程學報,1997,16(4): 368－374.TANG Chun-an. Numerical simulation of AE in rock failure[J]. Chinese Journal of Rock Mechanics and Engineering,1997,16(4): 368－374.

[4]馮增朝,趙陽升,段康廉. 巖石的細胞元特性及其非均質分布對巖石全曲線性態的影響[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 stressstrain[J]. Chinese Journal of Rock Mechanics and Engineering,2004,23(11): 1819－1823.

[5]張占榮,盛謙,楊艷霜,等. 基于現場試驗的巖體變形模量尺寸效應研究[J]. 巖土力學,2010,31(9): 2875－2881.ZHANG Zhan-rong,SHENG Qian,YANG Yan-shuang,et al.Study of size effect of rock mass deformation modulus based on in-situ test[J]. Rock and Soil Mechanics,2010,31(9): 2875－2881.

[6]張淵,趙陽升. 巖石非均質度與熱破裂的相關性分析[J].蘭州理工大學學報,2009,35(6): 135－137.ZHANG Yuan,ZHAO Yang-sheng. Analysis of correlation of rock thermal cracking with inhomogeneity[J]. Journal of Lanzhou University of Technology,2009,35(6): 135－137.

[7]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.

[8]NAPIER J A L,DEDE T. A comparison between random mesh schemes and explicit growth rules for rock fracture simulation[J]. International Journal of Rock Mechanics and Mining Sciences,1997,34(3－4): 217.

[9]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.

[10]TANG Chun-an,KAISER K P. 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 K P,TANG Chun-an. Numerical simulation of damage accumulation and seismic energy release during brittle rock failure-Part II: Rib pillar collapse[J].International Journal of Rock Mechanics and Mining Sciences,1998,35(2): 123－134.

[12]陳永強,鄭小平,姚振漢. 三維非均勻脆性材料破壞過程的數值模擬[J]. 力學學報,2002,34(3): 351－361.CHEN Yong-qiang,ZHENG Xiao-ping,YAO Zhen-han.Numerical simulation of failure processes in 3-D heterogeneous brittle material[J]. Acta Mechanica Sinica,2002,34(3): 351－361.

[13]王學濱. 峰后脆性對非均質巖石試樣破壞及全部變形的影響[J]. 中南大學學報(自然科學版),2008,39(5):2394－2399.WANG Xue-bin. Effects of post-peak brittleness on failure and overall deformational characteristics of rock specimen with random material imperfections[J]. Journal of Central South University(Science and Technology),2008,39(5): 2394－2399.

[14]康健,趙陽升,趙崢嶸,等. 隨機介質熱彈性平面軸對稱問題的解析解[J]. 巖土力學,2006,27(10): 1689－1692.KANG Jian,ZHAO Yang-sheng,ZHAO Zheng-rong,et al.Analytical solutions to thermoelastic plane axisymmetric problems for random non-homogeneous media[J]. Rock and Soil Mechanics,2006,27(10): 1689－1692

[15]康健,畢秀國,劉超. 非均質隨機分布對巖石熱破裂影響的數值試驗[J]. 巖土力學,2008,29(增刊1): 491－494.KANG Jian,BI Xiu-guo,LIU Chao. Numerical tests of influence of heterogeneous random distribution on thermal cracking of rock[J]. Rock and Soil Mechanics,2008,29(Supp.1): 491－494.

[16]許湘華,曲廣琇,方理剛. 基于節理幾何參數不確定性的邊坡可靠度分析[J]. 中南大學學報(自然科學版),2010,41(3): 1133－1139.XU Xiang-hua,QU Guang-xiu,FANG Li-gang.Reliability analysis of rock slope based on uncertainty of joint geometric parameters[J]. Journal of Central South University(Science and Technology),2010,41(3): 1133－1139.

[17]李世海,董大鵬,燕琳. 含節理巖塊單軸受壓試驗三維離散元數值模擬[J]. 巖土力學,2003,24(4): 648－652.LI Shi-hai,DONG Da-peng,YAN Lin. 3D-DEM numerical simulation for jointed reck under uniaxial press loading[J]. Rock and Soil Mechanics,2003,24(4): 648－652.

[18]張征,劉淑春,鞠碩華. 巖土參數空間變異性分析原理與最優估計模型[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.

[19]譚文輝,王家臣,周汝弟. 巖體強度參數空間變異性分析[J]. 巖石力學與工程學報,1999,18(5): 497－502.TAN Wen-hui,WANG Jia-chen,ZHOU Ru-di. Analysis of spatial variability of strength parameter of rock mass[J].Chinese Journal of Rock Mechanics and Engineering,1999,18(5): 497－502.

[20]胡小榮,唐春安. 巖土力學參數隨機場的空間變異性分析及單元體力學參數賦值研究[J]. 巖石力學與工程學報,2000,19(1): 59－63.HU Xiao-rong,TANG Chun-an. Spatial variation analysis of the random field of mechanical parameters for rock and soil and the parameter estimation of elements[J]. Chinese Journal of Rock Mechanics and Engineering,2000,19(1): 59－63.

[21]趙紅亮,馮夏庭,張東曉,等. 巖土力學參數空間變異性的集合卡爾曼濾波估值[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.