柱状节理岩体各向异性特性及尺寸效应研究_朱道建

第28卷第7期岩石力学与工程学报V ol.28 No.7 2009年7月Chinese Journal of Rock Mechanics and Engineering July，2009 柱状节理岩体各向异性特性及尺寸效应研究

朱道建1，2，杨林德1，2，蔡永昌1，2

(1. 同济大学地下建筑与工程系，上海 200092；2. 同济大学岩土及地下工程教育部重点实验室，上海 200092)

摘要：白鹤滩水电站坝基基岩及深部洞室围岩以柱状节理岩体为主，其力学特性极其复杂。为此，开发多弱面软

化本构模型，可以方便地定义和描述多组斜交节理面的力学属性和屈服后的软化特性；并考虑柱状节理岩体内部

介质的非均匀性，从细观角度建立整体力学属性满足Weibull分布的概率模型；在对概率分布参数进行敏感性分

析的基础上，结合刚性承压板试验结果，确定反映现场岩体结构的概率分布参数。同时，采用V oronoi算法构建

了无规则性和随机性排列的四面体、五面体和六面体柱面的随机模型，综合该三方面因素对柱状节理岩体的各向

异性特性展开研究，得出不同方向承载性能和变形特性的差异。对不同尺寸下多组试块进行分析计算，获取了不

同方向的弹性模量和抗压强度随试块尺寸的变化曲线，证明柱状节理岩体还具有明显的尺寸效应，并确定了其特

征尺度。

关键词：岩石力学；柱状节理岩体；各向异性；概率模型；V oronoi算法；尺寸效应

中图分类号：TU 45 文献标识码：A 文章编号：1000–6915(2009)07–1405–10 RESEARCH ON ANISOTROPIC CHARACTERISTICS AND SIZE EFFECT

OF COLUMNAR JOINTED ROCK MASS

ZHU Daojian1，2，YANG Linde1，2，CAI Yongchang1，2

(1. Department of Geotechnical Engineering，Tongji University，Shanghai200092，China；2. Key Laboratory of Geotechnical and

Underground Engineering of Ministry of Education，Tongji University，Shanghai200092，China)

Abstract：The mechanical properties of columnar jointed rock mass，which is the main composition of dam foundation and surrounding rock of deep cavern at Baihetan hydropower station，are very complex. In order to study the anisotropic characteristics of columnar jointed rock mass，three factors are considered. First，a mixed multi-weakness planes softening constitutive model is built. It can be used to define and describe the mechanical properties and the softening characteristics of multiple oblique joint planes after failure. Second，considering the heterogeneity of the columnar jointed rock mass，a probability model is established through assigning mechanical parameters throughout the numerical specimens by following Weibull distribution. On the basis of sensitivity analysis of probability distribution parameters，combining with the experimental results of rigid bearing plate，the reasonable distribution parameters that fit the structure of field rock mass are computed. Third，column section of jointed rock mass is composed by tetrahedron，pentahedron and hexahedron randomly and irregularly. V oronoi algorithm is used to achieve the randomness of geometric shape and irregularity of arraying order. Taking all these factors into account，the bearing capacity and deformation characteristics in different directions are obtained. The results turn out to be obviously distinct in different directions. At the same time，multiple test blocks with different sizes are analyzed；and the curves of equivalent elastic modulus and unconfined compressive strength in different

收稿日期：2008–12–03；修回日期：2009–02–20

基金项目：国家自然科学基金资助项目(50678135)

作者简介：朱道建(1981–)，男，1999年毕业于中国矿业大学建筑工程学院土木工程专业，现为博士研究生，主要从事隧道工程及岩石力学方面的研究工作。E-mail：zdj38@https://www.360docs.net/doc/f218759006.html,

? 1406 ? 岩石力学与工程学报 2009年

directions changing with size are acquired. It is shown that the columnar jointed rock mass has obvious size effect；meanwhile the characteristic scale is determined.

Key words：rock mechanics；columnar jointed rock mass；anisotropy；probability model；V oronoi algorithm；size effect

1 引言

岩石类材料是一种具有复杂力学性质的非均匀准脆性材料，本文所研究的柱状节理岩体更具有其自身的特殊性质。以往较多的研究将此类准脆性材料受力后的变形和断裂过程的非线性归结为弹塑性，仅用宏观上的弹塑性理论来表述。然而，这种基于经典力学理论的力学模型忽略了岩石材料内部细观结构的非均匀性，仅考虑材料力学特性的非线性特征，不足以表达岩石变形整个过程所表现的复杂性。事实上，在一个统一的变形场中，微破裂不断产生的原因除了荷载不均、形态不够光滑等结构因素形成的应力集中外，更主要的是细观单元体力学性质的不均匀性。所谓岩石介质的非均匀性，就是岩石力学性质(如弹性模量、强度等)在空间上分布的非连续性。如果将外力施加在岩石试样上，那么由于岩石结构的非均匀性，在岩石试样内出现的应力分布是相当复杂的，即岩石中的应力分布也相应表现出高度的非均匀性。就细观尺度(一般指组成岩石的晶粒、微裂缝等尺度)而言，岩石是一种典型的非均匀材料。因此，可以认为材料的非线性特征与其细观非均匀性有直接联系[1]。

白鹤滩水电站工程位于西南部金沙江之上，属于特大水电工程，其坝基基岩部分以柱状节理岩体为主，前期进行了大量的现场试验，耗费了大量的人力、物力和时间，但试验数据非常离散，例如承压板试验中，出现有两试验点相隔仅2 m，试验结果相差却有一倍之多的情况，这给试验结果的分析和参数选取带来了非常大的困难。同时，试验结果显示，同一地点岩体竖向变形模量明显小于横向变形模量，横向变形模量约为竖直向变形模量的2倍，表现出了明显的各向异性特性。因此，在稳定性分析过程中，如何合理有效地描述其各向异性特性是本工程的重中之重。如果将细观力学的思想应用到岩石破坏问题中，假定细观尺度的单元力学性质服从某种既定的分布，并沿用连续介质力学方法描述单元的行为，通过弹性模量和强度等力学性质的弱化或退化描述单元由于节理面和微裂隙引发的损伤，由不同属性的细观单元在空间上的不同分布形成整体模型的各向异性特性，再引入合理的本构模型来描述宏观上多组节理面的存在表现出的力学特性，并考虑节理分布的随机性，以此来分析柱状节理岩体这类特殊脆性材料的变形破坏情况效果将更好。

本文首先从细观角度来研究柱状节理岩体，成功实现了其材料力学参数在数量统计和空间上按照Weibull分布函数随机分布的概率模型；其次，在宏观上开发了合理的本构模型来描述多节理面的力学属性和屈服后的软化特性；再次，采用V oronoi算法构建了无规则性和随机性排列的四面体、五面体和六面体柱面的随机模型；最后，综合上述因素对柱状节理岩体的受力变形特性进行了计算分析，合理描述了其各向异性特性，并研究了柱状节理岩体的尺寸效应，更好地对柱状节理岩体力学特性进行了阐述，可以为相关实际工程提供参考并指导工作。

2 复合型多弱面软化模型的建立

本文所建立的复合型多弱面软化模型属于弹塑性本构模型，该模型从两方面对节理岩体进行了复合：(1) 将节理岩体材料视为由岩体和节理组成的一种广义宏观复合材料，岩体和节理材料分别服从自身的应力–应变本构关系，并设置了反映自身属性的评价指标，可分别进行相关属性的赋值。(2) 岩体材料采用弹塑性本构关系，其屈服准则是将Mohr-Coulomb和Hoek-Brown强度屈服准则[2，3]复合嵌入本构模型中，通过设定引入的屈服因子，可根据工程情况合理方便地选择相应的屈服准则。而硬化/软化阶段则根据不同的强度屈服准则，对相关的塑性参数指标进行软化曲线设置。

2.1 岩体强度屈服准则

本模型为弹塑性模型，其弹性阶段应力增量和应变增量均服从Hooke定律，即

e e e

111223

()

σαεαεε

?=?+?+?(1)

e e e

212213

()

σαεαεε

?=?+?+?(2)

e e e

313212

()

σαεαεε

?=?+?+?(3)

第28卷 第7期 朱道建，等. 柱状节理岩体各向异性特性及尺寸效应研究 ? 1407 ?

式中：1α，2α均为岩体材料常数，14/3K G α=+，22/3K G α=?，K 为岩体体积模量。

由于岩体的抗拉强度远不及其抗压强度，因此在应用Mohr-Coulomb 强度屈服准则的同时，引入了抗拉强度屈服准则[4]，即当岩体所承受的拉应力超过其抗拉强度时，岩体材料即发生拉破坏，同时还可以考虑材料的双线性特性。

与Mohr-Coulomb 和抗拉屈服准则相关联的流动法则表达式为

t p p 31()h a σσσσ=?+? (4)

其中，

p a N ?=+

p t 2N ?σσ=?

在引入Hoek-Brown 强度屈服准则时，分为两种情况：

(1) 当3σ≥ci b /s m σ?时，其屈服准则表达式为

f f 13ci b 3ci (/)0a F m s σσσσσ=??+= (5) (2) 当3σ＜ci b /s m σ?时，其屈服准则表达式为

f

f 1

3

ci b 3ci (/)0a

F m s σσσσσ=?++= (6)

为了反映岩石材料屈服阶段的体积行为，本文引入了流动参数γ来描述其流动法则，它是一个与应力相关的变量。

当cv

330σσ＜＜时，γ满足如下关系式：

1

cv af cv af 33

(1/)/γ

γγγσσ?=+? (7)

其中，

1af ci b 3ci b ci 1(/)(/)a a m s m γσσσσ?=++，cv 1γ=?

当cv

33σσ≥时，1γ=?。

2.2 岩体及节理面的软化

Mohr-Coulomb 强度准则中引入两个参数s k 和t k 来反映材料的硬化/软化特性，它是与单元的剪切

应变和拉应变相关的变量，那么塑性剪切硬化/软化参数增量也是与剪切应变相关的量。其中增量形式s k ?表达式为

s

ps ps ps a 1m ps ps ps

b 3m

ps

ps ps m 13()/3k εεεεεεεεε?=??=????

???=????

??=?+??

(8)

式中：ps 1ε?和ps 3ε?分别为第一和第三主应力方向的塑性剪切应变增量值。

t k 反映的是材料受拉应力作用产生变形后发

生软化效应的变量，它与拉应力产生的应变值有关，其增量的数学表达式为

t pt 3k ε?=? (9)

式中：pt 3ε?为单元所产生的拉应变增量。

Hoek-Brown 强度准则在反映材料软化的阶段性时，也是通过引入一个软化参数来进行的。该参数是一个与各方向剪切主应变量相关的变量，此处

定义为p

3e ，它是主程序执行过程中每步迭代后所产

生应变增量的总和，而应变增量的数学表达式为

1021

p

3210

F e F e e e F F ????=?=

? (10)

式中：0e ?为迭代前的应变初值，1e ?为迭代后的应变终值，0F 为应变为0e ?时由式(5)或(6)所求得的F 值，1F 为应变为1e ?时由式(5)或(6)所求得的F 值。

对于节理面来说，本模型中考虑了3组独立的节理弱面，每组节理面可以任意正交或斜交。可以按照工程的需要，设置所需的节理面位置和数目。各节理面均设置独立的局部坐标系[5]，参考基面为X -Y 坐标平面，节理面与基面形成的夹角为节理面

倾角，节理面法向向量在基面内的投影与Y 轴正向所形成的夹角为节理面倾向[5]。

各节理面强度屈服准则均采用Mohr-Coulomb 和受拉屈服的综合型强度屈服准则，采用该强度准则使节理面参数的确定简易化，可直接由现场试验数据获取。同时，对节理面也考虑了双线性特性，如图1中的AB 和BC 段，均满足复合型强度屈服准则，但两阶段与各自设置的塑性参数值有关，而CD 段则为受拉屈服段[6]。

图1 节理面强度屈服准则 Fig.1 Strength yield criteria of joint plane

节理面和岩体的软化特性均通过引入软化参数来进行设置，它主要与材料的塑性变形有关。其

? 1408 ? 岩石力学与工程学报 2009年

中，节理面参数s j k 和t j k 主要与节理面塑性剪切应变和拉应变相关，其增量形式的数学表达式分别为

s j k ?=

(11) t pt j 33k ε′′?=? (12)

式中：ps 3i ε′?(1 2 3i ′′′=，，)为节理面局部坐标系内各主方向应变增量；pt 33ε′′?为节理面拉应变增量。

采用VC++语言在FLAC 3D 的自写本构平台下进行二次开发，嵌入FLAC 3D 程序后则可进行节理岩体的计算研究及工程应用。本模型的特点是将节理岩体看成是由节理和岩体组成的一种广义复合材料，可对节理和岩体单独定义材料属性。两种材料均是弹塑性本构关系，但岩体材料可以通过设定其屈服因子选择Mohr-Coulomb 或Hoek-Brown 强度屈服准则。节理面则考虑了多组弱面的情况，同一单元内最多可设置3组相互正交或斜交的节理面，节理面强度屈服准则均服从Mohr-Coulomb 强度准则，同时可考虑节理材料的双线性特性，该特性可通过设定材料因子确定是否考虑其双线性特性。岩体材料和节理材料的软化特性均通过引入与应变相关的软化参数来进行描述，进入软化阶段后随着塑性变形的增大自动更新材料参数值。

3 柱状节理岩体各向异性特性

柱状节理岩体主要是岩浆喷发后，遇空气后冷凝面以特殊排列形式分布冷却而形成的。由于岩石非均质，导致该过程中发生了对流和冷却压缩的不均匀性，从而造成了柱面的不规则性，形成含有四面体、五面体和六面体等不规则形式的柱面。本文在研究柱状节理的各向异性特性[7

，8]

时，充分考虑

了3方面的因素：(1) 能合理描述多节理弱面力学特性的本构模型的开发；(2) 众多节理面的存在引起的岩体强度的弱化和非均匀性；(3) 柱面形状的不规则性和分布的随机性。前文已建立复合型多弱面软化模型来描述柱状节理岩体，此处主要针对节理岩体的非均匀性和节理面的无规则性展开讨论。 3.1 柱状节理岩体随机概率模型

3.1.1 Weibull 分布概率密度函数

为了描述柱状节理岩体材料性质的非均匀性，假定组成材料细观单元的力学性质满足Weibull 分

布，该分布可以按照如下分布密度函数来定义：

1

c c c 000()exp m m

m f σσσσσσ???

??

=???

????

??

(13)

式中：c σ为参数，可代表满足该分布的所有参数(例如强度、弹性模量、泊松比等)数值；0σ为与所有单元参数平均值有关的参数，但其数值并不等于参数平均值；m 为参数，定义了Weibull 分布密度函数的形状。

一般把0σ和m 称为材料的Weibull 分布参数，对于材料的每个力学参数，都必须在给定其Weibull 分布参数的条件下，按照式(3)给定的随机分布赋值。当m = 2.0，4.0和8.0时，Weibull 分布密度函数的曲线如图2所示。m 反映了参数的离散程度，当其由小到大变化时，材料细观单元强度分布密度函数由矮而宽变化到高而窄，细观单元强度分布变得较为集中，材料强度的均质性较为均匀，材料内部所包含的大部分细观单元近乎相同，接近于给定的参数[1]。显然，m 反映了数值模型中材料结构的均质性，本文称之为均质度，m 越大，组成材料的细观单元越趋于均匀。

图2 不同均质度时单元参数的分布密度函数 Fig.2 Distribution density functions of element parameter

corresponding to different homogeneities

3.1.2 Weibull 分布随机数的生成

利用[0，1]区间上均匀分布的随机数可以产生任意分布的随机数。主要的方法有反函数法、舍选法、离散逼近法、极限近似法和随机变量函数法等，常用的有反函数法和舍选法，这里不再赘述。

对于细观力学模型基元属性，首先需生成满足Weibull 分布的随机数。Weibull 分布概率密度函数的变换形式为

σ0

分布密度函数

第28卷 第7期 朱道建，等. 柱状节理岩体各向异性特性及尺寸效应研究 ? 1409 ?

1()exp x f x x ααααββ?????

??=??????????? (14)

对式(14)进行积分，即可得Weibull 分布函数为 1()()d {exp[(/)]}d F x f x x x x x ααααββ???==?=∫∫

1exp x α

β???

??

??????

????

??

(15) 令

1exp x u α

β???

=??????

(16)

如果u 为[0，1]区间上的随机数，可以得到Weibull 分布的随机数计算公式[7]为

1/[ln(1)]x u αβ=?? (17) 式中：α和β分别对应式(13)中的m 和0σ，此处进行变化是为了进行数学计算的方便和理解。

通过上述转换，利用FLAC 3D 程序中自带URAND 函数，生成[0，1]区间上的随机数，再通过FLAC 3D 中FISH 语言，

按上述方法完成转化过程，即可顺利生成满足Weibull 分布的随机数。 3.1.3 基元的随机分布与赋值

所谓的“基元”，即是细观模型的最小破裂单位，基元不但包含了介质物理力学性质特征，同时也是数值计算的单元。在成功生成满足Weibull 分布随机数的基础上，可实现基元属性的随机分布与赋值。此处的随机分布，不仅要求基元力学属性在数量统计上满足Weibull 随机分布，同时还包括空间上的随机分布，因为同样的概率分布在空间上的有序排列和无序排列，其力学破坏形态将截然不同[9]。由于柱状节理岩体节理面众多，且力学属性较弱，节理面周围的单元属性还要满足大于节理面所在单元的力学属性。总的来说，随机概率模型既从整体上考虑了大量节理面和微裂隙的存在对整体模型强度的弱化效应，即基元力学性质的不均匀性，又进一步考虑了节理面位置处和周边岩体弱化程度上的逻辑关系，使该概率模型能更好地描述其力学破坏行为。

基于上述几方面的要求，FISH 语言实现的基本思路为：(1) 采用循环函数，随机生成与单元个数相等的Weibull 随机数，循环次数为基元的总体个

数。(2) 随机选择基元的ID 号，逐个随机赋予各基

元的力学参数值。(3) 赋值过程中设置条件语句。条件一为：生成的ID 号若出现相等的情况，则再重新生成，直至随机生成的ID 号与已赋单元ID 不等；条件二为：如已赋值某条节理面位置处的单元属性，则赋予与该节理面单元相邻的单元时，该需判断两者的大小关系，若节理面周边单元属性小于该节理面单元属性时，则重新选择Weibull 随机数，直至满足节理面周边单元属性大于该节理面单元时为止。同样，如果先赋予节理面周边的单元属性，后赋予节理面位置处的单元属性，同样进行相应的判断。图3为细观模型基元属性随机分布的三维视图。

(a) 基元属性分布剖面图

(b) 基元属性三维空间分布图

图3 细观模型基元属性随机分布的三维视图 Fig.3 3D views of random distribution of meso-model

element properties

3.1.4 概率模型分布参数的确定

采用细观模型来模拟柱状节理岩体，其内部的微裂隙及节理面等结构面造成岩体的弱化，采用计算单元的力学属性按照Weibull 分布随机的方法赋值，从而从细观角度更好地描述了柱状节理岩体内部复杂多变的岩性。采用该方法来处理柱状节理岩体时，关键是确定合理的概率模型分布参数，使数值模型与现场的柱状节理岩体尽可能接近。

? 1410 ? 岩石力学与工程学报 2009年

本文采用的方法是通过建立细观模型来拟合现场承压板试验，对多个试验点的试验结果一一进行

拟合，得出一系列与各试验点吻合的概率模型分布参数。在此基础上，对概率分布参数进行统计分析及相关处理，确定出能代表多数试验点的概率模型分布参数的终值。

(1) 参数的敏感性分析

由于试验点个数众多，在对各试验点拟合计算其概率分布参数之前，本文对Weibull 分布的均质度m 和参数0σ进行了影响性分析，主要从两方面进行了讨论：一是研究同一种岩性时，不同的均质度对岩体受力变形的影响；二是研究相同的均质度

m ，不同的岩性强度对岩体受力变形的影响。在对这两个关键性参数分析的基础上，了解各自的敏感性，将对后续工作有很大的帮助。

基于上述分析目的，本文进行了大量不同概率分布模型的变形计算，图4为材料参数0σ= 22，不同均质度m 时，岩体竖向位移的变化曲线图。由图可知，随着均质度m 的不断增加，岩体整体强度逐渐提高，在相同荷载作用下，变形量呈逐渐减小趋势。图5为均质度m = 7.0，不同的材料参数0σ时，岩体竖向位移的变化曲线。由图可知，随着0σ取值的不断增大，均值度m 不变时，岩体整体强度逐渐提高，其竖向位移同样呈逐渐递减趋势。比较图4，5的计算结果可知，位移随两参数的变化趋势均服从指数型函数，但参数0σ的敏感性较大，均质度则相对偏小。

图4 岩体竖向位移随均质度变化曲线

Fig.4 Variable curve of rock mass vertical displacement

with homogeneity

拟合函数采用指数型二元一次函数：

12(/)(/)120e e x t x t y A A y ??=++ (18)

图5 岩体变形随参数σ0变化曲线

Fig.5 Variable curve of rock mass vertical displacement with

parameter σ0

图4，5拟合结果分别为

(/10.92)(/2.06)0.42e 1.410e 3.83m m y ??=++ (19) 00(/3.55)(/15.09)32.24e 10.864e 1.52y σσ??=++ (20) (2) 概率分布参数的确定

现场刚性承压板首先在试验部位清除爆破松动层，手工凿制成φ 60 cm 的圆形平面，用砂轮磨平，作为试验面；试验面周围1 m 范围亦大体凿制平整，以满足试验边界条件。清洗试验面及其边界岩体，进行地质描述、拍照。刚性承压板面积为2 000 cm 2；最大试验荷载为8或10 MPa ，分5级采用逐级一次循环法加压，采用4只千分表对称布置在承压板的4个方位，测量岩体变形，加载后立即测读变形值，以后每隔10 min 测读一次，当承压板上4个测表相邻2次读数差与同级压力下第一次变形读数和前一级压力下最后一次变形读数差之比小于5时，则认为变形稳定，再加(卸)下一级荷载[9]。

图6所示为现场刚性承压板试验的试验面。按上述要求处理完试验面后，即可按照试验要求施加荷载。图7为拟合刚性承压板试验所建立的三维细观模型，其半径为25.25 cm ，高100 cm 。整体模型划分网格单元66 080个，节点80 426个，每个单元

边长约1 cm 。单元的力学属性在数量和空间上均按Weibull 概率分布函数随机分配。荷载施加面为模型上顶面，垂直向下逐步施加荷载，节理面倾角与底面呈70°夹角。边界条件为：周边均约束法向变形，而底部设置为固定边界[10]。

本文以现场试验最大荷载为10 MPa 为例，进行分析。试验的初步结果为：节理岩体竖向压缩模量为7.85 GPa ，岩体最终变形量约为0.47 mm 。

m

岩体竖向位移/(10－

4 m )

2

468σ0

　岩体竖向位移/(10－

4 m )

第28卷 第7期 朱道建，等. 柱状节理岩体各向异性特性及尺寸效应研究 ? 1411 ?

图6 现场刚性承压板试验面 Fig.6 Surface of field rigid bearing plate test

图7 三维细观计算模型 Fig.7 3D meso-computing model

现场试验点中，III1级岩体试验点共14个，III2级岩体试验点共10个。通过对各试验点的试验结果进行数值拟合计算，得到了满足各试验点条件的概率模型参数。图8所示为现场试验和拟合计算后的压力–变形曲线对比图。

图8 压力–变形曲线试验与拟合结果对比 Fig.8 Comparison of testing and fitting results of pressure-

deformation curves

表1为各试验点概率参数拟合结果。由表1可知：

① III1级中较好岩体的弹性模量均值的平均值为28.26 GPa ，其概率模型参数可取为：0σ= 30，

m = 8.0；

表1 各试验点概率参数拟合结果

Table 1 Fitting results of probability parameters of testing

points

岩体级别

岩性σ0 m 弹性模量均值/GPa 较好

30 2.8，4.5，10， 12，15，20

26.69，27.38，28.54，28.74，28.99，29.21

22 1.9，2.2，8.0 19.55，19.56，20.75 一般

20

3.0，5.5 17.86，18.47 14

1.8，8.0 1

2.51，1

3.19

III1

较差

12 2.8

10.74

22 10，10，20 20.92，20.92，21.41 一般20 1.5，2.2 17.59，17.71 12 3.0，3.0，4.5 10.73，10.73，10.94 III2

较差10 2.5，4.2 8.87，9.10

② III1级中一般岩体的弹性模量均值的平均值为19.24 GPa ，其概率模型参数可取为：0σ= 21，

m = 4.0；

③ III1级中较差岩体的弹性模量均值的平均值为12.15 GPa ，其概率模型参数可取为：0σ= 13，均质度系数m = 5.5；

④ III2级中一般岩体的弹性模量均值的平均值为19.71 GPa ，其概率模型参数可取为：0σ= 21，

m = 5.0；

⑤ III2级中较差岩体的弹性模量均值的平均值为10.07 GPa ，其概率模型参数可取为：0σ= 11，

m = 4.2。

3.2 节理面分布的无规则性和随机性

节3.1通过建立柱状节理岩体的随机概率模型来反映岩体内部的非均匀性，并通过拟合现场刚性承压板试验获取了符合现场柱状节理岩体力学特性的概率分布参数。对于柱面形状的不规则性和分布的随机性，本文结合V oronoi 图的生成算法和编程软件成功地实现开发。V oronoi 图是一个关于空间划分的基础结构，在材料力学中经常用以研究晶体的力学特性[11]，在此不作详细介绍。应用和开发流程为：首先将V oronoi 图的生成算法在VC++中编程实现，可以在一定区域内生成满足设定要求的V oronoi 图；其次，将AutoCAD 作为图形输出界面，输出生成的V oronoi 图，在AutoCAD 内提取多边形顶点的几何坐标；再次，将已知各点的坐标编写成3DEC 命令文件，利用3DEC 软件能方便生成节理面和切割节理面的优点，实现柱状节理的三维化[12]；编写

2

4

6

8

10变形/(10－

3 mm)

压力/M P a

? 1412 ? 岩石力学与工程学报 2009年

接口程序将3DEC 生成的富含节理面的模型转入ANSYS 进行网格剖分；最后，建立与FLAC 3D 的程序接口，将ANSYS 中划分的单元转入FLAC 3D ，进行受力变形计算。为了消除模型在宽度和高度方向长度的不同而造成的影响，本文以III1级一般岩体为例，模型尺寸取为1 m×1 m×1 m 。

图9为随机生成的V oronoi 多边形柱状节理岩体模型，图10为节理面的分布形式。

(a)

(b)

图9 柱状节理岩体模型及三维网格划分 Fig.9 Model of columnar jointed rock mass and 3D grid

多弱面节理模型采用本文开发的复合型多弱面软化模型，可以合理描述众多节理面的分布形式及赋予相应的单元属性，而岩体的弱化和非均匀性则采用前文介绍的随机概率模型的建立方法和拟合确定的概率参数进行设置，从而综合地考虑了柱状节理岩体多方面因素的影响，进而可以更好地描述其各向异性的力学性质[13]。

本文主要从以下两个方面对柱状节理岩体的各向异性特性进行研究：(1) 恒定的加载速率下，水平向和竖向分别加载时，其最终承载力的差异；(2) 恒定的外荷载作用下，岩体受压变形水平向和竖向的等效弹性模量的差异。

(a)

(b)

图10 节理面的分布形式

Fig.10 Distribution forms of joint surface

如图11所示，当对计算模型竖向(z 向)和水平向(x 向)作用恒定速率时，两方向表现出了不同的承载力。x 向最大承载力为714 MPa ，z 向最大承载力

为621 MPa ，表现出明显的差异。达到最大承载力

后，随着变形的继续增加，模型整体进入软化，承载性能总体下降。

图11 恒定加载速率下的应力–应变曲线 Fig.11 Stress-strain curves under constant loading rate

当施加恒定的外荷载时，x 和z 向变形量也存在较大差异。如图

12所示，外荷载为6 MPa 时，模型在x 向的最终变形量为1.57×10－

4 m ，

而z 向最 0.000

0.008

0.016 0.024 0.032 0.040

应力/G P a

应变

第28卷 第7期 朱道建，等. 柱状节理岩体各向异性特性及尺寸效应研究 ? 1413 ?

图12 外荷载作用下的变形曲线

Fig.12 Deformation curves under external loads

终变形量则达2.82×10－

4 m ，较x 向变形增加近80%。

本文通过采用V oronoi 算法构造柱状节理岩体，并开发合理有效的本构模型，同时考虑节理岩体的非均匀性，综合该三方面的因素，对柱状节理岩体的各向异性特性进行相应的研究，计算结果符合现场的测试结果，证明了该研究方法的合理性。需要指出的是，水平向和竖向的应力差异不大，主要是变形差异非常显著。

4 柱状节理岩体的尺寸效应

依次选取柱状节理岩体模型尺寸为1 m×1 m ，2 m×2 m ，4 m×4 m ，6 m×6 m ，每一尺寸选取10组随机计算模型[12]，如图13所示。对模型不同方向进行单轴压缩计算，由计算获取的应力和应变值求解模型各方向的等效弹性模量。

(a) 1 m×1 m (b) 2 m×2 m

(c) 4 m×4 m (d) 6 m×6 m

图13 不同尺度下柱面形状及排列方式图 Fig.13 Shapes of cylinder section and arrangement ways

with different scales

图14为计算得到的不同尺寸试件x 向等效弹性模量散点图。模型尺寸为1 m×1 m 时，等效弹性模量波动区间较大，计算值分散；随着尺寸的逐渐增

图14 x 向等效弹性模量随模型尺寸变化的散点图 Fig.14 Scatter diagram of equivalent elastic modulus in

x -direction varying with model size

大，波动区间逐渐减小，当模型尺寸为4 m×4 m 时，

x 向等效弹性模量计算值已经较接近。同样，在y 和z 向也表现出同样的规律，但总体上y 与x 向差异较小。结果表明：柱状节理岩体表现出了明显的尺寸效应，随着模型尺寸增大，等效弹性模量趋于稳定。

对模型各方向施加恒定的速度，以考察不同尺寸下其单轴抗压强度的差异，如图15所示，随着模型尺寸的增大，模型各方向的承载性能也呈递减趋势，当模型尺寸为4 m×4 m 时，承载性能趋于稳定。而x 和y 向较接近，与z 向差异较大。

图15 单轴抗压强度随模型尺寸变化曲线

Fig.15 Variable curves of uniaxial compressive strength with

model size

上述两方面的计算结果表明，柱状节理岩体同时具有明显的尺寸效应[14]，其特征尺度可认定为4 m×4 m 。

5 结 论

研究柱状节理岩体的方法很多，对于众多节理

0.00.30.60.91.21.51.82.12.42.73.0变形量/(10－

4 m )

收敛时步

x 向等效弹性模量/G P a

模型边长/m

单轴抗压强度/G P a

模型边长/m

? 1414 ? 岩石力学与工程学报 2009年

面切割岩体的情况，一般来说，离散元可以更方便地来描述节理面分布情况。但本文采用连续介质来研究柱状节理岩体，主要是以下几方面的原因：(1) 考虑岩体材料的非均匀性，从细观角度将计算模型划分为满足基元尺寸的单元后，总体网格数量很大，采用离散元方法计算效率将大大下降；(2) 后期进行大型地下厂房开挖的稳定性计算时，整体模型达上千米之多，采用离散元计算该大型工程不太现实，而采用连续体来进行计算，已是很成熟的技术；

(3) 对于开发的复合型多弱面软化模型，是针对连续介质开发的，此处也是对该模型的应用效果和设置方法作检验和确认。

通过建立复合型多弱面软化模型，并考虑柱状节理岩体内部的非均匀性，建立了整体力学属性满足Weibull分布的概率模型，同时结合V oronoi算法综合考虑了节理面几何形态的无规则性和随机性，对柱状节理岩体的各向异性特性和尺寸效应进行了研究。结果表明，水平向和竖向整体强度和等效弹性模量都存在较大差异，表现出明显的各向异性，且随着模型尺寸的变化，其整体强度呈递减趋势，并逐渐趋于稳定，也具有明显的尺寸效应。这与现场勘察结果非常吻合，说明本文所建立的方法能够合理有效地描述柱状节理岩体的各向异性特性及尺寸效应[15]，这对后期细部关键部位的变形、受力分析及整体工程的稳定性研究均具有较高的指导意义和参考价值。

参考文献(References)：

[1] 唐春安，朱万成. 混凝土损伤与断裂——数值试验[M]. 北京：科学出

版社，2003.(TANG Chun′an，ZHU Wancheng. Concrete damage and fracture—numerical tests[M]. Beijing：Science Press，2003.(in Chinese)) [2] HART R，CUNDALL P，LEMOS J. Formulation of a three-dimensional

distinct element model—part II. mechanical calculations for motion and interaction of a system composed of many polyhedral blocks[J].

International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts，1988，25(3)：117–125.

[3] CUNDALL P，CARRANZA-TORRES C，HART R. A new constitutive

model based on the Hoek-Brown criterion[C]// Proceedings of the 3rd International Symposium on FLAC and FLAC3D Numerical Modeling in Geomechanic s.Sudbury：[s.n.]，2003：17–25.

[4] Itasca Consulting Group Inc.. FLAC3D(fast Lagrangian analysis of

continua in 3-dimension) user′s manual(version 2.10)[R]. Minneapolis：Itasca Consulting Group Inc.，1997.

[5] GERRARD C M. Elastic models of rock masses having one，two and

three sets of joints[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts，1982，19(1)：15–23.

[6] MARTI J，CUNDALL P. Mixed discretization procedure for accurate

modelling of plastic collapse[J]. International Journal for Numerical and Analytical Methods in Geomechanics，1982，6(1)：129–139. [7] 石安池，唐鸣发，周其健. 金沙江白鹤滩水电站柱状节理玄武岩岩

体变形特性研究[J]. 岩石力学与工程学报，2008，27(10)：2 079–2 086.

(SHI Anchi，TANG Mingfa，ZHOU Qijian. Research of deformation characterstics of columnar jointed basalt at Baihetan hydropower station on Jinsha river[J]. Chinese Journal of Rock Mechanics and Engineering，2008，27(10)：2 079–2 086.(in Chinese))

[8] 谢和平，陈忠辉. 岩石力学[M]. 北京：科学出版社，2004：60–

62.(XIE Heping，CHEN Zhonghui. Rock mechanics[M]. Beijing：

Science Press，2004：60–62.(in Chinese))

[9] 梁正召，唐春安，张永彬，等. 准脆性材料的物理力学参数随机概

率模型及破坏力学行为特征[J]. 岩石力学与工程学报，2008，27(4)：718–727.(LIANG Zhengzhao，TANG Chun′an，ZHANG Yongbin，et al. On probability model of physico-mechanical parameters of quasi-brittle materials and associated mechanical failure behaviors[J].

Chinese Journal of Rock Mechanics and Engineering，2008，27(4)：718–727.(in Chinese))

[10] 朱道建，杨林德，黄建勇. 厦门海底隧道地表沉降控制效果分析及

其预测[J]. 岩石力学与工程学报，2007，26(11)：2 356–2 362.(ZHU Daojian，YANG Linde，HUANG Jianyong. Analysis and prediction of effect of controlling ground settlement at Xiamen subsea tunnel[J].

Chinese Journal of Rock Mechanics and Engineering，2007，26(11)：

2 356–2 362.(in Chinese))

[11] 沈满德，余圣甫，王宁. 蒙特卡罗法模拟晶粒生长过程中的Voro-

noi模型[J]. 机械工程材料，2006，30(3)：11–13.(SHEN Mande，YU Shengfu，WANG Ning. Voronoi model in Monte Carlo simulation of grain growth[J]. Materials for Mechanical Engineering，2006，30(3)：11–13.(in Chinese))

[12] 宁宇，徐卫亚，郑文棠，等. 柱状节理岩体随机模拟及其表征单

元体尺度研究[J]. 岩石力学与工程学报，2008，27(6)：1 202–1 208.

(NING Yu，XU Weiya，ZHENG Wentang，et al. Study of random simulation of columnar jointed rock mass and its representative elementary volume scale[J]. Chinese Journal of Rock Mechanics and Engineerin g，2008，27(6)：1 202–1 208.(in Chinese))

[13] 胡大伟，朱其志，周辉，等. 脆性岩石各向异性损伤和渗透性演

化规律研究[J]. 岩石力学与工程学报，2008，27(9)：1 822–1 827.

(HU Dawei，ZHU Qizhi，ZHOU Hui，et al. Research on anisotropic damage and permeability evolutionary law for brittle rocks[J]. Chinese Journal of Rock Mechanics and Engineerin g，2008，27(8)：1 822–

1 827.(in Chinese))

[14] 杨圣奇，徐卫亚. 不同围压下岩石材料强度尺寸效应的数值模拟[J].

河海大学学报(自然科学版)，2004，32(5)：578–582.(YANG Shengqi，XU Weiya. Numerical simulation of strength-scale effect of rock materials under different confining pressures[J]. Journal of Hohai University(Natural Science)，2004，32(5)：578–582.(in Chinese)) [15] 尤明庆，邹友峰. 关于岩石非均质性与强度尺寸效应的讨论[J]. 岩

石力学与工程学报，2000，19(3)：391–395.(YOU Mingqing，ZOU Youfeng. Discussion on heterogeneity of rock material and scale effect on specimen strength[J]. Chinese Journal of Rock Mechanics and Engineering，2000，19(3)：391–395.(in Chinese))

–