考慮空間相關尺度特征的細觀力學模型及其應用

唐欣薇，周元德，張楚漢

（1.華南理工大學 亞熱帶建筑科學國家重點實驗室，廣州 510640；2.清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084）

1 引 言

巖石類材料內部存在著大量的微、細觀缺陷，是一種典型非均質材料。隨著試驗技術的發展以及現代數值分析技術和計算機技術的進步，學者們開始立足于材料的微、細觀結構，研究巖石的破損機制，提出了許多理論與數值模型，如格構模型[1]、梁-顆粒模型[2]、巖石破裂過程分析模型（RFPA）[3]、顆粒流模型（PFC）[4]和細胞自動機模型（EPCA）[5]等。但在以往的細觀分析模型中，普遍定義離散單元材料的力學屬性為隨機分布，而未考慮實際巖石材料成巖過程的影響，其各項物理力學指標的隨機分布應具有一定的空間相關尺度特征[6]，本文擬就這一個問題展開探討研究。

本文基于連續介質損傷力學分析細觀尺度下巖石的變形、損傷和破壞過程，假定該尺度下拉致開裂是導致裂紋發生與擴展的主要原因，采用彌散式損傷本構模型模擬各離散單元的力學響應。在所建模型中基于 Weibull統計分布函數描述材料的物理力學性質，假定材料的力學參數（如彈性模量、強度等）隨空間位置的不同而發生變化，以描述巖石材料的非均質特性。為了反映天然巖石材料自身的連續性，本文還對上述非均質隨機分布（Weibull分布）模型進行改進，通過引入空間相關尺度因子使得計算模型與天然巖石的非均值特征更加符合，并進一步研究了空間相關尺度因子對材料力學性能以及變形破壞過程的影響。

2 巖石細觀本構模型

本文模型基于連續介質損傷力學框架，描述巖石材料在細觀尺度下各離散單元的力學特性。按照Lemaitre應變等價原理，受損材料的本構關系可通過無損材料的名義應力-應變關系表示，即

式中：E0為初始彈性模量；為損傷后的彈性模量；D為損傷變量；σ為應力；ε為應變。

為簡化計算，本文采用雙折線本構關系。計算模型中采用最大拉應力準則，即當細觀單元的最大主拉應力達到給定的極限時，該單元開始發生拉伸損傷，計算單元滿足斷裂能守恒準則。

本文參考 Fenves模型中關于損傷定義由一維擴展至多維的方法[7]，將上述單軸應力狀態下應力-應變關系推廣至平面及三維應力狀態。

有關細觀本構模型的具體介紹可參見相關研究成果[8]。

3 考慮空間相關尺度特征的非均質隨機分布模型

3.1 Weibull統計分布的基本理論

在常規非均質隨機分布模型中，認為巖石材料的力學參數分布具有一定的隨機性，通常采用的一類假定是認為組成巖石材料單元的力學性質滿足Weibull概率統計分布，該分布依照如下分布密度函數來定義：

式中：u為滿足該分布參數（如強度、彈性模量、泊松比等）的數值；m為材料均質度；u0為一個與均值E(u)有關的參數。參數m定義了Weibull分布密度函數的形狀。如u0=30，均質度m分別取1.0、1.5、4.0和9.0時，Weibull分布密度函數曲線如圖1所示。可見，m值反映了統計模型中材料結構的均勻程度，當其由小到大變化時，材料細觀單元強度分布函數由扁而寬向高而窄變化，細觀單元強度趨于均勻，且材料內部所包含的大部分細觀單元接近給定的參數u0[9]。

圖1 不同均質度的分布密度函數曲線Fig.1 Curves of distribution density function of different homogeneity indices

此外，由式（3）所定義的均質度涵蓋了兩方面特征：①不同m值的Weibull分布參數代表均質度不同的非均勻材料；②對于相同的 Weibull分布參數m和u0，每次隨機產生的樣本，其材料參數的空間分布也是不同的。

3.2 考慮空間相關尺度特征的Weibull隨機分布模型

對于實際巖石材料，受天然成巖過程的影響，其物理力學性能指標并不是完全雜亂無章、隨機分布的，而是在局部范圍內存在一定的空間相關性，即材料單元的力學性能在空間分布上有一定的連續性和相關性。如圖2所示為天然巖石切口斷面的微缺陷具有局部的空間連續性就是這一相關性的例子。該分布特征在計算模型中可以通過結合隨機場和空間關聯函數來描述。在以往基于細觀力學對混凝土、巖石等材料的研究中，普遍僅給出了描述材料力學參數隨機分布的隨機場，并未對材料連續性對其力學性能的影響進行深入討論。本節引入空間相關尺度因子用以表征材料的空間局部連續程度，并對不同空間相關尺度因子對材料力學性能的影響進行了比較分析。

3.2.1 隨機場相關尺度的引入

為了描述隨機場內部各點之間的連續性，在常規的 Weibull隨機分布模型中引入空間關聯函數，用以表征材料的連續程度。以圖3為例，將二維隨機場Ω離散為n個有限單元，單元i的力學參數定義為

式中： fi、 fk為單元i、k的力學參數，其初值基于Monte Carlo法生成；ρi,k為Markov相關函數；di,k為單元i與單元k形心點之間的距離；為單元i、k的位置向量；Θ為空間相關尺度因子（spatial correlation length），用以描述空間兩點之間隨機變量的相關程度，Θ越大，隨機場變化越平滑；反之，變化越激烈[8,10]。

圖2 天然巖石切口[11]（灰度代表參量值）Fig.2 Cross-section of natural rock[11]

圖3 二維連續的隨機場Fig.3 Two-dimensional continuous random field

3.2.2 相關尺度的數值實現

如圖4所示，取一個隨機樣本，該樣本包括N個單元，首先生成符合給定 Weibull概率分布的初始隨機場。然后，根據研究對象的特性，引入空間相關尺度因子Θ，按照式（4）、（5）更新初始隨機場中各單元的力學參數，從而獲得一個力學參數滿足給定關聯性的隨機場。由于單元的力學參數已被更新，此時的分布并不符合給定的 Weibull統計形態，僅僅確定了單元力學參數的相關性，因而，還需對其進行修正。

將該隨機場中單元的力學參數從小到大排序，按順序記錄相應的單元號，即e1、e2，…，ei，…，en1+n2+…+nk。根據給定Weibull分布的形態，取前n1個單元e1～，將其力學參數映射至（u1，u2）。同理實現整個隨機場的修正，使其在保證各單元空間關聯性的同時，滿足既定的Weibull統計分布。

圖4 考慮空間相關尺度因子的隨機場實現過程Fig.4 Implementation process for random field with spatial correlation scale factor

4 算例分析

選取 Van Mier和 Shi[12]的巖石單軸拉伸試驗（見圖5），試件尺寸為90 mm×45 mm×10 mm，切口尺寸為2 mm×5 mm×10 mm，試件選自德國西部愷撒斯勞騰的黃砂巖，抗拉強度約為1.50～1.95 MPa，彈性模量為8～10 GPa，泊松比為0.2。

本文將巖石試件離散為1 mm×1 mm的平面有限元網格（見圖6），考慮材料力學性質的非均質性，假定各力學指標符合 Weibull分布，相關參數按試驗結果均值取為：抗拉強度ft=1.8 MPa，彈性模量E=8.8 GPa，泊松比ν=0.2，斷裂能Gf=60 N/m。基于上述改進的隨機分布模型，分別定義不同的空間相關尺度因子，對生成的細觀有限元離散模型進行重新定義，實現考慮巖石材料可能的不同空間相關尺度特征。限于篇幅，僅考慮均值度m = 8的情形，隨機生成一個樣本，然后引入空間相關函數，分別取Θ=2.0、1.0、0.5、0.2共4組，同時保證每一個樣本符合同一組Weibull分布，每組生成10個隨機樣本，其中，選取3個隨機樣本分布形態如圖7所示。圖中不同的灰度水平代表所處單元的強度值高低，灰度愈高，材料強度也越高。顯然，不同空間相關尺度因子Θ獲得的材料性能分布樣本具有顯著的差異。

基于生成的細觀隨機有限元模型，進行了系列數值仿真分析。計算得出的代表性結果，包括試件破壞形態特征、荷載-切口張開度關系曲線以及峰值荷載的變化比較等如圖8、9及表1所示。分析計算結果表明：

（1）隨著空間相關尺度因子的減小，試件力學性能的離散性增加，即各樣本的標準差（變異系數）增大。由于空間相關尺度因子考慮了隨機場的關聯性、連續性，當相關尺度較小時，隨機場各點力學性能的梯度較大，力學性能隨空間變化劇烈。反之，當空間相關尺度因子較大時，隨機場各點力學性能的梯度較小，隨機場越趨于平穩變化，因而在宏觀上呈現較小的離散性。

（2）圖8中以灰度水平代表損傷變量D值，灰度越高，代表D值愈大。當空間相關尺度因子Θ=2.0時，宏觀表現為試件均沿著切口開裂，裂紋擴展形態基本一致。隨著Θ的不斷減小，空間相關性下降，試件開裂的形態開始表現出更強的離散性，如圖8(d)所示，試件出現不同的裂紋擴展形態。

圖5 巖石單軸拉伸試驗示意圖[12]Fig.5 Sketch of uniaxial tensile test on rock[12]

圖6 計算采用的有限元網格Fig.6 Finite element meshes for numerical simulation

圖7 不同空間相關尺度因子的巖石樣本Fig.7 Numerical samples with different spatial correlation scale factors

圖9 考慮不同空間相關尺度因子的數值樣本荷載-切口張開度曲線Fig.9 Load-cut opening degreee curves of samples with different spatial correlation scale factors

（3）由圖9結果可知，引入空間相關尺度因子后，隨著空間相關尺度因子的變化，有可能得到一組數值仿真與試驗結果具有良好一致性的成果，圖中取Θ=0.2工況與試驗結果符合較好。還可以看出，在 0.2～2.0范圍中變化Θ值，對荷載-切口張開度曲線的上升段和峰值荷載影響很小，而對峰后軟化段的離散分布范圍有著較明顯的影響。

（4）由于各個數值樣本中細觀單元的力學分布符合同一種 Weibull統計分布，因而引入空間相關尺度因子后，各組樣本力學性能的均值變化并不大，如表1所列。

表1 考慮不同空間相關尺度因子的試件承載力數值仿真結果Table1 Numerical results of load bearing capacity of samples with different spatial correlation length factors

5 結 論

（1）在常規隨機分布模型基礎上，本文對常規Weibull隨機分布模型加以改進，引入空間相關尺度因子，從而實現考慮了關注材料在空間各點之間的關聯性，以表征材料的空間局部連續程度。算例結果表明，對于滿足相同統計分布的隨機場，變化空間相關尺度因子對試樣力學響應的均值影響并不顯著。

（2）變化空間相關尺度因子對其破壞形態特征以及荷載位移響應的離散性有著較明顯的影響。應當指出，本文只對單軸抗拉巖石試驗情形進行了系列分析探討，對于其他受載情形仍有待進一步深入研究。

[1]SCHLANGEN E,VAN MIER J G M. Simple lattice model for numerical simulation of fracture of concrete material and structures[J]. Materials and Structures,1992,25(9): 534－542.

[2]邢紀波,俞良群,王泳嘉. 砂巖類脆性無序介質連續破壞過程的細觀模擬[J]. 地質力學學報,1998,4(3): 28－35.XING Ji-bo,YU Liang-qun,WANG Yong-jia. Modelling of progressive failure of brittle disordered materials e. g.sandstone[J]. Journal of Geomechanics,1998,4(3): 28－35.

[3]TANG C A,LIU H,LEE P K K,et al. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—Part 2: Effect of heterogeneity[J].International Journal of Rock Mechanics and Mining Sciences,2000,37(4): 555－569.

[4]Itasca Consulting Group,Inc. PFC3D user’s guide[M].Mineapolis: Itasca Consulting Group Inc.,2000.

[5]潘鵬志,馮夏庭,周輝. 脆性巖石破裂演化過程的三維細胞自動機模擬. 巖土力學,2009,30 (5): 1471－1476.PAN Peng-zhi,FENG Xia-ting,ZHOU Hui. Failure evolution processes of brittle rocks using 3D cellular automaton method[J]. Rock and Soil Mechanics,2009,30(5): 1471－1476.

[6]FENTON G A,GRIFFITHS D V. Bearing-capacity prediction of spatially random c-φ soils[J]. Canadian Geotechnical Journal,2003,40(1): 54－65.

[7]LEE J,FENVES G L. Plastic-damage model for cyclic loading of concrete structures[J]. Journal of Engineering Mechanics,1998,124(8): 892－900.

[8]唐欣薇,張楚漢. 混凝土細觀力學模型研究: 非均質影響[J]. 水力發電學報,2009,28(4): 56－62.TANG Xin-wei,ZHANG Chu-han. Meso-scale modeling of concrete: Effects of heterogeneity[J]. Journal of Hydroelectric Engineering,2009,28(4): 56－62.

[9]唐春安,朱萬成. 混凝土損傷與斷裂數值試驗[M]. 北京: 科學出版社,2003.

[10]GRIFFITHS D V,GORDON F A. Probabilistic slope stability analysis by finite elements[J]. Journal of Geotechnical and Geoenvironmental Engineering,2004,130(5): 507－518.

[11]CHO N,MARTIN C D,SEGO D C. A clumped particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences,2007,44(7): 997－1010.

[12]VAN MIER J G M,SHI C. Stability issues in uniaxial tensile tests on brittle disordered materials[J].International Journal of Solids and Structures,2002,39(2): 3359－3372.