FLAC_3D_粘弹性模型的二次开发及其应用_徐平

t
+
Δt 时刻差分点的偏应力张量 ;e0ij 为 t 时刻差分点的
偏应变张量 ;e′ij 为 t +Δt 时 刻差分点 的偏应变 张
量 ;Δe′ij 为在 Δt 内差分点的偏应变张量增量 。
将式(6)整理并解得 t +Δt 时单元的偏应力张
量为
S ′ij
=(C
1S
0 ij
+C 2Δe′ij +C 3e0ij)· C4 ,
广义 Kelvin 模 型如 图 1 所示 , 其微分型的本 构关系可表示为[ 8]
图 1 广义 Kelvin 模 型 Fig .1 Kelvin-Voigt model
σ+p1 ﹒σ = q0 ε+q1﹒ε,
(1)
式中
:p1
=E
η1 +E1
;q0
=EE+EE1 1
;q1
=E
Eη1 +E
1
;η为材
尽管如此 , FLA C3D软件在应用过程中仍然暴露 出了不少的缺点和不足 。 如 F LA C3D的前处理部分 仍采用不太直观的文本命令格式 , 通过数据文件与 软件连接 , 虽然 F LAC3D 也提供了内嵌式的 F ISH 语 言 , 但在处理复杂网格时仍给使用者带来了很大的 不便 。文[ 6] 就其前处理部分的不足 , 借助于 Fortran 语言对此项进行了功能开发与研制 ;文[ 7] 为了 更细致地了解岩石变形局部化现象 , 研制了局部化 图形再现后的处理程序 。 由于岩体力学特性的复杂 性 , 其变形往往呈现出明显的时间性 , F LAC3D 程序 提供了 5 种适用于不同特性岩体介质的粘弹塑性分
F LAC3D粘弹性模型的二次开发及其应用
徐 平1 , 李云鹏2 , 丁秀丽1 , 王芝银2
(1 .长江科学院 岩基研 究所 , 湖北 武汉 430010 ;2 .石油大学 , 北京 102249)
摘要 :针对 F LAC3D(V ersion 2 .0)提供的粘弹性本构模型难以描述岩体变形随 时间发展趋于稳定的现象 , 对 F LAC3D 的粘弹性模拟功能进行了开发 , 导出广义 Kelvin 模型的 中心差分格 式 , 并 采用 F LAC3D内嵌 功能较强的 FISH 语 言 和 FLAC3D的相关命 令研制了广义 Kelvin 模 型的接口程序 。 通过检 验算例说 明了开 发思路 和接口 程序的 正确性 , 同时还指出了二次开发时应注意的一些问题 , 并 对基坑开挖进行了实例分析 , 结果符合基坑开挖变形规 律 。
stress iso chromatic zone pa ttern
图 5 x 向位移随时间变化规律及等色区图 Fig .5 Curv es of displacement in the x-direction vs
time and disp .isochromatic zone pattern
线为中心坐标点(7 .5 , 2 .5 , 0 .5)处单元 z 向正应力 值 , 其中 , 依次由上向下第一条曲线为解析曲线 , 第 二条曲线为开发的 F LAC 解 , 其相对误差为8 .5 %。 图 5 显示的是坐标点为(7 .0 , 3 .0 , 1 .0)处的 x 方向 的位移值 , 上方曲 线为解 析解 , 下方曲 线为 F LAC
模拟流变时间为8 .64 d , 图 4 中显示的应力曲
12
长江科学院院 报
2004 年
图 3 计算模型及载荷示意图 Fig .3 Illustration o f numerical model and
evenly load along x-axis direction
图 4 单元应力随时间变化规律及等色区图 Fig .4 Curves of element stress vs time and
挖深 度 为0 .8 m , 共 分 5 步 , 开挖区域为
1 m ×4 m ×2 m 。在
x (0), z(0), y (10) 面上均为 沿坐标 方
向单向约束 , y(0)顶
面为 自 由, 在 x (10), z (10)面上 作 用有沿 y 向线性变 化的 分布 载 荷 。 采
用广 义 Kelvin 模示相应物理量的增量 。
为了有效地与 F LAC3D 软件连接 , 采用中 心差
分格式将上式改写为
S
′ij
-S
0 ij
Δt
+
S
′ij +S 2 p′1
0 ij
=
2
q 1′Δe′ij p1′ Δt
+2
pq01′′(e′ij
+e0ij) ,
(6)
式中
:S
0 ij
为
t
时刻差分点的偏应力张量 ;S′ij 为
析模型 , 如 Maxw ell 模型描述的是变形随时间变化 而不断增大的情况 , 即属于非稳定流变模型 , 其对岩 体变形随时间的发展趋于稳定状态的情况就很难描 述 。F LA C3D中其它流变模型均有其特定的适用范 围 , 这就明显影响了 F LA C3D的广泛应用 。 故此 , 本 文以 F LAC3D(V ersion 2 .0)为基础 , 利用其 内嵌的 F ISH 语言工 具 , 开发研制了广义 Kelvin 模型接口 程序 , 从而扩大了 F LAC3D 程序的适用性 。 通过算 例分析表明 , 所研制的接口程序开发思路正确 , 可用 于分析工程实际问题 。
q1′=GG0 0+ηG1 1
;q0′=GG0 0+GG1 1
;p 1′=G 0
η1 +G1
;
(4)
S ij 为偏应力张量 ;eij 为偏应变张量 ;G0 , G1 分别为
弹性剪切模量和粘弹性剪切模量 。
2 .2 差分格式 将式(4)中 P″, Q″代入式(3), 改为增量形式并
化简得
p 1′ΔΔStij +Sij =2(q′ΔΔeti j +q0′ei j),
解 , 相对误差为9 .5 %。 由应力等值线和位移等值线
可见其规律性基本上与理论解分布规律接近 。 由图
5 可知 , 广义 Kelvin 模型所获结果达到稳定状态的 时间约为 2 d 左右 。
4 .2 工程应用
某深基坑 , 由对 称性取 1/ 4 区域(10
m ×10 m ×10 m)进
行分析研究 , 每步开
第2期
徐 平 等 F LAC3D粘弹性模型的二次开发及其应用
11
式中
:P′=(p 1
d dt
+1);Q′=(q 1
d dt
+q 0) 。
将式(2)扩充到三维 状态 , 得到广义 Kelvin 模
型的三维本构关系为
P″S ij =2Q″eij ,
(3)
式中 :P″=(p 1′ddt +1);Q″=(q1′ddt +q0′);
关 键 词 :粘弹性 ;二次开发 ;广义 K elvin 模型 ;FLAC3D软件 中图分类号 :T P311 .52;T U 452 文献标识码 :A
1 概 述
F LA C3D软件是一个用于分析和解决工 程力学 问题的三维显式差分法程序 , 可以有效地模拟岩土 体及其它材料结构的工程力学行为 , 与有限元法等 应用软件相比各自具有不同的优缺点 。 在 F LA C3D 软件中 , 具有能满足岩土工程结构分析的基本单元 模式 , 包含了进行弹性分析的各向同性 、异性和横观 各向 同 性 的 弹 性 本 构 模 型 , 塑 性 模 型 有 MohrCoulomb 弹塑性 、Drucker-P rager 弹 塑性模型等 ;描 述材料蠕变特性的有 Maxwell 、幂指数等 5 种蠕变 本构模型 。对于静态 、动态 、蠕变 、渗流 、温度效应及 其相互耦合的工程问题 , 通过选择适当的模式可以 实现对不同类型问 题的求解 。 目前 F LA C3D (Version 2 .0)软件在我国尤其是在岩土工程界得到了广 泛的认可和应用[ 1 ～ 5] , 并取得了具有工程指导意义 的应用成果 。
2 广义 Kelvin 模型及其差分格式
工程实践表明 , 许多岩体及支护结构在工程挠动
后产生一定变形 , 但随着时间的推移变形逐渐趋于稳
定 。为使 FLAC3D 软件能够有效地模拟这一变形特
性 , 本文给出了广义Kelvin模型的本构方程 , 并导出
适应于 F LAC3D软件的差
分格式 。
2 .1 广义 Kelvin 模型
荷载计算获得节点不平衡力 , 并以等效在节点上的 力和质量建立以节点位移速率为未知量的节点运动 方程 , 解得节点位移速率后并进一步求得某一时步 的单元应变增量 、单元应力增量 , 以及经累加得到当 前时步的单元总应力 。 然后再以单元应力等求得节 点不平衡力 。 如此循环上述过程 , 最终得到工程结 构的应力与变形结果 。 根据上述分析思路 , 并利用 F LAC3D内嵌的 F ISH 语 言开发的接口 程序主要流 程如图 2 所示 。
4 算例验证及应用
通过上述思路及 F LA C3D命令与 F ISH 语言 , 研 制并完成了广义 Kelvin 接口程序 。以 下给出检验 算例和一工程实例 , 以说明所开发程序的有效性 。 4 .1 验证算例
一粘弹性矩形厚板 , 尺寸为 20 m ×5 m ×2 m , 弹性剪切模量 G0 =5 ×108 P a , 体积模量 K =1 .5 × 109 Pa , 粘弹性剪切模量 G1 =3 .2 ×109 Pa , 粘性参 数 η1 =1 .574 ×1010 P a·d 。 由对称性取 1/4 模型进 行研究 , 在 x =0 .0面上 x , y 方向约束 , 在 z =0 .0 , z =2 .0面为 z 单向约束 。在 x =10 .0面上作用有 沿 x 方向的均匀边界载荷(σxx =20 MPa), 如 图 3 所示 。用开发的广义 Kelvin 模型进行计算 , 其结果 与理论解的比较列于图 4 和图 5 中 。