李世海 汪远年
（中国科学院力学研究所工程科学部 北京 100080）
摘要 针对土石混合体提出了一种随机计算模型。 在该模型中，可以模拟土石混合比、块石尺寸、块石形状等。给出了随机模型的实现方法，同时对该随机模型的算法可靠性进行了验证。通过统计分析的方法研究了单轴受压情况下土石混合体内部应力场分布与土石配比、岩石块度大小等因素的关系。并进一步给出了土石混合体的应力-应变特性和强度特性。
关键词 随机模型 土石混合体 统计分析 离散元
岩土工程中经常遇到非连续、不均匀介质，诸如由土石混合体组成的古滑坡体，由碎石堆积的面板堆石坝等。它是介于土体与岩体之间的一种特殊地质体，目前人们对于它的研究还处于探索之中。这类介质的主要特点是：①由松散的块体堆积而成；②块体的形状大小不均匀；③块体在空间随机分布；④块体之间经常有充填物。一般来说，块石和充填土的力学特性可以通过实验分别获得，而给出松散堆积体的力学特性比较困难，主要原因是：取样困难、难以进行小尺度的模型实验、大尺度的模型实验造价较高、实验离散度大难以掌握其规律性。因此通过数值模拟的方法研究随机块体构成的“岩土”混合体的力学特性是很有意义的。
D.R.Axelrad［1］对离散介质的统计模型及随机演化与可靠性方面做过比较系统的论述，综述了分子动力学模型、栅格模型、渗透模型三种离散介质统计模型，这三种模型皆是针对材料的微观结构提出的，并从能量的角度、运用拓扑几何以及随机过程的方法给出了比较详尽的理论分析。这类模型运用到土石混合体这类更加复杂和不确定的介质还有一定的困难。油新华等［2］对土石混合体做过大型野外水平推剪试验，在土石混合体的变形特点和抗剪强度方面有一定的认识，可以对数值模拟做一些参考。
本文在对现有3D-NURBM面-面接触弹性模型［3］块体单元划分的基础上引入了随机分布特征参数：块体密度、弹模、刚度；节理刚度、粘聚力、最大内摩擦角等。用统计分析的方法对单轴受压情况下土石混合体内部应力场分布进行了研究，分析结果表明土石混合比和块石尺寸是影响土石混合体内部应力场分布的主要因素。
1 三维离散单元法力学模型及计算方法
本文讨论的三维离散元力学模型是基于NURBM（Northwestern University Rigid Block Model）［3］，结合一种新的可变形块体模型［5］，基本假设如下：①作为单一介质的岩石块体或充填土体的力学特性是已知的；②每一个块石或充填土体都是由许多具有相同力学特性的小的岩石或土体块体单元组成的；③节理的力学参数分为3类：岩石块体单元之间，土体块体单元之间以及岩石块体单元和土体块体单元之间；④在低应力状态下，土体可看作弹性体；⑤每一块体单元的变形可通过应力状态和本构关系获得。
2 随机模型的实现方法
为了实现随机模型引入了4类参数：①土石混合比；②块体尺寸；③块体形状参数；④一定形状块体在空间的方位。
计算方法如下。
2.1 土石混合比的空间分布
现有3 D-NURBM块体单元划分是基于三组节理面切割而形成的许多相互接触的平行六面体单元。给定混合比，对土石块体单元进行分布的方法如下所述。
定义每一个块体单元为0到1之间的一个值。设F（x）＝x为0到1之间的均匀分布函数。如果岩块与土的混合比是ω，分离量A可定义为ω＝A/（1 -A），从而，A＝ω/（1＋ω）。对于每一块体单元可以赋值0 到1 之间的随机变量值x，如果F（x）小于A，块体单元标识为岩石，否则为土体单元。当标识完所有的块体单元后，岩石和土体块体单元的分配工作完成，并且整个研究区域的混合比接近1∶4。
2.2 块体尺寸大小选取
现场岩石块度的大小不均匀，形状也不相同。为了使岩石的尺寸接近真正的岩石块度大小，需要建立块体尺寸模型。
将整个研究区域划分成若干个子区域，子区域的三维尺寸与基于随机意义上的岩石块度的最大、最小长度相一致。如果子区域在三维方向上的层数分别为Li，Lj，Lk，那么每一子区域的块体单元总数（MT）为
2.3 块体形状参数
在每个子区域中，相同的MR有很多种集合方式。考虑到子区域单元的总数和石块所含单元数的个数，可以得到块体的形状参数。假定岩石块体在3个不同方向上的最大层数即形状参数为N1，N2，N3（分别不大于子区域参数Li，Lj，Lk），这些参数均是随机产生的，随机生成方法如下：
土石混合体
［］表示当［］内的数为整数时取其本身，否则取整后再加1；getstoc（m，n）是用来产生m到n之间的整数的随机生成函数。
2.4 一定形状块体在空间的方位
实际情况中，岩体的形状是千变万化的。为了更确切地描述岩体，随机选取平行六面体子区域的一个角点作为起始点进行集中排列岩石单元。当按照形状参数N1，N2，N3集中排列完所有岩石单元MR后，该子区域岩石单元的随机聚集排列工作就完成了。除了岩石单元，子区域内剩下的块体单元都被标识为土体单元。当循环完所有的子区域，整个研究区域的随机集中排列就结束了。有些情况下，会出现子区域内并非充满着实际的块体单元，即在某些研究区域边界的地方会出现部分子区域不完整现象，给出子区域的实际单元总数（MA）用来判断子区域是否完整。
2.5 块体单元和接触单元物理参数赋值
在三维离散元计算中，研究对象是以块体和接触两部分特理单元形式出现的，对应不同的块体和接触，有相应的物理参数：块体几何参数、密度、质量、转动惯量、弹模、泊松比、刚度；接触刚度、内聚力、最大内摩擦角等。
对应同一种介质这些参数都是固定的，而在不同介质情况下，相应不同介质块体单元的物理参数与不同接触单元（比如土-土、土-石、石-石接触）的物理参数均不同。
3 土石混合体空间分布模拟结果
3.1 混合比随机分布结果
验证按一定混合比随机生成的岩石块体单元数，算例中研究区域为20m×20m×20m、正交节理且三组节理间距皆为1m，混合比为1∶4，研究区域块体单元总数（NT）为8000。由表1可看出不同随机过程生成的岩石块体单元数基本符合预先给定的混合比要求。
表1 混合比随机分布计算结果（NT＝8000）
3.2 块石空间形态随机分布
图1是对4m×4m×4m的一个子区域、混合比为1∶4的块体单元不同次随机分布集中后的块体单元形心空间形态示意图。图2是对40m×60m×80m的研究区域、三组正交节理（节理间距均为1m）、4m×4m×4m的子区域、混合比为1∶4的块体单元随机分布集中后在x＝20.5m断面上的岩石块体单元型心分布情况。可以看出不同次随机过程产生的空间分布随机性比较好，对于土石混合体，可较接近真实地反映现场地质条件。
图1 子区域岩石块体单元随机集中空间形态示意图
图2 某一断面上岩石块体单元的分布情况（40%岩石）
4 土石混合体单向加载数值模拟试验
为研究土石混合体的力学特性，采用了单向加载数值模拟试验（轴向均布荷载：0.7MPa，侧向自由边界）。
4.1 计算参数选取
4.2 不同随机生成情况下的块体应力统计分析
参加的统计量均为各块体型心上的物理量，分别给出了各统计物理量的平均值和标准方差，并且有
表2 块体单元物理参数
表3 接触单元物理参数
算例中研究区域为20m×20m×20m、正交节理且三组节理间距皆为1m、子区域为4m×4m×4m、混合比为1∶4。
表4给出了相同混合比不同随机过程下的块体单元最大主应力和最大剪应力的统计分析结果，其中包括均值、最大最小值以及标准差。总共进行了9次随机模拟，结果显示均值、最大最小值、标准差的差别分别在1%，10%，5%左右，统计的最大值通常是均值的10倍左右。可以认为不同随机过程低于15%的差别是可以接受的。
表4 不同随机过程下块体单元应力统计分析结果（加载0.7MPa）
更为重要的是知道有多大区域的应力值高于强度值，也就是说在离散元计算中有多少块体单元的应力值超出强度范围。图3给出了在同一块度大小、不同随机过程（a）和不同混合比（b）情况下的块体单元最大剪应力的统计分布图，其中水平轴代表块体单元的最大剪应力值，竖直轴代表处在一定应力值下的块体单元数占总块体单元数的百分比。由图中我们可以明显看出：不同随机过程下的块体单元最大剪应力的统计分布差异很小；而不同混合比情况下则差异比较明显，当100%岩石时，即为均匀介质，块体单元的最大剪应力值基本上都稳定在0.4 MPa左右。
图3 不同随机过程（a）和不同混合比（b）情况下的块体单元最大剪应力τmax的统计分布图（相同块度大小）
由计算输出结果和统计分析可以看出，不同随机过程下的结果尽管有一些波动但基本上稳定在某一特定范围内。由于不同次随机过程块体单元的局部分布略有不同导致局部应力集中的差异，所以统计量中的最大值和最小值的最大误差（Max-error）相对于其他统计量显得更大。
4.3 不同混合比情况下的块体单元应力统计分析
算例中研究区域为20m×20m×20m、正交节理且三组节理间距皆为1m、子区域为4m×4m×4m。
统计结果显示（表5）：岩石块体含量在40%左右时应力空间分布的离散性最大，而在接近完全岩石或土体时离散性最小，可以忽略。当改变岩石块度大小时，结果有可能变化，但是我们可以确定混合比对应力的分布有很大的影响。
表5 不同配比情况下的块体单元应力统计分析（加载0.7MPa）
4.4 不同块度大小情况下的内部应力场分布
下面的算例是针对不同岩石块度大小的：研究区域为30m×30m×30m、正交节理且三组节理间距皆为1m、混合比1∶4。
图4给出了当块体尺寸分别为5m和15m（即子区域大小）时最大剪应力的分布。由于土体和岩石有着不同的强度，最大剪应力并不能直观反映混合体内部某一点处的破坏状态。在此有必要引入一个无量纲量——破坏度l，类似于安全系数：
图4 不同块度大小情况下断面x＝15.5m 上最大剪应力（τmax）等值线图
土石混合体
式中：c、φ分别为内聚力和最大内摩擦角。显然相同应力状态下对于岩石和土体来说破坏度是不同的。根据三维莫尔库仑准则，对于空间一点，当l＞1时处于稳定状态（不滑动），而当l＜1时块体会在节理面上产生滑动。
取岩体的内聚力和最大内摩擦角分别为2MPa，36°；土体的内聚力和最大内摩擦角分别为0.02MPa，27°。
图5给出了在不同块度下，断面x＝15.5 m上破坏度的等值线图，可以明显地分辨出哪些区域是破坏危险区域，哪些区域是相对稳定区域。图6给出了相应于图4，图5算例的断面x＝15.5 m上岩石块体单元的型心分布图，由结果我们可以看出岩石块体单元集中的地方往往是相对高应力状态区，在岩石和土体交界面的地方容易发生破坏。由于计算中没有考虑强度准则，即使大部分区域达到了破坏状态，但是，材料仍然没有失稳。
图5 不同块度大小情况下断面x＝15.5m 上破坏度l 的等值线图
图6 不同块度大小情况下断面x＝15.5m 上岩石块体单元的形心分布图
4.5 考虑破坏强度的计算结果
在离散元计算中，块体单元之间可以相互分离和滑动。在对土石混合体进行单向加载过程中，由于空间分布的不均匀性会造成局部应力集中，从而有可能使部分块体出现滑动和分离，但整体并未表现失稳状态，只有当较多块体发生滑动、分离并形成贯通面时，才认为整体发生了破坏。
由以上的统计分析结果可知，土石混合比和岩石块度大小是影响土石混合体宏观变形参数和强度参数的两个重要因素。对土石混合体进行单向数值模拟试验，分别改变混合比和块度大小，对其结果作宏观应力-应变分析。
从图7中可以看到：
（1）应力、应变关系有“塑性”段出现，塑性是由于滑移面增多的结果。
（2）当岩石含量达到80%时，出现脆性破坏。
（3）岩石含量越高，强度越大。
图7 相同块度大小（3m×3m×3m）不同混合比情况下应力-应变曲线图
（4）土石混合体在岩石含量较低的情况下表现为近似土体的性质，而当岩石含量较高时则表现为近似岩石的性质。
岩石含量为80%时其弹性模量近似为20GPa，与岩石的弹性模量Er＝30GPa比较接近。岩石含量为20%时其弹性模量近似为3GPa，与土体的弹性模量Es＝2GPa比较接近。
图8表明：
（1）岩石块度远小于研究区域时，应力、应变呈现线弹性，即可以认为是宏观均匀材料。
（2）岩石块度增大，应力、应变曲线呈现分段线性的现象，应力较大时，区域内滑移接触面增多，宏观表现为塑性过程。
（3）随着岩石块度增大，破坏强度明显降低。主要是由于岩石块度越大，土石混合体内部的弱面就越相对比较集中，从而容易导致整体的失稳。
（4）块度较小时，土体单元的变形相对充分，导致岩石块度很小时变形模量较小，而随着块度增大变形模量变大。
图8 相同混合比（40%岩石）不同块度情况下应力-应变曲线
5 结论
（1）该随机模型的可靠性较好，数值模拟可以给出土石混合体的随机分布状态。
（2）针对同一研究条件，不同随机过程稳定性较好。
（3）对于不同的土石混合比，在单轴加载下的内部应力场分布会有不同，根据统计分析得出：在土石比为3∶2时，应力场空间分布不均匀性最明显。
（4）岩石块度大小对内部应力场分布影响很大，岩石块体单元集中的地方一般是高应力区。
（5）土石混合体的混合比和岩石块度大小是影响其变形和破坏特性的两个重要因素。
参考文献
［1］ Axelrad D R.Stochastic mechanics of discrete media.Springer-Verlag Berlin Heidelberg，1993
［2］油新华，汤劲松.土石混合体野外水平推剪试验研究.岩石力学与工程学报，2002，21（10）：1537～1540
［3］ Dowding C H，O'Connor K M.Distinct element modeling and analysis of mining induced subsidence.Rock Mech Rock Engng，1992，25：1～24
［4］王泳嘉，邢纪波.离散单元法及其在岩土力学中的应用.沈阳：东北工学院出版社，1991
［5］董大鹏.三维可变形离散元改进算法及其应用.中国科学院力学研究所，2002
［6］ Meyer T，Einstein H H.Geologic stochastic modeling and connectivity assessment of fracture systems in the Boston Area.Rock Mech Rock Engng，2002，35（1）：23～44
［7］ Cundall P A.A computer model for simulating progressive，large scale movements in blocky rock systems.Proceedings of the International Symposium Rock Fracture，ISRM，Vol.Nancy，1971