工程地质数值法

某路基工程施工过程数值模拟

摘要

本文首先对FLAC3D软件进行了介绍，简明阐述了其特点、应用范围及不足；然后结合具体路堤工程，采用FLAC3D软件对施工过程进行了模拟，生成了初始竖向和水平应力云图、第一次填筑及填筑结束时的沉降云图及水平位移云图；最后生成了路基中心点和坡脚节点的沉降曲线。

关键词：FLAC3D；数值模拟；应力云图；沉降云图；位移云图

1 FLAC3D的功能与特性

自R.W.Clough 1965年首次将有限元引入土石坝的稳定性分析以来，数值模拟技术在岩土工程领域获得了巨大的进步，并成功解决了许多重大工程问题。特别是个人电脑的出现及其计算性能的不断提高，使得分析人员在室内进行岩土工程数值模拟成为可能，也使得数值模拟技术逐渐成为岩土工程研究和设计的主流方法之一。数值模拟技术的优势在于有效延伸和扩展了分析人员的认知范围，为分析人员洞悉岩体、土体内部的破坏机理提供了强有力的可视化手段。FLAC系列软件的出现，为岩土工程研究工作者提供了一款功能强大的数值模拟工具。

1.1 FLAC3D主要特点

FLAC(Fast Lagrangian Analysis of Continua)是由Itasca公司研发推出的连续介质力学分析，是该公司旗下最知名的软件系统之一，FLAC目前已在全球七十多个国家得到广泛应用，在国际土木工程（尤其是岩土工程）学术界和工业界享有盛誉。

FLAC3D界面简洁明了，特点鲜明。其使用特征主要表现为：命令驱动模式、专一性、开放性。作为有限差分软件，相对于其他有限元软件，在算法上，FLAC3D 有以下几个优点：采用“混合离散法”来模拟材料的塑性破坏和塑性流动，比有限元中通常采用的“离散集成法”更准确、合理；即使模拟静态系统，也采用动态运动方程进行求解，这使得FLAC3D模拟物理上的不稳定过程不存在数值上的障碍；采用显示差分法求解微分方程。采用FLAC3D进行数值模拟时，必须指定有限差分网格、本构关系和材料特性、边界和初始条件，这是FLAC3D求解的一般流程。

1.2 FLAC3D的应用范围

FLAC3D的应用范围已拓展到土木建筑、交通、水利、地质、核废料处理、石油及环境工程等领域，成为这些专业领域进行分析和设计不可或缺的工具。其研究范围主要集中在以下几个方面：

●岩体、土体的渐进破坏和崩塌现象的研究；

●岩体中断层结构的影响和加固系统（如喷锚支护、喷射混凝土等）的模

拟研究；

●岩体、土体材料固结过程的模拟研究；

●岩体、土体材料流变现象的研究；

●高放射性废料的地下存储效果的研究分析；

●岩体、土体材料的变形局部化剪切带的演化模拟研究；

●岩体、土体的动力稳定性分析、土与结构的相互作用分析以及液化现象

的研究等。

1.3 FLAC3D的不足

毋庸置疑，FLAC3D是十分优秀的岩土工程数值模拟软件，其实用性和专业性得到了广泛证实。但FLAC3D也存在着诸多不足，主要集中在以下几个方面：

●求解时间受网格尺寸影响很大；

●某些模式下的计算求解时间很长；

●前处理功能较弱。

●FLAC3D对于复杂三维模型的建立仍然十分困难。

2 某路基工程施工过程模拟

2.1工程概况

该路基工程地基计算深度为50m，分为两层，上部为回填土，厚度10m，下部为粘土，厚度40m；路基计算宽度为200m，填筑高度为5m，坡度为1:1.5。各土层物理、力学参数见表1。

图1 路堤施工几何模型

表1 各土层物理力学参数

土层名称ρ/kg·m-3 c /kPa φ / (°) E /MPa ν

回填土1500 10 15 8.0 0.33 粘土1800 20 20 4.0 0.33 2.2模型建立

由于几何模型具有对称性，可以采用1/2模型进行分析。首先建立坐标系，坐标系的原点O设置在地基表面与模型对称轴的交点，水平向右为X方向，竖直向上为Z方向，垂直于分析平面的方向为Y方向。

网格的建立按照分区域建模的思路进行，如图2所示。由于路基坡脚的位置存在一个关键点，所以将模型划分成5个矩形区域，对每个区域按照控制点利用brick单元建立网格，并进行分组后赋值。考虑到网格尺寸的一致性，Y方向只设置一个单元，该方向单元尺寸为5m。

图2 模型网格建立思路

网格建立以后，首先设置边界条件。对底部边界节点的X、Y、Z三个方向的速度进行约束，相当于固定支座，对X两侧的边界进行水平速度约束。由于Y 方向只设置一个单元长度，所以对模型中所有节点的Y方向速度均进行约束，相当于进行平面应变分析。

2.3初始应力计算

在路基施工前，需要将路基部分网格赋值为空模型，而将地基部分的网格赋值为Mohr模型。由于实例中存在null模型，不能采用solve elastic的求解方法获得初始应力，所以采用分阶段的弹塑性求解方法。先将Mohr模型的凝聚力c 值和抗拉强度σt赋值为无穷大进行求解，保证在重力作用下单元不至于发生屈服，然后再将Mohr模型参数赋值为真实值，再进行求解。

FLAC3D 3.00

Itasca Consulting Group, Inc.Minneapolis, MN USA

Step 1888 Model Perspective 17:55:28 Mon Apr 21 2014Center:

X: 5.000e+001 Y: 2.500e+000 Z: -2.500e+001Rotation: X: 0.000 Y: 0.000 Z: 0.000Dist: 2.790e+002

Mag.: 1Ang.: 22.500

Contour of SZZ

Magfac = 0.000e+000 Gradient Calculation

-8.5283e+005 to -8.0000e+005-8.0000e+005 to -7.0000e+005-7.0000e+005 to -6.0000e+005-6.0000e+005 to -5.0000e+005-5.0000e+005 to -4.0000e+005-4.0000e+005 to -3.0000e+005-3.0000e+005 to -2.0000e+005-2.0000e+005 to -1.0000e+005-1.0000e+005 to 0.0000e+000 0.0000e+000 to 9.7702e+000 Interval = 1.0e+005

图3 初始竖向应力云图

FLAC3D 3.00

Itasca Consulting Group, Inc.Minneapolis, MN USA

Step 1888 Model Perspective 17:55:56 Mon Apr 21 2014Center:

X: 5.000e+001 Y: 2.500e+000 Z: -2.500e+001Rotation: X: 0.000 Y: 0.000 Z: 0.000Dist: 2.790e+002

Mag.: 1Ang.: 22.500

Contour of SXX

Magfac = 0.000e+000 Gradient Calculation

-4.1992e+005 to -4.0000e+005-4.0000e+005 to -3.5000e+005-3.5000e+005 to -3.0000e+005-3.0000e+005 to -2.5000e+005-2.5000e+005 to -2.0000e+005-2.0000e+005 to -1.5000e+005-1.5000e+005 to -1.0000e+005-1.0000e+005 to -5.0000e+004-5.0000e+004 to 0.0000e+000 0.0000e+000 to 1.8693e+001 Interval = 5.0e+004

图4 初始水平应力云图