预应力锚索FLAC3D分析

G＝ E0 2(1 + u)
K = E0 3(1 − 2u)
（3－18）～（3－20）式中， u 为岩土体的泊松比。 各层所对应的组(Group)编号和各计算参数取值见表 5-1。
表 5-1 岩土体材料物理力学参数取值
FLAC 群组
岩土层
泊松比 重度 变形模量 剪切模量 体积模量 φ C 抗拉强度 （ µ ） KN/m3 （MPa） （MPa） （MPa） （o）（Kpa） (KPa)
5.3 计算模型
模型范围：计算模型前缘以公路中心线为界，后缘以公路中心线向后延伸
60m，沿坡面方向上取 52.5m，底面以路面以下 20m 为界。 地质概化模型：考虑地层单元参数取值的方便，综合 XX 省水文地质工程地
质勘察研究院的岩土工程勘察报告，将坡体结构概化为四层第四系坡积层、残积 层、风化花岗闪长岩层和新鲜花岗闪长岩层等四层。本次计算中，共划分了 5600 个单元，6290 个节点。
j
∆S
(c i
)
（5－5）
同理可求出 < ∂u& j > 值。 ∂xi
由几何方程可求得单元的平均应变增量：
< ∆eij
>=
1
⎡ ⎢<
2 ⎢⎣
∂u&i ∂x j
>+<
∂u& j ∂xi
⎤ >⎥∆t
⎥⎦
由广义虎克定律，各向同性材料的本构方程为：
（5－6）
σ ij = 2µε ij + λθ ⋅ δ ij
因此，下面分别对天然状态下以及两种加固措施下的边坡应力进行分析研 究。
由于垂直坡面方向上应力直接影响锚索框架（地梁）底部反力，同时也为了 便于与测试结果的对比研究，特别对计算结果利用下式（柯西公式）计算该方向 上的应力，并做重点分析。
σ N = l 2σ x + m2σ y + n2σ z + 2mnτ yz + 2nlτ zx + 2lmτ xy 式中，（l,m,n）坡面的方向余弦。
FLAC 程序的基本原理和算法与离散元相似，但它却象有限元那样适用于多 种材料模式与边界条件的非规则区域的连续问题求解；在求解过程中，FLAC 采 用了离散元的动态松驰法，不需要求解大型联立方程组（刚度矩阵）。同时，同 以往的差分分析方法相比，FLAC 不但可以对连续介质进行大变形分析，而且能 模拟岩体沿某一软弱面产生的滑动变形，FLAC 还能在同一计算模型中针对不同 的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为。此外， 该方法还可考虑锚杆、挡土墙等支护结构与围岩的相互作用。
计算方案：为系统研究各种坡面抑制结构（框架、地梁）对坡面及坡体内部 附加应力分布规律，本次数值模拟中，在边坡初始应力场模拟的基础上又计算了 两类边坡方案框架+地梁边坡和框架边坡在三级荷载（框架+地梁：300KN、 600KN、800KN，框架：300KN、500KN、650KN）的 6 种方案，共 7 种方案。
理网格上相应的结点的坐标 x，y 相对应，这一过程可以想象为数学网格是一张 橡皮做的网，拉扯以后可以变为物理网格的形状。
图 5-1 物理网格
图 5-2 数学网格
假定某一时刻各个节点的速度为已知，则根据高斯定理可求得单元的应变 率，进而根据材料的本构定律可求得单元的新应力。
根据高斯定理，对于函数 F 有：
边界条件：根据地貌特征及对称性，本计算模型中，在周边及底部施加约束 边界。
根据上述特征，建立如下计算模型（图 5-1~图 5-3）
图 5-1 初始 FLAC3D 计算模型
图 5-2 加锚索和工作框架计算模型 图 5-3 加锚索和破坏框架及地梁计算模型
5.4 模型参数
无论是常规计算还是数值计算，模型中地质力学参数的取值都是关系但计
5.1 FLAC 的基本原理
FLAC 是快速拉格郎日差分分析（Fast Lagrangian Analysis of Continua）的简 写。FLAC 是力学计算的数值方法之一，该名词渊源于流体动力学，它研究每个 流体质点随时间变化的情况，即着眼于某一个流体质点在不同时刻的运动轨迹、 速度及压力等。快速拉格郎日差分分析将计算域划分为若干单元，单元网格可以 随着材料的变形而变形，即所谓的拉格朗日算法，这种算法可以准确地模拟材料 的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以 及模拟施工过程等领域有其独到的优点。
FLAC 采用差分方法，每一步的计算结果与时间相对应。程序采用人机交互 式的批命令形式执行，在计算过程中可以根据施工过程对计算模型和参数取值等 进行实时地调整，达到对施工过程进行实时地仿真的目的。
具体地讲，FLAC 的基本原理如下： FLAC 用差分方法求解，因此首先要生成网格。将物理网格（图 5-3-1）y 映 射在数学网格（图 5-3-2）上，这样数学网格上的某个编号为 i，j 的结点就与物
后期的地应力场是初始天然应力场与施加结构荷载后附加应力的迭加。附加 应力的布辛内斯克解是利用弹性理论，假定地质体为均质线弹性体推导出来的， 实际上，根据一些学者的试验研究及测量结果，土体中附加应力的分布是受土体 的非线性、非均质和各向异性的影响的。也就是初始天然应力场在一定程度上影 响着施加结构荷载后的附加应力应力场，而附加应力场是结构设计的一个重要基 础。
新鲜花岗 1 闪长岩 0.18 29
5000 2118.6 2604.2 55 500 100
第四系坡 2 积层 0.3 19.05 27.6
10.6
23.0 17.56 16.2
1
3 残积层 0.29 18.62 36.5 14.2 29.0 29.04 25.6 1
风化花岗 4 闪长岩 0.2
27
∑ < ∂F ∂ xi
>= 1 V
N
Fi ni ∆Si
（5－3）
式中，∆Si 是多边形的边长；Fi 是 F 在∆Si 上的平均值。
假定以速度 u&i 代替式（3－3）中的 Fi，且 u&i 取边两端的结点（即差分网络的
角点）a 和 b 的速度平均值，则：
∑[( ) ] < ∂u&i >= 1 ∂x j 2V N
u&
a i
+ u&ib
n j ∆Si
≈ ∂u&i ∂x j
（5－4）
对于三角形单元（如图 3－3）：
[( ) ( ) ( ) ] < ∂u&i >= 1 ∂x j 2V
u&i(1)
+
u&
(2) i
n
j
∆S
(a i
)
+
u&
( i
2)
+ u&i(3)
n
j
∆S
(b) i
+
u&
(3) i
+ u&i(1)
n
锚索单元：为两节点的线单元，只允许轴向变形，一个自由度，只能承受拉 力或压力，而不能承受弯矩。
结构单元之间以及与岩土体之间的相互作用是联接来完成的，FLAC 中有两 种联接：节点与单元之间的联接和节点与节点之间的联接，其类型有三种：自由 联接（节点速度与和它的联接体无关）、刚性联接（节点与联接体连接在一起， 速度相同）、可变连接（节点与联接体之间为弹性）。力即通过联接关系来实现力 的传递，进而实现其相互作用。本次计算中采用的是刚性联接。
+
∆t
/
2) − u&i
∆t
(t
−
∆t
/
2)
（5－14）
（5－15）
u&i (t + ∆t / 2) = u&i (t − ∆t / 2) + u&&i (t)∆t
（5－16）
其中， u&i (t − ∆t / 2) 为结点上一时步的速度，而 u&&i (t)∆t 也已求出。 进一步得结点位移：
应变率
（5－7）
式中，λ，µ为拉梅常数；θ = εij = ε11 + ε22 + ε33 ，即体积应变；
δ ij
=
⎧1 ⎩⎨0
(i = j) (i ≠ j)
因此单元的平均应力增量可表达成：
<
∆σ ij
>=
λδ ij
<
∆θ
>
+
v E
I1δ ij
（5－8）
同时，若以应力表示应变，则其本构关系为：
<
∆eij
>= 1 + v E
（5－12）
∑ ρu&&i
=
1 V
< σ ij > n j ∆Si + ρgi
（5－13）
∑ 其中， < σ ij > n j ∆Si 为某结点周围单元作用在该结点的集中力。
u&&i
=
1 m
F
+
gi
式中，F 指作用在结点中的合力（净力）。 利用中心差分，得某结点加速度和速度：
u&&i
(t)
=
u&i
(t
本 次 计 算 中 ， 应 用 的 结 构 单 元 有 两 种 ： 梁 单 元 （ BeamSel ） 和 锚 索 单 元 （CableSel），其主要特点如下：
梁单元：假定为线性弹性材料，为线单元，具有两个节点，每个节点 6 个自 由度（3 个位移分量和 3 个转动分量），可以承受一定弯矩，允许一定的挠曲。
4200
1750 2333.3 42 300
50
5.5 计算结果分析
岩土体中的应力就其产生原因有两种：由岩土体本身自重引起的自重应力和 由外荷载引起的附加应力。一般来说，自重应力引起的地基压缩变形在地质历史 上绝大部分已经完成。所以，框架底部所承受的压力主要为由锚固力所引起的附 加应力。
锚索是一种主要的承受拉力的杆状构件，通过钻孔将钢绞线或高强度的钢丝 固定于深部稳定的地层中，并在被加固体表面通过张拉产生预应力，控制岩体的 变形和位移，改善岩体的应力状态，提高岩体强度，从而达到加固边坡稳定和限 制其变形的目的。而在预应力锚索框架中，是将锚索锚固到框架上，锚固力首先 作用框架上，通过框架传递给岩土体，从而产生附加应力。
式中， Nϕ = (1 + sinϕ )/(1 − sinϕ )。
（5－10） （5－11）
根据各单元 f 值的大小便可判断单元屈服与否（f<0 屈服；否则不屈服）。 上面已求出了各域（单元）的应力，下面来求各结点的平衡力。
由结点的运动方程：
∂σ ij ∂xi
+ ρgi
=
ρu&&i
式中， u&&i 为总加速度；gi 为重力加速度。 对（3－11）沿积分路径积分（见图 3－4）得：
5.2 FLAC3D 中结构单元计算原理
地质力学计算与设计的一个重要部分是结构的支护与抗滑稳定性设计。 FLAC3D 可以模拟任意形状、任意特性的结构体与岩土体的相互作用，以及力作 用在结构体上或岩土体时，结构体和岩土体的力学反应。FLAC3D 可以模拟 4 种结构单元（梁、锚索、桩及板壳单元）。FLAC3D 中，结构单元的计算模式与 其他计算一样，利用的依然是与隐式矩阵算法相对应的解全运动方程的显式拉格 郎日算法。结构单元力学性态的计算模拟可以通过几何大变形和小变形两种计算 模式计算，而且还设置动力状态来模拟结构与岩土体的动力反应。
∫ ∫ B Fni ds =
∂F dV V ∂xi
（5－1）
式中，V 是函数求解域（或单元）的体积；B 是 V 的边界；ni 是 V 的单位外法 线矢量。
定义梯度 ∂ F 的平均值为： ∂ xi
∫ < ∂F >= 1 ∂F dV ∂ xi V V ∂ xi
（5－2）
式中，< >表示求平均值。
对于一个具有 N 条边的多边形，上式可写成对 N 条边求和的形式：
算准确度和计算精度。所以各地层的参数取值是数值模拟中关键的一环。
在 FLAC 程序中，岩土体变形的变形参数采用的是剪切模量（G）和体积 模量（K），而不是直接采用压缩模量。因此，在具体计算时，需首先采用如下 公式将压缩模量 Es 转化为变形模量 E0：
E0
=
(1
−
2u 1−
2
u
)
E
s
然后再利用以下公式将变形模量转化为剪切模量和体积模量：
对于每个单元： 本构定律
新的应力
高斯定理理 对于每个节点：
速度
运动定律
Fi=σijnjL 节点力
图 3-5 FLAC 的计算循环
ui (t + ∆t) = u(t) + Fra Baidu bibliotek&i (t + ∆t)∆t
（5－ 17）
按照上述思路，通过迭代求解，便可求出各个时步边坡上各单元（或结点） 的应力、变形值，进而可模拟出整个边坡变形破坏过程。拉格朗日差分法计算循 环如图 5-5 所示。
5.预应力锚索框架(地梁)附加应力分布 FLAC3D 分析
传统的土压力的测试过程中，土压力盒的埋设对土压力有一定的扰动，而 且价值昂贵，只能在有限深度内测试有限个点，这对研究预应力锚索框架对边坡 加固后的附加应力分布是远远不够的。基于此，我们在对附加应力测试试验研究 的基础上，采用美国 ITASCA 咨询集团公司开发的 FLAC3D 快速拉格郎日差分 程序做了相应的数值模拟。
<
∆σ ij
>
+
v E
I1δ ij
式中，v 为泊松比；E 为弹模；I1 为应力第一不变量。
（5－9）
这样，通过上述各式的迭代求解，便可求出每一迭代时步相应各单元的应力
和应变值。
由莫尔库仑屈服准则：
τ n = −σ ntgϕ + C 将式（3－10）转换成用单元应力表示的形式：
( ) f = σ 3 − Nϕσ 1 + 2c Nϕ 1/ 2