5.预应力锚索框架(地梁)附加应力分布FLAC3D分析

FLAC 程序的基本原理和算法与离散元相似，但它却象有限元那样适用于多 种材料模式与边界条件的非规则区域的连续问题求解；在求解过程中，FLAC 采 用了离散元的动态松驰法，不需要求解大型联立方程组（刚度矩阵）。同时，同 以往的差分分析方法相比，FLAC 不但可以对连续介质进行大变形分析，而且能 模拟岩体沿某一软弱面产生的滑动变形，FLAC 还能在同一计算模型中针对不同 的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为。此外， 该方法还可考虑锚杆、挡土墙等支护结构与围岩的相互作用。
<
∆σ ij
>
+
v E
I1δ ij
式中，v 为泊松比；E 为弹模；I1 为应力第一不变量。
（5－9）
这样，通过上述各式的迭代求解，便可求出每一迭代时步相应各单元的应力
和应变值。
由莫尔库仑屈服准则：
τ n = −σ ntgϕ + C 将式（3－10）转换成用单元应力表示的形式：
( ) f = σ 3 − Nϕσ 1 + 2c Nϕ 1/ 2
∑ < ∂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
+
∆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 也已求出。 进一步得结点位移：
式中， Nϕ = (1 + sinϕ )/(1 − sinϕ )。
（5－10） （5－11）
根据各单元 f 值的大小便可判断单元屈服与否（f<0 屈服；否则不屈服）。 上面已求出了各域（单元）的应力，下面来求各结点的平衡力。
由结点的运动方程：
∂σ ij ∂xi
+ ρgi
=
ρu&&i
式中， u&&i 为总加速度；gi 为重力加速度。 对（3－11）沿积分路径积分（见图 3－4）得：
∫ ∫ 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 条边求和的形式：
5.1 FLAC 的基本原理
FLAC 是快速拉格郎日差分分析（Fast Lagrangian Analysis of Continua）的简 写。FLAC 是力学计算的数值方法之一，该名词渊源于流体动力学，它研究每个 流体质点随时间变化的情况，即着眼于某一个流体质点在不同时刻的运动轨迹、 速度及压力等。快速拉格郎日差分分析将计算域划分为若干单元，单元网格可以 随着材料的变形而变形，即所谓的拉格朗日算法，这种算法可以准确地模拟材料 的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以 及模拟施工过程等领域有其独到的优点。
（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
FLAC 采用差分方法，每一步的计算结果与时间相对应。程序采用人机交互 式的批命令形式执行，在计算过程中可以根据施工过程对计算模型和参数取值等 进行实时地调整，达到对施工过程进行实时地仿真的目的。
具体地讲，FLAC 的基本原理如下： FLAC 用差分方法求解，因此首先要生成网格。将物理网格（图 5-3-1）y 映 射在数学网格（图 5-3-2）上，这样数学网格上的某个编号为 i，j 的结点就与物
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&Fra bibliotek( i
2)
+ u&i(3)
n
j
∆S
(b) i
+
u&
(3) i
+ u&i(1)
n
（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.预应力锚索框架(地梁)附加应力分布 FLAC3D 分析
传统的土压力的测试过程中，土压力盒的埋设对土压力有一定的扰动，而 且价值昂贵，只能在有限深度内测试有限个点，这对研究预应力锚索框架对边坡 加固后的附加应力分布是远远不够的。基于此，我们在对附加应力测试试验研究 的基础上，采用美国 ITASCA 咨询集团公司开发的 FLAC3D 快速拉格郎日差分 程序做了相应的数值模拟。
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
理网格上相应的结点的坐标 x，y 相对应，这一过程可以想象为数学网格是一张 橡皮做的网，拉扯以后可以变为物理网格的形状。
图 5-1 物理网格
图 5-2 数学网格
假定某一时刻各个节点的速度为已知，则根据高斯定理可求得单元的应变 率，进而根据材料的本构定律可求得单元的新应力。
根据高斯定理，对于函数 F 有：