tceic.com
学霸学习网 这下你爽了
相关文章
当前位置:首页 >> >>

udec中文说明


通用离散元用户指导
(U D E C 3.1)

山东科技大学
2004.9 目
1 引



言 ............................................................. 1

1.1 总 论............................................................ 1 1.2 与其他方法的比较 ................................................. 2 1.3 一般特性 ......................................................... 2 1.4 应用领域 ......................................................... 3 2 开始启动 ........................................................... 4 2.1 安装和启动程序 ................................................... 4 2.1.7 内存赋值 ..................................................... 4 2.1.9 运行 UDEC .................................................... 5 2.1.10 安装测试程序 ................................................ 5 2.2 简单演示-通用命令的应用 ......................................... 5 2.3 概念与术语 ....................................................... 6 2.4 UDEC 模型:初始块体的划分 ......................................... 8 2.5 命令语法 ......................................................... 9 2.6 UDEC 应用基础 .................................................... 10 2.6.1 块体划分 ..................................................... 10 2.6.2 指定材料模型 ................................................. 16 2.6.2.1 块体模型 ................................................. 16 2.6.2.2 节理模型 ................................................. 17 2.6.3 施加边界条件和初始条件 ....................................... 19 2.6.4 迭代为初始平衡 ............................................... 21

2.6.5 进行改变和分析 ............................................... 24 2.6.6 保存或恢复计算状态 ........................................... 25 2.6.7 简单分析的总结 ............................................... 25 2.8 系统单位 ......................................................... 26 3 用 UDEC 求解问题 .................................................... 27 3.1 一般性研究 ....................................................... 27 3.1.1 第 1 步:定义分析模型的对象 ................................... 28 3.1.2 第 2 步:产生物理系统的概念图形................................ 28 3.1.3 第 3 步:建造和运行简单的理想模型.............................. 28 3.1.4 第 4 步:综合特定问题的数据 ................................... 29 3.1.5 第 5 步:准备一系列详细的运行模型.............................. 29 3.1.6 第 6 步:进行模型计算 ......................................... 29 3.1.7 第 7 步:提供结果和解释 ....................................... 30 3.2 产生模型 ........................................................ 30

3.2.1 确定 UDEC 模型合适的计算范围 .................................. 30 3.2.2 产生节理 ..................................................... 32 3.2.2.1 统计节理组生成器 ....................................... 32 3.2.2.2 VORONOI 多边形生成器 .................................... 34 3.2.2.3 例子 ................................................... 34 3.2.3 产生内部边界形状 ............................................. 35 3.3 3.4 变形块体和刚体的选择 ............................................ 38 边界条件 ........................................................ 42

3.4.1 应力边界 ..................................................... 42 3.4.1.1 施加应力梯度 ........................................... 43 3.4.1.2 改变边界应力 ........................................... 44 3.4.1.3 打印和绘图 ............................................. 44 3.4.1.4 提示和建议 ............................................. 45 3.4.2 位移边界 ..................................................... 46 3.4.3 真实边界-选择合理类型 ....................................... 46 3.4.4 人工边界 ..................................................... 46 3.4.4.1 对称轴 ................................................. 46 3.4.4.2 截取边界 ............................................... 46 3.4.4.3 边界元边界 ............................................. 49 3.5 初始条件 ........................................................ 50

3.5.1 在均匀介质中的均匀应力:无重力................................ 50

-1-

3.5.2 无节理介质中具有梯度变化的应力:均匀材料 ...................... 51 3.5.3 无节理介质中具有梯度变化的应力:非均匀材料 .................... 51 3.5.4 具有非均匀单元的密实模型 ..................................... 52 3.5.5 随模型变化的初始应力 ......................................... 53 3.5.6 节理化介质的应力 ............................................. 54 3.5.7 绘制应力等值线图 ............................................. 55 3.6 3.7 加载与施工模拟 .................................................. 57 选择本构模型 .................................................... 62

3.7.1 变形块体材料模型 ............................................. 63 3.7.2 节理材料模型 ................................................. 64 3.7.3 合理模型的选择 ............................................... 65 3.8 材料性质 ........................................................ 71

3.8.1 岩块性质 ..................................................... 71 3.8.1.1 质量密度 ............................................... 71 3.8.1.2 基本变形性质 ........................................... 71 3.8.1.3 基本强度性质 ........................................... 72 3.8.1.4 峰后效应 ............................................... 73 3.8.1.5 现场性质参数的外延 ..................................... 77 3.8.2 节理性质 ..................................................... 80 3.9 提示和建议 ...................................................... 81

3.9.1 节理几何形状的选择 ........................................... 81 3.9.2 设计模型 ..................................................... 81 3.9.3 检查模型运行时间 ............................................. 82 3.9.4 对允许时间的影响 ............................................. 82 3.9.5 单元密度的考虑 ............................................... 83 3.9.6 检查模型响应 ................................................. 83 3.9.7 检查块体接触 ................................................. 83 3.9.8 应用体积模量和剪切模量 ....................................... 83 3.9.9 选择阻尼 ..................................................... 84 3.9.10 给块体和节理模型指定模型和赋值............................... 84 3.9.11 避免圆角误差 ................................................ 85 3.9.12 接触嵌入 .................................................... 85 3.9.13 非联结块体 .................................................. 86 3.9.14 初始化变量 .................................................. 86 3.9.15 确定坍塌荷载 ................................................ 86

-2-

3.9.16 确定安全系数 ................................................ 86 3.10 解 释.......................................................... 88

3.10.1 不平衡力 .................................................... 88 3.10.2 块体/网格结点的速度 ......................................... 88 3.10.3 块体破坏的塑性指标 .......................................... 89 3.11 模拟方法 ....................................................... 89

3.11.1 有限数据系统模拟 ............................................ 89 3.11.2 混沌系统的模拟 .............................................. 90 3.11.3 局部化、物理的不稳定性和应力路径............................. 91

-3-

1 引 言
1.1 总 论
通用离散元程序(UDEC,Universal Distinct Element Code)是一个处理不连续介 质的二维离散元程序。UDEC 用于模拟非连续介质(如岩体中的节理裂隙等)承受静载 或动载作用下的响应。非连续介质是通过离散的块体集合体加以表示。不连续面处理为 块体间的边界面,允许块体沿不连续面发生较大位移和转动。块体可以是刚体或变形 体。变形块体被划分成有限个单元网格,且每一单元根据给定的“应力-应变”准则,表 现为线性或非线性特性。不连续面发生法向和切向的相对运动也由线性或非线性“力-位 移”的关系控制。在 UDEC 中,为完整块体和不连续面开发了几种材料特性模型,用来 模拟不连续地质界面可能显现的典型特性。UDEC 是基于“拉格朗日”算法很好地模拟块 体系统的变形和大位移。 UDEC 包含了功能强大的程序语言 FISH 函数。借助于 FISH 函数,用户可以编写自 己 的 功 能 函数 , 扩 展 UDEC 的 应 用 功 能 。 FISH 函 数 为 简 化 分 析 ,适 应 特 殊 要求 的 UDEC 的用户,提供了一个强有力的工具。 UDEC 采用的离散单元法理论由 Cundall(1971)首次提出,至今已经过了 20 多年 的发展。在 1985 年,Cundall 博士和 Itasca 公司在 IBM 系列兼容微机上开发了 UDEC 工 程计算应用程序。该软件为建立数以千块模型的高速计算而设计。基于浮点运算速度的 优势和低成本的内置 RAM,用 UDEC 程序可大大地提高了计算大规模问题的能力。例 如,在具有 4MB RAM 的微机上,UDEC 能够求解 2500 个刚体(或 1000 个具有 8 个自 由度变形体)的模型。该模型的求解速度大约为每分钟 200 次。在 RAM 确定的情况 下,其计算速度是与模型的块体数量成线性关系。 对于典型的模型,约 1500 个刚体(或 500 个变形体)或更少,在 UDEC 中采用的 显式解法,大约需要 2000~4000 计算步可以获得问题的解。例如,一个 500 个变形块 体的模型,计算 4000 步大约需要 6 min。因此,典型的工程问题用 UDEC 计算仅需几十 分钟或几个小时。 UDEC 是一个命令驱动(而不是菜单驱动)的计算程序。尽管菜单驱动程序易于初 次学习,但在 UDEC 中所提供的命令驱动结构具有如下优点: 1、输入的“语言”是基于可识别的文字命令,使你易于识别每一个命令的作用(例如 BOUNDARY 命令,是指施加模型的边界条件)。 2、工程模拟通常是按照系列施工顺序构成 ―― 即,构造原岩应力,施加作用的荷 载、开挖隧道、安装支护等。一系列(从文件或键盘上)输入命令完全对应于实际的施 工顺序。 3、根据文本编辑器,很容易对 UDEC 数据文件进行编辑和修改。几个数据文件能

-1-

相互连接,进行多个问题求解,这对于进行参数的灵敏度分析是十分有用的。 4、命令驱动结构允许用户开发前后处理程序,控制 UDEC 必要的输入/ 输出。用户 可以为一系列 UDEC 的模拟,编写节理模拟函数,产生特定的节理结构。可采用 FISH 程序语言,并插入到输入的文件中,使计算很容易实现。

1.2 与其他方法的比较
对于 UDEC 程序,一个共同的问题是,UDEC 是一个有限元程序还是离散元程序? 他们的主要区别是什么?UDEC 程序与其他程序有何关系?为回答上述问题,下面将给 予解释。 许多有限元、边界单元和拉格朗日有限差分程序都具有“界面单元”或“节理单元”, 使程序能够模拟问题中的不连续面,扩大程序的应用范围。然而,他们的公式在一个或 多个方面通常受到限制:首先,当考虑很多相互切割的节理就可能打乱系统的逻辑关 系;其次,不可能自动识别新的接触面进行自动考虑;第三,计算公式可能有小位移和 无转动条件限制,所以通常适用连续介质的程序。 术语“离散单元法”(Discrete element method)意味着: (a)允许离散块体发生有限的位移和转动,包括完全脱离; (b)在计算过程中,自动识别新的接触面。 在不连续介质中,如果没有第一个属性,程序不可能产生某些重要的机理。如果没 有第二个特性,程序将限制在事先已知的相互作用的有限块体数。离散元法(Distinct element method)是由 Cundall 和 Strack(1979)采用变形接触和显式、时间域的初始运 动方程(而不是变换,块体方程)提出的特殊的离散单元法程序。 离散单元法的计算机程序主要有以下四类: 1、Distinct Element Programs - 该类程序采用显式时间步直接进行运动方程的求解。 块体可以是刚体或变形体(通过细分成单元);接触面是可变形的。UDEC 就属此类。 2、Modal Methods - 该类方法类似于刚体离散单元法,但对于变形体采用模型叠 加技术。 3、Discontinuous Deformation Analysis - 接触是刚体,块体可以是刚体或变形体。 通过迭代算法可以获得非嵌入条件;块体变形性基于应变模型的叠加。 4、Momentum-Exchange Methods-接触面和块体都是刚体:块体接触面在瞬时碰 撞的过程中惯性矩发生交换,可以表征滑动和摩擦特性。

1.3 一般特性
UDEC 主要用于岩石边坡的渐进破坏研究及评价岩体的节理、裂隙、断层、层面对地 下工程和岩石基础的影响。UDEC 对研究不连续特征的潜在破坏模型是十分理想的工具。

-2-

当地质结构特征明显且易于明确描述的情况适宜使用该程序进行分析。UDEC 开发 了人工或自动节理生成器,用以模拟产生岩体中一组或多组不连续面。在模型中,可以 产生变化范围较大的节理模式。屏幕绘图工具允许用户随时观看节理模型。在最后确定 所选择的节理模型前,能容易进行调整与修改。 也可以获得不同的节理材料特性。基本模型是指定节理弹性刚度、摩擦角、粘结 力、张拉强度和剪胀特性的库仑滑动准则。对该模型的改进包括随着位移的发展而粘结 力和张拉强度的降低弱化。在此还可获得一个比较复杂的模拟连续屈服的节理模型,用 以模拟弱化为累积塑性剪切位移函数的连续变化特性。作为一个选择模型,还可获得 Barton-Bandis节理模型。节理模型和性质参数也可分别赋给单一节理或节理组。应当 注意,即使地质图上所显示的节理为直线段,节理的几何粗糙度也可以通过节理材料模 型加以表征。 UDEC 的块体可以是刚体或变形体。对于变形块体,开发了包括用于开挖模拟的空 模型(null)、应变硬化/软化的剪切屈服破坏模型以及非线性不可逆的剪切破坏和压缩 模型。因此,块体能被用来模拟回填、土体介质以及完整岩石。 UDEC 的基本公式假设为二维平面应变模型。此条件涉及断面保持为定值,并在平 行于该断面的平面上作用荷载的无限长结构。所以,非连续面也被假设为平面特性。另 外,UDEC 提供了一个平面应力问题的选择。对于平面应变分析,如果在垂直于平面方 向的应力 σ zz ,为最大或最小主应力,在垂直于平面方向, 块体可能出现塑性屈服, UDEC 的显式求解算法允许进行动态或静态分析。对于动态计算,用户指定的速度 或应力波可作为外部的边界条件或者内部激励直接输入到模型中。一个简单的动态波型 库也可以获取。UDEC 为动力分析设计了自由边界条件。 在静态分析中,包括了应力(力)和固定位移(速度为零)两种边界条件。边界条 件在不同的位置可以是不同的。同时,在 UDEC 中还可以获得边界元边界,用于模拟无 限弹性边界。也可以获得半平面解用来描述自由面效应。 UDEC 还能够模拟通过模型中的孔隙和不连续面的流体流动。在此认为块体是不可 渗透的。岩体的渗透率取决于节理的力学变形。也能够进行力学-流体全耦合分析。反 过来,节理水压也将影响力学特性。流体被处理为平行板的粘性流。 程序中的结构单元可用于模拟岩体加固和工程表面支护。加固包括端部锚固、全长 锚索和锚杆。表面支护模拟诸如喷射混凝土、混凝土衬砌和其他形式的隧道支护。 UDEC 包含一个强有力的程序语言,FISH,能够使用户定义新的变量和函数。FISH 是一个编辑器。通过 UDEC 数据文件进入程序被翻译并储存在内存中。

1.4 应用领域
UDEC 最初是为节理岩石边坡的稳定性分析开发的。对于块体不连续公式和运动方

-3-

程(包括惯性项)采用显式时间步求解方法???,便于块状岩体边坡的渐进破坏分析 和大变形运动研究。 UDEC 常用于采矿工程,已经进行了深部地下采矿洞室的静态与动态分析。洞室围 岩破坏诱发的断裂、滑移是用 UDEC 分析研究的实例之一。通过在模型的边界施加动应 力或速度波研究爆破影响。地震诱发的断层滑移也通过采用连续屈服节理模型进行了研 究。结构单元已经用于模拟全长岩锚和喷射混凝土的各种岩体加固系统。 UDEC 还应用于地下结构和深部高辐射废料的储存研究领域。通过应用热模型, UDEC 已经应用于模拟与核废料相关的热荷载效应。 UDEC 在作为一个计算设计工具,仍受到一定的限制。然而,程序较适用于研究节 理效应的潜在破坏机理。节理岩体特性是一个“有限数据系统”-即,在很大程度上内部 结构和应力状态是未知和不可知的。因此,建立一个完备的节理模型是不可能的。而 且,UDEC 是一个二维程序,除了特殊情况外,不可能表征具有三维结构的节理模型。 不过,应用 UDEC 程序,可以从现象学的角度研究节理岩体地下工程开挖响应。该方法 可加深岩石力学设计中各种不同现象的相互影响的理解。采用这种方法,工程师能够通 过识别地下工程可能产生不可接受的变形或加载导致的破坏机理,从而揭示工程所潜在 的诸多问题。 值得注意的是,UDEC 程序对于模拟颗粒流动或动态分析火山喷发是不适宜的。对 于该类研究,可以采用 PFC2D 程序。

2 开始启动
2.1 安装和启动程序
本节为首次使用 UDEC 的用户提供指导。如果你熟悉该程序仅仅是偶尔使用,你会 发现本节尤其是 2.6 节对于改变你原有印象是有帮助的。UDEC 程序共有 65 个主命令, 有接近 400 个关键词。 2.1.7 内存赋值 UDEC 自动调节内存大小达到 8MB。可以通过下列命令查询、改变内存值: Uedc m Uedc 14 Print mem 如果更多的内存可以获得,其内存能够通过应用环绕磁盘文件获得额外内存。表 2.2 给出了最大块体数与所需内存的关系。

-4-

表 2.年 2 RAM 与最大块体单元 RAM(MB) 2 4 8 16 最大刚性块体数 400 2500 7500 15000 最大变形块体数* 300 1000 3000 7000

* 假设每块体 8 个自由度。块体最大数随自由度的增加而减少。 2.1.9 运行 UDEC call file.dat 2.1.10 安装测试程序 有三个简单的数据文件,test1.dat 、test2.dat、test3.dat 用于程序测试。

2.2 简单演示-通用命令的应用
Block (0,0) (0,20) (20,20) (20,0) ; 产生一个块体 plot block 划分初始块体成小块体。 Crack (0,2) (20,8) Crack (5,3) (5,20) Crack (5,12) (20,18) 固定最下和最左块体,使之不可移动的命令如下: fix rang 0,20 0,5 fix rang 0,5 0,20 该命令固定形心处在 0<x<20,0<y<5 和 0<x<5,0<y<20 范围内的所有块体的当前 速度(即为零)。然后块体和节理所需的材料性质通过性质号予以赋值,即 prop mat=1 dens=2000 prop jmat=1 jkn=1.33e7 jks=1.33e7 jfric=20.0 对于该问题,所有的块体密度被指定为 2000kg/m3 。所有的节理切向刚度和法向刚 度分别被指定为 1.33e7 ,节理面的摩擦角为 20o 。下面将会发现,不同节理和块体可以 赋予不同性质参数。 其次,在 x 和 y 方向的重力加速度可以通过如下命令予以赋值: set gravity 0,-10.0 为了吸收振动能量,引入阻尼命令 damp local ; 显示该块体

-5-

上述命令是 UDEC 的缺省阻尼条件,因此,DAMP local 实际上并不需要。我们在 此仅仅是强调这是静态分析。 对于该点,问题很容易被执行。正如在后面能看到的,通过观测特定点的岩体运动 有助于进行工程特性判断。在该问题中,我们监测模型右角点 y 方向的速度,记录该运 动所采用的命令是: hist yvel (20,20) type 1 关键词 type 是在屏幕上以指定的间隔显示其值。 Step 100 ; 迭代次数

在计算过程中,当前的循环数,计算时间、最大不平衡力,在点(20,20)的 y 方 向速度以每间隔 10 次显示在屏幕上。 Plot hist 1 Title HEAD> ASIMPLE SLOPE STABILITY EXAMPLE EQUILIBRIUM STAGE Plot block Save slope.sav 通过最左边的块体来研究边坡的特性: delete range 0,5 0,20 命令 delete 将删除形心位于 0<x<5,0<y<20 范围内的所有块体。同时,采用 Step 或 cycle 命令继续进行计算。 Cycle 1000 Plot block velocity

2.3 概念与术语
UDEC 所涉及的一些术语大部分与其他应力分析程序类似。在 UDEC 模型中采用一 些特殊的术语来描述不连续面特征。按分类给出如下的基本定义。图 2.6 给出所给出的 术语定义。 UDEC MODEL -UDEC 模型:是用户为模拟实际的物理模型建立的。当称之为 UDEC 模型,就意味着为数值求解定义的求解条件的一系列命令。 BLOCK - 块体:是离散单元计算的基本单元体。通过切割一个块体成多个小的块 体产生 UDEC 模型。每一块可能与其他块体分离或通过界面力与其他块体相互作用的独 立块体。 CONTACT - 接触:每一块体通过点接触与相邻块体连接。接触可以认为是施加外 力到每一块体的边界条件。 DISCONTINUITY - 不连续面:是分离岩体成离散部分的地质特征。不连续面包 括岩体中的节理、裂隙、断层和其他不连续特征。

-6-

图 2.6 UDEC 模型的例子 DOMAIN - 区域:是指块体间的空洞或空间。DOMAIN 在 UDEC 模型中被处理 为实体。每一个 DOMAIN 是由两个或多个接触面确定的封闭区域。外 DOMAIN 是指围 绕 UDEC 模型的区域。 ZONE-单元:是由有限个单元组成的变形块体。在每一单元计算力学变化和温度 变化。在 UDEC 采用三角单元。 GRIDPOINT -结点(或节点):节点包括有限单元的角点。每一单元涉及三个节 点。一对 x 和 y 坐标定义每一个节点。因此确定了有限单元的精确位置。另一节点的术 语是 node。 MODEL BOUDARY – 模型的边界:是一个 UDEC 模型的周边。边界与模型的外 区域一致。内边界(即模型内的孔洞)也是模型的边界。每一内边界通过内区域定义。 BOUNDARY CONDITION - 边界条件:是约束或控制模型的边界(即对于力学 问题固定位移或外力)。 INITIAL CONDITION - 初始条件:模型受扰动(开挖)或加载(支护)前的原 岩应力状态。 NULL BLOCK-开挖块:表示模型中的空域(即材料不存在)。空块体可在后来 加以改变,例如,模拟回填(但一旦块体从模型中删除,就不可能恢复)。 STRUCTURAL ELEMENT - 结构单元:用来表征结构(如隧道衬砌、锚杆和锚 索)与岩体的相互作用的一维单元。结构单元也可以具有材料非线性。在大应变模型中 可以出现几何非线性。

-7-

STEP -求解(或迭代):尽管一个大的问题需要上万次计算才能达到稳定解,但 一般典型问题的求解需要 2000~4000 次循环,可以获得系统的平衡或稳态流。 STATIC SOLUTION - 静态解:当模型中动能的变化速率接近可以忽略的情况 时,UDEC 就认为达到了静态或拟静态解。 UNBALANCED FORCE - 不平衡力:表示当静力分析中的力所处于的不平衡状 态(即节理开始滑动或塑性流动)。 DYNAMIC SOLUTION - 动力解:尽管系统的缺省为静态求解过程,但可以进行 动态分析。对于动态分析,全运动方程(包括惯性项)可被求解。动能的消耗产生直接 影响。

2.4 UDEC 模型:初始块体的划分
UDEC 模型首先生成整个计算范围的单一块体。然后,通过用地质结构特征(如断 层、节理裂隙等)和工程结构(如地下洞室与隧道等)作为边界,切割该块体成小的块 体来考虑模型特征。 模型的所有块体都是通过块体质心和角点的坐标(x 和 y)确定。块体接触面以及变 形块体的节点也通过他们的坐标位置确定。产生模型包括由端点坐标(x,y 坐标)所定 义的线段(splits)切割模型块体。 UDEC 模型所有的条目(块体、角点、接触面、空区、节点和单元)都是通过位于 主数组中的地址编号,由 UDEC 自动的、唯一识别和确定。这些编码号也可以用作特殊 的单元。编码系统并不是顺序编码,所以用户必须通过绘图或打印加以识别。 例如,图 2.7 说明一个 UDEC 模型块体在 x 和 y 方向皆为 10 个单位(比如 10m)。模型通过一水平不连续面(x=0,y=5 to x=10,y=5)划分成两个块体。这两个块 体具有编号为 2 和 118。块体通过位于块体角点之间的接触面连接。接触号是 223 和 260。内部区域由两个接触面产生并由内部区域号 297 所识别。 显示在图 2.7 中的模型由列在下表中的命令产生。 block (0,0) (0,10) (10,10) (10,0) Crack (0,5) (10,5) plot hold block num cont num dnum 两块体的每一块可通过产生有限单元形成变形体。图 2.8 给出了上部块体划分为 8 个单 元和下部块体划分为 4 个单元的单元号。在两块体间产生了一个新的接触面(序号为 606)。位于块体棱上的任何节点总会产生。新的接触 606 对应于上部块体的棱产生的 新的角点。 Gen quad 11,6 range 0,10 0,5 Gen quad 10 Plot hold zone num cont num

-8-

图 2.7 UDEC 模型块体被划分成两个刚体

图 2.8 包含两个变形块体的 UDEC 模型

2.5 命令语法
UDEC 中所有命令都是面向单词,并由主要命令单词和随后的一个或多个关键词或

-9-

值构成。某些命令接受开关,即关键词修改命令的作用。每一命令都具有下列格式: COMMAND keyword value … < keyword value … > 在此,位于 < >内的参数为选择参数。而位于( )表示可以该参数为任意给定的值。 命令可依次写在命令行中。可能你已注意到,命令关键词仅前面几个字母为黑体。实际 输入时仅输入这些黑体字母就可由系统识别。

2.6 UDEC 应用基础
UDEC 是基于命令驱动格式。命令单词控制程序的运行。本节将提供给新用户一些基 本命令。为了建立一个 UDEC 模型进行模拟,必须考虑计算问题的基本成分: (1)由切割产生几何问题,由此建立离散单元块体模型; (2)本构特性和材料性质; (3)边界条件和初始条件。 块体模型定义问题的几何体。本构特性及所涉及的材料参数反映模型在受到干扰后 的力学响应。边界条件和初始应力定义了原岩状态,即在未受到扰动(开挖、支护、爆 破等)前的应力和位移状态。 在 UDEC 中定义了这些条件后,可以进行改变(即开挖材料或改变边界条件),从 而计算产生模型响应。像 UDEC 一类显式求解技术所获得问题的实际解与传统的隐式求 解方法的结果有所不同。UDEC 采用的是显式时间步求解代数方程,其解是在一系列计 算迭代后才获得。在 UDEC 中计算迭代步数可以通过用户控制。用户必须确定所进行的 求解迭代步数是否达到了实际问题的解。 图 2.9 给出了采用 UDCE 进行静态分析的求解一般过程。由于这求解程序符合实际 物理模型的生产工序和实际条件,因此其计算过程是方便的。采用上述过程进行简单的 应力分析的 UDEC 基本命令将叙述如下。

2.6.1 块体划分
UDEC 模型是通过切割初始的 UDEC 块体成小的块体代表模型的实际边界。采用下 述命令,建立模型块体。 Block x1,y1 x2,y2 x3,y3 在此,(x1,y1),(x2,y2),(x3,y3)… 是定义块体角点的坐标对。角点 必须按顺时针方向排列。角点应当与物理模型的边界条件一致。块体有很多角点,但通 常从 4 角点块体做起。 在 UDEC 中所有块体都有“圆角”,其目的在于避免块体悬挂在有棱角的节点上。由 于块体悬挂引起应力集中。然而,圆角值存在与模型有关的上限值。对于变形块体,最 大圆角长度应当不超过块体平均棱长的 1%。圆角长度可以如下命令加以改变: round d

-10-

在此,d 是圆角距离(缺省值是 d=0.5)。模型中的所有圆角长度都是相同的。 建议在 block 命令前指定圆角长度。在 block 命令后,键入 plot block 命令,就能够显 示圆角的效果。
建立模型 (1)生成模型块体,切割块体产生计算模型的几何体; (2)定义本构模型和材料参数; (3)指定边界条件和初始条件。

迭代计算使之平衡(模拟未扰动前状态)

检查模型响应

进行变化分析 例如: 地下开挖 改变边界条件

迭代计算

检查模型响应

重复其他变化进行响应分析

图 2.9 静态分析的一般过程 UDEC 有几个命令用于产生计算模型的几何体。生成地质结构(即节理)的两个主 要命令如下: Crack Jset Crack 命令用于产生块体中单一直线特征的裂缝。裂缝由端点坐标(x1,y1)和 (x2,y2)所确定。 Jset 命令则是自动节理组生成器。根据所给定的特征参数(即倾角、迹长、岩桥长 度、间距和空间位置)产生一组裂缝。 Crack 和 Jset 两个命令用于产生 UDEC 块体中的地质不连续面,即节理并不一定完

-11-

全将岩块切割成分离两个块体。然而,UDEC 需要连续断裂(即所有断裂都必须切割块 体)。由 crack 或 Jset 命令所产生的不连续面位置将被储存。为产生块体内连续断裂, 可采用 crack 命令产生的断裂。刚性块体在计算过程中,或变形块体在单元划分时,没 有连接形成完整块体的裂缝将被删除。 下面的例子说明用 Crack 和 Jset 命令切割块体。这两个命令的详细描述将在命令表 的第 1 节给出。节理生成器将在 3.2.2 节给予详细解释。 Example 2.4 产生简单的 UDEC 模型 round 0.1 block (0,0) (0,10) (10,10) (10,0) Crack (0,5) (10,5) 在这个最简单的模型中,切割块体涉及选择位置和指定裂缝。通过键入这些命令, 就可以产生 10×10 个单位的块体,然后劈裂成两个块体。CRACK 命令产生一个连续、 水平的贯通模型的裂缝。注意圆角的长度指定为 0.1。 通过键入如下命令,可以产生一个槽口: crack 2.5,10 5.0, 7.5 crack 5.0, 7.5 7.5,10 通过键入以下命令,就显示出包括块体地址号的块体图形。 Plot block num 通过采用 DELETE 命令,能从模型中删除一个块体。例如,为了删除槽口块体, 键入如下命令: delete range block 368 或 delete range 4.5 ,5.5 8,10 4.5< x <5.5 和 8 < y < 10 的范围必须包含被删除块体的形心。注意,当对模型进行 某些操作时,采用坐标范围是比较明智的。与问题相关的地址号有时发生变化。 对于JSET命令的参数需要 4 组数据对参数值。每一数据对中的第一个值是均值,而第二 个是对应于均值的最大均方差(相对于均匀概率分布形式)。第一组数据对是节理迹线 与x-坐标轴的正方向的夹角。第二对数据节理迹线长度;第三组数据是不连续节理的岩 桥长度;第四组数据是节理间距。还有一些选择参数,可以用于产生一组比较复杂的节 理模式。JSET命令的一般应用在 3.2.2 节加以讨论。采用JSET命令产生两组节理组的应 用在例 2.5 中得到说明。 例 2.5 两组连续节理组的产生 new round 0.01 block (0,0) (0,20) (20,20) (20,0)

-12-

jset (45,0) (5,0.5) (0.5,0) (2,0) jset (-10,0) (5,0.5) (0.2,0) (1.5,0) 在上例所生成的节理图如图 2.11 所示。第一个 JSET 命令产生一组与 x 轴方向夹角 为 45o 具有间距为 2 个单位的连续节理。第二个 JSET 命令产生与 x 方向夹角为 -10 o 、间 距为 1.5 个单位的连续节理。圆角长度的选择可能影响节理组的产状。节理的位置可能 由于当棱长小于 2 倍的圆角长度不能产生块体而可能改变节理的位置。 如果在 JSET 命令前,增大圆角的长度(比如说 0.1),模型中某些节理的位置将发 生改变。用 JSET 命令产生节理可能涉及某些试错法。第 3 节给出了进行产生节理过程 的建议。

图 2.11 产生的两组连续节理模型 当产生大小悬殊的块体时,建议从模型中删除较小块体,以提高模型的计算效率。 在例 2.5 中,块体尺寸的变化范围从 1.751×10 Print max 键入如下命令,删除极小块体: delete range area 3e-2 所有面积小于 3×10-2 的块体都从模型中删除。通常,将小于最大块体的 1%左右小 块被删除后对计算结果的影响并不显著。 最后,注意到 NEW 命令用在第二个例子,以便允许开始一个新的模型。当切割块 体(尤其当采用 JSET)时,一个重要问题是综合考虑块体数与计算速度的协调。计算速 度与模型的块体数(或变形体单元数)成函数关系。根据经验,模型具有大约 1200 刚
-3

到 3.679,可以由以下命令查找:

-13-

体(或具有 8 自由度的 500 变形体)进行 2000~4000 迭代步就能获得静态问题的解。 对于 90MHz 的微机,对于 500 个变形块体模型允许 4000 步大约需要 10 分钟。根据你 的计算机的计算速度,可以估算出一个模型所需的计算时间。 通过切割 UDEC 块体形成工程结构形状,这须在进行工程开挖前实施。通常采用三 个命令来产生形状: crack tunnel arc 前面已经给予介绍 CRACK 命令。TUNNEL 命令产生圆形形状。该圆由用户指定 的裂缝段数构成。ARC 命令由用户指定的角度,产生弧形断裂模型。可以结合这些命令 产生各种形状的 UDEC 块体。例 2.6 给出的命令产生断层切割一个圆形隧道的模型。 EXAMPLE 2.6 断层切割一个圆形隧道 New Round 0.1 Block -10,-10 -10,10 10,10 10,-10 Tunn 0,0 2 16 Crack -5,10 5,-10 Plo hold block num 所生成的模型如图 2.12 所示。圆形隧道的圆形坐标(0,0)、半径为 2 和划分成 16 个裂缝段。由于隧道全部处于块体内部,所以仅用 TUNNEL 命令不能产生独立的块 体 。 必 须 采 用 CRACK 切 割 模 型 块 体 的 边 从 而 产 生 新 的 块 体 。 如 果 用 户 运 行 仅 用 TUNNEL 命令所产生的模型,则隧道裂缝在运行前被删除。通过引入 CRACK 命令,连 接隧道裂缝延伸到模型外边界从而形成连续的裂缝,因此,形成有隧道和断层构成的块 体(如图 2.12 所示)。 应当注意,裂缝并不贯穿隧道的周边。如果 TUNNEL 命令先给出,随后的 CRACK 或 JSET 命令并不贯穿隧道。首先应用 TUNNEL 是较为方便的,因为隧道开挖仅涉及删 除一个块体,即 delete range block 1920 或 delete range -1,1 -1,1 将模拟圆形隧道的开挖。

-14-

图 2.12 断层切割圆形隧道 例 2.7 给出了一条断层切割一个马蹄形隧道: New Round 0.1 Block -10,-10 -10,15 10,15 10,-10 arc 0,5 2,5 180 8 Crack -2,0 -2,5 Crack -2,0 2,0 Crack 2,0 2,5 Crack -5,15 5,-10 隧道的形状如图 2.13 所示。隧道顶弧的圆心在(0,5),起始点在(2,5)和 180 的圆心角,逆时针画圆,划分成 8 段。前三个 CRACK 命令产生隧道的边墙和底 板,后一个 CRACK 产生一条切割隧道的断层。开挖隧道也可以通过删除块体编号或块 体形心位置来实现。另一选择就是指定包含该隧道块体的图象窗口,然后删除窗口中的 块体。这可以用下面命令实现: window -2,2 0,7 delete range window
o

-15-

图 2.13 断层与马蹄形隧道相交

2.6.2 指定材料模型
2.6.2.1 块体模型 一旦完成块体切割,必须对所有的块体和不连续面指定材料特性。缺省为所有的块 体皆为刚体。在多数分析中,块体应为变形体。仅仅在应力水平较低或岩块材料具有高 强度和低变形的情况才能够应用刚性块体的假设。 块体的变形特征通过以下命令定义: gen edge v 或 gen quad v GEN 命令激活三角形网格有限单元自动生成器。命令GEN edge v 将作用于任意形状 的块体。其v值定义三角形单元的最大边长,即v值越小,块体中的单元越小。应当注意 的是:具有高的边长比值的块体并不能产生单元,其极限的比重近似为 1:10。 通过 Plot zone 检查模型单元。 采用命令GEN quad v,指定模型为塑性材料模型的单元。该类型的单元提供了对于 塑性问题的精确解。然而,GEN quad 命令可能对某些形状的块体不起作用。在此情况 下,应当采用GEN edge 。 在 UDEC 中为变形块体(单元)开发了 7 种材料模型。对大部分用户,最常用的三 种模型如下: change cons=0 ; null model

-16-

change cons=1; elastic model change cons=3; Mohr-coulomb model CHANGE 命令改变块体为指定的变形块体。Cons=0 意味着模型块体材料被移出或 开挖。这允许用户改变块体在以后的某些阶段返回为弹性或弹塑性材料。如果块体被删 除,则以后计算阶段不可再恢复。 Cons=1 改变块体为各向同性弹性特性;而cons=3 则改变块体为摩尔-库仑模型,考 虑塑性特性。缺省值为所有变形体则自动改变为cons=1。 块体改变为cons=1 和cons=3 必须提供PROPERTY mat 命令给块体赋予材料参数 值。注意性质参数不要赋给特定的块体,而是赋给材料号。材料参数可以赋值给多达 50 种材料号。然后,材料号再赋给具有CHANGE mat 命令的块体。 对于弹性模型,需要的性质为: (1)密度 (2)体积模量 (3)剪切模量 注意:体积模量 K、剪切模量 G 与杨氏模量 E、泊松比 ν 之间的关系如下:

K=


E E , G= ( -2ν) 31 ( + ν) 21

E=

9 KG 3K-2G , ν= 3K+G (3K+G) 2

对于摩尔-库仑塑性模型,需要的性质为: (1)密度 (2)体积模量 (3)剪切模量 (4)内摩擦角 (5)粘聚力 (6)剪胀角; (7)抗拉强度 如果上述参数没有赋值,系统自动赋零值。 对于 UDEC 程序,对于上述两个模型,密度、体积模量和剪切模量必须赋予正值。 2.6.2.2 节理模型 除了给块体赋予材料模型外,还应对模型中的所有不连续面(即接触面)赋予材料 模型。对于不连续面,有四种本构模型。在 UDEC 中开发了四种节理本构模型。但对于 大部分模型分析,最适宜的模型有库仑滑动模型(完全弹塑性),可通过如下模量赋予

-17-

不连续结构面: change jcons=2 所有不连续结构面的缺省模型是Jcons=2。 节理材料模型也是通过 PROPERTY mat 命令赋予材料性质参数。如同块体,参数 并不是直接赋给不连续面,而是材料号。材料号是通过 CHANGE mat 命令赋给节理。 对于库仑滑动模型,所需要的参数是: (1)法向刚度 (2)切向刚度 (3)内摩擦角 (4)粘聚力 (5)剪胀角; (6)抗拉强度 如果所有的参数没有给赋值,他们的缺省值为零。在UDEC中,必须给节理的法向和切 向刚度赋值,并为正值。例 2.8 演示了材料模型的应用。 例 2.8 指定材料模型与性质参数 New Round 0.1 Block -10,-10 -10,10 10,10 10,-10 tunnel 0,0 2,16 jset -70,0 40,0 0,0 40,0 -1,-1 jset -50,0 40,0 0,0 3,0 gen edge 2.0 change jmat=2 range angle -51,-49 change jmat=5 range angle -71,-69 pro mat=1 d=2500 b=1.5e9 s=0.6e9 pro jmat=1 jkn=2e9 jks=2e9 jcoh=1e10 jten=1e10 pro jmat=2 jkn=2e9 jks=1e9 jfr=45 pro jmat=5 jkn=2e9 jks=1e9 jfr=10 change cons=0 range -1,1 -1,1 在上面的例子中,一条与 x 轴方向成-70o 单一断层切割圆形隧道。还有倾角为- 50o 和间距为 3m 的节理组。块体是变形的和具有最大边长是 2m 的三角形单元。图 2.14 显示了隧道、节理和单元。块体是弹性的,其性质参数是通过 PROP mat=1 赋值。50 o 的节理组是通过性质参数号 2 赋值。70o 的断层参数是由材料参数号 5 赋值。jmat=2, angle=-51, -49 仅将性质参数赋予-51o 和-49 o 之间。jmat=5,angle=-71, -69 仅将其值赋 0, 2

-18-

予范围角位于-71 o 和-69o 之间。所有其他不连续面被赋予材料号 1,这意味着将具有高 粘聚力和抗拉强度的粘结效应。这些不连续面是“虚拟的”节理和对应于后来将隧道开 挖。虚拟节理不发生滑动或张开。可用以下命令检查材料号。 Plot block mat 隧道块体被改变为 cons=0 是模拟开挖。 Plot zone 可得到下图。

图 2.14 圆形隧道与 70o 断层和 50o 节理组构成的模型

2.6.3 施加边界条件和初始条件
在完成所有块体切割(节理切割)和变形单元划分之后,应施加边界条件和初始条 件。施加力学边界条件通常采用BOUNDARY命令。该命令用来指定力、应力和速度 (位移)边界条件。边界力和应力能够施加到刚体和变形体的边界上,但速度(位移) 边界仅适用于变形块体(见施加刚性块体的边界命令FIX,FREE和LOAD)。表 2.4 提 供了边界条件命令的总结和效果。 BOUNDARY xload 和 yload 命令施加 x-和 y-方向的分力到边界角点。BOUNDARY stress 施加应力张量到边界上。BOUNDARY xvel 和 yvel 在选择的边界结点施加 x-和 y方向的速度。 注意的是:应用 BOUNDARY 命令所产生的条件或约束将不发生改变(除非用户再

-19-

次改变)。

表 2.4 边界条件命令总结 命 令 BOUNDARY Stress Xload Yload Xvel Yvel FIX FREE LOAD Xload Yload 效 果

施加总应力到刚体或变形体块体的边界上 施加刚体或变形体边界的 x 方向的荷载 施加刚体或变形体边界的 y 方向的荷载 施加变形体边界的 x 方向的速度(位移) 施加变形体边界的 y 方向的速度(位移) 固定刚体边界的速度(位移) 释放刚体的速度(位移) 施加 x 方向的荷载到刚体的边界 施加 y 方向的荷载到刚体的边界

初始应力条件能够被指定到所有的变形单元和所有刚体或变形体之间节理的法向应 力和剪切应力。INSITU 命令用来初始化应力。采用该命令,可以赋值初始应力。边界 的初始条件应用于例题 2.8,其命令应用见例 2.9。 例 2.9 施加的边界条件和初始条件 boundary stress -10e6,0,0 range -11,-9 -10,10 boundary stress -10e6,0,0 range 9, 11 -10,10 boundary stress -5e6,0,0 range -10,10 9,11 boundary yvel =0.0 range -10,10 -11, -9

insitu stress -10e6,0,-5e6 szz -4.8e6 10MPa压应力被作用到模型的左、右边界的x方向上。5MPa的压应力(负号为压) 被施加到上部边界的y方向。底部边界的y方向的运动被固定。固定边界的位移对于考虑 重力的情况下是尤其重要的。注意:应力边界影响所有的自由度。因此,在施加速度边 界条件前的同一边界,应当施加应力边界条件,否则,所指定的速度约束将不起作用。 而且,四个BOUND命令的每一个都应给出x和y的坐标范围。应当提醒的是为了确保通 过BOUND命令所产生的响应完全位于指定的范围。键入如下命令: Print bound 和 Plot bound xcond Plot bound ycond 来检查边界条件。

-20-

命令 INSITU 在 x 方向始化所有的应力为-10MPa 和在 y 方向初始化应力为-5MPa。 在平面之外的 z 方向也给予初值 σ zz =-4.8MPa。对于弹性块体分析,z 方向的应力未初始 化并不影响平面应变问题的解。然而,对于塑性分析,z 方向的应力可能影响到破坏状 态,因此,应当慎重选择 σ zz 的应力初始化。

2.6.4 迭代为初始平衡
UDEC 模型在进行开挖模拟前必须进行初始状态的平衡计算。施加合适的边界条件 和初始条件,使模型与初始的平衡状态相吻合。然而,对于复杂的几何形状和多介质材 料的情况,在给定的边界条件和初始条件下,进行计算获得平衡是十分必要的。其计算 可采用 STEP(或 CYCLE 或 SOLVE)命令。借助于 STEP 命令,为达到模型的平衡, 用户指定循环步进行计算。当每一刚体形心的结点力 或变形体的结点力接近零,模型就 处于平衡状态。当激活 STEP 命令后,最大的结点力矢量(称之为不平衡力)由 UDEC 进行监测,并在屏幕上显示。在此用户能够估计模型何时达到平衡状态。 对于任何模型的数值分析,不平衡力不可能完全达到零。但当最大的结点不平衡力 与初始所施加的总的力比较相对较小时,就可认为模型达到平衡状态。例如,如果最大 不平衡力从最初的 1MN 降低到 100N,此时(最大不平衡力与初始的不平衡力之比为 0.01%),则认为模型达到平衡。 采用 UDEC 进行数值分析判断模型平衡是一个重要的问题。用户必须确定模型在 何时达到平衡状态(即问题的解)。在 UDEC 中设置有一些特征,用于支持这种决策。 如记录最大不平衡力历史是其中之一: hist unbal 除此之外,还有速度历史(即某一结点的速度或位移),命令如下: hist xvel 5, 5 hist ydisp 0, 11 第一个是记录位移坐标(x=5, y=5)附近结点 x 方向的速度,而第二个是记录接近 坐标(x=0, y=11)位置处 y 方向的位移。在进行数百次或数千次迭代后,这些历史记录 将绘图和显示其平衡条件。 图 2.15 所示的模型为节理岩体中开挖一矩形区域。在此岩体中存在一组节理和一条断 层。模型受 10MPa 的静水应力场和重力作用。重力是用以下命令: set grav 0.0 , -9.81 第一个是 x 方向的加速度,第二个值为 y 方向的加速度为 9.81m/sec 2 (向下作用)。 当考虑重力所引起的应力变化不大可以忽略重力的作用。在例 2.10 中,尽管考虑重 力,有助于识别围绕洞室周围的松散块体受重力的作用,但是考虑重力引起的应力变化 小于 0.5MPa,相对于 10MPa 的原岩应力可以忽略不予考虑。 ro 0.1

-21-

bl -10,-10 -10,10 10,10 10,-10 cr -2,-2 -2,2 cr -2,2 2,2 cr 2,2 2,-2 cr 2,-2 -2,-2 jset 70,0 40,0 0,0 40,0 -2,0 jset -50,0 40,0 0,0 3,0 1,2 gen edge 2.0 chan jmat=2 range angle -51,-49 chan jmat=5 range angle 69,71 prop mat=1 d=2500 b=1.5e9 s=.6e9 prop jmat=1 jkn=2e9 jks=2e9 jcoh=1e10 jten=1e10 prop jmat=2 jkn=2e9 jks=1e9 jfr=45 prop jmat=5 jkn=2e9 jks=1e9 jfr=5 bound stress 0,0,-10e6 range -10,10 9,11 bound xvel=0.0 range -11,-9 -10,10 bound xvel=0.0 range 9,11 -10,10 bound yvel=0.0 range -10,10 -11,-9 insitu stress -10e6,0,-10e6 szz -10.0e6 set grav 0.0 -9.81 hist unbal hist ydis 0,2 ; solve for 10 step 700 save fall1.sav 如果位于开 挖体顶板的 块体被分离 ,则由于重 力作用将落 到洞室。这 将在后面的 2.6.5 节中详细说明。应当注意,当重力应力与原岩应力具有相同的量级,则用 INSITU 命令施加应力梯度(应力随高度的变化)以加速初始平衡状态的收敛。

-22-

图 2.15 节理岩体中的矩形开挖体 初始不平衡力近似为 2.0MPa。在进行 700 步后,降至约为 10N。通过绘制两个历 史可以发现,最大的不平衡力已接近零,位移接近 2.7×10 -3 m 常量(图 2.16 和 2.17)。

图 2.16 最大不平衡力历史

-23-

图 2.17 在(0,2)位置的 y 位移历史 模型的不平衡力在接触力和块体角点与结点的力略有差异,这与圆形角点有关。圆 形在块体角点产生很小的“洞”。为不平衡力的平衡增加某些迭代步是必要的。 设置不平衡力的值: solve force = f 在此,f 为用户定义的不平衡力值。 检查不连续面的破坏条件是十分重要的。对于这个问题,初始 x 发现和 y 方向的应 力分量是相同的,因此,模型中的节理不能滑动,这可以用以下命令加以验证: Plot bou slip 所显示的图形是由模型外边界和满足库仑滑动准则的节理组成。对于选择的模型条件, 有可能出现在原岩应力状态下节理发生滑动。例如,改变水平应力分量 σ zz 到-5MPa 并 回到例 2.10。则断层将沿着整个长度上发生滑移。如何识别初始应力在平衡计算过程中 节理发生的滑移,用户应当重新评估所选择的原岩应力参数和不连续面强度。沿着节理 长度滑移的模型表明该模型并非可靠。 在模拟开挖前确保模型处于平衡状态是十分重要的。通过记录几种历史以考察最大 的不平衡力的衰减。如果所进行的计算步超过模型达到平衡所需的计算步,并不会影响 计算结果。然而,如果不充分的计算步将影响模型的计算结果。 UDEC 计算可在任何时间通过按<Esc>被中断。更方便的是使用 STEP 命令进行高 次数的计算和周期的中断和再次分析,以确保达到平衡状态。

-24-

2.6.5 进行改变和分析
UDEC 允许在求解过程中的任意部位改变模型条件。这些变化可能具有以下形式: (1)开挖材料; (2)增加或删除边界荷载或应力; (3)固定或释放边界结点的速度(位移); (4)改变材料模型或块体和变形体的性质参数。 可以用 DELETE 命令或 CHANGE cons=0 命令模拟材料开挖。用 BOUNDARY xload,yload 或 stress 命令施加荷载和应力。通过采用 BOUNDARY xvel 或 yvel 命令固 定边界角点。通过 BOUNDARY xfree 和 yfree 命令移去边界约束。用 CHANGE 命令改 变变形块体和不连续面的材料模型。而用 PROPERTY 命令可改变材料性质参数。 很显然,几种命令可以重复应用,进行各种模型的改变。例如,从初始平衡状态, 应用这些命令继续例 2.10 获得例 2.11。 例 2.11 开挖隧道和监测其响应 rest fall1.sav delete -2,2 -2,2 reset disp reset hist hist unbal hist ydis 0,2 step 2000 plot blo stress disp save fall2.sav 由于采用 DELETE 命令,模拟开挖矩形洞室,导致模型应力的变化。结点位移与 历史记录被重新设定,仅由开挖所引起位移变化被监测。建议在块体被删除后重新设定 历史位置。 在矩形洞室开挖后产生很高的不平衡力,因此需要进行计算使之获得重新平衡。然 而,在此情况下,没有观测到不平衡力接近很小的值,而处于其值为 0.017MN 的常值。 进而 y 位移历史记录也显示在计算 2000 步后,在位置(0,2)的位移仍向下运动。开 挖顶板上的块体已经从围岩脱离和掉落到洞室内。这由图 2.18 清楚地看出。不平衡力不 可能接近于零,因为顶板块体应自由下落。 如果预计模型变化将导致破坏(即力的平衡条件不能获得)就不要用 SOLVE 命令求 解。

-25-

图 2.18 洞室顶板块体发生冒落

2.6.6 保存或恢复计算状态
当进行分步计算时,另外两个命令 SAVE 和 RESTORE 是有用的。在一个阶段的结 尾(即初始平衡),采用如下命令,可以保存模型状态。 Save file.sav 式中,file.sav 是一个用户定义的文件名。扩展名.sav 定义这个文件是一个保存文件。这 个文件可以采用如下命令进行恢复: rest file.sav

2.6.7 简单分析的总结
在表 2.5 中给出了本节所介绍的主要命令。更常见的情况求解问题需要从进行一个简 单的问题计算开始。

-26-

表 2.5 简单问题分析的基本命令 功 产生块体模型 能 ROUND BLOCK 切割块体 CRACK JSET TUNNEL ARC 块体和节理的材料模型和参数 GEN CHANGE PROPERTY 边界条件和初始条件 BOUNDARY INSITU 初始平衡(具有重力) DAMP local SET gravity STEP SOLVE 模型变化 DELETE CHANGE PROPERTY BOUNDARY CABLE 监测模型响应 HISTORY PLOT 保存或恢复当前状态 SAVE RESTORE 命 令

2.8 系统单位
表 2.6 系统单位-力学参数 SI Length Density Force Stress Gravity m kg/m N Pa m/sec
2 3 3

m 10 kg/m kN kPa m/sec
2 3 6

M 10 kg/m MN MPa m/sec
2 3

cm 10 g/m3 Mdynes Bar cm/sec 2
6

-27-

3 用 UDEC 求解问题
本章为应用 UDEC 求解岩石力学工程问题提供指导。在 3.1 节给出地质力学分析阶 段的建议。在 3.2 节通过例 3.10 明确了模型准则和求解过程必须考虑的问题,涉及的问 题如下: (1)产生模型(3.2 节); (2)刚体或变形块体的选择(3.3 节); (3)边界和初始条件(3.4 和 3.8 节) (4)加载和模拟顺序(3.6 节); (5)块体和节理模型和材料参数的选择(3.7 和 3.8 节); (6)改进模型效率的方法(3.9 节); (7)计算结果的解释(3.10 节)。 最后,地质力学领域的模拟原理参见 3.11 节。在该领域进行模型分析的新手可能希 望首先咨询该节。地质力学模拟方法与其他工程领域,如结构工程存在很大的不同。进 行地质力学分析始终记住这一点是十分重要的。

3.1 一般性研究
模拟地质工程的过程涉及一些特殊的考虑,其设计方法也与其他人工材料结构不 同。在岩土体上建造结构或在其中开挖分析与设计,必然面对相对少的现场数据以及材 料的变形和强度性质参数存在较大变化的情况。获得岩土工程现场完整的现场资料是不 可能的。例如,原岩应力、材料性质和不连续面特性等信息仅是部分的。 由于为设计预测所输入的必要信息是有限的,所以,地质力学数值模型主要用于理 解影响系统特征的力学机理。一旦掌握了系统的特性,然后,为工程设计过程探索一些 简单的计算。 面对地质工程研究总是缺少满意的数据以及缺乏对材料性质的充分理解,而在其他 领域,在具有充足的数据的情况下,应用 UDEC 可直接用于工程设计。当获得合理的数 据应用此程序时总能获得合理的结果。应当认识到,如图 3.1 所示的应用的过度阶段。

复杂的地质条件; 典型情况 不可获取的数据; 无试验经费 数 据 研 究 无 机理 研究 通过参数研究 现场特性分类

简单的地质条件 投入资金 进行现场研究 资料充分 预 测 (直接用以设计)

图 3.1 模型研究图谱

-28-

UDEC 程序可用于模型特性的预测(如图 3.1 的右边),或仅作为“数值试验”来测 试一些设想(图 3.1 的左边)。正是现场资料(和资金)而不是程序决定了应用情况。 如果具有高质量的和足够的数据,UDEC 就能够给出好的预测。 由于大部分的 UDEC 分析是处于较少数据的情况下进行的,所以,本节将探讨类似 于试验研究的数值模型研究技术。数值模型决不应当被认为是一个一端接受信息而另一 端输出结果的“黑箱”。为了获得可以合理的解释,必须十分注意准备数值“样本”和多样 本“试验”。表 3.1 列出了进行成功的数值模拟试验的建议步骤。下面分别讨论: 表 3.1 地质力学问题的数值分析步骤
第 1步 第2步 第3步 第4步 第5步 第6步 第7步 定义模型分析的对象 产生模型系统的概念图形 建造和允许简单的理想模型 搜集模型所需的计算数据 准备一系列用于分析的详细模型 进行模型计算 提供结果和解释

3.1.1 第 1 步:定义分析模型的对象
一个分析模型所了解内容与深入程度常常取决于分析的目的。例如,如果是为解释系 统的特性所提出的两种相互冲突机理的决策,此时可建造一个较粗糙的模型,用于两种 机理的研究。如果试图涉及存在于实际模型中的复杂条件,然而,如果可能对模型的响 应产生微不足道的影响或与模型计算的目的毫不相关的计算特征可以被忽略。

3.1.2 第 2 步:产生物理系统的概念图形
重要的是构思出实际问题的图形,便于初步估计在所施加的条件下,预测系统的基 本特性。当准备这个图形时,应当回答几个问题。例如,该系统是否稳定?主要力学响 应是线性还是非线性?是否存在可能影响特性的不连续面?是否存在地下水的影响?实 际的系统物理结构是否还存在其他几何问题? 这些考虑将表征了诸如模型的几何形状、块体材料模型、边界条件以及初始平衡条 件等数值模型的总体特征。这将决定是否采用三维模型或二维模型?

3.1.3 第 3 步:建造和运行简单的理想模型
当为数值分析理想化一个物理系统时,较有效的方法是在构筑详细的模型之前,首 先建造和运行一个简单的测试模型。为产生数据和对问题的理解,应在尽可能早的阶段

-29-

产生这样的一个简单模型。其结果可用于进一步系统的概念图形。在简单模型运行后可 能需要重复第 2 步。 简单的模型能揭示一些问题,以便在进行深入分析之前加以修正。例如,所选择的 材料模型是否能够代表所期望的特性?边界条件是否影响模型的响应?基于简单模型的 计算结果能够有助于指导对分析起重要影响作用的数据研究方案。

3.1.4 第 4 步:综合特定问题的数据
对于一个模型分析所需的数据类型包括: 详细的几何参数(即地下洞室形状、地表形态、坝形状,岩石或土体结构); 地质结构的位置(即断层、层理,节理组等); 材料特性(即弹性或塑性性质,峰后特性); 初始条件(即原岩应力状态,孔隙压力,饱和度); 外部加载(即爆破荷载、洞壁压力)。 由于分析所涉及的条件(尤其应力状态、变形和强度性质)存在很大程度的不确定 性,为研究必须选择参数的合理变化范围。基于简单模型的计算结果(第 3 步)常常能 够有助于确定变化范围。

3.1.5 第 5 步:准备一系列详细的运行模型
通常数值分析用于研究不同的破坏机理、研究一定范围变化的计算参数的系列分析 。当为计算准备一系列计算模型时,应考虑如下一些方面的问题: (1) 每一个计算需要花费多少时间?如果模型运行的时间过长,为达到有用的结 论所需要获得足够的信息可能是困难的。为缩短计算时间,可以考虑在多个 计算机上,运行参数变化的计算。 (2) 应考虑保存所需要的模型在计算过程中的中间状态,以便每一参数的变化不 必重复计算。例如,如果分析几个加载或卸载阶段,用户应当能够返回运行 已经进行的任一阶段,以便改变一个参数后从那一阶段继续计算。 (3) 在模型中是否设置足够的监测位置(历史记录),为进行清楚地解释模型计 算结果和不同计算参数的比较分析提供足够的信息?在模型中设置几个参数 变化的监测点,对计算过程监测是有帮助的。尤其模型中的最大不平衡力应 当被监测,以便检查在分析的每一阶段的平衡或破坏。

3.1.6 第 6 步:进行模型计算
在进行一系列模型分析之前,最好首先选择一个或两个模型进行详细地分析。这些 运行应当随时被中断,确保达到预期的效果。一旦能够确信模型的计算是正确的,几个 模型数据文件被联系在一起,进行一系列模型的连续分析。在连续运行的任何时间,应 有可能中断计算,查看结果,然后继续或修改模型。

-30-

3.1.7 第 7 步:提供结果和解释
求解的最后一步是为进行清楚地解释分析提供计算结果。最好是通过直接在屏幕上 显示或输出的图形结果。图形结果应当提供便于进行计算与现场观测结果的比较方式。 图形应当能够清楚地分析所感兴趣的区域,例如应力集中位置、模型中稳定与不稳定区 域。模型中任何变量的数值也能够容易获得,为详细解释模型的响应。 为有效地进行地质工程问题分析,我们建议了上节介绍应遵循的七个步骤。下面章 节将论述在 UDEC 的应用中为模型研究中涉及到的这些步骤中的每一方面所涉及的特定 问题。

3.2 产生模型
UDEC 程序在产生几何模型的方式与传统的数值分析程序有所不同。首先产生计算 范围的单一块体。然后,这个块体被切割成小的块体。模型中块体的边界是地质结构面 或工程结构(如开挖体边界)。这种切割处理被称之为节理生产的几何体。然而,“节 理”代表物理模型中的实际地质结构和人造结构边界或将被移去或在以后连续的计算步中 改变材料。对于后者,节理是虚拟的,其存在不应影响模型的计算结果。虚拟节理的表 征将在 3.2.3 节中讨论。

3.2.1 确定 UDEC 模型合适的计算范围
UDEC 几何模型必须具有足够大的范围,在感兴趣的区域内,包含主要的地质结构 特征,由此代表真实的实际的物理问题。考虑的方面如下: (1) 处于何处的地质结构(即断层、节理和层面)应详细描述? (2) 模型边界的位置对模型的影响程度如何? (3) 如果应用变形块体,在关心的区域,何种密度的单元可满足问题的精度? 上述三个方面决定了实际分析 UDEC 模型的规模。 如果仅有很少的地质结构(即两个或三个切割断层或遍布空间的节理组),这些可 能通过 CRACK 或 JSET 命令单独输入。记住 UDEC 是一个二维程序。除了特殊情况, 忽略三维效应。如果地质结构不能用垂直于分析平面的二维特征加以表征,则可能需要 采用三维分析(例如 ITASCA 程序 3DEC)。 单个独立特征的断裂可用一种或两种方式进行输入。一种是用 CRACK 命令,他给 出断裂的两个端点;另一种是用 JSET 命令,给出断裂的倾角和断裂通过的位置。例 如,或者 Crack (0,0) (10,10) 或 Jset (45,0) (20,0) (0,0) (100,0) (5,5)

-31-

都可用来定义一个倾角为 45 o ,且通过坐标点(x=5, y=5)的节理。为了在 UDEC 计算 中得到认可,节理必须是连续的(即完全劈裂块体成两个部分)。然而,节理可以由分 段连续、具有不同角点的分段构成。CRACK 和 JSET 产生不连续节理段(SPLIT 命令 与 CRACK 具有相同的形式,但他不能产生不连续节理)。例如,燕尾状节理可以用例 3.1 中的命令生成。 例 3.1 燕尾状节理 Round 0.1 Block 0,0 0,10 10,10 10,0 cr 0,5 2.5,6 cr 2.5,6 5,5 cr 5,5 7.5,6 cr 7.5,6 10,5 在上述例子中,几个CRACK命令的顺序没有限制。如果被后来产生的裂缝交切,内 部裂缝被保存在一个临时性的文件中和后来被应用。在节理生成过程中,任何内部或部 分贯穿裂缝,在模型运行过程中被删除。当块体被赋予可变形的(GEN),内部或部分 裂缝也自动被删除。也可采用JDELETE命令实施人工删除。 对于块体作指定的圆角长度可能局部影响节理的产生。最小块体棱长定义为圆角长 度的两倍。因此,节理段可能背离这个准则。例如,例 3.2 显示了指定一个裂缝的端点 位置处于角点圆角长度的两倍,该裂缝通过角点的位置坐标发生错位。UDEC 并对任何 这样错位不会给出任何警告,所以用户应当通过采用 PLOT block 命令随时进行检查。 例 3.2 圆角长度对产生裂缝的影响 New Round 0.2 ; Rou 0.1 Block 0,0 0,10 10,10 10,0 cr 0.3,0 9.7,10 如果圆角长度小到 0.1,则裂缝将处在指定的位置。 SET edge 命令运行用户人为定义的最小块体棱长。基于这个命令,用户为精确求解 而避免块体具有小的棱长,或反之具有很大的边长比,能够设定一个较小的圆角长度。 例如,如果设置命令 SET edge 0.4 和 ROUND 0.1,则不会产生小于 0.4 的块体棱长和 块体的圆角长度将为 0.1。这些命令必须在 BLOCK 命令之前给出。 模型中用于描述地质特征的节理数(即块体数)存在一个限制。这涉及模型的范围 和块体的单元数(如果采用变形单元)。实际的限制依赖于如同表 2.2 所示的可利用计 算机内存。在进行节理生成时必须考虑这个限制。根据经验,最后的方式总是从较少节

-32-

理开始,然后,如果有必要再逐渐增加节理来达到预期的效果。应当避免试图建立一个 复杂节理模型的诱惑,关于此问题将在 3.11 节中进一步讨论。 节理自动生成器可以在 UDEC 中得到,可根据实际测量的参数(即节理倾角、间 距、长度和岩桥长度等)产生节理组。JSET 命令也起到节理生成器的作用。如在第 3.2.2 节的例子所描述的那样。在 UDEC 模型还可获得一种特殊的生成 Voroni 形状多边 形的节理生成器。在 3.2.2 节给出该程序的应用实例。高级用户还可编写自己的节理生 成器。这可通过将节理生成器所定义节理段端点坐标(x,y)对组成的表列来实现。 FISH 函数可以完成这个自动处理过程。这个列表能够被 UDEC 直接读取。另外,图形 数字化仪也可以用来产生节理端点的坐标对。 请记住,在 UDEC 内的节理是作为直线段显示。许多线段可能需要用不规则的节理 形态来拟合。模拟者必须决定用 UDEC 模拟实际节理模式的水平。几何的不规则对节理 的影响也可以通过节理材料模型性质(即变化节理面的性质参数)加以考虑。 如上所述,模型边界必须具有足够的远,以致模型对边界不产生影响。第 3.4 节论 述了边界的影响结果。一般地,对于单一地下开挖工程,边界离开挖边界的距离应当大 于开挖跨度的 5 倍左右。 然而,合适的距离取决于分析的目的。如果分析目的主要用于 考虑破坏,然而,模型边界可以靠近一些。如果关注的是位移(变形),则距离边界的 距离需要增加。 借助于经验估计边界对模型的影响是重要的。从一个粗糙的模型和矩形边界开始, 分别研究固定边界和自由边界条件研究改变边界的距离的影响。考察感兴趣的模型区域 内的应力或位移,随边界距离变化所发生的变化,来评价边界对结果的影响。参见第 3.4 节所给出的边界影响研究实例。 一旦完成块体切割和确定模型的边界位置,下一步就是考虑应采用块体单元的大小 与网格密度。较密的网格单元应当处在高应力区或高梯度变形区(即在开挖区附近)。 为了高精度,单元形状尺寸之比(即三角边与高之比)也应尽可能接近于 1。对于 5:1 的情况可能是不精确的。同时建议相邻块体单元的大小不应有较大的突变。合理的精度 是相邻两单元面积之比不应当超过 4:1。

3.2.2 产生节理
UDEC 提供了两个节理生成器:一个统计节理生成器,由传统的岩石力学参数所定 义的参数产生节理;另一个是 Voronoi 分块式节理生成器,用于产生随机尺寸的多边形 块 体 。 统 计节 理 生 成 器是 采 用 JSET 命 令 和 岩体 不 连 续 面的 几 何 参 数统 计 特 征 值 。 Voronoi 生成器是采用 VORONOI 命令和划分块体为随机大小和形状的子块。该生成器 将在下节加以描述。 3.2.2.1 统计节理组生成器 JSET 节理生成器是根据所选定的统计参数生成节理模式。基于这样的节理模式,

-33-

特殊的几何参数对力学特性产生的影响可以进行定量描述。同时,在获得现场观测的节 理模式情况下,采用人工模式与观测模式相匹配方法决定生成的参数。 一个节理组可以通过 8 个生成参数表述:4 个几何参数的均值和 4 个随机参数的均 方差。由 JSET 命令给出的参数如下: JSET a m、a d 、t m、t d 、g m、g d 、s m、s d < x0、y 0 >< ad 0 > 在此, a -节理与 x 轴的夹角;

t -节理段迹线长度;

g -两节理段间的长度(即岩桥长度);
s -垂直于节理迹线的间距。
对于上述的每一对值,前一个带有下标 m 的是均值;第二个带有下标 d 的是均方 差。图 3.2 给出了参数的说明。

图 3.2 节理组参数 最后的三个参数是选择参数:

x0、y 0 - 节理起始坐标(整体坐标轴)。节理将从点( x0、y 0 )开始产生节
理,将充满由选择参数 range 所定义的整个范围。

ad 0 - 所有节理与所给定的节理迹线方向的偏差。
通过定义的一个限制区域(range),生成的节理能够限制在所选择的模型区域。见

-34-

理论与背景第 1.1.3 节中的 range 关键词的描述。在大部分情况下,通过 jregion n 关键 词定义各种区域,在此 n 是指定 JREGION 命令的参考序号 id。 JREGION 定义了一个凸多边形区域来限制节理组的生成范围。该命令参数如下: JREGION id n x1 y1 x2 y2 x3 y3 x4 y4 <delete> 每一个节理区域是通过 id 序号识别。区域的坐标按顺时针方向定义了节理产生的边 界。如果给出选择的关键词 delete ,在此之前由 JSET、VORONOI 或 CRACK 命令生成 的节理将全部被删除。这可避免当指定多个节理区域情况下,被相邻区域产生的节理所 切割。 其他相对于 JSET 命令(即 mat n)的关键词 range 也能用来限制块体被切割。如果 对 JSET 命令没有指定范围,将产生遍及整个区域的节理。 JSET 生成器也可能通过设定很大的 t m 和 s m 值用来产生贯穿UDEC块体的一条节理。 请记住,JSET能够产生不连续节理。不连续节理完全处在块体里边不能看见。当一个变 形块体开始生成单元或当一个刚体模型开始运行时,所有的节理段将被删除。 3.2.2.2 VORONOI 多边形生成器 VORONOI 生成器产生随机大小的多边形块体。UDEC 模型中的一个或多个块体能 够在细分成任意大小的 Voronoi 子块。该节理生成器对模拟裂缝扩展是有用的。 当 Voronoi 块体间的节理强度被超过时将发生断裂。 VORONOI 命令具有下列形式: VORONOI edge l < iterations n> <round v> < range … > 对于 Voronoi 多边形指定平均棱长。多边形具有随机尺寸,但具有平均棱长 l 。 Voronoi 块体的大小尺寸能够通过增加迭代次数变得均匀缺省值 n=5。也能指定圆角长 度。圆角长度,v 必须至少小于块体棱长 l 的 20 倍。 当应用 VORONOI 命令时,多边形区域应当略大于被细分的块体区域。这将抑制边 界影响。Range 关键词在产生 Voronoi 块体区域的用法,与 JSET 命令具有相同的方式。 Voronoi 算法根据随机分布点从多边形区域开始。然后允许内部点移动。迭代过程 运动到这些点。迭代步越高,点间距越均匀。接下来,所有点产生三角形。最后,通过 做具有公共边所有三角形的垂直平分线,生成 Voronoi 多边形。多边形在所指定区域被 其边界所截取。 3.2.2.3 例子 为了说明节理的生成过程,下面给出了产生节理的几个例子。 例 3.3 四组规则节理组 new ro 0.01 blo 0 0 0 20 20 20 20 0

-35-

jset 0 0 50 0 0 0 3 0 jset 90 0 50 0 0 0 3.5 0 jset 30 0 50 0 0 0 4 0 jset 50 0 50 0 0 0 6 0 ret

3.2.3 产生内部边界形状
当完成一个计算区域的 UDEC 建模后,还必须定义与实际物理模型的边界条件相一 致的块体边界条件。这些可能是代表地下开挖,或孔洞或诸如土坝的人工结构或天然边 坡表面特征的外部边界。如果实际物理问题具有复杂的边界,重要的问题是评估简化模 型对对计算结果的影响程度,即简化后的简单模型是否能够重新给出重要的模型机理。 在模拟所有的物理边界(包括后来的模型中增加或开挖的区域)在求解之前必须加 以定义。在以后的分析中增加结构形状必须在事先定义,并通过 CHANGE cons=0 或 ZONE model null 移动,只有当在出现的阶段(如分析筑坝阶段)再出现。 用下列命令实施边界形状的产生: CRACK JSET TUNNEL ARC 每 一 个 命 令 切 割 块 体 成 一 条 或 多 条 裂 缝 , 集 合 组 成 所 要 求 的 形 状 。CRACK 和 JSET 命令产生节理的直线段。TUNNEL 和 ARC 命令分别形成圆形线段和弧段。 在大部分情况下,在人造特征之前应当产生模型的自然特征(如节理、层面、地表 形状、坡面等)。例如,在例 3.9 中首先产生节理结构,然后形成圆形开挖体。 这个问题的模型显示在图 3.11 中。有一组倾角为 110o 和间距 0.5 的节理组和三组三 条独立节理。每一条的倾角为 5o 。圆形隧道的圆心为(0,0)、半径为 2 和边界由 32 段组成。由 JSET 产生的节理贯穿隧道周长。 注意,在UDEC中,正倾角是从x轴的正轴方向顺时针到节理进行测量的角度。 因为隧道是圆形的,所以隧道的开挖可以直接用 DELETE annulus 命令删除。 Delete range annulus 0,0 0,2 另外,如果 TUNNEL 命令在 JSET 命令前给出,则节理并不贯穿隧道周边(如图 3.12)。因此给出坐标范围也可以删除隧道,即 delete range -0.1,-0.1 -0.1,-0.1 该命令仅适用于 TUNNEL 命令。 如果 TUNNEL 命令被首先给出,应当检查模型以确保隧道内部无节理化区域的存 在并不影响初始应力状态。比较在模型中,无隧道和有隧道边界两种情况的应力分布,

-36-

确信应力状态不受边界的影响。

图 3.11 在生成节理组后产生隧道

图 3.12 在节理组生成前产生隧道 不考虑天然节理与内部边界节理(虚拟的)的产生顺序。但以后必须考虑沿隧道周 边的虚拟节理并不影响模型响应。周围块体间的虚拟节理必须被粘合在一起,以阻止滑

-37-

动或张开。这可通过设定较高的粘聚力和抗拉强度值来实现。剪切刚度和法向刚度也必 须赋予足够高的值以阻止两块体产生不同的节理位移。 正是给节理赋予很高的刚度以阻止虚拟节理发生运动。然而,由于UDEC的计算步 取决于刚度,因此,当给予很好的刚度值,系统的响应和收敛速度将十分缓慢。所以, 建议最低的刚度应与小裂隙变形分析所采用的刚度相一致。一个好的经验是节理刚度, k n 和 k s ,取 10 倍于相邻单元的最高刚度等效值。单元在法线方向的等效刚度(产生单 位位移所需要的应力)是:

? (K + 4 G ? 3 max ? ? Δz min ? ?
在式中, K 、 G - 分别是体积模量和剪切模量;

(3.1)

max [ ] 符号表示对于与节理相邻的所有单元的最大值(即相邻节理可能存在几种材 料)。

Δz min - 相邻单元在垂直方向的最小宽度(图 3.13)。

图 3.13 单元尺寸用于刚度计算 通过采用 CHANGE 命令能够改变虚拟节理的材料性质(见 2.6.2 节)。例如,在 上述的例子中,隧道节理性质由下列命令进行改变: change jmat=2 range ann 0,0 1.99, 2.01

在圆心为(0,0)、半径 r1 =1.99 和 r2 =2.01 之间的环状节理将具有节理性质号为 2。 该变化将由以下命令进行查看: plot mat joint

-38-

虚拟节理的性质参数通过 PROPERTY jmat=2 命令予以赋值。 现论述建立模型的最后一个问题。对于连续介质分析时,通常利用开挖体的对称条 件,以利于减小模型尺寸。但对称条件并不易应用于不连续介质,因为,不连续特征的 存在,限制模型的对称性,除非在特殊情况。

3.3 变形块体和刚体的选择
不连续介质分析的一个重要问题,是采用刚性块体还是采用表征完整岩石特征的变 形块体的决策。本节将探讨选择刚体或变形体的问题。如果需要进行变形分析,由几个 不同的模型可以用于模拟块体的变形特征。这些在 3.7 节讨论。 如在理论与背景的第 1 节所述,早期的离散元程序假设块体是刚性的。然而,对于 考虑块体的变形特性的重要性已经得到共识,尤其研究地下工程的稳定性分析和深埋结 构的动力响应更为重要。在离散元分析中考虑块体的变形特征的一个最明显的理由是表 征受侧向约束的“泊松比效应”。 岩石力学问题通常对岩体泊松比的选择是十分敏感的。这是因为节理和完整岩石对 压力是敏感的:其破坏准则是侧向压力的函数(即摩尔-库仑准则)。获得节理岩体的 真实泊松特性,对于具有实际意义的数值模拟是十分关键的。 岩体有效泊松比是由两个部分组成:(1)由节理位移产生;(2)由完整岩石的弹 性性质引起。除非在浅埋条件或低侧向应力水平,完整岩石的压缩性对整个岩体的压缩 性起到重要作用。因此,完整岩石泊松比对节理岩体泊松比产生重要影响。 严格地讲,单一的泊松比 ν ,仅仅定义弹性材料。最多仅适用于极少节理产生各向 弹性特性。所以,为讨论各向异性的岩体材料,定义“泊松效应”是十分必要的。 泊松效应将定义为:当在垂直方向施加荷载的条件,不允许在水平方向产生应变 (位移)时水平应力与垂直应力之比。 假设平面应变条件。作为一个例子,对于各向同性弹性材料的泊松效应为:

σ xx ν = σ yy 1 ? ν

(3.2)

推导垂直节理模型的泊松效应如图 3.14 所示。如果该节理用刚性块体模拟,施加的 垂直应力将根本不产生水平应力。这显然不实际的,这是因为忽略了由完整岩石的泊松 比产生的水平应力。 节理和岩块以串联方式相互作用。换句话说,作用在节理上和岩块上的应力是相等 的。节理岩体的总应变是节理应变和岩块应变之和。岩体弹性性质可以通过叠加节理和 岩块的变形得到:

?ε xx ? rock + C jo int ing ?ε ? = C ? yy ?

(

)?σ ?
?

σ xx ?
yy

? ?

(3.3)

-39-

图 3.14 岩体内具有水平和垂直节理的泊松效应模型 如果以各向同性弹性材料模拟完整岩块,其变形矩阵为:

C rock =

1 +ν E

?1 ? ν ?? ν ?

-ν ? 1 -ν ? ?

(3.4)

节理变形矩阵为:

C jo int ing

? 1 ? sk n =? ? ?0 ?

? 0? ? 1 ? sk n ? ?

(3.5)

式中, s - 节理间距, k n - 节理的法线刚度。 如果方程式(3.3)中的 ε xx =0,则
total σ xx C12 = total σ yy C11

(3.6)

式中, C

total

= C rock + C jo int ing 。

因此,节理岩体的总泊松效应为:

σ xx ν (1 + ν ) = σ yy E /( sk n ) + (1 + ν )(1 ? ν )

(3.7)

方程(3.7)给出了以 E /( sk n ) 作为变量的函数(如图 3.15 所示)。而且此图是为了检 验此公式所进行 UDEC 模拟的结果。其比值 E /( sk n ) 是与节理刚度相关的完整岩块刚度

-40-

的度量。对于 E /( sk n ) 的低值,岩体泊松效应主要受控于完整岩块的弹性性质。对于

E /( sk n ) 的高值,泊松效应受控于节理。
现在考虑变化的不同倾角岩体的泊松效应。泊松效应是节理产状和弹性参数的函 数。考虑如图 3.16 所示的特例。岩体包含两组等间距,倾角为与 x 轴的夹角 ? 的节理 组。节理弹性性质是由法线刚度和切向刚度组成。完整岩块假设为完全刚性的。

图 3.15 垂直节理岩体的泊松效应( 完整岩石的泊松比 ν = 0.3)

图 3.16 与 x 方向夹角为 ? 、间距为 s 的节理岩体泊松效应模型 对于这种节理模型的泊松效应为

-41-

σ xx cos 2 ? [(k n / k s ) ? 1] = σ yy sin 2 ? + cos 2 ? (k n / k s )

(3.8)

上述公式中几个不同的 ? 值的图形如图 3.17 所示。在此图也给出了 UDEC 数值模 拟结果。UDEC 模拟的结果与方程(3.8)十分接近。

图 3.17 具有不同节理角的节理岩体(块体是刚体)的泊松效应 方程(3.8)表明了在数值分析中选用合理的节理剪切刚度是重要的。剪切刚度与法 线刚度之比值戏剧性的影响岩体泊松效应。如果剪切刚度等于法线刚度,则泊松效应为 零。对于从 2.0 到 10.0 比较合理的 k n / k s 比值,泊松效应是相当的高,达到 0.9。 其次,完整岩块的弹性性质对其影响将在 ? =45 o 的情况下加以研究。在下面的分 析中,节理为垂直分布模式,处理完整岩石为各向同性弹性材料。岩体总弹性性质将由 叠加节理和岩石的变形获得: 具有倾角为 45o ,等间距的两组节理的为:

C jo int ing =

1 2 sk n k s

?k s + k n ?k ? k n ? s

ks ? kn ? ks + kn ? ?

因此,岩体总泊松效应为:

σ xx [ν (1 + ν )] / E + (k n ? k s ) /(2 sk n k s ) = σ yy [ν (1 + ν )(1 ? ν )] / E + (k n + k s ) /(2sk n k s )

(3.9)

-42-

方程(3.9)中对于泊松比 ν =0.2 和几个不同的 E /( sk n ) 比值的情况下,给出的图形 如图 3.18 所示。对于较低的 E /( sk n ) 比值,岩体的泊松效应受控于完整岩石的弹性性 质。对于较高的 E /( sk n ) 比值,泊松效应受控于节理。

图 3.18 倾角 ?=45 、等间距两组节理岩体的泊松效应(泊松比 ν =0.2)
o

3.4 边界条件
在数值模型中,边界条件是由表征模型边界的场变量(即应力或位移)组成。边界 分 为 两 类 :实 际 的 和 人工 的 。 实 际边 界 是 模 型中 存 在 的 物理 对 象 ( 即隧 道 界 面 或地 面)。人工界面并不真正存在,但为了封闭所选择的单元数(即块体)必须引入。能够 施加的边界类型是类似的。这些条件将首先讨论。然后(在第 3.4.4 节),对涉及到边 界的位置和人工边界的选择以及对计算结果的影响提出某些建议。 力学边界具有两类:位移和应力。自由边界是应力边界的特例。两类力学边界在 3.4.1 和 3.4.2 中讨论。第三种类型,“边界元边界在 3.4.4 中论述。用于动力分析的粘滞 边界和自由场边界在理论和背景中的第 4 节中讨论。

3.4.1 应力边界
UDEC 模型的缺省边界是无约束的自由边界。力或应力可以通过 BOUNDARY 命 令 施 加 到 任 意 整 个 边 界 或 部 分 边 界 上 。 用 stress 关 键 词 可 以 指 定 平 面 应 力 张 量 ( σ xx、σ xy、σ yy )的每一个单独应力分量。例如,命令 Boundary stress 0,-1e6,-2e6 range 0,10 -1,1 将施加 σ xx = 0、σ xy = ?10 、和σ yy = ?2 × 10 到位于坐标窗口 0< x <10,-1<y<1 范
6 6

-43-

围内模型边界上。用户应检查窗体周围的所有边界角点所指定的边界条件,这可用下面 的命令实现: Print boundary 每一外边界的角点将以表的形式列出指定的边界值。边界在模型计算过程中可能运 动,所以,用户必须检查坐标窗口是否足够大,使窗口包含所有的边界角点。 另外,边界条件也可在两边界的角点间施加。作用的边界从第一个边界角点地址顺 时针到第二个边界角点地址施加。角点的地址可以用下面的命令进行查看: Print boundary state 角点地址以标题 COR 和列的方式显示。最好给出坐标范围而不是角点地址,因 为,应用的地址与模型相关,即如果改变命令排序,地址可能发生改变。 在 UDEC 中按照一般符合惯例,压应力为负号。而且,UDEC 实际上施加应力分 量作为力或者是产生作用到给定的边界平面上的应力张量的摩擦力。摩擦力被分成两个 分量:永久的和瞬时的。动力分析采用随时间变化的永久摩擦力是常荷载和瞬时荷载。 采用 xload 和 yload 关键词也可将单个 x 方向和 y 方向的力的分量施加边界上。用 LOAD 命令也可指定荷载于刚性块体上。刚性块体是施加到块体的形心上。 3.4.1.1 施加应力梯度 BOUNDARY 命令可以增加关键词 xgrad 和 ygrad,此关键词允许应力或力在指定 的边界上按线性变化。在上面的每一个关键词后面的关键词,给出了在 x 或 y 方向变化 的应力分量: Xgrad sxxx axyx syyx Ygrad sxxy sxyy syyy 应力从坐标原点(x=0,y=0)沿边界按距离线性变化:
0 σ xx = σ xx + ( sxxx) x + ( sxxy) y
0 σ xy = σ xy + ( sxyx) x + ( sxyy ) y 0 σ yy = σ yy + ( syyx) x + ( syyy ) y

式中, σ xx , σ xy , σ yy 是原点的应力分量。
0 0 0

上式运算过程由下面的例子加以解释: boundary stress 0,0,-10e6 ygrad 0,0, 1e5 range -0.1,0.1 -100,0 原点的应力是 σ xx =0, σ xy =0, σ yy =-10×106 。
0 0 0

在 y 方向的应力变化 σ yy 的计算式为

σ yy = -10×106+(105) y

-44-

σ yy 在 y=100 的值为-20×106。在边界上,y 的变量是从原点开始线性变化。
3.4.1.2 改变边界应力 正如上述讨论,瞬时荷载可用 history 关键词,为动力分析施加边界条件。对于静 力分析,在 UDEC 模拟过程中,改变边界的应力值也是必要的。例如,作用到基础的荷 载可能发生改变。为了考虑应力或荷载突然变化的影响,需要给出新的 BOUNDARY 命 令,对原来边界角点施加改变后的应力或荷载。 在这种情况下,新的边界应力或荷载值将被叠加到原有的值上。如果应力被移动, 当前的值应当以相反的符号给出。如果一瞬时荷载被改变(即用history 关键词指定荷 载),用相同的history 所给出新的荷载被叠加到存在的荷载上。然而,用不同的history 所给出的一个新的瞬时荷载,将取代老的荷载。 3.4.1.3 打印和绘图 边界应力和荷载可以通过命令 PRINT bound 和 PLOT bound xcond 或 PLOT bound ycond 进行验证。PRINT bound 命令列出每一角点所指定的边界角点地址和对应的当前 的值。一旦 BOUNDARY 命令被发出,整个模型的外边界生成边界角点列表。选择的关 键词可以用 PRINT bound 命令来检查边界的不同条件。例如 print bound force 列出永久力(fx,fy)和在当前荷载步增加的荷载增量(fxi,fyi)的表列。如果施加 瞬时荷载(用 BOUND … hist 命令),总荷载为永久荷载与当前循环步的瞬时荷载之和。 需要施加边界条件的边界角点地址能够通过以下命令输出显示。 Print bound state PLOT bound xcond 或 ycond 命令通过符号表示施加在 x-或 y-方向边界条件的类型。 3.4.1.4 提示和建议 本节将探讨施加应力边界条件所存在的一些琢磨不定的困难。应用 UDEC 程序,可 将应力施加于一个没有位移约束物体的边界上(不同于很多的有限元程序,在此需要某 些约束)。物体将完全按照实际物体所作用的方式产生反作用,即如果边界应力没有平 衡,则整个物体将产生运动。例 3.10 说明这个效应。 图 3.19 给出了所生成的图形。施加的 σ 11 引起水平应力作用到物体上。由于该物体 是倾斜的,所以力导致运动引起物体的转动。相似但影响有所不同的是从一个物体中开 挖应力边界附近的材料。物体最初在重力作用下平衡,但移去材料减小重量,整个物体 开始向上运动(例 3.11)。 在运行中这个数据文件(例 3.11)所遇到的困难能够通过固定底部边界予以消除。 第 3.4.4 节包含涉及这样的人工边界的位置的信息。 最后,应力边界影响所有的自由度。所以速度边界条件必须在应力边界条件之后给 出。如果应力边界在速度边界之后施加,给出速度的响应将会消失。例 3.12 说明此问题。

-45-

图 3.19 力作用于倾斜物体产生的旋转位移

3.4.2 位移边界
UDEC模型不能直接控制位移。事实上,他们对计算过程并不起作用。为了施加已 知的位移到边界上,有必要固定边界或对于给定的迭代步,给出边界的位移速率(采用 COMMARY命令)。如果要求的位移是D,对应于时间为增量T(计算步),施加的速 度为V(即D=VT),在此, T = ΔtN , Δt 是时间步,N是迭代步数(循环数)。实际 上,为减小模型系统的波动,V应当赋予小值而N应大值。 BOUNDARY 命令用来沿着边界(而不是与 x 或 y 坐标轴一致)固定在 x 或 y 方向 (BOUND xvel 或 yvel)或法线或切线方向(BOUND nvel 或 svel)的变形块体的结点 速度。刚性块体的速度可以用 FIX 命令。如果 FIX 命令放在 INITIAL 命令之前,速度 能够用用户所选定值加以固定。也可用 FISH 函数加以改变。 随时间变化的速度历史能够通过 BOUND… hist 命令对刚体或变形体予以施加。这 个 history 关键词必须出现在指定速度历史边界 BOUND xvel 或 BOUND yvel 的同一 行。历史也可通过 FISH 函数施加。正如在 3.4.1.4 节所讨论的,速度边界应当在应力边 界之后给出。变形块体的固定速度边界条件能够通过 BOUND xfree 或 BOUND yfree 命 令得以释放。对应刚体,则用 FREE 命令。 通过在指定的速度关键词之后增加关键词 gvel,速度也允许在指定的边界范围按线 性变化。Gvel 关键词后有六个参数,描述了速度分量在 x 或 y 方向的变化: gvel vx0 vy0 vxx vxy vyx vyy 从坐标原点(0,0),速度随距离线性变化:

-46-

vx=vx0+(vxx)x+(vxy)y vy=vy0+(vyx)x+(vyy)y 在式中,vx0 和 vy0 是在原点的速度分量。 (3.11)

3.4.3

真实边界-选择合理类型

准确地识别计算模型中某些特殊表面边界条件的类型有时是很困难的。例如,在一 个三轴试验模型中,通过台板施加的荷载是应力边界?或台板作为刚体而视为位移边 界?当然,包括台板的整个试验机能够一并模拟,但可能是非常耗时的。记住如果存在 很大差异的刚度,UDEC 花很长时间才能收敛。一般情况,如果施加荷载的物体与样本 相比非常“刚”(比如说,是弱者的 20 倍)。则可以作为应力控制边界加以模拟。很清 楚,作用在物体表面上的流体压力属于后一类。位于节理岩体上的基础,常常可移去常 值的速度而作为刚性边界加以处理,用于研究岩体的破坏荷载。这个研究还有另一优点 -它非常容易地控制试验和获得较理想的荷载或位移图形。人所共知的刚性试验机远比 一般试验机稳定。

3.4.4 人工边界
人工边界分为两类:对称轴和截取边界线。 3.4.4.1 对称轴 当几何形状和荷载是关于一轴或多轴对称时,就可能考虑其特点进行模型的简化。 例如,如果一结构的对称轴是垂直的,则对称轴上的水平位移均为零。所以,我们能够 取其垂线作为边界,边界上的所有结点,用 BOUND xvel=0 设置 x 方向的位移为零。如 果在对称轴上的速度为零,也可采用此命令进行设置为零。在此情况下,y 方向的位移 并不受到影响。类似的情况是水平轴为结构的对称轴。命令 BOUND nvel=0 能够用来设 置与坐标轴成一定角度的对称轴。 正如在 3.2.3 节所讨论,不连续面的存在使得对称性的利用变得困难。当采用对称 性,应当十分小心地考虑节理产状的影响。 3.4.4.2 截取边界 当模拟无限大物体(即地下隧道)或非常大物体时,由于内存或计算时间的限制, 其计算模型不可能涉及整个物体。人工边界应置于所关心区域的足够远,以致不影响模 型效应。知道这些人工边界应置多远以及对模型的计算应力和位移产生多大的误差是有 用的。基于包含两个隧道的弹性介质模型进行一系列数值试验。该模型仅包含用于产生 两个隧道的虚拟节理。图 3.20 和图 3.21 分别显示最小和最小的模型。 在所有的计算中,块体中的单元尺寸是相同的,以致消除单元离散的影响。考虑两 个物理量:在较大跨度隧道顶板中心点的垂直位移, u y ,和两个隧道间柱中点的垂直应 力, σ yy 。原岩应力比是 2:1(垂直与水平之比)。对应每一矩形边界的几何体,进行两

-47-

中边界条件的计算:长应力边界和零位移边界。另外,还进行了 UDEC 边界元边界的计 算。例 3.13 包含上述三种边界的数据文件。

图 3.20 包含两个隧道的小模型

图 3.21 包含两个隧道的大模型 图 3.22 用无量纲的形式显示了计算结果。观测的位移和应力相对于他们的无影响值

-48-

进行正则化处理。“边界尺寸”的值是平均宽度和高度,“物体尺寸”是贯穿两隧道的平均 距离。从图 3.22 中给予我们几点提示:

图 3.22 包含两隧道的截断边界对模型影响的 UDEC 数值试验结果 1、固定边界低估了应力和位移(即计算的应力和位移小于实际工程值),而应力 边界则相反。 2、两种类型的边界条件“包含”了真实解,因此,有可能采用小的边界进行两种试 验,通过取两种结果的均值,可以给出真实解的合理估计。 3、作为初步指导,对于边界距离与跨度比为 5 的情况,对于固定和应力边界的应

-49-

力和位移误差约为 6%。 4、对于边界元边界的误差约为 0.5%。 从上述结果表明,模拟无限介质,边界单元是提供人工边界的最好的方式。然而, 存在一个困难。边界元边界求解需要计算和储存刚度矩阵。对于较大模型,这种储存能 利用额外内存。甚至存在这种限制,边界单元对于大部分问题仍是最有效的方式。在从 事一系列计算前,应当进行不同边界类型的比较测试试验。 上述报道的数值试验是弹性体。这可能表征的最差情况。因为当存在塑性特性,位 移和应力变化受到较大限制。人工边界可能被置于没有严重误差的边界稍远些。然而, 任何人工边界不必足够地近。 3.4.4.3 边界元边界

边界元边界条件是一种模拟各向同性、线弹性材料的无限边界效应(半无限边界) 的人工边界。其公式是一种在外边界与离散单元法耦合的直接边界元求解算法。应用的 耦合公式的算法,是由 Lorig 和 Brady(1983)提出。离散元区域是置于弹性区域。该 区域首先承受原岩应力。如果在离散元区域发生改变(如开挖),在离散元和边界元耦 合的边界,应满足位移和应力的连续。 在目前的公式中,在离散元区域的周围,边界元结点与块体角点(网格结点)必须 一致。在每一边界单元结点的开挖位移,通过乘以结点位移矢量刚度矩阵直接确定结点 力:

f = Ku
式中,K - 是边界单元刚度矩阵;

(3.12)

f - 次生结点力矢量
u - 结点位移矢量。
次生结点力与总结点力的关系为:

ft = f0+ f
在此, f 0 -代表初始结点力矢量。 因此,方程 3.13 可被改写为:

(3.13)

f t = Ku + f 0

(3.14)

在迭代计算中,用这种方式确定的结点力,能被施加位于离散元区域的周边块体角点。 因此,由于继续这种释放过程,在接下来的计算迭代中,结点位移和力被修改和应用。 两个区域的耦合,必须确保所有结点的位移和力在其界面是连续的。这些条件需要 非线性材料特性被限制在离散单元区域内。节理滑动面或分离也不能在界面内嵌。 边界单元公式采用直线单元,离散单元由直线定义。所以,能够确保物理的相容 性。离散元区域的每一次循环导致块体位移和转动。通过识别块体角点位移和边界元界 面的结点位移,检查在界面是否满足动态连续性。因为结点是两个块体的接触面,所

-50-

以,结点位移取两个相邻角点位移的平均值。在界面上两个块体的接触必须是弹性的, 因此,位移差是较小和平均算法是近似合理的。因为,位移变化假设为线性的,所以, 能够满足在边界元区域与块体的棱相接触界面的位移连续。 每一次迭代产生界面位移,采用边界元刚度矩阵和该位移计算界面力。如果在离散 元区域内考虑重量作用,为模拟开挖离散元区域的单元运动,将导致作用在界面力的不 平衡。这种不平衡对模型产生影响,以致导致刚性物体向上运动。迭代过程的收敛要求 固定边界元区域的某些结点界面位移。这需要修改刚度矩阵,对于每一组界面位移,使 得所选择的内部点,P 的位移为零。 两个位置的点需要固定:一个是 x 方向,另一是 y 方向的运动。固定点通常应位于 离散元区域的外边界 1~3 倍 UDEC 模型之间 用户应当知道,对应边界元边界施加的应力单位必须是 MPa。这样为避免在 PC 机 上计算刚度矩阵的截断误差。参考表 2.6 的统一单位。同时,边界单元也仅用在平面应 变条件;同时,边界元结点的最大数限制在 300 个以内。 对于边界元边界的应用例子见计算实例的第 7 节和在验证例子中的第 5 和 7 节。

3.5 初始条件
在土木和采矿工程中,任何开挖或建造之前,地下存在一原始应力状态。在 UDEC 模型中需要通过设定初始条件,模拟原岩应力状态。因为,此状态对后来模型特性产生 影响。理想地,初始应力状态的信息应来自于现场实际测量结果,但当这些资料无法获 取时,模型应当在一个可能条件的范围内运行。尽管该范围是很大的,但有很多约束因 素加以控制(即系统必须是平衡的,选择的屈服和滑移准则在任何位置都不能违背)。 对 于 具 有 自由 表 面 的 均匀 分 布 的 岩土 层 , 垂 直应 力 通 常 等于 岩 土 体 自重 应力,

gρ z ,在此, g 是重力加速度, ρ 是材料的体积密度和 z 是距地表的距离。然而,现
场的水平应力很难估计。存在共同的但又是错误的是,确信存在由 ν /(1 ? ν ) 所给定的水 平与垂直应力的“天然”比值。这个公式是基于重力突然施加到侧向位移被阻止的弹性 材料的假设推导的。该条件很难在实际工程中应用,其原因是重复的构造运动,材料破 坏、地层运动和由于断裂和局部的应力封闭(见 3.11.3 节)。当然,如果具有足够特定 区域的历史知识,我们可以模拟整个过程,为工程研究规划获得初始条件。这种研究通 常是行不通的。典型地,我们做出这样的让步:在模型中采用一系列应力,然后运行 UDEC,直到获得平衡状态。对于给定的系统,认识到存在很多平衡状态的是重要的。 在下一节,我们逐步考察较为复杂的情况以及可能模拟初始应力的方式。用户要鼓 足勇气去试验所提供的各种数据文件。

3.5.1 在均匀介质中的均匀应力:无重力
对于深部开挖,由于与作用到模型上的岩石条件相比,其应力较小,可以忽略重力

-51-

应力的作用。当 SET gravity 命令被省去,而引起重力加速度的缺省值为零。初始应力 可以用 INSITU 命令施加,即 insitu stress -5e6 0.0 -1e7 szz = -5e6

应力分量 σ 11 (或 σ xx ), σ 22 (或 σ yy )和 σ 33 (或 σ zz )被赋给整个模型的压应 力分别为 5×106 、107 和 5×106 。由于所有的塑性本构模型都涉及到 σ 33 应力,所以,当 进行UDEC模型计算时,记住给 σ 33 赋值是重要的。如果忽略,其缺省值是零,这可能 引起在垂直于平面方向上的破坏。INSITU命令设置所有的应力为给定的值,但并不保 证应力处于平衡状态。至少存在两个方面的问题:首先,变形块体的应力可能违背所给 定的非线性本构模型的屈服准则。在这种情况下,在给出STEP命令后,块体单元立刻 出现塑性流动以及应力重新调整。这可通过试算一次进行检验和考察其响应(即PLOT plas)。第二,在网格边界指定的应力不可能等于给定的初始应力;在此情况下,只要 给出STEP命令,边界结点就发生运动;同时,对出现这种可能性应当通过输出,进行 检查(即PLOT vel)。在例 3.14 中的命令产生初始应力,并在给定的边界条件下处于 平衡状态的单一块体。

3.5.2 无节理介质中具有梯度变化的应力:均匀材料
在接近地表,不能忽略应力随深度的变化的影响,SET grav 命令被用来为模型形成 重力加速度运算。重要的是要理解 SET grav :他不能直接引起模型的应力,仅产生体 积力作用在变形块体的所有结点(或刚性块体的形心)。这体积力对应着围绕每一结点 材料的重量。如果存在初始应力,这重力将引起材料在重力方向运动(在迭代中),直 到由单元应力产生等于或相反的力。事实上,当给出合适的边界条件(即固定底部), 模型将产生与作用的重力相容的重力场。然而,这个计算过程的效率不高,需要运行几 百次才达到平衡。最好初始化内部应力,以致满足平衡和梯度变化条件。INSITU 命令 必须包括 xgrad,ygrad 和 zgrad 参数以致应力梯度匹配重力梯度 gρ 。内部应力在应力 边界也必须与边界条件相匹配。 例如,考虑一个 20m×20m 无节理均质材料的块体模型位移地下 200m,固定底部边 界,其他三边为应力边界。例 3.15 对这个条件产生平衡系统。 在这个例子中,水平应力和梯度等于一半的垂直应力和梯度。但可设置任意不违背 屈服准则(在此为摩尔-库仑)的值。在准备一系列数据文件后(如上述例子),进行 一个循环步计算并监测不平衡力。当不平衡力大约达到施加荷载的量级时,将显示内部 应力和边界应力相匹配的任何破坏。注意,在上述例子中,如果忽略 σ zz ,则在垂直于 平面方向材料发生破坏。

3.5.3 无节理介质中具有梯度变化的应力:非均匀材料
当存在不同密度的材料时,给出合理的初始应力是困难的。考虑具有自由面的层状

-52-

系统,具有滚轴侧向边界和固定底部边界的闭合块体。假设材料具有如下密度分布: 从 0 到 10m 深度:1600kg/m3 从 10 到 15m 深度:2000kg/m3 从 15 到 25m 深度:2200kg/m3 通过在例 3.16 中的数据文件,产生平衡状态。 对应每一种材料密度产生一个独立的块体。虚拟节理分割每一块体。内部应力可由 每一块体已知埋深计算。注意,这个例子被简化。在实际情况下,弹性模量是变化的, 并存在水平应力。如果在一层内存在水平应力,这也可用 INSITU 命令产生。 此例在进行一次迭代并不平衡,大约需要 500 次迭代。当初始应力与边界应力匹配 时,虚拟节理的存在阻止了模型的自动平衡。这是圆形角点的逻辑问题,即沿着块体棱 计算的接触力与角点力略有不同。圆角长度越大,则初始不平衡力也越大。对于实际和 虚拟节理都存在这种情况。例如,如果在上述数据文件中的棱长从 0.1 改变到 0.0001, 则模型在一次计算后将处于平衡状态。不平衡力接近于零。然而,事实上,圆角长度不 应当设置到这样极其小的值,因为当沿着不连续面发生运动时,该长度也影响新的接触 产生。当剪切位移大于 2 倍的圆角长度时,就产生新的接触。因此,如果这个长度极 小,就可能产生很多新的接触。这可能大大降低计算速度。对应小的 ROUND 值也可能 发生超过界面而产生接触嵌入。考虑实际运行时间,圆角长度应设置均值棱长的 1%。 无论真实节理还是虚拟节理,甚至当内部应力与边界应力相匹配的情况,节理模型仍需 要进行一定的迭代计算方可获得平衡状态。

3.5.4 具有非均匀单元的密实模型
含有非均匀单元模型在重力作用下达到平衡状态时,可能观察到令人费解的结果。由 不同尺寸的变形体组成的模型,通常将产生非均匀单元。当设定块体为摩尔-考虑准则 或其他非线性本构模型时,即使边界为平面和水平自由面,最终的应力状态和位移模式 并非均匀。在例 3.17 说明这个事实。图 3.23 显示了位移矢量和垂直应力等值线。对于 重力加载于弹性系统,采用目标阻尼(DAMP auto)以提供快速平衡求解。 由于在模型的两侧采用了滚轴边界,可能期望材料模型在两边向下平行运动。然 而,临近边界和垂直节理的单元小于块体中间的单元。对于静力分析,UDEC 对所有单 元试图保持同步,以致在块体的边界结点增加惯性质量,从而对小块体单元进行补偿。 这些结点以与比里边块体产生较慢的加速度。这对最终线性材料的状态不产生影响,但 在与应力路径相关材料产生非均匀性。对于无粘结强度的摩尔-库仑材料,这种情况类 似于砂粒从某一高度下落到容器中,期望他最终状态是均匀的。实际上,因为侧向应力 并没有立刻集合,将产生大量的塑性流动。甚至对于均匀单元的模型,这个研究也并非 理想。因为水平应力依赖于动力过程。

-53-

图 3.23 非均匀应力和位移 最合理的求解是采用 INSITU stress 命令设置初始应力以形成要求的 λ 值(水平与 垂直应力之比)。例如,在前面的数据文件中的 STEP 1000 命令可用下一行替代: insitu stress (-1.5e5, 0,-2.0e5) szz (-1.5e5) & ygrad (1.5e4,0,2.0e4) zgrad (0,1.5e4) 获得 λ =0.75 的稳定应力场,仅需要很少的计算步(STEP 100)就能获得平衡,且 应力状态是均匀的。 另外,为最终状态,首先采用弹性模型计算,然后改为非线性模型。用下一行命令 替代 STEP 1000: zone tens=1e10 coh=1e10 step 750 zone tens=0 coh =0 step 250 图 3.24 显示了最后状态的位移矢量和垂直应力等值线图。在压缩过程中阻止材料 发生屈服。当达到平衡时,保存初始性质参数。

3.5.5 随模型变化的初始应力
可能存在这种情况:采用一种材料的变形块体模型计算,使其达到要求的应力分 布;但在另一种模型用来连续模拟。能够对整个块体模型(通过 CHANGE)或在块体内 的单元区域(通过 ZONE 命令)发生改变。在另一情况,如果模型用一非零模型取代, 在影响单元的应力仍保持。如在例 3.18 中所述:

-54-

图 3.24 均匀应力和位移 基于此进行计算,由初始弹性模型所产生的应力仍存在,并作为含有摩尔-库仑新 模型的初始应力。 应当记住两点。首先,如果在模型的任何部分产生一个开挖块体(cons=0),即使 在后来由另一非开挖单元替代,则所有的应力将从零单元被移出;第二,如果一种材料 模型被零所替代,同时,在新模型中实际的应力为零,然而,必须对这个区域用 INSITU 命令设置为零。如果块体矿石被挖出,然后用充填体充填将会出现这种情况。 充填应当以无应力状态开始计算。

3.5.6 节理化介质的应力
在节理和裂隙化介质中,能够模拟初始应力状态空间的不均匀性。在涉及地质不连 续面在不同的施工阶段所出现的滑移、张开和断裂的地质构造和模型变化过程中,必将 产生与路径有关的应力。应力状态空间不均匀性是进行地下开挖设计的重要因素,尤其 如果十分应力集中不利于开挖结构的稳定性情况。 确定节理模型中所形成的应力状态是否代表初始状态是十分困难的。正如在后面的 3.11.2 节所讨论的,统计分析可为研究模型表征的可信度提供了一种有用的工具。Brady 等(1986)报道了采用 UDEC 模型进行研究的一些例子。 当生成一个节理模型使之平衡时,可以借鉴的几个问题:首先,INSITU 命令应在 产生所有的节理后给出。然后,将初始化节理的切向应力和法向应力,分解到每一节理 面所对应的初始应力。 如上所述,一个节理模型开始并非平衡,甚至当内部应力被设置与边界应力相匹配

-55-

的情况。因此,需要某些循环步进行计算,并对不平衡力进行监测。除此之外,应对模 型中的各个位置的位移或速度历史进行记录。用户在进行下一阶段的分析前,应确信模 型的平衡应力状态确实达到应停止计算的平衡状态。 节理刚度也将影响趋于平衡状态计算的迭代步。为节理刚度值的选择提供一些经验 性建议:对于虚拟节理,(方程 3.1)给出可作为节理刚度的极限条件的一般性建议。 如果法向或剪切刚度被限制在大于 10 倍方程 3.1 所定义的值,系统的特性不会产生有意 义的影响。该方面将在 3.6 节进一步讨论。 对于给定的初始应力状态和节理强度性质,也有可能使某些节理在模型平衡状态的 迭代过程中产生滑动、张开。限制模型内的节理滑移是合理的;在节理端部将会导致 “锁闭”应力。然而,用户应避免节理破坏扩展到模型边界的情况。这表明模型条件没 有较好地模拟。可能需要重新评价给出应力状态、节理性质、节理产状与位置。如果这 样的条件仍使节理破坏扩展到边界,则应当考虑一个固定边界条件。这意味着节理在边 界被截断。 采用如下命令: plot block slip lmag open yell 来识别已经发生滑动(通过鲜艳的色彩识别)或张开(通过黄色识别)的节理长度。 例 3.19 中的数据文件给出了一节理倾角为 60o,限制在 20o 倾角的范围内。对于给 定的初始应力,60o 的节理发生滑移,而 20o 的节理没有发生。在此,所有节理的内摩擦 角都是 30o 。 例 3.19 以受约束节理的滑移 图 3.25 显示了节理滑动区域由粗线表示的块体图形(关键词 lmag 产生一粗线)。

σ xy 的等值线也在图中给出,显示了 60o 节理端部附近的锁闭应力区域。该图形由以下
命令给出: Plot pen block slip lmag sxy max=0 zero

3.5.7 绘制应力等值线图
通过在 UDEC 模型区域,覆盖网格等值点便形成应力等值线图。单元应力被移植到 每一单元的结点。然后,对储存在结点的值进行线性插值,计算其等值线的值。缺省的 网格尺寸在 x 方向和 y 方向均为 20 点。 如果网格间隔小于单元尺寸,可能观察到应力等值线图的某种非均匀性。所有位于 一个单元的结点具有相同的应力值。这可能产生等值线图的变形(例如图 3.24)。最好 的效果,网格尺寸应当近似与单元尺寸具有相同的量级。网格尺寸可以用 grid 绘图命令 和调整。例如,为调整图 3.24 所示的等值线网格尺寸,键入如下命令: plot block syy grid 10 10 disp

-56-

由此获得的等值线图较为均匀(图 3.26),但对于渐变单元的情况,可能不允许调整等 值线网格来匹配单元尺寸。如果允许,用 WINDOW 命令调整绘图窗口来定义近似等于 单元尺寸的区域。

图 3.25 限制节理的滑移,图形显示剪切应力等值线

图 3.26 对于 10×10 的等值线网格的应力等值线

-57-

3.6 加载与施工模拟
通过在不同的计算阶段施加不同的荷载条件,可以模拟诸如开挖和建造的施工顺序 等模型荷载的变化。有多种方式进行荷载变化:或者施加新的应力或位移边界,或者通 过改变块体材料模型为零模型或不同材料的模型,或通过改变材料的性质等。 重要的是认识模拟的顺序应对应于工程实际工程的施工阶段。在大部分分析中,每 一施工步骤对应于不同的静态加载的变化,即物理时间并不是一个参数。UDEC 能够进 行节理瞬时流动、热传导和动力分析的模拟。对于所有的情况,必须满足平衡应力状态 的静力解。例如,对于受爆破作用的动力计算或通过节理流动瞬态计算。 在另一方面,不能被直接模拟时间参数。必须采用某些工程判断来估计时间效应。 例如,在预测的位移或应变已经出现的情况后改变模型参数。该位移可估计已经超出给 定的时间周期。 为了考察模型对变化的响应,必须进行荷载变化引所起不平衡力的研究。所以,改 变弹性参数将没有影响;如果变化引起当前应力状态超过破坏极限,改变强度性质将会 产生影响。 下面的例子中研究了模拟顺序。这个问题涉及节理岩体地下开挖的稳定性分析,包 括动力和静力计算阶段。进行分析的阶段如下: (1)原岩应力状态的平衡; (2)隧道开挖; (3)动力荷载的应用。 研究的目的在于探讨原岩应力条件下的开挖结构在受动力波作用的稳定。 位于岩体中的隧道包含具有 ± 15 o 倾角和 1m 的间距的两组节理。其计算模型由例 3.20 所列出的一系列命令: 例 3.20 一个地下开挖的稳定性分析-初始模型 图 3.27 和图 3.28 显示其模拟结果图。隧道是一水平洞室,所以垂直对称轴线(即 模型的左右边界)是位于洞室的中线。顶部边界位于地表面。节理仅在包围隧道足够大 的有限范围内产生。为加速计算,面积小于 0.1 的非常小的块体被删除。 在开挖体附近产生较小的网格。根据动力分析的特性选取单元的棱长。动力模拟的 单元划分要求在理论与背景的第 4.4 节中描述。单元划分命令是: gen edge 3.5 range -10.33 10.33 15 45 gen quad 4.0 为模型的初始平衡给定的边界和初始条件如下: insitu stress -3.5e6 0 -6.8796e6 grav 0, -9.81 bound xvel 0 range -19.4 -17.9 -40.1 300.1 ygrad 11666.7 0 23932

-58-

bound xvel 0 range 17.9 bound yvel 0 range -19.3

19.3 -40.1 300.1 19.3 -40.1 39.9

在这个分析中,假设破坏仅限于节理。块体具有弹性特性。对于块体性质的命令是: prop mat=1 den=2340 bu=8.39e9 sh=6.29e9

图 3.27 隧道区域模型的全图

图 3.28 隧道附近区域的局部图形

-59-

节理具有内摩擦角 20o ,粘聚力为 1MPa 和零抗拉强度。如果节理发生剪切或张拉 破坏,粘聚力便消失(jcons=5)。命令是 change jcons=5 prop jmat=1 jkn=1.0e11 jks=1.0e11 jcoh=1.0e6 jfric=20.0 jtens=0.0 set jcondf=5 jmatdf=1 采用如下命令,将节理性质参数赋给围绕开挖体边界的虚拟节理: change jmat=2 range ang -1 1 change jmat=2 range ang 89 91 prop jmat=2 jkn=1.0e11 jks=1.0e11 jcoh=1e20 jten=1e20 用如下命令进行求解,并监测初始阶段: hist unb ydisp 0,35 damp auto step 3000 注意 DAMP auto 被用来使模型快速趋向于弹性平衡(见在 3.9 中的注意点 9)。 保存这一阶段的模型状态: save stepq.sav 通过绘制不平衡力和垂直位移历史图来证实模型的平衡状态。图 3.29 显示了点 (x=0,y=35)的垂直位移已经收敛趋于一常数。 下一步,通过如下命令立刻开挖隧道:

图 3.29 监测在点(x=0,y=35)的垂直位移历史

-60-

del -2.5 2.5 30,35 del 2.0 2.3 35,35.5 del -2.3 -2 35,35.5 del -1.55 1.55 35,36.2 这是用 UDEC 模拟开挖的最常用的方式。在此假设开挖是突然进行的(即通过爆破 掘进)。产生模型的非线性响应取决于卸载速率,所以,模拟者必须决定这种模拟开挖 是否适用于实际情况。另一方面,也可通过逐步减小开挖边界的应力来实施地下开挖。 可产生不同的响应。对于与荷载路径的影响将在 3.11.3 节中进一步讨论。 通过监测围绕隧道位置的位移来判断第二计算步的解: hist ydis 0,36 36.2 hist xdis 2.5 32.5 记录隧道拱顶 y 方向的位移和隧道两帮中点 x 方向的位移。求解由以下命令开始: step 5000 当完成计算步,绘制位移历史。图 3.30 显示位移已经收敛于拱顶为 7.5mm 和两帮 为 3.8mm 的常值位移。

图 3.30 监测隧道拱顶(历史 1)和边墙中点(历史 2)的位移历史 在图 3.31 绘制的模型图形显示了在开挖体每一边墙的一个块体已经脱离岩体。这些 块体在保存这个阶段的模型之前被删除。 Delete -2.9 3.2 31.9 32.6 Save step2.sav

-61-

图 3.31 隧道开挖的应力平衡

图 3.32 地震波周期为 1s 时不稳定块体 动力波是地震荷载的一种简单描述。10Hz 正弦剪切波被施加到模型基底,并允许 向上传播。在模型底部的最大速度是 2.5m/sec。对于动力分析阶段的命令见例 3.21 中: 这些命令施加一种正弦剪切波和粘性边界条件到模型基底,模型侧边施加自由边 界。对于解释这些边界条件参见理论与背景中的第 4.5 节。注意,为了补偿粘滞性边 界,作用的应力值是实际值的两倍(见理论与背景中的第 4.6 节)。

-62-

这个例子在动力计算阶段,启动了 Reyleigh 阻尼来模拟实际材料的阻尼。命令 DAMP 0.0001 10 施加阻尼为 10Hz 频率的 0.01%。参见理论与背景的第 4.3.2 节可指导 选择动力分析的阻尼参数。 正弦剪切波按 1s 的周期施加。对于 2s 的周期按 0.5s 增量进行计算。图 3.32 显示在 施加 1s 的动力荷载的隧道。由此显示了顶底板出现不稳定。 从地面产生的地震波使顶板块体失稳是合理的。正如在图 3.33 所示,在施加正弦波 荷载后,当反映波到达开挖体时,在隧道边墙中点(history 4)的位移值是 0.4s 时的约 2 倍。

图 3.33 地震荷载周期为 1s 时隧道拱顶和边墙位移历史图 为了研究支护对隧道的稳定性影响,通过安装锚索单元或结构单元来重复模拟顺序 (即 RESTORE step2.sav)。如果对块体和节理采用不同的材料性质,必须再次进行隧 道静态分析(RESTORY step1.sav)。如果评价隧道不同位置和节理组产状,有必要重 新产生模型和使其再次达到平衡状态。不要忘记,当改变荷载时,模型必须达到平衡状 态。

3.7 选择本构模型
本节概述 UDEC 所提供的块体和节理本构模型,并对模型的应用提出一些建议。在 理论与背景中的第 2 节中,为块体本构模型的计算公式提供了背景信息。在理论与背景 中的第 1.2.4.4 和 3.2 节和特殊特征的第 3 节介绍了节理模型。

-63-

3.7.1 变形块体材料模型
在 UDEC 中开发了 7 种块体材料模型: (1)开挖模型(null)(CHANGE cons=0 或 ZONE model null); (2)各向同性弹性模型(CHANGE cons=1 或 ZONE model elastic); (3)Drucker-Prager 塑性模型(CHANGE cons=6); (4)Mohr-Coulomb 塑性模型(CHANGE cons=3 或 ZONE model mohr); (5)堆砌节理模型(ZONE model ubiquitous); (6)应变软化 / 硬化模型(ZONE model ss); (7)双屈服模型(ZONE model dy); 注意,对于开挖模型,弹性模型和摩尔库仑模型,可采用两种方式的任何一种方式 来定义模型。可采用 CHANGE cons 命令给一个多个块体指定模型。模型性质参数可以 为材料性质号,而合适的材料号用 CHANGE mat 命令赋给块体,用 PROPERTY mat 命 令指定性质参数。可选择地,采用 ZONE model 命令,能够给一部分块体或块体组进行 赋值。在这种情况下,应用 ZONE 命令直接给单元指定性质。 进行每一块体模型的设计,表征地质材料所特有的本构特性。开挖模型用来代表从 模型中被移去的材料。各向同性模型对表现为线性应变特征的均质、各向同性的连续材 料是有效的。D-P 塑性模型是一种简单的破坏准则,在此,材料屈服是独立应力的函 数。M-C 模型假定材料受剪切屈服破坏,但屈服应力仅依赖最大和最小主应力。堆砌 节理模型与 M-C 模型一致,但适用于具有较强的各向异性特性的材料。应变软化模型 的基础是 M-C 模型,但适合于当剪切加载超越其极限,表现出剪切弱化的材料。双屈 服模型是应变弱化模型的扩展模型,用于模拟不可逆压缩以及剪切屈服破坏的材料。 UDEC 的材料模型主要应用于地质工程,即地下开挖、建造、采矿、边坡稳定性、 基础、土石坝。当为特殊工程分析选择本构模型时,应当注意以下两个问题: (1)模拟的材料具有何种特征? (2)模型分析的目的是什么? 表 3.2 提供了 UDEC 中具有代表性材料和可能应用模型的总结。弹性块体模型通常 用于沿不连续面滑移是模型失稳控制性机理的情况。当应力水平达到足以使完整岩体产 生破坏时,应考虑采用 M-C 准则。M-C 库仑模型的粘聚力和内摩擦角参数通常比其 他模型中所涉及的地质参数更容易获取。堆砌节理模型、应变硬化或软化模型和双屈服 模型实际上是 M-C 模型变形形式。如果这些模型中增加的材料参数赋予很高的值,则 由此所获得与 M-C 模型具有相同的结果。D-P 模型是比 M-C 模型较为简单的模 型,但一般并不适用描述地质材料的破坏。在 UDEC 中包含有 D-P 模型主要是与其他 数值分析程序进行比较。注意到,当摩擦角等于零时,M-C 模型将变成 Tresca 模型, 而 D-P 模型则变为 Mises 模型。

-64-

表 3.2 UDEC 块体本构模型
模 型 代表性材料 空 洞 应用实例 钻 孔、开 挖、待回填的空区等 荷载低于极限强度的人造材料(即钢铁), 安全系数计算 与有限元程序比较的通用模型 一般土或岩石力学问题(即边坡稳定性和地 下开挖) 峰后效应研究(即 渐进坍塌,矿柱屈服,地 下塌陷) 封闭的层状地层中开挖

开挖模型 弹性模型

均质、各向同性、连续、线性

D-P 塑性模型 M-C 塑性模型

低摩擦角软粘土,应用范围有限 松 散 和 粘 结 颗粒 材 料 , 土 、岩 石和混凝土 具 有 明 显 的 非线 性 硬 化 或 软化 的颗粒材料 材 料 强 度 具 有显 著 各 向 异 性的 薄层状材料 压力引起孔隙永久性减小的低粘 结性的颗粒材料

应变软化/硬化 M-C 模型 堆砌节理模型

双屈服模型

水力装置充填

D-P 模型和 M-C 模型是最有效的塑性模型。而其他的塑性模型需要为计算扩展 内 存 和 消 耗时 间 。 例 如, M -C 模型 并 不 直 接计 算 塑 性 应变 ( 见 理 论与 背 景 的第 2 节)。如果需要塑性应变,则必须采用应变软化或双屈服模型。这两个模型主要是用于 研究峰后破坏效应显著的材料,即矿柱屈服、崩落或充填研究。 抗拉破坏准则与 M-C 模型、堆砌节理模型、应变软化模型和双屈服塑性模型是一 致的。定义抗拉破坏准则与剪切和体积强度准则无关。抗拉强度准则由岩体的抗拉极限 强度确定,峰后破坏由相关联塑性流动准则控制。当出现抗拉破坏时,对于 M-C 模型 和堆砌节理模型,给定的抗拉强度为零(瞬时弱化)。张拉软化可以用应变硬化或软化 或双屈服模型控制。在这些模型中,如果没有赋值,则缺省的抗拉强度是零。

3.7.2 节理材料模型
有四个模型和一个选择模型可用于表征不连续性特征: (1)点接触-库仑滑移(CHANGE jcons=1 或 JOINT model point); (2)节理面接触-库仑滑移(CHANGE jcons=2 或 JOINT model area); (3)节理面接触-具有残余强度库仑滑移(CHANGE jcons=5 或 JOINT model residual); (4)连续屈服(CHANGE jcons=3 或 JOINT model cy); (5)Barton-Bandis 节理(CHANGE jcons=7 或 JOINT model bb,选择模型)。 可采用两种中的任一种方式给出节理模型。通过应用 CHANGE jcons 命令,指定模 型 一 个 或 多 个 接 触 面 。 然 后 , 用 PROPERTY jmat 命 令 为 材 料 性 质 号 , 并 通 过 用 CHANGE jmat 命令指定材料号来给节理模型性质赋值。另外,采用 JOINT model 命令

-65-

可以给节理模型的单一或一组接触面赋值。在这种情况下,JOINT 命令可直接赋给接触 面赋性质参数。 节理本构模型是用来表征实际工程岩体节理而开发的。点接触模型描述两块体间的 接触面相对于块体的尺寸是非常小的接触的情况。节理面接触模型用于具有面接触的封 闭块体的表征。这个模型提供了节理刚度和屈服极限的线性描述,建立在岩体的节理弹 性刚度、摩擦性质、粘聚力、抗拉强度和剪胀特性的基础上开发的。该模型的残余强度 版本采用摩擦力、粘聚力和(或)抗拉强度在开始出现剪切或张拉破坏参数减小或消失 来模拟节理的位移软化特性。连续屈服节理模型是一个比较复杂的模型,他考虑了性质 参数为累积塑性剪切位移连续的函数关心,模拟节理面连续弱化的特性。Barton-Bandis 模型也是一非线性模型,他直接利用由挪威地质研究所 Barton 和 Brandis 博士推导的室 内节理试验性质指标参数(Bandis 等,1983)。 表 3.3 对 UDEC 节理模型、典型材料以及可能应用范围所进行概括。面接触库仑滑 移模型是一般工程研究最常用的模型。库仑摩擦和粘聚性质比其他节理性质常常更容易 获得。如果预测模型块体间所发生的点接触和面接触,仅指定面接触模型是可以接受 的。当最小接触面积等于模型圆角长度的两倍,则模型自动赋予点接触。如果点接触是 控制性接触,即松软的无规则形状颗粒封闭系统,应当指定点接触模型。 连续屈服模型和 B-B 节理模型属于经验性描述,在此需要比较详细的节理特性参 数。连续屈服模型性质来自节理剪切应力和法向位移的室内节理试验结果。如上所述, B-B 节理性质参数由节理指标试验确定。通常建议,为了探索对节理效应的本质理 解,在应用较为复杂的节理模型前,开始首先应采用库仑滑移模型。关于此问题在下节 进一步讨论。在验证问题中第 1 节,给出了库仑滑移模型与连续屈服模型分析比较的例 子。在理论与背景的第 3 节给出了连续屈服模型的响应示范和需要的性质参数。 表 3.3 UDEC 节理本构模型
模 型 点接触 典型材料 应用有限,颗粒材料,无规则形状的 松散挤压块体 岩体中的节理,断层、层面 显现明显的峰值/残余强度特性 表现渐进损伤和滞后特征的岩体节理 应用实例 破碎和断裂岩体,受强扰动边坡的稳 定性 一般岩石力学问题(即地下开挖) 一般岩石力学问题 具有显著的滞后循环加载和反向加载; 动力分析 评价节理岩体的渗透特性

面接触 位移弱化的面接触 连续屈服

Barton-Brandis

由 Barton-Bandis 指标性质定义的岩 体节理

3.7.3 合理模型的选择
任何问题的分析应当从简单的块体和简单的节理模型开始。在大部分情况下,应首

-66-

先考虑弹性 块体模型( cons=1 或 ZONE model elastic)和节理面接触 库仑滑动模 型 (jcons=2 或 JOINT model area)。弹性块体模型仅需要体积密度、体积模量和剪切模量 3 个参数(见 3.8.1.2 节)。库仑滑移模型需要 6 个参数:法向和剪切刚度、内摩擦角、 粘聚力、抗拉强度和剪胀角。在 3.8.2 节给出了这些参数的估算和参考资料。这些材料 模型提供了应力变形特性的简单透视。其分析结果有助于用户是否采用复杂的还是简单 的本构模型来描述块体或节理特性。例如,如果块体的应力和变形与节理位移相比较是 小的,则简单的刚性块体模型可能就满足了。 在采用 UDEC 求解全空间的边界值问题前,选择材料模型进行简单的试算常常是有 帮助的。这能够洞察模型响应,并与已知实际材料响应比较。 下面的例子说明简单测试模型的应用。该例子是地下开挖围岩节理滑移分析。为评 价库仑滑移模型和连续屈服模型表征承受剪切荷载的节理的响应的合理性,为此建立了 一个简单的模型。该模型是对单一水平节理构成的直剪试验的模拟:它先是受到法向常 应力作用,然后产生单向剪切位移。图 3.34 显示这个模型。该节理由 4 个接触面定义。

图 3.34 直剪试验模型 首先,施加具有代表性的 10MPa 法向应力。然后,在块体端部施加水平应力产生剪 切位移。此值也具有代表性的。为了说明这个目的,仅施加小于 1mm 的剪切位移。 采用 FISH 函数(av_atr)记录了节理的平均法向和剪切应力、法向和剪切位移。基 于此记录结果,我们能够确定不同模型的峰值和残余剪切强度以及剪胀特性。例 3.22 包 含了应用库仑滑移模型进行该试验模拟的数据文件。 图 3.35 和图 3.36 分别给出了节理的平均剪应力与位移的关心曲线和平均法向位移

-67-

对剪切位移的图形。这两个图形表明:对于给定的模型性质和边界条件,节理已经发生 滑移。在峰值剪切强度接近 6MPa 前,图 3.35 中的加载斜率为线性。正如在图 3.36 所显 示的,当节理剪切位移约为 0.15mm 时,不能抗剪,节理开始剪胀。其剪胀一致延续到 其极限剪切位移(zdilation=0.4mm)。最大的平均剪胀接近于 0.027mm。

图 3.35 平均剪切应力与剪切位移的关系

图 3.36 平均法向位移与剪切位移的关系

-68-

正如上述结果所示,库仑滑移模型(jcons=2)仅限于节理的极限剪切强度,在节理 开始滑移后才发生剪胀。在剪胀极限范围内,剪切位移是剪胀角成线性函数。这些函数 在理论与背景的第 1.2.4 节中加以描述。 在库仑滑移模型中也可以获得节理特性的其他调整。例如,当节理滑移时,通过包 含有效内摩擦角的剪胀特性,近似逼近位移弱化特性。该特性通过在 jcons=2 的模型中 应用 SET add_dil on 命令进行激活(见例 3.22)。通过将剪胀角加上内摩擦角的输入参 数给出有效内摩擦角。当剪胀角为零达到极限剪切位移时即为残余强度。图 3.37 和图 3.38 给出了如何考虑这些选择特性。注意,在模型中仅调整内摩擦角。

图 3.37 平均剪切应力与剪切位移的关系

还有其他方式修正库仑模型用于模拟节理剪切强度的峰值和残余强度。一种方式就 是应用 jcons=5 或 JOINT model residual 的库仑模型。在例 3.22 中,在对节理进行剪切 试验前,添加下一行(注意,对于初始法向加载,我们仍应用 jcons=2 模型,这阻止在 初始平衡计算过程中节理强度的降低)。 change jcons=5 set jcondf=5 &

prop jmat jkn=40000 jks=40000 jfric=35 jrfric=10 jcoh=2.0 jrescoh=0 jdil=6 zdil=4e-4 & &

-69-

图 3.38 平均法向位移与剪切位移的关系 图 3.39 显示的结果说明了模型中峰值和残余强度的变化。能够给模型赋予内摩擦 角、粘聚力、抗拉强度的峰值和残余强度。在此例减小内摩擦角和粘聚力。注意,强度 下降是突然发生的。如图 3.40 所示,对于相同的极限剪切位移,用 jcons=5 模型计算的 最大剪胀角是较低的,这是因为在出现较大的剪切位移后,节理开始破坏(比较图 3.40 和图 3.36)。

图 3.39 平均剪切应力与剪切位移的关系

-70-

图 3.40 平均法向位移与剪切位移的关系 另外,用户定义的函数能用于降低节理强度。由此采用 FISH 来处理节理性质(即 节理摩擦角和粘聚力)。在例 3.23 中,其性质参数作为剪切位移的线性函数予以减小。 注意,JOINT model area 命令用来指定局部接触的节理性质。然后,对于每一接触面, 通过 c_jex(ic)FISH 变量给节理模型和性质参数赋值,在此,ic 是接触面指标号。对 于储存剪切位移的分支索引,采用符号名$KSC,获取剪切位移。接触面数组指标在文 件 CONTACT.FIN 中定义。节理性质索引指标符号在文件 JMAT.FIN 中定义。例如,在 指定的接触面时,$AC_PHI 是节理摩擦角,$AC_COH 是节理粘聚力(见 FISH 函数 的第 4 节)。 节理弱化特性可用连续屈服节理模型(jcons=3 或 JOINT model cy)自动考虑。该 模型模拟了节理受剪切作用过程的渐进损伤过程。例 3.24 包含了应用连续屈服模型的直 剪试验的数据文件。 在这些例子中,库仑滑移模型和连续屈服节理模型材料性质的选择,是考虑了节理 剪切强度和节理承受单向剪切的剪胀具有大致相同的响应。对于实际的应用,选取参数 应能模拟在期望的加载条件下节理的响应。在大部分情况下,库仑模型参数容易获取; 同时,也可对库仑模型可进行简单的修正,足以接近节理的特性。 对于其他节理模型,如连续屈服模型和 B-B 节理模型,性质参数较为复杂的。为 了使用连续屈服模型,有必要运行一系列剪切试验来拟合模型参数。如果可以获得 B- B 模型的节理指标,则 B-B 的节理模型可能是合适的。建议忽略节理模型的选择,进 行简单的剪切试验,在考虑的问题条件,具有期望的节理特性。

-71-

如果有必要模拟复杂的节理响应,则可能需要一个比较复杂的节理模型。然而,在 进入较复杂的模型之前,首先应用一个简单的模型,为评价复杂节理特性的影响建立基 础。例如,几个研究者描述了在室内剪切试验,节理承受循环剪切作用的复杂节理特 性。这些试验显示了两个共同的特性: (1)剪切强度的变化依赖逆向荷载; (2)在逆向剪切期间,节理剪胀反弹。 应用库仑模型可以获得对上述效应的近似模拟。如果剪切位移增量与总剪切位移方 向相同或相反,则在 SET add_dil on 命令(前面描述)中,在输入的内摩擦角加上或减 去剪胀角。在例 3.25 所给出的数据文件,能够应用 SET add_dil on 命令,进行循环剪切 试验。 通过这个模型模拟的节理特性可能足以评价循环加载的影响。然而,如果加载导致 节理的渐进损伤,则库仑模型就不合适了。因为该模型不能涉及这一特性。连续屈服模 型确实可以模拟渐进损伤,但不能模拟节理剪胀在逆向剪切过程中的反弹。也可以考虑 B-B 模型。模拟者必须决策哪一模型是最合适的,或是否必须开发新的公式。

3. 8 材料性质
UDEC 需要完整岩块和不连续面的材料性质。本节对用于描述节理岩体特性的典型 参数进行概述以及对选择合适参数给予指导。还有一些特殊问题,例如峰后破坏特性的 定义以及室内试验性质向现场参数的外延。这些问题也给予讨论。 参数的选择常常是产生模型中最困难的问题。因为性质参数存在很大的不确定性。 应当记住,当进行分析时,尤其涉及地质力学问题,总是涉及到有限数据系统。现场信 息决不可能完全已知。然而,基于获得的信息,洞察实际的物理问题,仍可获得合理的 性质参数。该研究在 3.11 节中进一步讨论。

3.8.1 岩块性质
岩块性质参数通常由实验室试验确定。下列 4 节将描述岩块的基本性质(实验室结 果),并给出各种岩石的常用参考值。 3.8.1.1 质量密度 UDEC 模型的所有非空介质都需要质量密度。该性质具有质量除以体积的单位,并 不涉及重力加速度。在很多情况下,给出材料单位重量。如果用力的单位除以体积给出 的单位重量,则作为输入 UDEC 的数据之前,该值必须除以重力加速度。 3.8.1.2 基本变形性质 UDEC 的所有变形块体的材料模型,有两个弹性常数,体积模量 K 和剪切模量 G。 在此描述的弹性材料,均假设为各向同性材料特性。UDEC 中采用的是两个常数 K 和 G 而不是杨氏模量 E 和泊松比 ν 。因为可以确信,描述材料最基本特性的是体积模量和剪

-72-

切模量而不是 E 和 ν 。 (E、 ν )和(K、G)的变换公式为:

K=

E E , G= ( -2ν) 31 ( + ν) 21

(3.15)

当 ν 接近于 0.5 时,方程(3.15)就失效。由于 K 的计算值将达到不符合实际的高 值,而解的收敛速度也变得十分缓慢。如果能较好的确定实际的 K 值,则可从 E 和 ν 计 算 G 值。某些典型的弹性常数值列于表 3.4 中。

3.8.1.3 基本强度性质 UDEC 的岩石材料破坏基本准则是摩尔-库仑关系,在此对应的是剪切破坏的线性 破坏面:

f s = σ 1 ? σ 3 N φ + 2c N φ
式中, N φ = (1 + sin φ ) /(1 ? sin φ )

(3.16)

σ 1 - 最大主应力;

σ 3 - 最小主应力;

φ - 内摩擦角;
c - 粘聚力。 如果 f s < 0,就发生剪切屈服。两个强度常数 φ 和 c 容易从实验室的三轴试验进行
导出。 当法向应力变为拉应力时,摩尔-库仑准则就失去实际意义。但为简化,屈服面扩 展到 σ 3 等于其抗拉强度张拉强度 σ 的区域。最小主应力不能超过抗拉强度:
t

ft = σ 3 ? σ t

(3.17)

如果 f t > 0,发生张拉屈服。岩石和混凝土的抗拉强度通常由巴西试验获得。注 意,抗拉强度不能超过 σ 3 的值,该值对应摩尔-库仑关系的上限。最大值由下式确定:

-73-

t σ max =

c tan φ

(3.18)

对于一组具有代表性岩石(样本),其粘结力、摩擦角和抗拉强度的代表性值列于表 3.5 中。

D-P 模型强度参数可从粘聚力和内摩擦角推导。例如,假设 D-P 破坏包络线限制 了 M-C 包络线,则 D-P 参数 qφ 和 kφ 与 φ 和 c 的关系为:

qφ = kφ =

6 3 (3 ? sin φ ) 6 3 (3 ? sin φ )

sin φ

(3.19)

c cos φ

(3.20)

两个参数的关系进一步解释,可见理论与背景的第 2 节。 对平面弱化的堆砌节理模型也需要强度性质参数。节理性质在下面的 3.8.2 节中讨 论。节理粘聚力和内摩擦角也在堆砌节理模型中应用。 3.8.1.4 峰后效应 在很多情况,尤其是采矿工程,材料在开始破坏后的响应是工程设计的重要因素。 因此,在模型中必须考虑材料的峰后特性。UDEC 定义了如下三种类型的峰后效应: (1)剪胀效应; (2)剪切硬化或软化效应; (3)抗拉软化效应。 如同 M-C 关系所定义,这些性质只有材料开始破坏后才出现。在 M-C 模型中考 虑了剪胀效应;在堆砌节理模型和双屈服模型中,考虑了应变硬化或软化效应。在 M-

-74-

C 模型和堆砌节理模型中,还给出了瞬时张拉软化(当抗拉破坏出现时,抗拉强度即为 零)效应。 (1)剪 胀 剪胀效应是指随着材料发生剪切而出现体积变化的特性。剪胀效应由剪胀角, Ψ 表 征,它涉及塑性体积变化与塑性剪切应变的比值。UDEC 的塑性材料模型也可给出其 角。确定剪胀角典型的方法是三轴试验或剪切盒试验。例如,根据 M-C 破坏剪切面, 基于三轴试验所获得剪胀角的理想关系如图 3.47 所示。由体积应变与轴应变的关系曲线 可求出剪胀角。注意,该图的初始斜率对应弹性区域,而用于计算剪胀角的斜率对应于 塑性区域。

图 3.47 取自三轴试验的剪胀角的理想关系 对于土、岩石和混凝土一类材料,剪胀角通常明显小于材料的内摩擦角。Verneer 和 Borst(1984)报道了下列典型的剪胀角的值: 密度砂 15o 松散砂 < 10o 正常固结粘土 0o 砾石和完整大理岩 12o ~20o 混凝土 12o Verneer 和 Borst 发现,无论是土、岩石还是混凝土的剪胀角,其值大致位于 0 o ~ 20o 的 范围内。UDEC 所有的本构模型,剪胀角的缺省值是零。 剪胀角也可在堆砌节理模型中赋给节理。该参数典型的获取方法采用节理的直剪试 验。常用的值可参考 3.8.2 节。

-75-

(2)剪切硬化或软化 一旦出现塑性屈服,材料的塑性硬化或软化就开始渐进发展。在破坏开始时刻,混 凝土和岩石内产生微裂缝或土颗粒的滑移,变形特性越来越非弹性化。这也导致材料强 度的蜕化以及初始剪切带的产生。这种涉及局部弱化现象在 3.11 节中进一步讨论。 在 UDEC 中,通过构造 M-C 模型材料参数(随剪胀过程的粘聚力和摩擦角)作为 塑性应变的函数来模拟剪切硬化和软化特性。这些函数可以从应变软化和双屈服模型中 获取,并采用 TABLE 给予赋值。 对于每一特定的分析,通常应用室内三轴试验结果反演硬化和软化参数值,并进行 标 定 。 这 通常 是 一 个 反复 的 过 程 。研 究 者 已 经对 硬 化 或 软化 特 性 进 行改 进 。 例 如, Verneer 和 Borst(1984)提出了摩擦硬化关系:

sin φ m = 2

epe f epe f

sin φ

对于e p ≤ e f

sin φ m = sin φ
式中, φ - 峰值摩擦角;

对于e p > e f

(3.21)

φ m - 滑动摩擦角;
e p - 塑性应变; e f - 硬化常数。
数值试验条件也影响模型剪切硬化或软化特性效果。加载速率产生惯性效应,可通 过监测不平衡力和减小加载速率加以控制。FISH 函数可以用来自动控制加载速率。计算 结果与网格也有关。因此,无论在何种情况下,考虑剪切硬化或软化特性进行分析时, 考虑不同单元尺寸和网格方位来评价模型特性是十分重要的。 (3)体积硬化或软化 体积硬化对应不可逆压缩。增加压力可引起体积永久性减小。该特性在诸如胶结 砂、砾石和水力冲积物等材料是共有的。 体积的应变硬化或软化也可在双屈服模型中考虑。该模型假设硬化仅依赖于塑性体 积应变。 在双屈服模型中,采用 TABLE 考虑材料的硬化特性。尽管这些代表性的数据不能 给出,但理论与背景的第 2.4.5.5 节给出了利用三轴试验获取这些参数的试验方法。该三 轴试验的侧向压力保持不变,轴向应力采取常均应力加载。另一可考虑的方法是从单轴 应变试验反演参数。 Clark(1991)为矿山充填进行了一系列单轴试验,提出了体积屈服面的硬化准 则。用一般的形式表示如下:

-76-

? ev ? p o cp = W ? ? + ep H ? ev ? ? p ? ?
式中, c p - 当前屈服面面压力;

α

(3.22)

c 0 - 初始屈服面压力 p c v - 塑性体积应变; p H,W,α - 参数。
根 据 Clark ( 1991 ) 进 行 的 单 轴 应 变 试 验 的 结 果 进 行 拟 合 , 得 到 H=0.28 , W=1.15×107 和 α =1.5。由此上述公式为:

? ev ? p c p = 1.15 × 10 ? v ? ? 0.28 ? e p ? ? ?
7

1.5

+ 10 4 Pa

(3.23)

(4)张拉软化 张拉破坏一出现,材料的抗拉强度通常降至零。因此,在 M-C 模型中,对于出现 张拉破坏的单元,抗拉强度就赋零(即瞬时软化)。抗拉强度降低的速率,或出现张拉 软化也可通过塑性张拉加以控制。该函数可从应变软化模型中获取,并可应用 TABLE 模量进行定义。 一个简单的张拉试验说明了脆性张拉破坏。这个模型是采用矩形块体和 M-C 材料 构成的一个张拉试验。样本的端部以常速度张拉。 轴应力与轴应变的关系如图 3.48 所示。在模型顶板,平均轴向应力降低到零。减小 模型的宽度导致开始出现张拉破坏,然后,随着出现张拉软化而张开(图 3.49)。

图 3.48 对张拉弱化材料进行张拉试验中的轴向应力与轴向应变的关系

-77-

图 3.49 对张拉弱化材料进行张拉试验中的正交应变与轴向应变的关系 注意,采用缺省阻尼(局部阻尼)对减小当张拉破坏突然断裂时出现的波动是合适 的。另外,如果通过 DAMP auto 命令,给出适当的总体阻尼,则在应力应变图上可以 观察到这种振动。用合适的总阻尼,阻尼参数随着模型的拉伸而连续减小。当出现张拉 破坏,总阻尼参数是低的,可能产生影响求解结果的波动。借助于局部阻尼,从结点到 结点的阻尼变量值与不平衡力成正比。该阻尼减小了在突然发生抗拉破坏的波动。 张拉软化的脆性可以通过用应变软化模型代替 M-C 模型,由塑性张拉应变函数控 制。当采用剪切软化模型时,张拉软化模型必须对每一特定问题和网格尺寸进行标定, 由于其结果与网格划分相关。 3.8.1.5 现场性质参数的外延 UDEC 模型中的材料性质应尽可能接近物理问题的实际参数。室内试验参数通常不 能直接用于 UDEC 的计算模型。模型中的不连续面的存在,说明了尺寸效应对参数影 响。然而,某些块体性质仍可能需要判断,评价裂隙、层理以及其他地质不连续面等不 均质体对岩体性质的影响。 为获得现场的岩体性质,目前已经提出几个经验研究结果。现对某些公认的方法进行 讨论。 岩体的变形通常由变形模量, E m 定义。如果包含有一组相对平行、连续和等间距 的节理,那么可以将岩体处理为横观各向同性的连续体来估算 E m 值。3.8.2 节中的关系

-78-

可用于计算在垂直于节理组方向的 E m 值。对于多于一组节理组的情况,也能够估算其 变形模量。参考在 3.8.2 节中所提供的多组节理的解。 实际上,岩体结构常常是无规则或没有足够的数据用于上述研究。通常可从现场压 缩试验所获得的力-位移曲线来确定 E m 。这样试验包括承压板试验、千斤顶试验和膨 胀仪试验。 Bieniawski(1978)根据遍及世界各地的现场试验结果,提出了计算 E m 的经验公 式。该公式是根据岩体分类指标(RMR)。对于指标值大于 55,试验数据拟合可给出 以下的近似公式:

E m =2(RMR)-100
式中, E m 的单位是 GPa。

(3.24)

对于 E m 的值位于 1~10GPa,Serafim 和 Pereira(1983)给出较好的拟合,即

E m = 10
验岩体强度准则为:

RMR ?10 40

(3.25)

预测岩体强度为人们所接受的方法是由 Hoek 和 Brown(1980)提出的。提出的经

σ 1s = σ 3 + (mσ cσ 3 + sσ c2 )1 2
式中, σ 1s - 峰值强度的最大主应力;

(3.26)

σ 3 - 最小主应力;
m和s - 常数;

σ c - 完整岩石的抗压强度。
岩体无侧限抗压强度由下式给出:

qm = σ c s1 2
岩体的单轴抗拉强度为

(3.27)

σ t = σ c [ m ? ( m 2 + 4 s )1 2 ]
表 3.6(Hoek 和 Brown,1988)提供了受扰动和无扰动岩体的 m和s 的代表值。

1 2

(3.28)

由 H-B 准则(1990)也可以计算 M-C 模型的摩擦角和粘聚力。对于给定的 σ 3 ,函 数方程 3.26 的切线将给出了等效的 M-C 屈服准则具有下列形式:

σ 1 = N φ σ 3 + σ cM
式中, N φ =
M
1+ sin φ 1? sin φ

(3.29)

= tan 2 ( φ + 45 o ) 。 2

代入, σ c 为

σ cM = σ 1 ? N φ σ 3 = σ 3 + (mσ cσ 3 + sσ c2 )1 2 ? N φ σ 3 = σ 3 (1 ? N φ ) + (mσ cσ 3 + sσ c2 )1 2
式中, σ c - 对应于 σ 3 的岩体单轴抗压强度。
M

-79-

方程 3.26 的切线由下式定义:

N φ (σ 3 ) =
M

σ cm ?σ 1 = 1+ ?σ 3 2 σ 3σ c m + sσ c2

(3.30)

从 N φ 和 σ c ,可以获得粘聚力(c)和内摩擦角( φ ):

φ = 2 tan ?1 N φ ? 90 o
c=

(3.31) (3.32)

σ cM
2 Nφ

-80-

3.8.2 节理性质
节理性质通常从室内试验获得(即三轴或直剪试验)。这些试验能够给出节理诸如 内摩擦角、粘聚力、剪胀角、抗拉强度以及节理切向和法向刚度等力学参数。节理粘聚 力和内摩擦角对应于 M-C 强度准则中的参数。 对 于 夹 有 软 土 和 淤 泥 的 岩 石 节 理 , 其 切 向 与 法 向 刚 度 的 代 表 值 大 约 从 10 到 100MPa/m。而在花岗岩和玄武岩中的闭合节理,其值超过 100GPa/m。关于岩石节理刚 度性质发表的文章是有限的。 近似刚度值可从节理岩体变形特性、节理结构以及岩石的变形特征推求。如果节理 岩体假设基于等效弹性介质的变形响应,则可以推求出节理岩体性质的等效连续岩体性 质的参数。 对于含有一组产状垂直于加载方向,且等间距的节理组,具有下面关系:

1 1 1 = + Em Er k n s


kn =
式中, E m - 岩体的杨氏模量;

Em Er s( E r ? E m )

(3.33)

E r - 岩石的杨氏模量;

k n - 节理的法向刚度;
s - 节理间距。
对于节理剪切刚度,类似表达式也能够给出:

ks =
式中, G m - 岩体的剪切模量;

Gm Gr s (Gr ? Gm )

(3.34)

Gr - 岩石的剪切模量;

k s - 节理的切向刚度。
当扩展到三组正交节理组时,基于等效连续假设可给出下列关系;

? 1 1 Ei = ? ?E + sk i ni ? r

? ? ? ?

(i = 1, 2, 3)
? ? ? ? (i,j = 1, 2, 3)

(3.35)

? 1 1 1 + + Gi j = ? ?G ? r s i k si s j k sj

(3.36)

限制最大节理刚度对 UDEC 模型的应用是合理的。如果实际法向和剪切刚度小于 10 倍相邻单元的等效刚度(在 3.9 节中的式 3.37),则实际问题中应用不存在问题。如 果该比值大于 10,求解时间将显著增长,对其比值限制到 10,则系统的特征没有显著 影响。应指出的是,为了提高求解效率,应减小切向和法向刚度。如果法向刚度, k n 非

-81-

常低,还可能产生块体嵌套问题。应对节理的法向位移做出大致的估算,此位移取决于 施加于系统的代表性应力( u = σ / k n ) 。该位移与典型的单元尺寸相比较小,如果与相 邻单元的尺寸比较大,比如说 10%,那么单元序号存在误差或者刚度应当被提高。 已经发表的节理强度参数远比节理刚度更容易获得。摩擦角的变化范围可从小于 10 到超过 50 o ,分别对应于软岩(如凝灰岩)的光滑节理和硬岩(如花岗岩)的粗糙节 理。节理粘结力也可能从 0 到接近于围岩的抗压强度。 室内典型试验所测试的节理性质并不能代表现场实际节理参数。与尺寸相关的节理 性质是岩石力学的主要问题。选择合理参数的仅有的方式是比较与现场试验类似的节理 性质。然而,现场试验结果是极其有限的,某些结果由 Kulhawy(1975)报道。
o

3.9 提示和建议
当用 UDEC 求解问题时,重要的事情是为有效地进行问题求解而进行模型优化设 计。本节为改进运行模型提出一些建议。而且,还给出了进行 UDEC 计算时应避免可能 出现的一些陷井。

3.9.1 节理几何形状的选择
对于选择输入到 UDEC 模型中的节理几何形状是分析中的关键一步。对具有 10 至 数百条节理的多裂隙岩体,考虑所有的节理结构是不现实的。通常,为实际分析产生具 有合理大小的模型和运行速度,仅很小比例的节理能够被输入模型中。因此,为模拟力 学响应,模拟者必须对节理的几何数据进行过滤,选择关键节理。通常从钻孔和节理图 开始,需要一个反复过程。选择关键节理的困难集中在确定性的 UDEC 模型中如何表征 这些统计数据。对于多组节理组,应当假设给定的一组节理,其产状是相同的,且变化 较小(比如说小于 10o 到 15o)。节理岩体范围可由特定的节理被划分成子区域,且在子 域内的节理被认为是连续的。这是对有限的长度和不连续性特性影响的一个上限估计。 下一步是应用一个滤波准则来识别那些在给定的荷载作用下,容易滑移或张开的节 理。该准则可能涉及范围从:(1)识别是否提供充分的自由度允许滑移,到(2)利用 块体理论来检查节理位移的可动性,到(3)与现场观察和记录进行比较 - 例如,从微 震监测识别关键节理。

3.9. 2 设计模型
模拟者总是试图建立尽可能考虑详细地质结构的 UDEC 模型。反对这种研究的主要 理由在于:(1)期望采用充分的数据模拟节理岩体的每一细节是无用的;(2)对于工 程问题,计算机硬件需要一个快速模型,接近典型的工程实际;(3)最重要的是,考 虑得较为详细节理对模拟结果的理解与控制的效率会更低。 产生 UDEC 模型时应考虑两个方面:首先是不连续分析是否真正的需要。这在大部

-82-

分情况下,取决于研究实际模型尺寸与节理的平均间距之比。例如,如果包含一组平均 间距不大于 1m 的岩体中开挖,开挖的最小尺寸是 10m,然而,用堆砌节理材料模型进 行连续分析可能是合理的。在这种情况下,连续分析产生的响应与模拟节理的计算结果 在总体上是等效的。用 UDEC 模型分析给出破坏机理较为详细分析,但比连续分析花费 更多的计算时间。在实际模型尺寸与节理间距之比值大于 10:1 的情况下,连续介质分析 是完善的。当存在连续分析是否能代表不连续效应时的问题时,应采用连续和非连续介 质计算。 第二个考虑是包含详细地质结构的计算模型的范围。最关键节理结构的详细描述 (见上述第 1 点)通常仅需要包含感兴趣的有限区域,即在几倍隧道半径的围岩范围 内。一般地,考虑详细节理的范围从感兴趣的区域扩展到足够的距离,包围可能产生破 坏的区域。详细的地质结构应在延伸到节理滑移和张开的范围之外。

3.9.3

检查模型运行时间

UDEC 的求解时间是模型中刚性块体或变形块体结点和接触面数的函数。如果在模 型中几乎没有接触面,则计算时间与 N
32

成正比(在此, N 是刚性块体数或变形块体

结点数)。该公式对弹性问题也成立。对于塑性问题,运行时间略有变化,但没有实质 上的差异。 求解时间随模型中接触面的增加而增加。重要的是检查你的计算机在求解特定问题 的计算速度。最好的方式就是运行在 5.1 节所给出的测定试验。然后,采用这个速度, 基于结点数和接触面数进行插值计算,估算特定模型的计算速度。

3.9.4 对允许时间的影响
如果存在下列情况,UDEC 将需要运行很长时间才收敛: (1)块体材料与节理材料的刚度或性质反差较大; (2)块体或单元尺寸存在较大差异。 当上述反差变得较大时,程序的计算效率很低。刚度反差的影响应在详细分析前进 行研究。例如,节理法向和切向刚度的力学计算应小于 10 倍节理毗邻单元的等效刚度。 即

? ? K + (4 / 3)G ? ? k n 和 k s ≤ 10.0 ?max ? ?? ? Δz min ?? ?

(3.37)

式中, K 和 G 分别是块体材料的体积模量和剪切模量; Δz min 是毗邻节理单元在法线方 向的最小宽度。如果节理刚度大于 10 倍等效刚度,则模型求解时间将明显大于其刚度比 限制在 10 以内的情况,而在系统特征上不会发生有意义的变化。 在另一方面,如果法向刚度, k n 很低,还可能存在问题。应对模型代表性应力所产 生的节理法向位移有一个大致的估计( u = σ / k n )。该位移与代表性的单元尺寸相比是

-83-

小的。如果该值大于毗邻单元尺寸的 10%左右,则一个单元存在错误或刚度应提高。

3.9.5 单元密度的考虑
UDEC 模型的变形块体采用的是常应变单元。如果具有很高梯度的应力或应变,描 述 变 量 分 布则 需 要 很 多单 元 。 可 采用 不 同 单 元密 度 运 行 同一 问 题 来 检查 其 效 果 。在 UDEC 中采用常应变单元是因为当模拟塑性流动时,较多低阶单元比较少高阶单元具有 更高的精度(见理论与背景的第 1.2.5 节)。 期望使单元尽可能的均匀,尤其在兴趣的区域,尽量避免单元的边长比大于 5:1 的 长瘦单元。

3.9.6 检查模型响应
UDEC 显示系统的特性,应进行简单试验,检查你所做的和你所期望得到的结果。 例如,如果加载条件和几何形状是对称的,检查其响应对称性。在对模型进行改变后, 运行很少几步(比如说 5 或 10 步)来验证初始响应在适当的位置显现出正确的符号。 如果施加给模型一剧烈振动,应获得一剧烈响应。如果给予模型一种非常合理的的 事情,必将期望一个奇怪的结果。如果在给定的分析阶段获得意料之外的结果,重新查 看所施加该阶段的计算步。 在用模型模拟之前,关键的是考察输出。例如,对于施加块体角点较大速度,如果 仍显示合理的预测结果,则不要继续进行,直至对此做出合理的解释。在此情况下,可 能没有适当地固定边界角点。

3.9.7 检查块体接触
当产生复杂的节理模型时,STEP 0 命令可用在 PLOT block contacts 之后,显示接 触位置。确信所有的接触存在于块体之间,并且块体本身并不接触。例如,下列数据文 件将一个圆形隧道与模型边界连接,并由此与接触面产生一个块体: rou 0.1 block 0,0 0,10 10,10 10,0 tunnel 5,5 2 16 crack 0,5 3,5 plot block contact 环绕隧道的块体是与本身接触。在此都不产生接触。CRACK 命令应当完全延伸到 模型以阻止这种情况。

3.9.8 应用体积模量和剪切模量
对于 UDEC 的弹性性质,采用体积模量 K 和剪切模量 G 比杨氏模量 E 和泊松比 ν 更有好处。

-84-

对于不违背热动力学原理的所有弹性材料,一对( K , G )参数都是可行的。但对 于某些可接受的材料,一对( E , ν )值就失去意义。一个极端,我们使材料发生体积 变化但没有剪切。另一极端材料发生剪切但没有体积变化。第一种类型材料对应有限的

K 值和零 G 值;第二种材料对应零 K 值和有限的 G 值。然而,一对( E , ν )并不能
够表征上述两种类型材料的任何一种。如果我们排除两种极限情况(通常, ν =0.5 和

ν =-1),则方程

3K (1 ? 2ν ) = E 2G (1 + ν ) = E

(3.38)

涉及两组常数。然而,只有当我们所研究的参数接近(但没有达到)这种极限情况下, 方程才成立。对于物理试验,我们并不需要涉及可能或不可能行不通的问题。该方程是 两种定义比例系数方式的简化结果。假设我们有一种逐渐较小抵抗变形的力,但仍保持 常体积变化材料。在此情况下, ν 接近于 0.5。方程 3K (1 ? 2ν ) = E 仍必须被满足。有 两种可能性(关于代数基础,而不是物理的争论):或者 E 保持为有限值(非零)和 K 趋向于任意大值,或 K 保持为有限值和 E 趋向于零。我们排除第一种可能,因为,对于 所有的材料,都存在压缩的极限值(即对于水是 2GPa,在此,泊松比为 ν =0.5)。这只 可能是第二种情况,即使我们假设材料的主要弹性抗力模式是卸载,E 值也戏剧性变 化。我们推论,用( E , ν )表达材料特性是不合适的。

3.9.9 选择阻尼
在大部分情况下,建议在静力分析中采用 DAMP local 。这是在 UDEC 中缺省的阻 尼模式,通常对于静力分析是合适的。在理论与背景的 1.2.6 节给出解释。而且,正如 在例 3.26 中所看到的,局部阻尼是适宜极小化由于模型突然发生破坏所出现的波动。 在某些情况下,尤其当计算初始平衡状态时,应用 DAMP auto 具有较好的计算效 率。正如在理论与背景的第 1 节所讨论的那样,当结点速度分量周期性的通过零点,局 部阻尼是最有效的。因为,体积调整过程取决于速度的符号变化。换句话说,采取总体 阻尼施加不受速度符号变化影响的常量阻尼引子。如果在一个方向的速度起主导性作用 (即由于重力加载),则采用局部阻尼比总体阻尼可能需要消耗更长的时间。正是由于 该原因,例 3.17 和例 3.20 采用 DAMP auto 。当存在疑虑时,最好分别采用 DAMP local 和 DAMP auto 进行计算,比较达到收敛所需要的计算步。

3.9.10 给块体和节理模型指定模型和赋值
UDEC 中可以获得 ZONE 和 JOINT 命令,使得用户可以给单元和接触面直接指定 模型参数而不是性质编号。如果模型性质参数发生改变或者如果用户希望用 FISH 控制 性质参数的变化,则应当采用 ZONE 和 JOINT 命令。如果性质上很少没有变化,则 CHANGE 和 PROPERTY 命令满足要求,且所需内存较少。

-85-

3.9.11 避免圆角误差
在 UDEC 中大部分计算都是单精度(即坐标优于位移)。通常,单精度可以获得有 意义的图形,32 比特计算机计算能够避免圆角误差对计算的影响。然而,可能存在圆角 误差分析结果的影响。例如,圆角误差使模型收敛速度很低,模型计算需要数以万次。 为了极小化圆角误差问题,产生模型时避免采用大坐标值。例如,当为采矿应用建 议模型时,采用与采矿规划和设计图形具有相同的坐标值,容易参照图形和结果返回到 采矿图上,但在计算中放大了圆角误差。 UDEC 总是随变形的发展而修改坐标。如果初始坐标是非常大(即坐标范围 10000 到 10100),当较小位移增量叠加到坐标上,可能被失去精度。而且,包括坐标值间的 差异块体搜索,也可能变得不现实的。通过改变坐标从 0 到 100,能够避免这个问题。

3.9.12 接触嵌入
如果一个块体插入另一块体太深,则可能出现“接触嵌入太大”的错误信息。程序 允许的最大嵌入量是圆角长度的一半。如果出现这样的错误,通常有必要从早期状态重 新开始这个问题。然而,在重新开始以前,重要的是识别引起错误的原因和校正方法。 涉及接触嵌入位置的有用信息在接触嵌入信息前显示。而且,采用 PLOT overlap 命令来 识别所涉及的块体。下面可能是块体嵌入的原因: 1、节理法向刚度太低 如果节理法向刚度相对于施加的荷载是不切实际的低,则块 体将插入到其他块体太深。这个原因常常通过如下识别:(a)绘制影响区域闭合图;或 (b)打印在同一区域的接触面。对于修正该问题,见上述注意 4。 2、数值不稳定性 数值不稳定性,通过增加波动幅度表征导致时间步太大。显示不 可控制的波动历史记录图是数值不稳定的标志。若不改变其他计算参数,校正数值不稳 定的唯一方式是使用 FRAC 命令减小时间步。增加体积阻尼参数能隐藏不稳定性,但不 可能消除他。UDEC 对大部分情况自动地确定一个稳定的时间步。然而,当这个时间步 太大可能出现情况。已经识别引起数值不稳定的情况如下: (1)与高频率阻尼成比例地采用大值刚度(见理论与背景的第 4.2 节); (2)采用了较高的节理剪胀角; (3)采用几何模型中有一个块体在一个边与很多(多于 3 个)其他块体接触; (4)采用了非反映(即粘结)边界,在此边界材料的刚度远远高于模型内部材料 的刚度。 如果接触嵌入的原因不能找到,有必要采用 SET cscan 命令。该命令引导围绕角点 的中心位置被快速修改,使得接触与嵌入获得较精确地计算。

-86-

消除该 误差 的另一 方法 是在计 算开 始,采 用 SET ovtol 命令,增 大嵌 入“容许 值”。如果初始圆角长度非常小,和(或)模型包含有锐角的块体,该方法是有用的。 然而,如果应用时没有认真考虑,这个方法可能是非常危险的。这是因为,如果 ovtol 太大,可能出现错误的结果。

3.9.13 非联结块体
UDEC 的数据结构被专门设计用于模拟挤压岩体(即紧密堆积块体)。该程序能够 通过“模拟接触面”处理与他们的相邻块体松散接触的离散情况,在此情况,块体之间 的连接是虚构的。然而,程序并没有考虑多面的虚拟接触。在此情况就变得无效。如果 没 有 任 何 兴趣 , 建 议 删除 这 些 不 连接 块 体 。 对于 包 含 许 多非 连 接 块 体的 分 析 , 采用 CONFIG cell 选择命令。

3.9.14 初始化变量
为帮助解释不同开挖阶段的数值模拟结果,初始化结点位移和节理面的接触位移 (RESET disp jdisp)是共同的习惯。因为程序并不需要计算过程的位移,这对结点施加 即可。为方便用户,从结点速度确定他们。对于接触面,如果节理剪胀,重新设置剪切 位移可能影响结果。接触位移用以副本储存,以致这些值能被打印和绘图而重新设置, 但并不影响模型结果。

3.9.15 确定坍塌荷载
为了确定坍塌荷载,常常采用“应变控制”边界条件而不是“应力控制”- 即,施 加常速度和监测边界反力,而不是施加力和监测位移。施加荷载的系统很难控制坍塌荷 载,这符号于一个真实系统。

3.9.16 确定安全系数
UDEC 不能直接计算安全系数。如果你需要安全系数,在给定的边界条件下,通过 产生破坏时的荷载与设计荷载之比或实际参数与破坏时参数之比来确定。例如,

FL = Fφ =

产生破坏的荷载 设计荷载

tan(实际摩擦角) tan(破坏时摩擦角)

注意,大值总是除以小值(假设系统在实际条件下并不破坏)。用户必须给出破坏 的定义。

-87-

也可以绘出变形块体单元的强度与应力之比所定义的“强度安全系数”图形。可以 采用两个强度准则:M-C 准则和 H-B 准则。根据 M-C 准则,PLOT mohr 命令可以画 出强度与应力比值的等值线。在理论与背景中方程 2.25 给出了 M-C 准则。第 3.8.1.5 节中的方程 3.26 给出了 H-B 准则。 任意单元的应力状态能够用主应力 σ 1 和 σ 3 表达。通常,该应力状态在摩尔图(图 3.50)中是作为半径为 ra 的圆“a”给出。如果该圆刚接触包络线,破坏就发生。由圆“a” 表示的应力状态的强度可通过保持 σ 3 为常数而增加或减小 σ 1 直到半径为 rb 的与包络线 相切圆“b”。两个圆半径的比值( F = rb / ra )是强度/应力比。F 也被称之为“破坏指 标”或“安全系数”。注意, F < 1 时,对于所有圆“a”的点位于包络线的外边。在

σ 3 大于抗拉极限的任何情况,F 都等于零。对于 M-C 准则,这是采用较小的抗拉强度
和粘聚力/tan(摩擦角)。对于 H-H 准则,则极限值取 sσ c / m (如方程 3.26 所定义)。

M-C 准则

H-B 准则

图 3.50 对于 M-C 和 H-B 破坏准则的强度/应力之比 方程 3.39 给出了 M-C 准则,而方程 3.40 则为 H-H 准则。方程 3.41 是强度/应力比。

σ1f = ?

1 + sin φ ? 1 + sin φ ? ?σ 3 ? 2c 1 ? sin ι ? 1 ? sin ι ?

(3.39)

σ 1 f = σ 3 ? ? mσ cσ 3 + sσ c2
rb / ra =

(3.40)

σ 3 ?σ1f σ 3 ? σ1

(3.41)

3.10 解 释

-88-

因为 UDEC 模型是一个非线性系统,在完成计算阶段,解释分析结果比传统的有限 元 程 序 所 给出 的 “解 ” 更为 困 难 。 有几 个 指 标 能够 用 作 估 计静 力 分 析 的数 值 模 型 状态 (即,系统是稳定的,不稳定的,或处于稳态塑性流)。不同指标叙述以下:

3.10.1 不平衡力
力被集中到每一刚性块体的形心和变形块体的每一结点上。在平衡状态-或变形块 体的稳态塑性流状态 - 这些力的代数和几乎是零(即,作用在形心块体一条边或结点 上的力几乎与作用在另一条边平衡)。在迭代求解过程中,对整个模型,计算最大的不 平衡力。这个力连续在屏幕上显示。也可以作为历史被保存和图形显示。不平衡力是静 力分析评价模型状态的重要指标,但其值必须与作用到模型上的内力比较。换句话说, 有必要知道“小力”是多少。采用模型中感兴趣区域的典型值,对于变形块体的典型内结 点力可以通过应力乘以与力相垂直的单元长度求得。用 R 表示最大的不平衡力与代表性 内力之比值作为百分比,R 值将永远不会降低到零。然而,1%或 0.1%可被认为是平衡 状态的可接受的值,但此值还取决于精度(即如果在报告或论文中需要最终的应力和位 移值,则可以采用 R=0.1%)。注意,一个较低的 R 值仅表示所有结点的力处于平衡, 然而,稳定塑性流可能没有加速度的出现。为了区别这个条件和真正的平衡,诸如将下 面描述的其他指标,应当考虑。

3.10.2 块体/网格结点的速度
刚性块体的速度和变形块体的结点可通过绘制速度场图(采用 PLOT vel 命令)或 选择模型中的某些关键点的历史记录(HIS xvel 或 HIS yvel)加以评估。这两种类型的 图形都是有用的。如果速度历史记录在最后阶段显示出水平曲线,就可认为达到稳态条 件。如果都已经收敛接近于零(与他们的起始值比较),然而,表明已经达到绝对平衡 状态。如果历史记录显示收敛于一个非零值,则所对应的历史的块体冒落或结点出现稳 态塑性流。如果一个或多个速度历史图显示波动速度,则该系统很可能出现瞬变条件。 速度场矢量图的解释较为困难,由于该模式的大小和特性是重要的。由于结点力、速 度决不会精确地降至零。如果计算步达到较大值(即 1000),速度的大小应与将出现的 位移有关。例如,如果系统当前位移是 1cm,则在速度图上的最大速度是 10-3/sec 和时间 步是 10-5,则 1000 步将增加 10 5m 或 10 3cm 的位移,该值是当前位移的 0.1%。在此情
- -

况下,可以说系统是平衡的,即使速度似乎是在一个方向“流动”。更为常事,其矢量在 方向和(可能)在大小的出现是随机的(或者几乎是随机的)。当结点力低于计算机精度 限度时就出现这种条件。低振幅的随机速度场是块体稳定性和无塑性流的可靠指标。 如果速度场矢量是连续的(即,具有某些系统模式)和他们的值是相当大(采用上

-89-

述描述的准则),则块体正在出现块体冒落、滑移或塑性流,或者系统仍在弹性调节 (正在发生阻尼弹性波动)。为证实正在出现连续塑性流动,应参考如下面所述的塑性 指标图。然而,如果运动包含弹性波动,则应当观测到该值的大小为表明这样的一个运 动是否有有意义的,似乎也可能发现平均波动模式。然而,如果幅度是小的,则该运动 没有物理意义。

3.10.3 块体破坏的塑性指标
对于 UDEC 中的大部分非线性块体模型,PLOT plas 命令显示应力满足屈服条件的 单元。这样的一个指标通常表明塑性流动正在出现,但块体可能是位于屈服面上,没有 发生有意义的流动。重要的是寻找塑性指标的全部模式来发现一种机理是否出现。如果 出现连接两个面的塑性接触线,就表明一个破坏机理。如果速度图也显示运动对应相同 的机理,则证实了该诊断。注意,初始塑性流动常出现在模拟的开始,但连续的应力重 新分布对屈服单元卸载以致他们的应力不再满足屈服准则(“在过去屈服”)。仅活动的 屈服单元(“正在屈服”)对识别破坏机理是重要的。如果在边界单元间存在正在屈服单 元的接触线或带,在应当 500 步计算值之前和后进行两种模式的比较。正在屈服的区域 是增加还是减少?如果是减少,则系统可能正在区域平衡。如果在增加,则到达破坏是 可能的。 如果连续塑性流动的条件已经确定,进一步的问题是:屈服流动带是否毗邻人工边界 的单元?术语“人工边界”意味着不是实际的物理边界,但为限制模型所进行的简化所形成 的边界。如果塑性流动沿着这样的边界出现,则其解是不实际的,因为破坏机理受非物理 因素的影响。此解释仅适用于最后的稳态解,中间阶段可能显现沿边界的塑性流动。 在很多问题中,存在特别感兴趣的变量(即,在一个问题中可能是位移,而在另一 个问题中就可能是应力)。应当采用 HIST 命令,以便在感兴趣的区域对所感兴趣的变 量进行追踪。在进行一些迭代计算后,这些历史提供发现系统正在做什么的一种途径。

3.11 模拟方法
3.11.1 有限数据系统模拟
在诸如地质力学领域,数据并不总是能够获得,因此,数值模拟所采用的方法应当 不同于工程力学领域。Starfield 和 Cundall(1988)为研究给出了一些建议,为进行有限 数据系统的合理模拟。在用 UDEC 进行任何数据分析之前,拜读该文是值得的。实质 上,本研究认识到场数据(例如原岩应力、材料性质和地质特征)不可能完全获得。所 以,在输入数据存在大量的不确定性的情况下,期望计算模型为设计提供数据(例如位 移)是没有用的。然而,数值模型在认识可能出现在某些特殊工程问题中的破坏机理图

-90-

形仍是有价值的。通过给出一系列“因-果”例子,模型起到培养设计工程师直感的作 用。模型可能是简单的,具有与已知场地和工程判断结果一致的数据。努力去建立很大 和复杂的模型是徒然的,对模型的理解如同对实际模型一样的困难。 当然,如果获得现场的充分数据,也可以构造可以直接产生设计信息的详细模型。 然而,更为常见的是有限的模型不能直接给出这样的信息,但提供洞察可能出现的破坏 机理的关键信息。然后设计者根据破坏机理,进行一些简单的计算,给出感兴趣的参数 判断稳定性条件。

3.11.2 混沌系统的模拟
在某些计算中,尤其在求解不连续材料,其结果可能初始条件的微小变化或加载顺 序的微小调整极为敏感。咋一看,此情况似乎不能令人满意的,以致于怀疑计算机的问 题。然而,敏感性存在于被模拟的实际物理系统。对应这种似乎是飘忽不定的特性的出 现,至少有两个方面的原因: 1、存在一定的几何不连续模型迫使系统在两个不同结果的随机选择;后来的发展 取决于实现了哪一中选择。例如,图 3.51 显示了节理岩体的一小部分。如果迫使块体 A 相对于块体 B 向下运动,它可能通过 B 的左侧,或者右侧。此选择取决于几何形状、性 质参数或动能的微观的不规则性。 2、在系统中存在称为“弱化”过程,或者更为普遍地称为正反馈。在相当均匀的应 力场,小的干扰在后来的发展中被放大,因为在一次正反馈循环,在张拉较紧的区域, 弱化较多,所以,拉力较紧等。 两 种 现 象 给 出 了 以 极 端 形 式 的 混 沌 ( Gleick,1987 和 Thompson 和 Stewart , 1986)。混沌系统的研究揭示了这样的系统的发展是不可预测的,甚至是原理。在计算 模型所观察到的对初始条件或模拟因素微小变化的敏感性,是对真实世界对微小不规则 性敏感的类似反映。没有任何意义追求任何较“精确”的计算,因为产生的模型并不代表 真实世界,在此条件是不完善的。在面对混沌系统我们的模拟策略是什么?似乎我们从 这样的一个模型中最好的期望预测特性是一有限频谱;混沌系统的统计分析是最好的定 义。我们需要建造包含初始不规则性分布的模型,即通过对节理特征采用带有随机方差 的 JSET 命令。每一模型应当采用不同的节理分布运行几次。在这些条件下,我们就可 能预测通过所施加的非常规所产生特性的波动。我们能够以统计形式表达其结果。

-91-

图 3.51 节理岩体的一小部分

3.11.3

局部化、物理的不稳定性和应力路径

在能够用 UDEC 模拟的很多系统中,可能存在其解依赖于初始条件微小变化的几条 路径。该现象被称之为“分叉”。例如,基于弹性或塑性材料的剪切试验可能变形是均匀 变形或可能显示剪切带,在此剪切应变是局部的而不是均匀分布的。似乎是如果数值模 型有足够的自由度(即,足够的单元),则可预测到这个位置。的确,关于分叉过程的 理论研究表明,即使材料没有应变软化,在剪胀角低于摩擦角的情况下也形成剪切带。 如果在求解的一个或多个局部化带存在足够的单元,简单的 M-C 材料应当始终表现出 局部化。应变软化材料是比较倾向于产生剪切带。 某些计算机程序似乎不能产生“带”的形式,尽管该现象在实际物理世界是存在的。 然而,UDEC 能够允许“带”的产生与发展,部分因为它模拟动力运动方程(即,伴随剪 切带形成的动能以实际有形的方式被释放和消失)。几篇论文报道了采用 FLAC 模拟剪 切带形式。为详细了解求解过程应当查阅这些论文。用 UDEC 没有很好处理的问题就是 剪切带的厚度。实际上,剪切带的厚度是通过如颗粒尺寸的材料内部特征确定。这一特 征在 UDEC 的本构模型中考虑。因此,剪切带从坍塌减小到根据网格重新求解的最小宽 度:如果该带平行于网格,则其值为一倍网格块度;如果带以任意角度切割网格体,其 值大约为 3 倍网格宽度。尽管剪切带形式的总体特性可用 UDEC 正确地模拟,但带的厚 度和带的间距与网格有关。进而,如果应变软化模型应用于软弱材料,对强烈受网格影

-92-

响的模拟试验,则荷载或位移关系由 UDEC 产生。这是因为集中于剪切带的应变取决于 与单元尺寸有关的剪切带的宽度(长度单位)。 涉及混沌、物理的不稳定性和分叉的主题是与路径有关的。在大部分非线性、非弹 性系统,存在无数多个满足平衡、相容和本构关系的解。对物理问题不存在正确的解, 除非指定路径。如果没有指定路径,所有可能的解都是正确的。这种情况可能引起模拟 者和用户无休止的争论,尤其求解中似乎不相关参数是否影响最终结果。所有的解在数 值上是无效的。例如,用低阻尼模拟一个采矿开挖,可能显示很大的过量和大的位移; 而高阻尼将消除振动和给出很小的最终位移。哪一个是比较符合实际?他取决于路径。 如果开挖通过边坡(即突然地),则用过量可能是合适的结果。如果开挖采用人工开挖 (即渐进地),则第二种情况可能是比较合理的。在取决于路径的情况是一种因素,模 拟应当按照系统模拟实际进展的方式加以模拟。
block 0 0 140 0 140 -99 0 -99 round 0.5 crack 0 -50 140 -50 crack 0 -60 140 -60 crack 0 -70 140 -70 crack 0 -74 140 -74 crack 0 -79 140 -79 crack 0 -99 140 -99 gen edge 5.0 bound stress -1.02e5 0 -1.02e5 ygrad 1.34e4 0 1.34e4 insit stres -1.02e5 0 -1.02e5 szz -1.02e5 ygra 1.34e4 0 1.34e4 zgra 0 1.34e4 gravity 0 -10 change cons=3 prop mat=1 d=1000 s=5.0e7 co=6.00e6 dil=4 bul=1.5e7 fric=15 change mat=1 range 0 140 -50 0 prop mat=2 d=2300 s=8.0e8 co=3.5e7 dil=9 bul=4.5e7 fric=45 change mat=2 range 0 140 -60 -50 prop mat=3 d=2500 s=1.52e9 co=5e7 dil=10 bul=5.5e7 fric=70 change mat=3 range 0 140 -70 -60 prop mat=4 d=2200 s=4.00e8 co=2.0e7 dil=7 bul=3.0e7 fric=40 change mat=4 range 0 140 -74 -70 prop mat=5 d=1600 s=7.00e7 co=1.00e7 dil=8 bul=2.7e7 fric=27 change mat=5 range 0 140 -79 -74 prop mat=6 d=2000 s=8.00e8 co=3.9e7 dil=9 bul=4.5e7 fric=40 change mat=6 range 0 140 -99 -79 change jcons=2 prop jmat=1 jkn=1e9 jks=1e9 jcoh=1e10 jtens=1e10 solve del 50 90 -79 -74 solve

材料性质 *material 1 粗砂岩-砾岩 unmark i 1 96 j 1 33 mark i 1 96j 1 4 prop d 2500 s 1.52e9 co 5e7 dil 10 bul 5.5e7 fric 70 & ctab 1 dtab 2 ftab 3 i 1 95 j 1 4 table 1 0 5.000e7 0.001 3e7 0.1 3e7

-93-

table 2 0 10 0.001 5 0.1 5 table 3 0 70.000 0.001 40 0.1 40 *material 2 砂质泥岩 unmark i 1 96 j 1 33 mark i 1 96 j 4 5 prop d 2200 s 4.00e8 co 2.0e7 dil 7 bul 3.0e7 fric 40 & ctab 4 dtab 5 ftab 6 i 1 95 j 4 5 table 4 0 2.00e7 0.001 3e7 0.1 3e7 table 5 0 7 0.001 4.0 0.1 4.0 table 6 0 40.000 0.001 16 0.1 16 *1,2,3,4,5,6 material 3 2-3 煤 unmark i 1 96 j 1 33 mark i 1 96 j 5 10 prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & ctab 7 dtab 8 ftab 9 i 1 95 j 5 11 table 7 0 1.000e7 0.001 5e6 0.1 5e6 table 8 0 8 0.001 4.0 0.1 4.0 table 9 0 27.000 0.001 14 0.1 14 *4material 3 2-3 煤 *unmark i 1 95 j 1 32 * mark i 1 95 j 5 6 *prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & * ctab 7 dtab 8 ftab 9 i 1 95 j 5 6 *table 7 0 1.000e7 0.001 5e6 0.1 5e6 *table 8 0 8 0.001 4.0 0.1 4.0 *table 9 0 27.000 0.001 14 0.1 14 *3material 3 2-3 煤 *unmark i 1 95 j 1 32 * mark i 1 95 j 5 6 *prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & * ctab 7 dtab 8 ftab 9 i 1 95 j 6 7 *table 7 0 1.000e7 0.001 5e6 0.1 5e6 *table 8 0 8 0.001 4.0 0.1 4.0 *table 9 0 27.000 0.001 14 0.1 14 *2material 3 2-3 煤 *unmark i 1 95 j 1 32 * mark i 1 95 j 6 7 *prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & * ctab 10 dtab 11 ftab 12 i 1 95 j 7 8

*table 10 0 1.000e7 0.001 5e6 0.1 5e6 *table 11 0 8 0.001 4.0 0.1 4.0 0.1 14

*table 12 0 27.000 *1material 3 2-3 煤

0.001 14

*unmark i 1 95 j 1 32 * mark i 1 95 j 7 8 *prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 &

-94-

* ctab 13 dtab 14 ftab 15 i 1 95 j 8 9 *table 13 0 1.000e7 0.001 5e6 0.1 5e6 *table 14 0 8 0.001 4.0 0.1 4.0 0.1 14

*table 15 0 27.000 *material 4 细砂岩

0.001 14

unmark i 1 96 j 1 33 mark i 1 96 j 10 14 prop d 2000 s 8.00e8 co 3.9e7 dil 9 bul 4.5e7 fric 40 & ctab 16 dtab 17 ftab 18 i 1 95 j 11 15 table 16 0 3.900e7 0.001 1e7 0.1 1e7 table 17 0 9.0 table 18 0 40.000 0.001 4.0 0.001 22 0.1 4.0 0.1 22

*material 5 2-1 煤层 unmark i 1 96 j 1 33 mark i 1 96 j 14 16 prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & ctab 19 dtab 20 ftab 21 i 1 95 j 15 17 table 19 0 1.00e7 0.001 5e6 0.1 5e6 table 20 0 8 0.001 4.0 0.001 14 0.1 4.0 0.1 14

table 21 0 27.00 *material 6 泥岩

unmark i 1 96 j 1 33 mark i 1 96 j 16 19 prop d 1600 s 7e7 co 1.20e7 dil 8 bul 2.7e7 fric 25 & ctab 22 dta 23 ftab 24 i 1 95 j 17 20 table 22 0 1.20e7 0.001 5e6 0.1 5e6 table 23 0 8 0.001 4.0 0.001 18 0.1 4.0 0.1 18

table 24 0 25.00

*material 7 1-1 1-2 煤 unmark i 1 96 j 1 33 mark i 1 96 j 19 21 prop d 1600 s 7.00e7 co 1.000e7 dil 8 bul 2.7e7 fric 27 & ctab 25 dta 26 ftab 27 i 1 95 j 20 22 table 25 0 1.00e7 0.001 5e6 0.1 5e6

-95-

table 26 0 8

0.001 4.0 0.001 13

0.1 4.0 0.1 13

table 27 0 27.00

*material 8 砂质泥岩 unmark i 1 96 j 1 33 mark i 1 96 j 21 27 prop d 2200 s 4.00e8 co 2.0e7 dil 7 bul 3.0e7 fric 35 & ctab 28 dtab 29 ftab 30 i 1 95 j 22 28 table 28 0 2.00e6 0.001 1e6 0.1 1e6 table 29 0 7 0.001 4.0 0.1 4.0 0.1 22

table 30 0 35.000 *material 9 砂岩

0.001 22

unmark i 1 96 j 1 33 mark i 1 96 j 27 28 prop d 2300 s 8.00e8 co 3.5e7 dil 9 bul 4.5e7 fric 45 & ctab 31 dtab 32 ftab 33 i 1 95 j 28 29 table 31 0 3.5000e7 0.001 10e6 0.1 10e6 table 32 0 9 0.001 4.0 0.1 4.0 0.1 22

table 33 0 45.000

0.001 22

*material 10 黄土层 unmark i 1 96 j 1 33 mark i 1 96 j 28 33 prop d 1000 s 5.0e7 co 6.00e6 dil 4 bul 1.5e7 fric 15 & ctab 34 dtab 35 ftab 36 i 1 95 j 29 34 table 34 0 6.00e6 0.001 3e6 0.1 3e6 table 35 0 4 0.001 2.0 0.001 7 0.1 2.0 0.1 7

table 36 0 15.000

Tunnel Support Loading 3 - 7 ; tunnel support loading ;block 0 -30 60 -30 60 -90 0 -90 round 0.1 crack 0 -70 60 -70 crack 30 0 30 -90 crack 42 0 42 -90 tun 30 -70 4.11 12

-96-

tun 42 -70 2.62 8 tun 30 -70 5.5 12 tun 42 -70 5.5 12 gen edge 2.0 ; initial stress state bound stress 1.02e5 0 1.02e5 ygrad 1.34e4 0 1.34e4 insit stres 1.02e5 0 1.02e5 szz 1.02e5 ygra 1.34e4 0 1.34e4 zgra 0 1.34e4 gravity 0 -10 ; rock properties prop mat=1 d=1340 zone model mohr zone shear=.33e9 bulk=.99e9 coh=1e6 fric=30.0 ;prop mat=1 d=1340 g=.33e9 k=.99e9 coh=1e6 fric=30.0 ;change con=3 ; elastic joint properties prop jmat=1 jkn=1e9 jks=1e9 jcoh=1e10 jtens=1e10 ; cycle to initial equilibrium hist solve type 1 solve rat 1e-5 save tun1.sav ; excavate service tunnel del 40 44 -72 -68 bound -1 1 -91 0 xvel=0.0 bound -1 90 -91 -89 yvel=0.0 bound 59 61 -91 0 xvel=0.0 ; histories around tunnel 1 hist ydis 42 -67 sxx 42 -67 hist ydis 42 -73 sxx 42 -73 hist xdis 39 -70 syy 39 -70 hist xdis 45 -70 syy 39 -70 ; histories around tunnel 2 hist ydis 30 -65.0 sxx 30 -65.0 hist ydis 30 -75.0 sxx 30 -75.0 hist xdis 25.0 -70 syy 25.0 -70

-97-

hist xdis 35.0 -70 syy 35.0 -70 reset disp jdisp solve rat 1e-5 save tun2.sav ; line service tunnel struct gen xc=42 yc=-70 npoint=16 mat=5 thick=0.37 fang=-11.25 theta 360 prop mat=5 st_d=2400 st_ymod=24.0e9 st_prat=0.20 st_yield=1e10 prop mat=5 if_kn=1e8 if_ks=1e7 if_coh=1.0e10 ; excavate main tunnel del 28 32 -72 -68 reset disp solve rat 1e-5 save tun3.sav ; line main tunnel struc gen xc=30 yc=-70 npoint=8 mat=5 thick=0.46 fang=22.5 theta -360 ; add additional load representing raised water level bound stress 0.0 0.0 -1.0e6 range -1 91 -31 -29 ; add hydrostatic loads to tunnel liners struct apply press 0.0 2.06e6 reset disp solve rat 1e-5 save tun4.sav ret 开挖程序 round 0.1 block 0 0 0 100 150 100 150 0 crack 0 30 150 30 crack 0 40 150 40 crack 0 60 150 60 crack 70 30 70 33 crack 70 33 74 33 crack 74 33 74 30 crack 74 30 70 30 change cons=3

-98-

jregion id=1 0 0 0 30 150 30 150 0 jregion id=2 0 30 0 40 150 40 150 30 jregion id=3 0 40 0 60 150 60 150 40 jregion id=4 0 60 0 100 150 100 150 60 jset 0,0 150,0 0,0 6,0 jset 90,0 6,0 6,0 15,0 (15 0) jset 90,0 6,0 6,0 15,0 (7.5 6) block being split jset 0,0 150,0 0,0 2,0 jset 90,0 jset 90,0 2,0 2,0 4,0 (4 30) 2,0 2,0 4,0 (2 32) range jreg 2 range jreg 2 range jreg 2 range jreg 3 range jreg 1 range jreg 1 range jreg 1

jset 0,0 150,0 0,0 5,0

jset 90,0 5,0 5,0 10,0 (10 40) range jreg 3 jset 90,0 5,0 5,0 10,0 ( 5 45) range jreg 3 jset 0,0 150,0 0,0 10,0 range jreg 4

jset 90,0 10,0 10,0 15,0 ( 15 60) range jreg 4 jset 90,0 10,0 10,0 15,0 (7.5 70) range jreg 4 change change change change mat=1 mat=2 mat=3 mat=4 range 0 150 0 30

range 0 150 30 40 range 0 150 40 60 range 0 150 60 100 0 30

generate edge 5 range 0 150

generate edge 2 range 0 150 30 40 generate edge 4 range 0 150 40 60 generate edge 10 range 0 150 60 100 delete range 70 74 30 33 bound stress -10 0 -10 insitu stress -10e6 0 -10e6 szz-10e6 ; rock blocks ; joints (no cohesion) prop m=1 dens=2500 k=5e9 g=3e9 prop jm=1 jkn=10e10 jks=10e10 jfric=30 prop m=2 dens=1000 k=5e8 g=3e8 prop jm=2 jkn=10e9 jks=10e9 jfric=20 prop m=3 dens=1500 k=1e9 g=1e9

-99-

prop jm=3 jkn=10e10 jks=10e10 jfric=25 prop m=4 dens=2500 k=5e9 g=3e9 prop jm=4 jkn=10e10 jks=10e10 jfric=32 change jcons=5 set jcondf=5 jks=40000 jfric=35 jrfric=10 & &

prop jmat jkn=40000 jcoh=2.0 jrescoh=0 jdil=6 zdil=4e-4

bound xvel 0 bound xvel 0 bound yvel 0

range -0.01 0.01

-0.01 100.01

range 149.99 150.01 -0.01 100.01 range -0.

推荐相关:

udec中文说明.pdf

udec中文说明_计算机软件及应用_IT/计算机_专业资料 暂无评价|0人阅读|0次下载 | 举报文档 udec中文说明_计算机软件及应用_IT/计算机_专业资料。对UDEC4.0的操作...

UDEC中文指导说明.doc

UDEC中文指导说明 - 通用离散元用户指导 (U D E C 3.1) 目 1

BYSOFT 中文说明.doc

控制中文说明 1页 1下载券 udec中文说明 106页 1下载券 说明表一(中

udec中文说明(已打印).pdf

udec中文说明(已打印)_电力/水利_工程科技_专业资料 暂无评价|0人阅读|0次下载 | 举报文档 udec中文说明(已打印)_电力/水利_工程科技_专业资料。udec 模拟...

UDEC4.0使用说明_图文.pdf

UDEC4.0使用说明 - 菜单驱动模式运行离散元 1、菜单驱动模式运行离散元

UDEC3.1中文说明.doc

UDEC3.1中文说明_电力/水利_工程科技_专业资料 暂无评价|0人阅读|0次下载|举报文档 UDEC3.1中文说明_电力/水利_工程科技_专业资料。udec相关...

DETEX230D中文说明.doc

更新说明(中文) 1页 免费 控制中文说明 1页 1下载券 udec中文说明 1

UDEC3.0中文手册33-63.pdf

UDEC3.0中文手册33-63 - 理开始,然后,如果有必要再逐渐增加节理来达

UDEC3.0中文手册64-101.pdf

UDEC3.0中文手册64-101 - 3.7 选择本构模型 本节概述 UDEC 所提供的块体和节理本构模型,并对模型的应用提出一些建议。在 理论与背景中的第 2 节中,为块体...

ITASCA文档资料集_图文.pdf

14.ITASCA_FLAC中文说明.rar 15.ITASCA_cad导入cfd.rar 16.ITASCA_UDEC一个例子.zip 17.ITASCA_汶川地震波.rar 18.ITASCA_Fishtank函数库.rar 19.ITASCA_3DEC...

3DEC用户手册(中文版).pdf

3DEC用户手册(中文版) - 译者在使用该软件的过程中,利用业余时间翻译说明,其中 有很多不妥的地方,请谅解。该说明为 3dec4.0 说明,4.1 和 4.0 变化不大, 4.1 ...

网站首页 | 网站地图
All rights reserved Powered by 学霸学习网 www.tceic.com
copyright ©right 2010-2021。
文档资料库内容来自网络,如有侵犯请联系客服。zhit325@126.com