考虑空间相关尺度特征的细观力学模型及其应用

唐欣薇，周元德，张楚汉

（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.