岩土工程数值计算方法.docx

《岩土工程数值计算方法.docx》由会员分享，可在线阅读，更多相关《岩土工程数值计算方法.docx（15页珍藏版）》请在三一办公上搜索。

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

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

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

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

5、一库仑模型中相同的相关塑性修正，然后对更新的应力在弱面上产生的破 坏进行分析，同时对这些应力分别进行进一步的校正。弱面内的破坏准则存在于 包含了拉应力路径的摩尔一库仑屈服条件的局部形式中，与局部剪切流动法则不 相关联而与局部拉应力流动法则相关联。图1. 2显示了整体坐标系（x，y）和局部坐标系（x，y）下存在于摩尔一 库仑体内的软弱面。y图1.2整体坐标系中沿0角方向的软弱面为了简化这部分的符号，定义。订它对应于由各阶段总体破坏的塑性修正 的应用引起的应力分量，这种总应力分解成局部应力后可表示为：b =b cos20 + 2b sin0 cos6+b sin20b b sin20 - 2b s

6、in0 cos0 +b cos20b，b3333b，=-（b -b ）sin0 cos6+b （cos2。一sin20）式中：0节理夹角（从x坐标轴逆时针方向算起）。依照约定，由T表示弱面上切向引力分量的大小，相应的应变变量为Y,可 以看到：b，12e，12有了这两个符号，弹性增量的法则的局部表达式可以表示为:b，=以&。+ 以(&e +&e)1111122233b，=以&。+ 以(&e + ef e )2213321122b，=以& e + 以(Aee + ee )3313321122T = 2丫 e式中：a：K+4G/3，气二K2G/3，上标e代表“弹性部分”。软弱面的破坏准则可以在图1.

7、3所示的b，t平面内表示出来。22Cj / tan%图1.3 FLAC中的软弱面破坏准则将摩尔一库仑破坏准则定义为fs =0,从点A到点B的局部破坏包络线定义为:将拉应力破坏准则定义为ft =0,从点B到点C的局部破坏包络线定义为:ft =bt 一。 22式中：七，七，。j 分别为弱面的摩擦角、粘聚力和抗拉强度。对于摩擦角不为零的弱面，抗拉强度的最大值定义如下：Cb t = jj ,max tan。j剪切和张拉势函数穿和g对应于不相关联的流动法则，剪胀角* j对应于相关联的流动法则，它们分别表示为：gs = 一r 一 b tan*22 jgt =b22破坏准则边界附近的流动法则用FLAC中摩尔

8、一库仑模型中所述的方法定 义。这里，由函数 b 22，T ）=0表示（b 22，T ）平面中f = 0和f =0曲线的对 角线（图1.4），此函数形式为：h =r r p a P （b b t）图1.4节理化模型中用于定义弱面流动法则的区域其中T P和a P是两个常量，定义如下:r p = C tan 8 b t弹性假设和破坏准则不一致，分别由L 22, T )平面中位于1区或2区(分 别对应于h=0区域内-或+区域)。如果位于1区，说明平面内是剪切破坏，应用 由势函数gs确定的流动法则，应力点回归到fs = 0的曲线上。如果位于2区，说明是局部拉应力破坏，应用由势函数M确定的流动法则，应力点

9、回归到ft =0 的曲线上。首先考虑平面内的剪切破坏，流动法则如下：Aefp = Ms 虫ii 8b11ep =Ms 空228b22kep = Ms 虫338b33其中上标p表示弱面上同破坏有关的塑性部分，M s是待定参数，应用式gs =-T-b tanv，进行偏微分法后，上式可写为：22 je p = 0 iie p = Ms tanw22je p = 0 33y p = Ms最终，由于局部应力修正而分解到总体坐标轴下得到的弱面上剪切破坏的总应力修正式为：b = 2br cos0 sin0 + b cos2 0 + b sin2 0b = 2Aor cos0 sin0 + b sin2 0

10、+ b cos2 0b = 2Ab，b =Ab (cos20 sin20) + (b b )sin0 cos0 12121122这些修正被添加到应力分量b上，bj包括了总体破坏的应力修正，即这些应力修正会计算出该时步新的应力状态。现在考虑弱面上的拉应力破坏，这种情况下，流动法则的形式为:* p我冬118b11e p d 室228b22e p d 也338b33 p孔也 8t这里人t是待定的参数，应用式g，=f ，进行偏微分后，上式可写为:22e e p 011e p 人 t22e p = 033y p = 0应用如上的推理，可以得到:b n =b +人t a11112b n =b +Xt a2

11、2221b n =b +人 ta33332.ft (b)Ki =23a1分解到整体坐标轴下以后，应力修正变成了:b= (bt b )(%cos20 + sin20)1122 a1b= (b t b )(_sin2 0 + cos2 0)2222 a1Ab = (b t b ) 23322 a1，a八八Ab =(。t S )(1 t)cosO sin20 1222 a1在大应变模式下，为考虑由于刚体转动和变形引起的转动，对软弱面的方位 角0进行调整，修正值A0由一个区域内所有角度的平均值得出，其表达式如下:ao = e +12其中:e = (e 一 e )sin0 cos0 + e (cos2

12、0 sin2 0)12112212=2 (& 汲)且A0以弧度形式表示。三、采场开采模拟分析鞍千矿业有限公司许东沟采场设计长1740m、宽280490m，最高海拔高度 为+242.3m。设计边坡阶段高度12m，并段高度24m，阶段坡面角上盘(东帮)65。， 下盘(西帮)55。，运输平台、清扫平台、安全平台宽度分别取18m、8m、8m。最 终边坡角：上盘(东帮)45。55。、下盘(西帮)40。50。矿层走向145165， 倾角很陡，一般都超过80。露天采场现状最低开采标高为+96m，+144m水平以上的平台已经靠帮到界， 形成边坡垂直高度为98m的高陡边坡，开采现状图见图1.5。局部边坡+168

13、m至 +240m设计边坡角为55，由于开采初期没有设计，实际边坡境界+204m水平以 上的部分已经超挖，所以现边坡角度为50。，小于设计的边坡55。的角度。表1.1 FLAC数值模拟计算参数表岩性名称密度(g/cm3粘聚力(kPa)内摩擦角( )体积模量(GPa)剪切模量(GPa)抗拉强度(kPa)绿泥石英片岩2.80400329.005.40200绿泥石英片岩片理20030辽河群千枚岩3.35450349.505.70500含铁石英岩2.701000408.336.25500FLA C3D 3.00Step 2930 Model Perspective 17:36:09 Fri May 16

14、 2014Center:X: 5.762e+002Y: 6.787e+001Z: 2.412e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 1.56Ang.: 22.500Block GroupItasca Consulting Group, Inc.图1.5开采现状图FLAC3D 3.00Step 9573 Model Perspective 18:57:47 Thu May 15 2014Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500

15、Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Block Group，89Itasca Consulting Group, Inc.Minneapolis, MN USA图1.6开采最终设计图采场边坡分7步FLAC模拟开挖，模拟开采至72米、48米、24米、0米、-24米和-48米等四个水平的最大主应力。等值线图、最小主应力。3等值线图、 最大剪应变Ymax增量等值线图、水平向位移七等值线图、垂直向位移u,等值线 图。四X7Itasca Consulting Group, Minneapolis, MN USAMin

16、neapolis, MN USAa开挖到72m时塑性区b开挖到48m时塑性区FLA C3D 3.00Step 6757 Model Perspective 18:33:50 Thu May 15 2014Center:Rotation:X: 3.460e+002X: 280.000Y: 2.728e+001Y: 340.000Z: 1.561e+002Z: 170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500Block State None shear-n shear-p shear-pshear-p tension-ptension-pFLAC3D 3.

17、00Step 7836 Model Perspective18:39:54 Thu May 15 2014Center:Rotation:X: 3.460e+002 X: 280.000Y: 2.728e+001 Y: 340.000Z: 1.561e+002 Z: 170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500ock State None shear-n shear-p shear-pshear-p tension-ptension-pMinneapolis, MN USAMinneapolis, MN USAc开挖到24m时塑性区d开挖到0m时塑

18、性区FLA C3D 3.00Step 8779 Model Perspective 18:45:34 Thu May 15 2014FLAC3D 3.00Step 9573 Model Perspective 18:55:42 Thu May 15 2014Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500ock StateNone shear-n shear-p shear-p shear-p ten

19、sion-p tension-pCenter:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500Minneapolis, MN USAe开挖到-24m时塑性区ock StateNone shear-n shear-p shear-p shear-p tension-p tension-pMinneapolis, MN USAf开挖到-48m时塑性区图1.7矿区各开挖步塑性区变化情况(1)矿山分6步模拟开采后，上盘边坡

20、192米96米水平间的岩体在台阶开挖的瞬间出现了岩体单元的塑性屈服和弱面破坏，开采扰动塑性区宽度约30 米，边坡应力场重新调整后，屈服单元重回弹性状态；下盘边坡岩体则仅在台阶 表层出现塑性屈服状态。FLA C3D 3,00Step 4296 Model Perspective 17:09:55 Thu May 15 2014FLA C3D 3.00Step 5594 Model Perspective 18:21:32 Thu May 15 2014Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Disf 2 501e+003Rotation:X:

21、290.000Y: 340.000Z: 170.000Mag :1Ang.: 22.500Contour of SMax(agfac = 0.000e+000 radient Calculation1 -2.5656e+006 to-2.5000e+006 to-2.0000e+006 to-1.5000e+006 to-1.0000e+006 to-5.0000e+005 to0.0000e+000 to2.5000e+0062.0000e+0061.5000e+0061.0000e+0065.0000e+0050.0000e+0004.1132e+005Interval = 5.0e+00

22、5Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Disf 2 501e+003Rotation: X: 280.000 Y: 340.000Z: 170.000Mag:0 8Ang.: 22.500Contour of SMaxagfac = 0.000e+000 radient Calculation -2.5071e+006 to -2.5000e+006 -2.5000e+006 to -2.0000e+006 -2.0000e+006 to -1.5000e+006 -1.5000e+006 to -1.0000e+006 -1.0000e

23、+006 to -5.0000e+005 -5.0000e+005 to 0.0000e+0000.0000e+000 to 4.1635e+005Interval = 5.0e+005Minneapolis, MN USAMinneapolis, MN USAFLA C3D 3.00Step 6757 Model Perspective 18:35:11 Thu May 15 2014FLA C3D 3.00Step 7836 Model Perspective 18:40:33 Thu May 15 2014Center:X: 3.460e+002Y: 2.728e+001Z: 1.561

24、e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500Contour of SMaxI agfac = 0.000e+000 radient Calculation -2.4638e+006 to ，-2.2500e+006 to -2.0000e+006 to -1.7500e+006 to -1.5000e+006 to -1.2500e+006 to -1.0000e+006 to -7.5000e+005 to -5.0000e+005 to -2.5000e+005 to 0.

25、0000e+000 to 2.5000e+005 to2.2500e+0062.0000e+0061.7500e+0061.5000e+0061.2500e+0061.0000e+0067.5000e+0055.0000e+0052.5000e+0050.0000e+0002.5000e+0054.1205e+005Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500Contour of SMaxI ag

26、fac = 0.000e+000 radient Calculation -2.4334e+006 to ，-2.2500e+006 to -2.0000e+006 to -1.7500e+006 to -1.5000e+006 to -1.2500e+006 to -1.0000e+006 to -7.5000e+005 to -5.0000e+005 to -2.5000e+005 to 0.0000e+000 to 2.5000e+005 to2.2500e+0062.0000e+0061.7500e+0061.5000e+0061.2500e+0061.0000e+0067.5000e

27、+0055.0000e+0052.5000e+0050.0000e+0002.5000e+0053.8717e+005Interval = 2.5e+005Interval = 2.5e+005Minneapolis, MN USAMinneapolis, MN USAFLA C3D 3.00Step 8779 Model Perspective 18:46:14 Thu May 15 2014Center:Rotation:X: 3.460e+002 X: 280.000Y: 2.728e+001 Y: 340.000Z: 1.561e+002 Z: 170.000Dist: 2.501e+

28、003 Mag.: 0.8Ang.: 22.500Contour of SMax.agfac = 0.000e+000 Iradient Calculation. -2.4181e+006 to -2.0000e+006-2.0000e+006 to -1.5000e+006-1.5000e+006 to -1.0000e+006-1.0000e+006 to -5.0000e+005-5.0000e+005 to 0.0000e+0000.0000e+000 to 5.0000e+0055.0000e+005 to 6.0766e+005Interval = 5.0e+005FLAC3D 3

29、.00Step 9573 Model Perspective 18:52:22 Thu May 15 2014Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500Contour of SMaxI agfac = 0.000e+000 radient Calculation -2.4086e+006 to -2.2500e+006 to -2.0000e+006 to -1.7500e+006 to -1.

30、5000e+006 to -1.2500e+006 to -1.0000e+006 to -7.5000e+005 to -5.0000e+005 to -2.5000e+005 to 0.0000e+000 to 2.5000e+005 to2.2500e+0062.0000e+0061.7500e+0061.5000e+0061.2500e+0061.0000e+0067.5000e+0055.0000e+0052.5000e+0050.0000e+0002.5000e+0054.7489e+005Interval = 2.5e+005Itasca Consulting Group, Mi

31、nneapolis, MN USAMinneapolis, MN USA图1.8各开挖步主应力变化情况FLA C3D 3.00Step 4296 Model Perspective 17:12:45 Thu May 15 2014Center:Rotation:X: 3.460e+002 X: 290.000Y: 2.728e+001 Y: 340.000Z: 1.561e+002 Z: 170.000Dist: 2.501e+003 Mag.: 1FLA C3D 3.00Step 5594 Model Perspective 18:23:00 Thu May 15 2014Ang.: 22.

32、500Contour of SMinI agfac = 0.000e+000 radient Calculation -8.0813e+006 to -8.0000e+0068.0000e+006 to7.0000e+006 to6.0000e+006 to5.0000e+006 to4.0000e+006 to3.0000e+006 to2.0000e+006 to7.0000e+0066.0000e+0065.0000e+0064.0000e+0063.0000e+0062.0000e+0061.0000e+0061.0000e+006 to -6.5882e+004Interval =

33、1.0e+006Center:Rotation:X: 3.460e+002X: 280.000Y: 2.728e+001Y: 340.000Z: 1.561e+002Z: 170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500Contour of SMinagfac = 0.000e+000 radient Calculation -8.0966e+006 to -8.0000e+006 -8.0000e+006 to -7.0000e+006 -7.0000e+006 to -6.0000e+006 -6.0000e+006 to -5.0000e+006

34、 -5.0000e+006 to -4.0000e+006 -4.0000e+006 to -3.0000e+006 -3.0000e+006 to -2.0000e+006 -2.0000e+006 to -1.0000e+006 -1.0000e+006 to -8.2339e+004Interval = 1.0e+006FLA C3D 3.00Step 6757 Model Perspective18:36:09 Thu May 15 2014Center:Rotation:X: 3.460e+002 X: 280.000Y: 2.728e+001 Y: 340.000Z: 1.561e

35、+002 Z: 170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500Contour of SMinI agfac = 0.000e+000 radient Calculation -8.1201e+006 to -8.0000e+006 to -7.0000e+006 to -6.0000e+006 to -5.0000e+006 to -4.0000e+006 to -3.0000e+006 to -2.0000e+006 to -1.0000e+006 to8.0000e+0067.0000e+0066.0000e+0065.0000e+0064.00

36、00e+0063.0000e+0062.0000e+0061.0000e+0061.3624e+005Interval = 1.0e+006FLAC3D 3.00Step 7836 Model Perspective 18:41:15 Thu May 15 2014Center:Rotation:X: 3.460e+002 X: 280.000Y: 2.728e+001 Y: 340.000Z: 1.561e+002 Z: 170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500Contour of SMinS agfac = 0.000e+000 radie

37、nt Calculation -8.1452e+006 to -8.0000e+006 to -7.0000e+006 to -6.0000e+006 to -5.0000e+006 to -4.0000e+006 to -3.0000e+006 to -2.0000e+006 to -1.0000e+006 to8.0000e+0067.0000e+0066.0000e+0065.0000e+0064.0000e+0063.0000e+0062.0000e+0061.0000e+0061.2770e+005Interval = 1.0e+006Minneapolis, MN USAMinne

38、apolis, MN USAFLA C3D 3.00Step 8779 Model Perspective 18:47:15 Thu May 15 2014FLAC3D 3.00Step 9573 Model Perspective 18:53:01 Thu May 15 2014Center:Rotation:X: 3.460e+002X:280.000Y: 2.728e+001Y:340.000Z: 1.561e+002Z:170.000Dist: 2.501e+003 Mag.: 0.8Ang.: 22.500Contour of SMinagfac = 0.000e+000 radie

39、nt Calculation -8.1732e+006 to -8.0000e+006 -8.0000e+006 to -7.0000e+006 -7.0000e+006 to -6.0000e+006 -6.0000e+006 to -5.0000e+006 -5.0000e+006 to -4.0000e+006 -4.0000e+006 to -3.0000e+006 -3.0000e+006 to -2.0000e+006 -2.0000e+006 to -1.0000e+006 -1.0000e+006 to -1.1393e+005Interval = 1.0e+006Center

40、:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.: 0.8Ang.: 22.500Contour of SMinS agfac = 0.000e+000 radient Calculation -8.1984e+006 to -8.0000e+006 to -7.0000e+006 to -6.0000e+006 to -5.0000e+006 to -4.0000e+006 to -3.0000e+006 to -2.0000e+006 to

41、-1.0000e+006 to8.0000e+0067.0000e+0066.0000e+0065.0000e+0064.0000e+0063.0000e+0062.0000e+0061.0000e+0069.2931e+004Interval = 1.0e+006Minneapolis, MN USAMinneapolis, MN USA图1.9各开挖步最小主应力变化图(2)由图1.8和图1.9各开挖部最大和最小主应力变化云图可知，受矿山开采影响边坡岩体进行应力场重新调整，坡面处应力完全释放，各剖面的最大主应力。和最小主应力。3从坡面向坡体内逐渐增加，并过渡到原岩应力状态。坡面附近最小主应力

42、。总体上为00.5MPa,个别平台坡顶处出现了拉应力，如168米水平、48米水平、24米水平和一24米水平。Itasca Consulting Group, Inc.Minneapolis, MN USAFLA C3D 3.00Step 5594 Model Perspective 18:26:55 Thu May 15 2014Center:X: 3.460e+002Y: 2.728e+001Z: 1.561e+002Dist: 2.501e+003Rotation:X: 280.000Y: 340.000Z: 170.000Mag.:0.8Ang.: 22.500Jontour of X-

43、DisplacementI agfac = 0.000e+0001 -6.9277e-003 to -5.0000e-003-5.0000e-003 to -2.5000e-003-2.5000e-003 to 0.0000e+000-0.0000e+000 to 2.5000e-003-2.5000e-003 to-5.0000e-003 to-7.5000e-003 to-1.0000e-002 to1.2500e-002 to1.5000e-002 to5.0000e-0037.5000e-0031.0000e-0021.2500e-0021.5000e-0021.5479e-002Interval = 2.5e-003Itasca Consulting Group, Inc. Minneapolis, MN USAFLA C3D 3.00Step 6757 Model Perspective 18:37:12 Thu May 15 2014FLAC3D 3.00Step 7836 Model Perspective 18:41:58 Thu May 15 2014Center:X: 3.460e+002Y: 2