唐欣薇，周元德，張楚漢

（1.華南理工大學(xué) 亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，廣州 510640；2.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100084）

1 引 言

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

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

2 巖石細(xì)觀本構(gòu)模型

本文模型基于連續(xù)介質(zhì)損傷力學(xué)框架，描述巖石材料在細(xì)觀尺度下各離散單元的力學(xué)特性。按照Lemaitre應(yīng)變等價原理，受損材料的本構(gòu)關(guān)系可通過無損材料的名義應(yīng)力-應(yīng)變關(guān)系表示，即

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

為簡化計算，本文采用雙折線本構(gòu)關(guān)系。計算模型中采用最大拉應(yīng)力準(zhǔn)則，即當(dāng)細(xì)觀單元的最大主拉應(yīng)力達(dá)到給定的極限時，該單元開始發(fā)生拉伸損傷，計算單元滿足斷裂能守恒準(zhǔn)則。

本文參考 Fenves模型中關(guān)于損傷定義由一維擴(kuò)展至多維的方法[7]，將上述單軸應(yīng)力狀態(tài)下應(yīng)力-應(yīng)變關(guān)系推廣至平面及三維應(yīng)力狀態(tài)。

有關(guān)細(xì)觀本構(gòu)模型的具體介紹可參見相關(guān)研究成果[8]。

3 考慮空間相關(guān)尺度特征的非均質(zhì)隨機(jī)分布模型

3.1 Weibull統(tǒng)計分布的基本理論

在常規(guī)非均質(zhì)隨機(jī)分布模型中，認(rèn)為巖石材料的力學(xué)參數(shù)分布具有一定的隨機(jī)性，通常采用的一類假定是認(rèn)為組成巖石材料單元的力學(xué)性質(zhì)滿足Weibull概率統(tǒng)計分布，該分布依照如下分布密度函數(shù)來定義：

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

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

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

3.2 考慮空間相關(guān)尺度特征的Weibull隨機(jī)分布模型

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

3.2.1 隨機(jī)場相關(guān)尺度的引入

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

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

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

圖3 二維連續(xù)的隨機(jī)場Fig.3 Two-dimensional continuous random field

3.2.2 相關(guān)尺度的數(shù)值實(shí)現(xiàn)

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

將該隨機(jī)場中單元的力學(xué)參數(shù)從小到大排序，按順序記錄相應(yīng)的單元號，即e1、e2，…，ei，…，en1+n2+…+nk。根據(jù)給定Weibull分布的形態(tài)，取前n1個單元e1～，將其力學(xué)參數(shù)映射至（u1，u2）。同理實(shí)現(xiàn)整個隨機(jī)場的修正，使其在保證各單元空間關(guān)聯(lián)性的同時，滿足既定的Weibull統(tǒng)計分布。

圖4 考慮空間相關(guān)尺度因子的隨機(jī)場實(shí)現(xiàn)過程Fig.4 Implementation process for random field with spatial correlation scale factor

4 算例分析

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

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

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

（1）隨著空間相關(guān)尺度因子的減小，試件力學(xué)性能的離散性增加，即各樣本的標(biāo)準(zhǔn)差（變異系數(shù)）增大。由于空間相關(guān)尺度因子考慮了隨機(jī)場的關(guān)聯(lián)性、連續(xù)性，當(dāng)相關(guān)尺度較小時，隨機(jī)場各點(diǎn)力學(xué)性能的梯度較大，力學(xué)性能隨空間變化劇烈。反之，當(dāng)空間相關(guān)尺度因子較大時，隨機(jī)場各點(diǎn)力學(xué)性能的梯度較小，隨機(jī)場越趨于平穩(wěn)變化，因而在宏觀上呈現(xiàn)較小的離散性。

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

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

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

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

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

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

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

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

5 結(jié) 論

（1）在常規(guī)隨機(jī)分布模型基礎(chǔ)上，本文對常規(guī)Weibull隨機(jī)分布模型加以改進(jìn)，引入空間相關(guān)尺度因子，從而實(shí)現(xiàn)考慮了關(guān)注材料在空間各點(diǎn)之間的關(guān)聯(lián)性，以表征材料的空間局部連續(xù)程度。算例結(jié)果表明，對于滿足相同統(tǒng)計分布的隨機(jī)場，變化空間相關(guān)尺度因子對試樣力學(xué)響應(yīng)的均值影響并不顯著。

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

[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ì)波,俞良群,王泳嘉. 砂巖類脆性無序介質(zhì)連續(xù)破壞過程的細(xì)觀模擬[J]. 地質(zhì)力學(xué)學(xué)報,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]潘鵬志,馮夏庭,周輝. 脆性巖石破裂演化過程的三維細(xì)胞自動機(jī)模擬. 巖土力學(xué),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]唐欣薇,張楚漢. 混凝土細(xì)觀力學(xué)模型研究: 非均質(zhì)影響[J]. 水力發(fā)電學(xué)報,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]唐春安,朱萬成. 混凝土損傷與斷裂數(shù)值試驗(yàn)[M]. 北京: 科學(xué)出版社,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.