当前位置:首页 >> 能源/化工 >>

FLAC第一讲(很好)


第一讲 FLAC 技术的基本原理和应用范围 1、FLAC 基本简介与本构关系
1.1 FLAC 程序简介

FLAC(Fast Lagrangian Analysis of Continua,连续介质快速拉格朗日分析)是 由 Cundall 和美国 ITASCA 公司开发出的有限差分数值计算程序,主要适用地质 和岩土工程的力学分析。该程序自 1986

年问世后,经不断改版,已经日趋完善。 前国际岩石力学学会主席 C. Fairhurst 评价它: “现在它是国际上广泛应用的可靠 程序”(1994)。 根据计算对象的形状用单元和区域构成相应的网格。 每个单元在外载和边界 约束条件下, 按照约定的线性或非线性应力—应变关系产生力学响应,特别适合 分析材料达到屈服极限后产生的塑性流动。 。由于 FLAC 程序主要是为岩土工程 应用而开发的岩石力学计算程序, 程序中包括了反映岩土材料力学效应的特殊计 算功能,可解算岩土类材料的高度非线性(包括应变硬化/软化)、不可逆剪切破坏 和压密、粘弹(蠕变)、孔隙介质的固—流耦合、热—力耦合以及动力学行为等, 另外,程序设有界面单元,可以模拟断层、节理和摩擦边界的滑动、张开和 闭合行为。支护结构,如砌衬、锚杆、可缩性支架或板壳等与围岩的相互作用也 可以在 FLAC 中进行模拟。此外,程序允许输入多种材料类型,亦可在计算过程 中改变某个局部的材料参数, 增强了程序使用的灵活性,极大地方便了在计算上 的处理。同时,用户可根据需要在 FLAC 中创建自己的本构模型,进行各种特殊 修正和补充。 FLAC 程序建立在拉格朗日算法基础上,特别适合模拟大变形和扭曲。FLAC 采 用显式算法来获得模型全部运动方程(包括内变量)的时间步长解,从而可以追踪 材料的渐进破坏和垮落,这对研究工程地质问题非常重要。FLAC 程序具有强大 的后处理功能, 用户可以直接在屏幕上绘制或以文件形式创建和输出打印多种形 式的图形。 使用者还可根据需要,将若干个变量合并在同一副图形中进行研究分 析。
1.2 本构模型

FLAC程序中提供了由空模型、弹性模型和塑性模型组成的十种基本的本构 关系模型, 所有模型都能通过相同的迭代数值计算格式得到解决:给定前一步的 应力条件和当前步的整体应变增量,能够计算出对应的应变增量和新的应力条 件。注意,所有的模型都是在有效应力的基础上进行计算的,在本构关系调入程 序之前, 将孔隙压力把整体应力转化成有效应力。下面将简要的介绍各种模型的 理论基础。 1、空单元模型

空单元用来描述被剥落或开挖的材料,其中应力为0,这些单元上没有质量 力(重力)的作用。在模拟过程中,空单元可以在任何阶段转化成具有不同材料特 性的单元,例如开挖后回填。 2、各向同性弹性模型 弹性模型的特点是卸载时材料的变形可逆, 应力—应变之间保持线性的变化 规律,和加载路径无关。 各向同性弹性模型对材料的力学行为提供了最简单的描述。 它适用于只产生 线性应力—应变关系, 在卸载条件下没有产生滞后效应的均质材料、各向同性材 料和连续介质。 在平面应变条件下,模型的应力—应变增量表达式由胡克定律得:
?? 11 ? ?1?e11 ? ? 2 ?e22 ?? 22 ? ? 2 ?e11 ? ?1?e22 ?? 12 ? 2G?e12

(6-44) ( ?? 12 ? ?? 21 )

?? 33 ? ? 2 (?e11 ? ?e22 )

式中: ?1 ? K ? ?4 3?G ;

? 2 ? K ? ?2 3?G ;
K——体积模量; G——剪切模量。
? ? ? ? 1 ?? ui ? u j ? ?eij ? ? ?t 2 ? ?x j ?xi ? ? ?

(6-45)

式中: ?eij ——应变张量增量;
u i ——位移量;
?

?t ——时步。

在平面应力条件下,模型的应力—应变增量表达式为:
?? 11 ? ?1?e11 ? ? 2 ?e22
?? 22 ? ? 2 ?e11 ? ?1?e22 ?? 12 ? 2G?e12

(6-46) ( ?? 12 ? ?? 21 )

式中: ?1 ? ?1 ? ??

? 2 ? ? 2 ? ??

2 2 2 2

?1 ?;

?? 33 ? 0

?1 ? 。

3、横观各向同性弹性模型 横观同性模型能够模拟层状的弹性介质, 模型中垂直于层状介质的弹性模量 和平行于介质的弹性模量明显不同。弹性模型的各向同性平面位于X-Z平面内,

弹性模量的相关定义如下: E1 ( Ex )——各向同性平面内的弹性模量; E2 ( Ey ) ——垂直于各向同性平面内的弹性模量; G12 ( Gxy) ——X-Y平面或者Y-Z平面内的剪切模量; G13 ( Gxz)——各向同性平面内的剪切模量;

? 21 (? yx )——由单轴Y方向上的应力引起的各向同性平面内X方向上的法向
应力和垂直于各向同性平面上的Y方向上的法向应力组成平面内的泊松比;

? 31 (? zx )——由于单轴Z方向上的应力引起的各向同性平面内X方向上的法
向应力和垂直于各向同性平面上的Z方向上的法向应力组成平面内的泊松比。 横观各向同性弹性材料的特征是包含五个独立的常量。 各弹性模量参数之间 的具体关系如下:
E 3 ? E1

(或 E z ? E x ) (或? zx ? ? xz ) (或? yz ? ? yx ) (或 G yz ? G xy ) (或 G xz ?
Ex ) 2?1 ? ? zx ?
Ex ) Ey

? 31 ? ? 13
? 23 ? ? 21
G23 ? G12

G13 ?

E1 2?1 ? ? 31 ? E1 E2

? 12 ? ? 21

(或? xy ? ? yx

图6-6 横观各向同性坐标系定义规则(X-Z方向为各向同性面)

模型中关于材料的弹性模量变化有一定的限制(Amadei,1982),其限制条件 为:

Ex ? 0
Ey ? 0

G xy ? 0
2 ? xy ? 1
2 ? xz ? 1

(6-47)

?1 ? ? xz ? ?

2 2 E x? yx

Ey

?0

4、Drucker-Prager(德鲁克-布拉格)塑性模型 Drucker-Prager塑性模型适用于低摩擦角的软质粘土。 (1)增量弹性定律 Drucker-Prager模型根据两种广义应力(剪切应力τ和平均法向应力σ)定义如 下:

? ? J2

??

1 ?? 11 ? ? 22 ? ? 33 ? 3

(6-48)

公式中J2是应力偏量第二不变量,可以表示为:
J2 ? 1 2 ?? 11 ? ? 22 ?2 ? ?? 22 ? ? 33 ?2 ? ?? 11 ? ? 33 ?2 ? ? 12 6
' ?? ? 2 ?J 2

?

?

(6-49)

剪切应变增量 ?? 、体积应变增量 ?e 和τ以及σ的关系如下: (6-50)

?e ? ?e11 ? ?e22 ? ?e33
' 公式中, ?J 2 为应变偏增量第二不变量:

' ?J 2 ?

1 2 ??e11 ? ?e22 ?2 ? ??e22 ? ?e33 ?2 ? ??e11 ? ?e33 ?2 ? ?e12 6
?? ? ?? e ? ??
p

?

?

(6-51)

应变增量可以分解如下: (6-52)
?e ? ?e e ? ?e p

公式中上标e和p分别指弹性部分和塑性部分,在弹性变形阶段,塑性变形不为0, 胡克定律的增量表达式根据广义应力和应变表述如下:
?? ? G?? e

(6-53)

?? ? K?e e

公式中,K和G分别为剪切模量和体积模量。 (2)屈服函数

图6-7描述了(σ,τ)平面内的破坏判断准则,由Drucker-Prager模型的屈服函数 可以得到点A到点B的破坏包络线:
f s ? ? ? q? ? ? k?

(6-54)

B点到C点的拉破坏屈服函数如下:
f t ? ? ?? t

(6-55)

公式中, q? 和 k ? 是材料力学性质常量; ? t 是Drucker-Prager模型中的抗拉强度, 被定义为研究材料平均法向应力的最大值。对于 q? ? 0 的材料,抗拉强度不能超 过最大值 ? max t :
t ? max ?

k? q?

(6-56)

图6-7

Drucker-Prager模型破坏准则

(3)材料力学参数的注解 Drucker-Prager模型的剪切破坏准则 f s ? 0 在主应力空间 ?? 1 , ? 2 , ? 3 ? 中, 通过 以 ? 1 ? ? 2 ? ? 3 为轴, ?? 1 , ? 2 , ? 3 ? ? ?a, a, a ? ( a ? k? q? ,图6-8)为顶点的圆锥体进 行了说明。莫尔-库伦剪切破坏准则中有两个特征参数:粘聚力c和内摩擦角 ? , 可以通过一个同轴的不规则六棱锥进行描述,六棱锥有三个向外的边缘,有三个 向内的边缘(图6-9)。参数 q? 和 k ? 根据Drucker-Prage圆锥体通过莫尔-库伦六棱锥 的外部边缘或内部边缘而变化,对于通过六棱锥的外部边缘有:
q? ? 3 ?3 ? sin ? ? 6 6 sin ? c cos?

(6-57)

k? ?

3 ?3 ? sin ? ?

对于通过六棱锥的内部边缘有:

q? ?

3 ?3 ? sin ? ? 6

6

sin ? c cos?

(6-58)

k? ?

3 ?3 ? sin ? ?

在 q? ? 0 的特殊情况下,Drucker-Prager屈服准则可简化为Mises屈服准则, Mises屈服准则和主应力空间中的圆柱体相对应。Tresca屈服准则是 Mohr-Coulomb准则中 ? ? 0 的一种特殊情况,在主应力空间中可以通过规则的六 棱柱进行描述。Mises圆柱已包容了棱柱,因此有:
q? ? 0

(6-59)
c

k? ?

2 3

图6-8 主应力空间中Drucker-Prage屈服准则和Mises屈服面

图6-9

主应力空间中Mohr-Coulomb屈服准则和Tresca屈服面

5、Mohr-Coulomb(莫尔-库伦)塑性模型

Mohr-Coulomb模型通常用于描述土体和岩石的剪切破坏。模型的破坏包络 线和Mohr-Coulomb强度准则(剪切屈服函数)以及拉破坏准则(拉屈服函数)相对 应。 (1)增量弹性定律 FLAC程序运行Mohr-Coulomb模型的过程中,用到了主应力 ? 1 、 ? 2 和 ? 3 , 以及平面外应力 ? zz 。主应力和主应力的方向可以通过应力张量分量得出,且排 序如下(压应力为负):

?1 ? ? 2 ? ? 3
对应的主应变增量 ?e1 、 ?e2 和 ?e3 分解如下:
?ei ? ?eie ? ?eip i ? 1,3

(6-60)

(6-61)

公式中,上标e和p分别指代弹性部分和塑性部分,且在弹性变形阶段,塑性应变 不为零。根据主应力和主应变,胡克定律的增量表达式如下:
e e ?? 1 ? ? 1 ?e1e ? ? 2 ?e2 ? ?e3 e e ?? 2 ? ? 1 ?e2 ? ? 2 ?e1e ? ?e3

?

?

e e ?? 3 ? ? 1 ?e3 ? ? 2 ?e1e ? ?e2

?

?

?

(6-62)

?

公式中, ?1 ? K ? 4G 3 ; ? 2 ? K ? 2G 3 。 (2)屈服函数 根据公式(6-60)的排序,破坏准则在平面 ?? 1 , ? 3 ? 中进行了描述,如图6-10。 由Mohr-Coulomb屈服函数可以得到点A到点B的破坏包络线为:
f s ? ? 1 ? ? 3 N ? ? 2c N ?

(6-63)

B点到C点的拉破坏函数如下:
f t ? ? t ??3

(6-64)

公式中, ? ——内摩擦角; c——粘聚力;

? t ——抗拉强度。
N? ? 1 ? sin ? 1 ? sin ?

(6-65)

注意到在剪切屈服函数中只有最大主应力和最小主应力起作用; 中间主应力
t 不起作用。对于内摩擦角 ? ? 0 的材料,它的抗拉强度不能超过 ? max ,公式如下:

t ? max ?

c tan?

(6-66)

图6-10 Mohr-Coulomb强度准则

6、统一节理塑性模型 此模型是一种各向异性塑性模型,用于模拟符合Mohr-Coulomb破坏准则的 材料中含有特定方位节理面的岩体。 统一节理模型中可能骨架发生破坏,也可能节理面发生破坏,或者两者都发 生破坏,这主要依赖于应力条件、节理面的方向以及骨架和节理面的力学参数。
' 理面的塑性变形。图6-11给出了 ?? 22 ,? ? 平面内,节理面的破坏判断准则:

FLAC程序中关于Mohr-Coulomb模型在上面作了介绍, 本部分将主要介绍节 A点到B点的局部破坏包络线由Mohr-Coulomb破坏准则定义为 f s ? 0 和:
' f s ? ?? ? ? 22 tan? j ? c j

(6-67)

B点到C点的拉破坏判断准则定义为 f t ? 0 和:
' f t ? ? tj ? ? 22

(6-68)

公式中,? j 、cj和 ? tj 分别为节理面的内摩擦角、粘聚力和抗拉强度。注意,对于 摩擦角不为零的节理面,其抗拉强度的最大值为:
cj tan ? j

? tj ,max ?

(6-69)

图6-11

统一节理模型中节理面的破坏准则

7、应变硬化/软化塑性模型 应变硬化/应变软化模型用于模拟非线性材料的软化和硬化力学行为,模型 的实现基于Mohr-Coulomb模型中力学参数(粘聚力、内摩擦角、膨胀角、抗拉强 度)的变化,这些力学参数都是塑性偏应变的函数。 此模型依赖于上面提到的含有剪切判断准则和拉破坏判断准则的 Mohr-Coulomb模型。它们的不同点在于应变硬化/软化模型在塑性屈服开始时, 粘聚力、内摩擦角、膨胀角和抗拉强度具有变硬或者变软的可能性,而在 Mohr-Coulomb模型中,这些力学参数假设是恒定的。在这种模型中,使用者能 够把粘聚力、内摩擦角以及膨胀角定义为测量塑性剪切应变参数的分段线性函 数, 抗拉强度的分段线性软化规律可以根据另外一个测量塑性拉应变的硬化参数 定义。 模型的屈服函数、塑性变形法则和Mohr-Coulomb模型中相同。 (1)硬化/软化参数 塑性剪切应变由剪切硬化参数eps测量,其增量形式定义如下:
1

?e Ps
公式中:

?1 Ps ? ? ?e1Ps ? ?s m ?2

?

?

2

1 Ps ? ?em 2

?

?

2

2 1 Ps Ps 2 ? ? ?e3 ? ?em ? 2 ?

?

?

(6-70)

Ps ?em ?

1 Ps ?e1Ps ? ?e3 3

?

?

?e Ps , j ? 1,3 是主塑性剪切应变增量。 j

拉硬化参数ept测量了累积的拉塑性应变,其增量形式定义如下:
Pt ?e Pt ? ?e3

(6-71)

公式中, ?e3Pt 是主应力方向上的拉应变增量,(拉应力为正)。 附加的上标s指和剪切屈服面相关的塑性应变(而不是拉屈服面)。 ?eiPs 是塑

形主应变增量,而不是剪切应变增量。 ?e3Pt 的上标t指和拉屈服面相关的塑性应 变。 (2)自定义粘聚力、内摩擦角、膨胀角和抗拉强度函数 一维的应力-应变关系曲线σ-e,在屈服阶段发生软化,并且保留一定的残余 强度:

屈服点

图6-12

应力-应变曲线

在到达屈服点以前,曲线保持线性变化,在这个阶段,仅仅存在弹性应变, 即 e ? e e 。超过屈服点以后,全应变由弹性应变和塑性应变两部分组成,即:
e ? e e ? e p 。在应变软化/硬化模型中,使用者把粘聚力、内摩擦角、膨胀角和抗

拉强度的变化定义为全应变中的塑性应变部分 e p 的函数,这些函数的例子见图 6-13,在FLAC程序中调整为线性变化的线段(图6-14)。

图6-13

和塑性应变对应的粘聚力(a)和内摩擦角(b)的变化曲线

图6-14

和线段近似的曲线

8、双线性应变硬化/软化一致节理模型 双线性应变硬化/软化一致节理模型能够模拟地质材料以及节理面的软化和 硬化行为,模型通过一致节理模型力学性质(内聚力、内摩擦角、膨胀率、抗拉 强度)的变化来实现,这些力学性质参数是偏塑性应变和拉塑性应变的函数。双 线性模型还能模拟材料强度随平均应力的变化规律。 双线性应变硬化/软化一致节理模型是上面描述的统一节理模型的一种普遍 化,在模型中,骨架以及节理的破坏包络线是两个Mohr-Coulomb强度准则和根 据特定的规律可以硬化或软化的拉破坏判断准则组成的组合。 骨架和节理的软化行为通过四个独立的硬化参数确定, 其中两个参数用于描 述骨架,另外两个参数用于描述节理,它们分别对应着塑性剪切应变和拉应变。 在此模型的计算过程中,如果塑性变形已经发生,硬化参数增加,粘聚力、内摩 擦角、膨胀角和抗拉强度的参数通过自定义数据对骨架和节理进行参数调整。 (1)骨架的失效准则和变形规律 模型中关于骨架的破坏准则在图6-15的主应力平面(σ1,σ3)内进行了描述(压 力为负, ? 1 ? ? 2 ? ? 3 )。 且 破坏包络线含有A-B段和B-C段内的两个Mohr-Coulomb 破坏判断准则: f 2s ? 0 和 f1s ? 0 ;以及C-D段的拉破坏判断准则: f t ? 0 。 剪切破坏判断准则的普通形式为: f s ? 0 ,在于A-B段内,粘聚力 c ? c2 、 内摩擦角 ? ? ? 2 ;在B-C段内,粘聚力 c ? c1 、内摩擦角 ? ? ?1 。拉破坏准则通过 抗拉强度 ? t (正值)进行描述:

图6-15

双线性骨架破坏准则

f s ? ? 1 ? ? 3 N ? ? 2c N ?

(6-72) (6-73) (6-74)

f t ? ?3 ?? t

公式中:

N? ?

1? s i n ? 1? s i n ?

对应于 f 2s ? 0 和 f1s ? 0 交叉点处的σ3的数值如下:

? 3I ?
(2)节理面的破坏准则

2c2 N ?2 ? 2c1 N ?1 N ?2 ? N ?1

(6-75)

图6-16描述了节理面的破坏准则,破坏准则和Mohr-Coulomb准则的两种判 断形式相符,在A-B段内: f 2s ? 0 ;在B-C段内: f1s ? 0 ;以及一个拉破坏判断 准则,在C-D段内: f t ? 0 。每种剪切破坏判断准则都有普通形式 f s ? 0 ,其特 点是在A-B段内, 粘聚力 c j ? c j 2 , 内摩擦角 ? j ? ? j 2 ; 在B-C段内, 粘聚力 c j ? c j1 , 内摩擦角 ? j ? ? j1 。拉破坏的判断准则通过抗拉强度 ? tj (正值)进行了说明。因此 有:

图6-16
s

节理的破坏判断准则

' f ? ? ? ? 22 tan ? j ? c j ' f t ? ? 22 ? ? tj

(6-76) (6-77)

注意,对于摩擦角 ? j1 ? 0 的节理面,抗拉强度的最大值为:

? tj max ?
(3)硬化参数

c j1 tan ? j1

(6-78)

双线性应变硬化/软化一致节理模型发生塑性变形后,一些区域或者所有区 域内骨架以及节理的参数(粘聚力、内摩擦角、膨胀角和抗拉强度)根据硬化函数 对开始的输入值,进行分段线性调整。 模型中用到了四个独立的硬化参数: ① ? s ,对应骨架的塑性剪切应变,用于更新骨架的粘聚力、内摩擦角和膨 胀角; ② ? t ,对应骨架的塑性体积拉应变,用于更新骨架的抗拉强度; ③ ? sj ,对应节理的塑性剪切应变,控制着节理的粘聚力、内摩擦角和膨胀 角的更新; ④ ? tj ,对应节理的的塑性体积拉应变,用于更新节理的抗拉强度。 9、双屈服模型 双屈服模型能够模拟一些材料的力学行为,这些材料除了发生剪切屈服外, 还会发生一种和剪切屈服明显相互作用的屈服模式, 例如液压回填材料或者由混 凝土胶结形成的小颗粒材料。 模型中考虑了由于各向同性压力引起的永久性体积变化,还考虑了FLAC程 序中应变软化/硬化模型中剪切破坏和拉破坏判断准则以及体积屈服面(帽盖)。 为 了简便起见,体积屈服面由体积压力来定义,即 pc ? 0 ,和剪切应力无关;它由 剪切应力和平均应力图上的一条垂线组成。 体积屈服面的硬化行为被体积塑性应 变激活,并且符合分段线性规则。当根据因子R定义的特殊规律发生塑性体积应

变时,切线体积模量和剪切体积模量得到了发展,其中R假设为连续的,并且定 义为弹性体积模量和塑性体积模量的比值。 除了和应变软化模型相关的参数外, 此模型还需要另外两个参数和一个切平 面: ① pc的初始值,它等于材料曾经经历的最大平均压力值; ② R值,大于1,在体积应力卸载时,控制了应力-应变曲线的斜率; ③ 硬化曲线的切平面表达式,把体积屈服面的压力pc和塑性体积应变epv联 系在一起。 (1)弹性定律增量形式 此模型中,用到了主应力 ? 1 、 ? 2 和 ? 3 以及平面外应力 ? zz ,这些主应力以 及主应力的方向通过应力张量分量计算得出,且有以下排序(压力为负):

?1 ? ? 2 ? ? 3
对应的主应变增量形式分解如下:
?ei ? ?eie ? ?eip i ? 1,3

(6-79)

(6-80)

公式中上标e和p分别指代弹性部分和塑性部分,在弹性屈服阶段,塑性应变不是 零(向外的应变为正), 假设剪切屈服、 拉屈服以及体积屈服的塑性和都是附加的, 则有:
?eip ? ?eips ? ?eipt ? ?eipv

(6-81)

公式中,上标ps、pt和pv分别代表塑性剪切应变、塑性拉应变和塑性体积应变。 根据定义,这部分中的符号 ?e 用于指负的体积应变增量( ?e1 ? ?e2 ? ?e3 )、塑性 应 变 增 量 ?e p 和 弹 性 应 变 增 量 ?e e 。 符 号 ?e pv 指 负 的 塑 性 体 积 应 变 增 量 ( ?e1pv ? ?e2pv ? ?e3pv )。 根据主应力以及主应变的胡克定律表达式可以得到增量表达式的形式:
e e ?? 1 ? ? 1 ?e1e ? ? 2 ?e2 ? ?e3 e e ?? 2 ? ? 1 ?e2 ? ? 2 ?e1e ? ?e3

?

?

e e ?? 3 ? ? 1 ?e3 ? ? 2 ?e1e ? ?e2

?

?

?

(6-82)

?

公式中, ? 1 ? K c ? 4Gc 3 ; ? 2 ? K c ? 2Gc 3 ,Kc、Gc是根据下面的考虑定义的 切向体积模量和剪切模量。 它的塑性刚度 ?dpc de pv ? 通常是增加的,既然物体的颗粒距离更近,弹性刚度增 的弹性刚度Kc是常数因子R和目前逐渐增加的塑性刚度的乘积。体积模量K和剪 切模量G由使用者提供,被认为是Kc和Gc的上限,并且假设 K c / Gc ? K / G 的比 率为常数。使用增量符号,定义如下: 考虑用逐渐增加的压力pc对物体进行各向同性压缩,当材料变得更密实时,

加似乎也很合理。 在一般的加载条件下模型采用了一个非常简单的准则,即增加

Kc ? R Gc ? G Kc K

?pc ?e pv

? K c: m i nK c , K ? ?
(6-83)

公式中,比率R已知,

?p c 是pc切平面的当前斜率。 ?e pv

图6-17展示了双屈服破坏模型的力学行为, 图中给出了非线性体积加载曲线 和几条卸载曲线。卸载曲线是弹性的,其斜率和卸载点的塑性刚度以及R相关。 (2)屈服函数 剪切屈服函数 f s 和拉屈服函数 f t 有以下的形式:
f s ? ? 1 ? ? 3 N ? ? 2c N ?

(6-84) (6-85)

f t ? ? t ??3

公式中, N ? ? ?1 ? sin ? ? ?1 ? sin ? ? ;

? ——内摩擦角;
c——粘聚力;

? t ——抗拉强度。
体积屈服函数 f v 定义如下:
fv ? 1 ?? 1 ? ? 2 ? ? 3 ? ? pc 3

(6-86)

公式中pc为体积屈服面压力。

径 路 载 加

卸 载 路 径

压力

体积应变

图6-17 弹性体积加载/卸载路径

10、改进的剑桥(Cam-clay)模型 需要考虑变形模量引起材料的体积变化以及对剪切破坏的阻碍时, 例如软质 粘土,可以采用改进的剑桥模型。该模型是一种逐渐增加的硬化/软化弹塑性模

型。它的特点包括特殊格式的非线性弹性变化以及由体积塑性应变(由密度决定) 决定的硬化/软化行为。破坏包络线形状上自我相似,并且在主应力空间中和围 绕着平均应力轴旋转的椭圆体相对应。(以上所有模型都是根据有效应力进行计 算的,这部分的压力也指有效压力。) (1)弹性定律的增量形式 改进的剑桥模型用三个变量进行描述:平均有效压力p,偏应力q和特定的体 积v。在FLAC程序运行此模型的过程中,使用到了主应力 ? 1 、 ? 2 、? 3 以及平面 外应力 ? zz (牵引力和膨胀角为正)。 广义应力分量p和q根据主应力表达式如下:
p??
q? 1 2

1 ?? 1 ? ? 2 ? ? 3 ? 3

?? 1 ? ? 2 ?2 ? ?? 2 ? ? 3 ?2 ? ?? 1 ? ? 3 ?2

(6-87)

(注意, q ? 3J 2 ,其中J2是有效偏应力张量第二不变量。) 和-p 、q相关的应变增量是体积应变增量 ?e 和扭曲应变增量 ?eq ,因此有:
?e ? ?e1 ? ?e2 ? ?e3

?eq ?

2 3

??e1 ? ?e2 ?2 ? ??e2 ? ?e3 ?2 ? ??e1 ? ?e3 ?2

(6-88)

公式中, ?e j , j ? 1,3 是主应变增量,假设主应变增量可以分解为弹性部分和塑性 部分,则:
?ei ? ?eie ? ?eip i ? 1,3

(6-89) (6-90)

特殊的体积v定义为:

v?

V Vs

公式中,Vs是固体颗粒体积,假设不可压缩,包含在体积为V的土体中。体积应 变e和特定体积之间的增量关系如下:
?e ? ?v v

(6-91)

假设初始特定体积为v0,对于小体积应变增量有:
v ? v0 ?1 ? e ?

(6-92)

公式中,e是目前累积的体积应变。 主坐标轴上的胡克定律增量表达形式为:
e e ?? 1 ? ? 1 ?e1e ? ? 2 ?e2 ? ?e3 e e ?? 2 ? ? 1 ?e2 ? ? 2 ?e1e ? ?e3

?

?

e e ?? 3 ? ? 1 ?e3 ? ? 2 ?e1e ? ?e2

?

?

?
(6-93)

?

公式中: ?1 ? K ? 4G 3 ; ? 2 ? K ? 2G 3 。

另外,使用应力张量、应变张量增量的偏张量部分可以得出:
?si ? 2G? ?ie i=1,3

? ?p ? K?e e

(6-94)

公式中:

?si ? ?? i ? ?p
? ?ie ? ?eie ? ?e e 3
e e ?e e ? ?e1e ? ?e2 ? ?e3

在剑桥模型中,体积关系式(6-94)中的切线体积模量K被更新,用来反映各 向同性压缩试验中的非线性规律。图6-18给出了半对数坐标系中,典型的各向同 性压缩试验的结果图。

法向固 结线

膨胀线

图6-18 各向同性压缩试验中的法向固结线及加载-卸载曲线

由于法向固结力p增加,材料的特定体积v减小。描述材料沿着垂直固结线运 动的点可以定义为:
v ? v? ? ? ln p p1

(6-95)

公式中, ?? 和 v ? 是两个材料参数,p1是参考压力(注意,在参考压力下, v ? 就是 特定体积)。 图中A-B之间的卸载-重新加载曲线将沿着斜率 ? 的弹性膨胀线移动, 回到法 向固结线,然后重新继续移动。膨胀线的方式如下:
v ? v? ? ? ln p p1

(6-96)

公式中, ? 是材料常量,特殊线上的 v? 依赖于卸载法向固结线上点的位置。 特殊体积 ?v e 可重获的变化由公式(6-96)微分后求得:
?v e ? ?? ?p p

(6-97)

用v区分所有的参数,和公式(6-91)比较,可得: (6-98) ?e e ? 在改进的剑桥模型中, 假设平均压力的任何变化对应着上面的表达式中的体 积弹性变化,和公式(6-94)对比,该模型材料的切向体积模量为:
K? vp ? ?p ? vp

?

(6-99)

(2)屈服函数 和固结压力的特殊值pc对应的屈服函数为:
f ? q 2 ? M 2 p? p ? pc ?

(6-100)

公式中M是材料常数,屈服函数 f ? 0 通过平面(q,p)中的一个椭圆、横轴pc以及 纵轴Mpc进行描述(图6-19)。注意到椭圆形穿过原点,因此这个模型中的材料不 能承受多方面的拉应力。

临界状态线 塑性膨胀角 塑性压实状态

图6-19

改进的剑桥模型破坏准则

基本应力空间中的一个以平均应力为轴转动的椭圆体描述了破坏的判断准 则(恒定平均有效应力p点通过屈服面的任何一部分都是一个圆形)。 (3)初始应力状态 FLAC程序中的改进剑桥模型仅仅用于模拟和压缩平均有效应力对应的应力 状态的材料,介质中初始有效压力p0必须为正。 (4)超固结率 超固结率R被定义为预压密压力和初始压力的比值,即:
R? pc0 p0

(6-101)

在描述改进的剑桥模型材料的力学行为时,这个比率非常有用。

2 、FLAC 建模的基本方法和程序

2.1 FLAC 程序建模方法 通过建立数值计算模型求解不同的工程地质问题。 下面给出 FLAC 程序了具 体的解题步骤以及应注意的相关问题。 1、根据实际工程规划计算模型,主要包括以下五个方面的内容: (1)设计模型尺寸:计算模型范围的选取直接关系到计算结果的正确与否, 模型范围太大,白白耗费了计算机能源,模型范围太小,计算结果失真,不能给 实际工程指导性的意见,因此合理的选择计算模型的范围至关重要; (2)规划计算网格数目和分布:计算模型的尺寸一旦确定,计算网格的数目 也相应确定,程序中所能容纳的计算网格数目和计算机的 CPU 以及内存有重要 的关系, 因此一台配置较好的计算机是非常重要的。程序中为了减少因网格划分 引起的误差,网格的长宽比应不大于 5,对于重点研究区域可以进行网格加密处 理; (3)安排工程对象(开挖、支护等):对于需要开挖或者支护的工程,应在建模 过程中进行规划,调整网格结点,安排开挖以及支护的位置等; (4)给出材料的力学参数:在建模时,应根据实际工程确定本构关系,给模 型赋以相应的力学参数,力学参数往往来源于现场或者试验; (5)确定边界条件:模型的边界条件包括位移边界和力边界两种(包括模型内 部出适应力和位移),在计算前应确定模型的边界状况。 2、做好以上的规划准备后,须在计算机上建立模型,建模过程中具体操作 步骤和常用的语句介绍如下: (1)网格生成:grid i,j 模型开始建立时, 首先应该给出模型在 X 方向上总的单元格数目 i 和模型在 Y 方向上总的单元格数目 j,例如: grid 30,20 mod e

Y

j=20

X i=30

图 6-20

模型网格图

(2)网格规划:gen x1,y1 x2,y2 x3,y3 x4,y4 模型总的网格数目给定后, 需要对模型的整体区域进行圈定,因此也就指定 了模型的尺寸,(x1,y1)、(x2,y2)、(x3,y3)和(x4,y4)分别为区域从左下角起按顺时针 旋转的四点坐标,例如: grid 30 20 mod e gen 0,0 0,20 30,20 30,0
2(x2,y2) 3(x3,y3)

1(x1,y1) 图 6-21 模型网格规划示意图

4(x4,y4)

(3)分区规划网格:考虑到实际工程中地质条件的影响以及施工过程中的需 要,程序中需对整体模型进行区块划分,命令如下: gen xI1,yI1 xI2,yI2 xI3,yI3 xI4,yI4 i=1,n1 j=1,m1 (I 区) gen xII1,yII1 xII2,yII2 xII3,yII3 xII4,yII4 i=n1,n2 j=1,m1 (II 区) 其中(xi1,yi1)、(xi2,yi2)、(xi3,yi3)和(xi4,yi4)为第 i 个区域从左下角起按顺时针旋转的 四点坐标,例如: grid 30 20 mod e gen 0,0 0,10 10,20 10,0 i=1,11 j=1,21(i、j 为区域沿 X、Y 方向的结点号) gen 10,0 10,20 30,20 30,0 i=11,31 j=1,21
i=11,j=21 i=31,j=21

i=1,j=21 Ⅱ Ⅰ

i=1,j=1

i=11,j=1 图 6-22 模型分区网格图

i=31,j=1

(4)特殊形状网格的形成:在实际工程中涉及到地质状况非常复杂,仅仅掌 握这些规则区域形成的命令语句还是不够的,为了能形象地刻画出洞室、复杂的 矿体埋藏以及地表分布状况等,还需要用到以下特殊形状网格形成的语句。 ① 圆形: circle xc,yc rad, gen 其中(xc,yc)为圆心坐标, 为圆形半径, rad 例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0 gen circle 15,10 5 (定义了以(15,10)为圆心,5 为半径的圆形)

rad

(xc,yc)

图 6-23

圆形网格图

② 弧线:gen arc xc,yc xb,yb theta,其中(xc,yc)为圆心坐标,其中(xb,yb)为弧线 起点坐标,theta 为弧线按照逆时针旋转的角度,例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0

gen arc 15,10 20,15 90 (定义了以(15,10)为圆心,(20,15)为弧线起点,逆时针 旋转 90? 的弧线)

(xb,yb) theta

(xc,yc)

图 6-24

弧线网格图

③ 线段:gen line x1,y1 x2,y2,其中(x1,y1)和(x2,y2)分别为线段的两个端点坐 标,例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0 gen line 10,5 20,15(形成了端点分别为(10,5)、(20,15)的线段)

(x2,y2)

(x1,y1)

图 6-25

线段网格图

④ 任意形状:table n x1,y1, x2,y2, ?,xm,ym, x1,y1 gen table n 其中 table n 代表圈定的第 n 个区域,(x1,y1)、(x2,y2)和(xm,ym)分别代表区域 典型点的位置, 可以根据习惯逆时针圈定或者顺时针圈定,但是圈定区域必须为 闭合区域,也就是第一个坐标和最后一个坐标应该相同,gen table n 则定义了这 个区域,例如: grid 30,20

mod e gen 0,0 0,20 30,20 30,0 table 1 5,5 8,12 15,15 25,12 22,5 12,3 5,5 gen table 1

(x3,y3)

(x2,y2)

(x4,y4)

(x1,y1) (x6,y6) 图 6-26 任意形状网格图

(x5,y5)

3、模型建立以后,在计算前需要赋给单元材料性质,下面介绍常用的两种 力学模型中需要赋值的参数。 (1)弹性模型 模型采用弹性本构关系进行计算时,单元材料只需要三个参数值,具体语句 如下: mod e (弹性) prop d=? bu=? sh=? i= i1,i2 号) 其中:d——材料的容重; bu——材料的弹性模量, K ? sh——材料的剪切模量, G ?
E 3(1 ? 2 ? ) E 2(1 ? ? ) (bu lk modulus) ; (sh ear modulus) 。

j= j1,j2(i,j 分别为沿模型 X 方向和 Y 方向的单元

注意:在赋予材料力学性质时,应注意量纲的统一,表 6-1 给出了量纲之间的关 系。
表 6-1 Length Density Force m kg/m3 N SI m 3 10 kg/m3 kN 力学参数单位对应表 m 10 kg/m3 MN
6

cm 10 g/cm3 Mdynes
6

Imperial ft slugs/ft3 Ibf

In snails/in3 Ibf

Stress Gravity Stiffness

Pa m/sec2 Pa/m

kPa m/sec2 kPa/m

MPa m/sec2 MPa/m

bar cm/s2 Bar/cm

Ibf/ft2 ft/sec2 Ibf/ft3

psi in/sec2 Ib/in3

*刚度指节理面的法向和剪切刚度。

(2)Mohr-Coulumb 模型 模型采用莫尔-库伦本构关系进行计算时,单元材料需要六个参数值,语句 如下: mod m (弹塑性 Mohr-Coulumb 准则) prop d=? bu=? sh=? c=? fri=? ten=?i=i1,i2 j=j1,j2 其中:d——材料的容重; bu——材料的弹性模量; sh——材料的剪切模量; c——材料的粘聚力; fri——材料的内摩擦角; ten——材料的抗拉强度。 4、模型边界条件 (1)位移边界(结点) 通常由于选取模型的范围比模拟对象的范围大, 因此对模型无明显影响的边 界条件可以采用限制边界位移方法处理,具体语句如下: fix x i=i1,i2 j=j1,j2 fix y i=i3,i4 j=j3,j4 例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0 gen circle 15,10 5 mod nu reg i=15 j=11 fix x i=1(限制了模型左边界 X 方向上的位移) fix y j =1(限制了模型下边界 Y 方向上的位移) (限制了 x 方向上的位移) (限制了 y 方向上的位移)

fix x

i=1

fix y 图 6-27

j=1 位移边界网格图

(2)力边界条件 (结点) 有些模型建立以后, 根据实际情况需要对模型施加外部载荷, 也就是边界力, 具体语句如下: apply yf=? i=i1,i2 j=j1,j2 或 apply syy=? i=i1,i2 j= j1,j2 (施加 Y 方向上的外力) apply xf=? i=i3,i4 j=j3,j4 或 apply sxx=? i=i3,i4 j=j3,j4 (施加 X 方向上的外力) 其中 yf(syy)为 Y 方向上的外力,方向与坐标轴相同;xf(sxx)为 X 方向上的外力, 方向与坐标轴相同。例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0 gen circle 15,10 5 mod nu reg i=15 j=11 apply yf=-10 i=1,31 j=21 apply xf=-5 i=31 j=1,21 5、初始应力场的形成 地应力始终是工程地质界的研究重点,特别是在有构造应力的条件下,水平 地应力往往大于垂直地应力, 应力环境对地质工程的受力状态起重要的作用,在 程序中正确的输入地应力成为不可缺少的一项,FLAC 程序很好的解决了地应力 输入的问题,具体语句如下: ini sxx=? var n1 n2 i=i1,i2 j=j1,j2 ini syy=? var m1 m2 i=i3 i4 j=j3 j4 其中:sxx——X 方向地应力; syy——Y 方向地应力; 或 apply syy=-10 i=1,31 j=21 或 apply sxx=-5 i=31 j=1,21

var n1 n2——地应力的变化量从 n1 逐渐过渡到 n2; i,j——单元。
yf=-10 j=21

图 6-28

力边界网格图

例如: grid 30,20 mod e gen 0,0 0,20 30,20 30,0 gen circle 15,10 5 mod nu reg i=15 j=11 ini sxx=-10 i=1,30 j=1,20 ini syy=-6.5 i=1,30 j=1,20

图 6-29

地应力场

6、开始计算 程序在计算前,还应该设置以下参数:

xf=-5

i=31

(1) 根据实际情况施加重力加速度,施加重力加速度的语句如下: set grav 9.81 (重力加速度的数值可根据实际情况进行选取) (2) 假如程序中允许大变形发生,则使用下面语句: set large (值得注意的是采用大变形后, 在计算过程中网格由于变形过大易出 问题) (3) 根据实际工程设置程序计算的步数: step n (n 为计算的步数) (4) 将计算结果存盘,以便以后调用: save test.sav (test 为存储的文件名) 7、结果显示 模型计算完成后, 需要根据计算结果的应力场、破坏场以及位移场等对工程 进行评价和分析,主要用到以下语句: (1) plot grid:显示计算模型的网格; (2) plot bo:显示模型的边界; (3) plot plas:显示塑性区; (4) plot sig1 fill:显示最大主应力 ?1; (5) plot sig2 fill:显示最小主应力 ?2; (6) plot sdif fill:显示主应力差(?1-?2); (7) plot str:显示主应力矢量场; (8) plot xdis fil:显示 X 方向位移; (9) plot ydis fil:显示 Y 方向位移; (10) plot disp:显示位移矢量场。 8、调用结果命令 (1)call test.txt(或 ca test.dat):调用名为 test 的数据段进行计算; (2)new:重新开始新的建模过程; (3)rest test.sav:调用已经存盘的名为 test 的结果; (4)quit:退出程序,停止使用。 6.6.2 FLAC 程序解题技巧 除了以上介绍的基本使用命令外,FLAC 程序中还有许多其它的使用命令, 这里就不一一介绍了,下面主要介绍一些常用的解题技巧。 1、关于开挖的模拟 涉及到开挖模拟的工程以边坡、洞室最为常见,开挖模拟的语句如下: mod nu i=i1,i2 j=j1,j2 (开挖规则的形状) 例如: grid 30,20 mod e

gen 0,0 0,20 30,20 30,0 mod nu i=20,30 j=15,20

i=20,30 j=15,20

图 6-30

开挖规则形状网格图

或 mod nu region i,j(一个封闭区内任一点的单元号被给定,这个区域将被挖掉, 主要用于不规则形状的开挖) grid 30,20 mod e gen 0,0 0,20 30,20 30,0 table 1 5,5 8,12 15,15 25,12 22,5 12,3 5,5 gen table 1 mod nu region i=15 j=10

region i=15 j=10

图 6-31

开挖不规则形状网格图

2、关于锚杆支护的模拟 图 6-32~6-34 给出了锚杆支护在程序模拟中的相关定义和解释。在工程应用 中,锚杆有三种类型,即端锚、全长锚固和预应力锚固,这三种锚固方式在程序 中的模拟方式不同,具体表述如下。

张拉力

u2[b] u1[b] y u2[a]
端点 b

屈服极限 E 压缩 1 张拉 轴向应变

端点 a

u1[a] x
压缩力

图 6-32

锚杆单元结点位移定义

图 6-33

锚杆材料力学特性

图 6-34

理想的砂浆锚杆体系

(1)端锚 锚杆端锚的具体模拟语句如下: struct cable begin grid i1,j1 end grid i2,j2 seg n prop n1 stru prop n1 e=? yield=? a=? sbond=? kbond=? sfri=? perimeter=? den=? 其中: 锚杆的两端均为结点号; seg n——表示将锚杆划分的区段; prop n1——表示锚杆的力学参数为第 n1 种; e——锚杆的弹性模量; yield——锚杆的屈服强度; a——锚杆的横截面积; sbond——锚杆的变形模量; kbond——锚杆的剪切模量; sfri——锚杆和砂浆之间的内摩擦角; perimeter——钻孔的周长; den——锚杆的密度。 例如:

grid 30,20 mod m gen 0,0 0,20 30,20 30,0 prop d 2.5e-3 bu 2.5e3 sh 1.2e3 c 1.5 fri 42 ten 0.1 mod nu i=20,30 j=15,20 struct cable begin grid 20,16 end grid 10,14 seg 10 prop 1 stru prop 1 e=2e5 yield=0.5 a=0.235e-3 sbond=0.0001 kbond=0.0001 stru prop 1 sfri=0 perimeter=0.2723 den=7.5e-3 fix x i=1 fix x i=31 fix y j=1 set grav 9.81 step 1000

i1,j1 i2,j2

图 6-35

锚杆端锚模拟

(2)全长锚固 全长锚固的力学特性模拟体系见图 6-36,模拟语句如下: struct cable begin grid i1,j1 end x,y seg n prop n1 stru prop n1 e=? yield=? a=? sbond=? kbond=? sfri=? perimeter=? den=? 其中: 锚杆的起始端为结点号,终端为坐标; seg n——表示将锚杆划分的区段; prop n1——表示锚杆的性质为第 n1 种; e——锚杆的弹性模量; yield——锚杆的屈服强度; a——锚杆的横截面积; sbond——锚杆的变形模量; kbond——锚杆的剪切模量; sfri——锚杆和砂浆之间的内摩擦角; perimeter——钻孔的周长;

den——锚杆的密度。

L

锚杆

砂浆体

Kb m

Kb m Ks cs,? Ks

Kb m cs, ?

锚杆轴向刚度 砂浆剪切刚度

锚固结点

砂浆粘聚力,摩擦力

图 6-36

锚杆全长锚固砂浆材料的力学特性模拟体系

例如: grid 30,20 mod m gen 0,0 0,20 30,20 30,0 prop d 2.5e-3 bu 2.5e3 sh 1.2e3 c 1.5 fri 42 ten 0.1 mod nu i=20,30 j=15,20 struct cable begin grid 20,16 end 10,14 seg 10 prop 1 stru prop 1 e=2e5 yield=0.5 a=0.235e-3 sbond=0.42 kbond=4.37 stru prop 1 sfri=30 perimeter=0.2723 den=7.5e-3 fix x i=1 fix x i=31 fix y j=1 set grav 9.81 step 1000

i1,j1 x,y

图 6-37

锚杆全长锚固模拟

3 FLAC 应用举例分析
以上内容详细地介绍了使用 FLAC 程序建模的过程, 以及建模过程中应注意的 问题,下面将用实例进一步说明 FLAC 二维程序以及三维程序的具体应用。 1、北京市引进陕甘宁天然气市内扩建工程圆明园调压站进出线下穿京包铁 路及城市铁路暗挖工程。 (1)工程简介 本工程为北京市引进陕甘宁天然气市内扩建工程圆明园调压站进出线土建 工程的一部分。该工程位于海淀圆明园前八家西侧、下穿城市铁路及京包铁路。 暗挖隧道呈东西走向,全程 61.455m,隧道中线与京包线相交处铁路里程为 K19+579.6,与城铁相交处城铁里程为 K7+897.4。为了确保北京市引进陕甘宁天 然气市内扩建工程圆明园调压站进出线下穿京包铁路及城市铁路暗挖工程施工 期间的安全,本着稳妥可靠、安全第一、经济合理、技术先进的原则,将对隧道 在开挖过程中引起的上覆岩土层和暗挖隧道移动变形进行数值计算, 以确定隧道 施工引起的铁路地基沉降量和隧道变形量。 (2)数值计算模型 由于隧道范围较长且左右对称,假设整个隧道在轴线方向变形很小,用平面 应变模型假设,即垂直于计算剖面方向的变形为零,因此本工程选择其中一个剖 面的右半部分进行数值模拟和力学分析。

?y=2.7 kPa

高填土 7.35m 素填土

注浆区

4.2m

粉质粘土 24.15m 14.25m 城铁暗挖工程岩性分布图

砼衬砌

12.6m

图 6-41

隧道宽度 5.7m,模型宽度取半个隧道宽度的 5 倍,即 14.25m,模型总高度 为 24.15m,其中:自地表至隧道顶部深度 7.35m,隧道高度 4.2m,隧道底部向 下取隧道高度的 3 倍,即 12.6m,图 6-41 是数值计算模型岩性分布图。根据模 型的尺寸, 模型共划分为 15295 个平面单元,构成计算模型单元网格尺寸平均为 0.15× 0.15m(图 6-42), 模型两侧限制水平方向移动, 模型底面限制垂直方向移动。

14.25m 图 6-42 城铁暗挖工程数值计算网格图

(3)数值模拟计算结果分析 图 6-43、6-44 给出了城铁暗挖工程开挖后的位移场,从图中可以看出,开 挖引起的地面最大位移为 17.47mm, 在隧道一侧的地表下沉影响宽度为 8m 左右, 隧道顶部沉降量为 10.52mm。 计算结果表明: 整个隧道的垂直位移场和水平位移场的分布是合理的,引起 的路面沉降量也在合理的范围以内,因此采用标号 C25 的混凝土衬砌以及注浆 加固部分土体可以满足工程的稳定性需要。

24.15m

图 6-43 开挖后隧道垂直位移场

图 6-44 开挖后隧道水平位移场

2、首钢水厂铁矿高陡边坡优化设计 (1)工程简介 首钢水厂铁矿位于河北迁安,是北京地区最大的露天开采铁矿,也是首钢重 要的钢铁原材料基地之一。 由于首钢水厂铁矿边坡最终坡角设计值较小,造成了 剥岩量增加,采剥比减小,很大程度的降低了矿山的经济效益,因此在保证边坡

稳定性的基础上需要对铁矿最终边坡角进行优化。 本研究选择具有代表性的边坡 剖面对边坡进行破坏机理以及稳定性分析,并给出边坡角的具体优化方案。 (2)数值模拟计算模型 典型剖面原设计坡角整体 41? ,根据边坡岩体力学特性,优化后的边坡最终 境界在标高+50m 处分为上、下两段,上部边坡角 43?,下部边坡角 45?。取模型 宽度 1093.85m,垂直方向从水平-600m 起,一直模拟到地表,高度为 925m。模 型共有 41140 个平面单元,单元网格尺寸平均为 5× 5m,其中包括 F36、F11 和 F8 三条断层。图 6-45、图 6-46 显示典型剖面计算模型开挖前的岩性分布和边坡 最终开挖形态。
火山岩溶 岩、角砾岩 长石石英砂 岩、角砾岩 第四系人工堆积物

砾岩

铁矿

925m

铁矿

铁矿

黑云母斜长片麻岩 铁矿

F36

F11

F8

1093.85m
图 6-45 典型剖面边坡数值计算岩性分布图(开挖前)

925m

1093.85m

F36

F11

F8

图 6-46

典型剖面边坡最终开挖形态

(3)数值计算结果分析 计算结果表明,受矿体分布形态和浅部人工堆积物的影响,典型地质剖面边 坡开挖过程中岩体的破坏区比较发育,塑性变形区深入岩体内部较深,大部分塑 性变形区在边坡形成后重新恢复到弹性应力状态, 在断层之间和边坡底部浅层岩 体出现拉破坏和弧形剪切破坏。 当典型剖面边坡采取优化方案(43?—45?)时,边坡的破坏性质基本没有发生 变化,边坡的破坏范围减小,边坡坡脚处没有出现滑移迹线,不会发生局部垮落 (图 6-47)。边坡的最大水平位移量减小到 25cm(图 6-48)。 通过以上的研究表明,水厂铁矿边坡的开挖采用优化方案(43?—45?)是最佳 方案。采取优化方案(43?—45?)时,边坡角有所提高,边坡破坏场的范围有所减 小,边坡坡脚处不出现剪切滑移带,边坡总体产生的水平位移量减小。

图 6-47 典型剖面边坡优化方案(43?—45?)中岩体破坏场

图 6-48 典型剖面边坡优化方案(43?—45?)中岩体水平位移场

4、抚顺老虎台矿开采引发矿震的研究 (1)工程简介 老虎台矿位于辽宁省抚顺市市区南部,抚顺煤田中部。井田范围:西部以矿 区坐标 E3450 米为界,与西露天矿、胜利矿井田相邻;东部以矿区坐标 E8400 米为界,与龙凤矿相邻;南至煤层露头,北至 F1 断层和 F18 断层。老虎台井田 东西走向 4.95 公里,南北宽 2 公里左右,井田面积约 10 平方公里(图 6-55)。

图 6-55

老虎台矿平面位置图

老虎台矿井田地层自下而上分别为:太古界鞍山群、中生界下白垩系、新生 界下第三系抚顺群和第四系。 煤系地层赋存于下第三系抚顺群古城子组,该组为 特厚复合煤层,由 2?38 个分层组成,厚度 8.57?110.5m,平均厚度为 59.58m。 煤层由南向北倾斜,平均倾角 30?左右。顶板岩层主要由油母页岩(平均厚度 194m), 泥岩(平均厚度 422m)和页岩(平均厚度 225m)构成,底板和周围岩层有凝 灰岩、凝灰砂岩、玄武岩和片麻花岗岩等。在井田范围内,有大型断裂 14 条, 其中对开采影响较大的断层有 F25,F16,F18,F26,F7-1,F1 和 F1-A 断层(图 6-56)。

F1-A F1 F16 F18

F26 F25

1100m

F7-1

z y x 3800m

3000m

图 6-56

老虎台矿煤层与断层分布图

老虎台煤矿已有近 100 年的开采历史(表 6-2),目前开拓水平-880m,开采水 平-780m。 近年来, 随着开采深度的增加以及开采位置逐步逼近断层和向斜构造, 开采动力现象日渐频繁, 冲击地压和开采引起的矿震时有发生,冲击强度和震级 急剧增大(表 6-3,图 6-57),并导致井下伤亡事故。
表 6-2 老虎台矿开采深度、时间与开采方式关系对照图表

表 6-3 1998.01?2000.05 抚顺老虎台矿矿震级别与次数统计 矿震级别(ML) 矿震次数 1?1.5 5721 1.6?2.0 1675 2.1?2.5 298 2.6?3.0 63 3.1?3.6 19 合计 7776

图 6-57 1988-2000 年抚顺大于 2.0ML 矿震事件统计

(2)数值模拟计算模型

三维计算模型容括了老虎台矿大部分开采区域。 计算模型沿走向长 3800m(矿 区坐标 E4200?E8000), 沿倾斜宽 3000m(大地坐标 4500?7500), 高度约 1100m(自 地表至-1000m)(图 6-56)。三维模型共划分有 119340 个单元,135298 个结点,其 中包括七条对开采影响较大的断层,有 F25,F16,F18,F26,F7-1,F1 和 F1-A 断层。图 6-58 是三维模型网格图,图 6-59 是典型位置处的剖面图。

图 6-58

三维计算模型网格图
F25 F1 F1-A

F26 (a)

F18

F16

F1

F7-1

(b) 图 6-59 典型位置处的剖面图:(a) E5200,(b) E7550

(3)数值模拟结果分析 通过对漫长的开采历史过程的模拟反演,对应力场有三点基本认识: ① 随着开采深度的增加,煤岩体中的应力集中程度逐渐加大,具体表现在 煤岩体集中应力值的增加和集中范围的局部化; ② 在逐步接近断层区域开采时,断层附近的应力集中分布模式发生本质性 的变化,由跨断层对称分布转变为以断层为核心的非对称集中分布; ③ 在井田范围内, 由于周围开采形成的“孤岛”煤层存在较高的应力集中区, 这给本区煤层开采将造成影响。 通过对漫长的开采历史过程的模拟反演,对破坏场也有三点基本认识: ④ 浅部开采引起的地表破坏角度小于深部开采,岩体破坏角变化范围 50?65?; ⑤ 断层带对岩体破坏具有“屏蔽”效应,当地下开采引起的岩体破坏波及到 断层带附近时,岩体破坏扩展可在断层处截止; ⑥ 开采引起的覆岩破坏在断层附近表现出局部化特征,当该处岩体受开采 扰动并具有运动空间时,局部化破坏带将产生较大规模的运动并伴随能量的释 放。这是构成开采诱发矿震和冲击地压的力学机理。

~1992

图 6-60(a) E5200 剖面剪应力分布

1999-2000
图 6-60(b) E5200 剖面剪应力分布

~1992

图 6-61(a) E5200 剖面岩体破坏区

1999-2000

图 6-61(b) E5200 剖面岩体破坏区
0.9 0.8 0.7 (1992-2000)

Shear Disp., mm

0.6 0.5 0.4 0.3 0.2 0.1 0 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2

78001-1

9.4

9.6

Step

(a)

30 (1992-2000) 25

Shear Disp., mm

20 15 10 5 0 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2

78001-1

9.4

9.6

Step

(b) 图 6-62 1992?2000 年井田西部剪切动力学特征 (a) F25 断层露头移动特征,(b) F25 断层与煤层交界处


相关文章:
FLAC3D(Fast Lagrangian Analysis of Continua)问题讲...
FLAC3D 是一个三维有限差分程序,目前已发展到 V2。1 版本。 FLAC3D 的输入和一般的数值分析程序不同,它可以用交互的方式,从键盘输入各种 命令,也可以写成命令(...
1FLAC 3D基本介绍
任何的本构模型而不用调整解决的算法程序, 在这一点上 FLAC 3D 是非常好的...1.3.2. 可选部分这里主要讲了动力分析、热力学分析、蠕变材料的性能和用户自...
FLAC3D接触面建立方法
impgrid 1.flac3d 得到的就是已经建立好接触面的模型 下面进行简单的计算,以验证...j Kn 法向刚度通常设置一个较大的值,但是这个值不能太大,否则又会存在收敛...
FLAC3D学习笔记
经过这个过程后,FLAC 学习一阶段完成,达成目标:熟练掌握全部基本操作、基本命 令,能独立处理简单问题,养成好的操作习惯,形成规范的 FLAC 思维方式。 以后的事情,...
FLAC讲义
FLAC讲义_电力/水利_工程科技_专业资料。第一讲 FLAC 技术的基本原理和应用范围...2、做好以上的规划准备后,须在计算机上建立模型,建模过程中具体操作 步骤和常用...
FLAC3D中一些问题的讨论
FLAC3D 中显式的时步对应于最小区域中信息从一个...对于这些方法,本人也曾做过一些试验,效果并不好,...而且对这么多单元来讲每个计算只 有 500 步好象太...
FLAC3D基础知识介绍1
FLAC3D基础知识介绍1_建筑/土木_工程科技_专业资料。...第一主应力 2.xdisp、ydisp、zdisp、disp 用 ...相比局部阻尼可以使稳定状态达到更好的收敛,这时网 ...
FLAC_3D快速入门(手册翻译版——一米)
FLAC_3D快速入门(手册翻译版——一米)_计算机软件及应用_IT/计算机_专业资料。FIAC3D入门文档,欢迎下载快 速 入 门 (GETTING STARTED) 版本:flac3d 3.0版(FTD...
FLAC3D 实例命令流1
FLAC3D 实例命令流1_建筑/土木_工程科技_专业资料。FLAC3D 实例命令流第1 部分 命令流按照顺序进行 2-1 定义一个 FISH 函数 new def abc abc = 25 * 3 ...
FLAC3D学习笔记(自己总结版)
的面,采用相同的集合变化 率,网格才会重合,ratio 对于减少单元的数目有很好的...第一:相对收敛准则: 一般而言,大多数问题可以采用FLAC3D默认的收敛标准(或称相对...
更多相关标签: