岩土工程数值计算方法

1、岩土工程数值计算方法报告学院：土木与环境工程学院姓名：xxxxxx学号：xxxxxxxx三维有限差分稳定性分析一、 FLAC3D 基本原理FLAC(Fast Lagrangian Analysis of Continua)是由美国Itasca 咨询公司研究开发的显式有限差分程序，可用于工程力学计算，模拟岩石、土等材料的力学行为。由于其采用了显式拉格朗日算法及混合离散划分单元技术，使得该程序能较好地模拟地质材料在达到强度极限或屈服极限时发生的破坏和塑性流动，分析渐进破坏和失稳，特别适用于模拟大变形。材料通过单元和区域表示，根据计算对象的形状构成相应的网格。每个单元在外载和边界约束条件下，按照给出

2、的本构关系产生力学响应。FLAC 软件主要是为岩土工程稳定性分析开发的岩石力学计算程序，它包括了反映地质材料力学效应的特殊计算功能，能够计算地质类材料的高度非线性(包括应变硬化/软化)、不可逆剪切破坏和压密、粘弹(蠕变)、空隙介质的应力渗流耦合及动力学行为等。FLAC 提供了多种材料本构模型：各向同性弹性模型、横观各向同性弹性模型、摩尔-库仑塑性模型、应变硬化/软化塑性模型、德鲁克-普拉格塑性模型、遍布节理模型、双屈服塑性模型、霍克-布朗模型、空单元模型等。另外，程序设有界面单元，可以模拟断层、节理和摩擦边界的滑动、张开和闭和行为。支护结构，如砌衬、锚杆、支架等与围岩的相互作用也可以在FLAC

3、 中进行模拟。同时，用户可根据自己的需要在FLAC 中创建自己的本构模型，进行各种特殊修正和补充。FLAC 采用显式算法来获得模型全部运动方程的时间步长解，从而可以追踪材料的渐进破坏和跨落，这对研究开采的时间效应和空间效应是非常重要的。此外，程序允许输入多种材料类型，亦可在计算过程中改变某个局部的材料参数，增强了程序使用的灵活性，用来提供采动区域的跨落过程和开采中的充填过程。FLAC 具有强大的后处理功能，用户可以直接在屏幕上绘制图形，或以文件形式创建和输出打印多种形式的图形。使用者还可以根据需要，将若干个变量合并在同一幅图形中进行研究分析。FLAC3D 的计算过程如图1.1 所示，在每个时步

4、内，首先，根据高斯定律，由节点速度求出新的应变速度；根据本构方程计算出各时步内单元新的应力；然后，由运动方程计算新的结点速度和位移。在各个时步内进行循环计算，当计算达到平衡状态时，不平衡力趋近于0，如果有塑性流动产生，则不平衡力趋近某一数值。图1.1 FLAC3D循环求解过程二、节理化模型节理化模型是各向异性塑性模型，它包括了包含在摩尔库仑体内特殊方向上的弱面。根据应力状态、弱面走向以及模型体和弱面的材料特性的不同，屈服可能发生在模型体内，或者发生在弱面上，或者在两个部位同时发生。这种模型在FLAC中的实现方法是首先判别总体破坏，同时应用到和FLAC中摩尔库仑模型中相同的相关塑性修正，然后对更

5、新的应力在弱面上产生的破坏进行分析，同时对这些应力分别进行进一步的校正。弱面内的破坏准则存在于包含了拉应力路径的摩尔库仑屈服条件的局部形式中，与局部剪切流动法则不相关联而与局部拉应力流动法则相关联。图1.2显示了整体坐标系（x，y）和局部坐标系(x，y)下存在于摩尔库仑体内的软弱面。图1.2 整体坐标系中沿角方向的软弱面为了简化这部分的符号，定义ij，它对应于由各阶段总体破坏的塑性修正的应用引起的应力分量，这种总应力分解成局部应力后可表示为：式中：节理夹角（从x坐标轴逆时针方向算起）。依照约定，由表示弱面上切向引力分量的大小，相应的应变变量为，可以看到：有了这两个符号，弹性增量的法则的局部表达

6、式可以表示为：式中：1=K4G3 ，2=K2G3 ，上标e代表“弹性部分”。软弱面的破坏准则可以在图1.3所示的，平面内表示出来。图1.3 FLAC中的软弱面破坏准则将摩尔库仑破坏准则定义为0，从点A到点B的局部破坏包络线定义为： 将拉应力破坏准则定义为0，从点B到点C的局部破坏包络线定义为： 式中：，分别为弱面的摩擦角、粘聚力和抗拉强度。对于摩擦角不为零的弱面，抗拉强度的最大值定义如下：剪切和张拉势函数和对应于不相关联的流动法则，剪胀角对应于相关联的流动法则，它们分别表示为：破坏准则边界附近的流动法则用FLAC中摩尔库仑模型中所述的方法定义。这里，由函数h(，)=0表示(，)平面中0和0曲线

7、的对角线（图1.4），此函数形式为：图1.4 节理化模型中用于定义弱面流动法则的区域其中和是两个常量，定义如下：弹性假设和破坏准则不一致，分别由(，)平面中位于1区或2区（分别对应于h=0区域内-或+区域）。如果位于1区，说明平面内是剪切破坏，应用由势函数确定的流动法则，应力点回归到0的曲线上。如果位于2区，说明是局部拉应力破坏，应用由势函数确定的流动法则，应力点回归到0的曲线上。首先考虑平面内的剪切破坏，流动法则如下：其中上标表示弱面上同破坏有关的塑性部分，是待定参数，应用式，进行偏微分法后，上式可写为：最终，由于局部应力修正而分解到总体坐标轴下得到的弱面上剪切破坏的总应力修正式为：这些修正

8、被添加到应力分量上， 包括了总体破坏的应力修正，即这些应力修正会计算出该时步新的应力状态。现在考虑弱面上的拉应力破坏，这种情况下，流动法则的形式为：这里是待定的参数，应用式，进行偏微分后，上式可写为：应用如上的推理，可以得到：分解到整体坐标轴下以后，应力修正变成了： 在大应变模式下，为考虑由于刚体转动和变形引起的转动，对软弱面的方位角进行调整，修正值由一个区域内所有角度的平均值得出，其表达式如下：其中： 且以弧度形式表示。三、采场开采模拟分析鞍千矿业有限公司许东沟采场设计长1740m、宽280490m，最高海拔高度为+242.3m。设计边坡阶段高度12m，并段高度24m，阶段坡面角上盘（东帮）

9、65o，下盘（西帮）55o，运输平台、清扫平台、安全平台宽度分别取18m、8m、8m。最终边坡角：上盘（东帮）45o55o、下盘（西帮）40o50o。矿层走向145°165°，倾角很陡，一般都超过80°。露天采场现状最低开采标高为+96m，+144m水平以上的平台已经靠帮到界，形成边坡垂直高度为98m的高陡边坡，开采现状图见图1.5。局部边坡+168m至+240m设计边坡角为55°，由于开采初期没有设计，实际边坡境界+204m水平以上的部分已经超挖，所以现边坡角度为50°，小于设计的边坡55°的角度。表1.1 FLAC数值模拟计算参数

10、表岩性名称密度(g/cm3)粘聚力(kPa)内摩擦角(°)体积模量(GPa)剪切模量(GPa)抗拉强度(kPa)绿泥石英片岩2.80400329.005.40200绿泥石英片岩片理20030辽河群千枚岩3.35450349.505.70500含铁石英岩2.701000408.336.25500图1.5 开采现状图图1.6 开采最终设计图采场边坡分7步FLAC模拟开挖，模拟开采至72米、48米、 24米、 0米、-24米和-48米等四个水平的最大主应力1等值线图、最小主应力3等值线图、最大剪应变max增量等值线图、水平向位移ux等值线图、垂直向位移uy等值线图。 a开挖到72m时塑性区

11、 b开挖到48m时塑性区 c开挖到24m时塑性区 d开挖到0m时塑性区 e开挖到-24m时塑性区 f开挖到-48m时塑性区图1.7 矿区各开挖步塑性区变化情况（1）矿山分6步模拟开采后，上盘边坡192米96米水平间的岩体在台阶开挖的瞬间出现了岩体单元的塑性屈服和弱面破坏，开采扰动塑性区宽度约30米，边坡应力场重新调整后，屈服单元重回弹性状态；下盘边坡岩体则仅在台阶表层出现塑性屈服状态。 a b c d e f图1.8 各开挖步主应力变化情况 a b c d e f图1.9 各开挖步最小主应力变化图（2）由图1.8和图1.9各开挖部最大和最小主应力变化云图可知，受矿山开采影响，边坡岩体进行应力场

12、重新调整，坡面处应力完全释放，各剖面的最大主应力1和最小主应力3从坡面向坡体内逐渐增加，并过渡到原岩应力状态。坡面附近最小主应力3总体上为00.5MPa，个别平台坡顶处出现了拉应力，如168米水平、48米水平、24米水平和24米水平。 a b c d e f图1.10 各开挖步水平向位移变化图 a b c d e f图1.11 各开挖步垂直位移变化情况（3）图1.10和图1.11为各开挖步水平和垂直位移云图，分析水平向位移ux等值线图和垂直向位移uy等值线图表明，边坡岩体质点主要为弹性卸荷回弹变形，量级较小，位移较大区域位于边坡中下部，边坡中下部比较危险需要重点考虑一些防护措施。 a b c d e f图1.12 各开挖步剪应变增量（4）图1.12为各开挖部剪应变增量云图，从最大剪应力max（即主应力差值13）等值线