FLAC采矿工程数值计算方法

采矿工程数值计算方法——FLAC建模技巧与工程应用采矿工程数值计算方法——FLAC建模技巧与工程应用学术讲座主讲：王金安教授北京科技大学数值方法分类数值方法分类ƒ连续介质方法¾有限元、有限差分、边界元、无单元ƒ非连续介质方法¾离散元、颗粒元、流形元提纲1.FLAC简介2.FLAC理论基础3.FLAC应用实例4.FLAC建模方法5.FLAC运用技巧1.FLAC简介1.1FLAC特点ƒFLAC——FastLagrangianAnalysisofContinuaƒFLAC建立在拉格朗日算法基础上，采用有限差分显式算法来获得模型全部运动方程(包括内变量)的时间步长解，从而可以追踪材料的渐进破坏和垮落，这对研究采矿工程设计是非常重要的。ƒFLAC适用模拟计算岩土材料力学行为，特别适合模拟大变形和扭曲，包括材料的高度非线性(应变硬化/软化)、不可逆剪切破坏和压密、粘弹(蠕变)、孔隙介质的应力—渗流耦合、热—力耦合以及动力学问题等。(1)各向同性弹性材料模型；(2)横观各向同性弹性材料模型；(3)莫尔-库仑弹塑材料模型；(4)应变软化/硬化塑性材料模型；(5)双屈服塑性材料模型；(6)遍布节理材料模型；(7)空单元模型，可用来模拟岩体开挖和开采。1.2岩体力学本构模型εσ1-σ3莫尔-库仑准则与应变软化加载路径卸载路径压力P体积应变eV岩石充填体岩石应变软化/硬化模型LNNTkkM单元SPA侧B侧sn断层断层模型1.3支护构件与本构模型(1)杆单元模型；(2)梁模型；(3)桩模型；(4)群支架模型。预应力砂浆锚杆的模拟体系cs,ϕmmm锚杆砂浆体锚固结点砂浆剪切刚度砂浆粘结刚度砂浆粘聚力,摩擦力KbSbKbKbSbcs,ϕ自由段锚固段土体砂浆体锚杆体σm=p’(a)(b)(c)FtutFLAC举例(a)巴西试验(b)边坡滑移（间接拉伸）P2.FLAC理论基础2.1有限差分法02220032⎟⎟⎠⎞⎜⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂−=xfhxfhff02220012⎟⎟⎠⎞⎜⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂+=xfhxfhffhffxf2310−=⎟⎠⎞⎜⎝⎛∂∂20310222hfffxf−+=⎟⎟⎠⎞⎜⎜⎝⎛∂∂hffyf2420−=⎟⎟⎠⎞⎜⎜⎝⎛∂∂20420222hfffyf−+=⎟⎟⎠⎞⎜⎜⎝⎛∂∂xyh0123456789101112h2.2简单力学问题分析mF(t)•牛顿运动定律dtudmamF&⋅=⋅=•对于一个连续体，可写成广义形式：ijijigxdtudρ+∂σ∂=ρ&式中,ρ——体密度,xj——坐标矢量(x,y)σij——应力张量分量gi——重力加速度u,u,u&&&2.3应力—应变关系除了运动定律，连续介质必须服从本构关系——应力与应变之间的关系•对于弹性材料：•一般地，本构关系可表述为：式中teGeGKijkkijijijΔ+−+=}2)32({:..δσσ式中，δ——Kronecker记号；ijΔt——时间步；G，K——分别是剪切模量和体积模量。),,(:.keMijijijσσ=][21...ijjiijxuxue∂∂+∂∂=2.4有限差分基本格式•在有限差分法中,每个在前一个公式中导出的量（运动量或应力应变）由所在网格特定位置处的相关变量代数式所代替。•代数式是完全显式的；在代数式右侧的所有量都是已知的。在一个计算时步，FLAC网格每个元素（域或结点）与相邻单元表现出在物理上是隔绝的。选取的时步非常小，乃至在此时步间隔内实际信息不能从一个单元传递到另一个单元（事实上，所有材料都有传播信息的某种昀大速度）。•计算循环的基本要求输入应变输出应力非线性定律应变应力2.5基本显式计算循环平衡方程(运动方程)应力-应变关系(本构方程)对于所有结点LnFjiji对于所有单元σ=新的应力应变速率结点力高斯定理速度ijijigxdtudρσρ+∂∂=&例如,弹性teGeGKijkkijijijΔ+−+=}2)32({:..δσσFLAC的网格内部由三角形构成，三角形组合成四边形单元。推导差分方程的方式如下：单元叠加三角形单元结点力向量速度向量abababui(b)ui(a)ΔSFini(1)ni(2)s(1)s(2)高斯定理：∫∫∂∂=SAiidAxffdSn用于推导任意形状单元的有限差分公式。)b(iu&结点速度ba)a(iu&结点速度ΔS对于多边形，公式为：∑≈∂∂SiiSnfAxfΔ1下式用于计算应变增量,Δeij:()txuxueSnuuAxuijjiijSjbiaijiΔΔΔ⎥⎥⎦⎤⎢⎢⎣⎡∂∂+∂∂=+≈∂∂∑&&&&&2121)()(一旦计算出全部应力，可以从作用每个三角形边界上产生的牵引力计算得到结点力。例如：然后，用“传统”的中间差分公式获得新的速度和位移：(…大变形模式))(21)2()2()1()1(SnSnFjjiji+=σ∑+=−+mtFuutittittiΔΔΔ)()2/(.)(.tuxxttitittiΔΔΔ)2/(.)(.)(.+++=Fini(1)ni(2)s(1)s(2)在时间域上的解算方法所有单元:{}{}()σΔ=Δ,ufF(非线性定律)所有结点:{}tmFuΔ⎭⎬⎫⎩⎨⎧Σ=Δ&重复n个时步时步内不作迭代信息在每个时步内不在单元之间传播假设(u)固定假设(F)固定改正，若pminCxtΔΔp-波速度位移u力FΔxF应力σu数值结点显式隐式{}[]{}uKF=单元[]{}[]{}{}FuKumΣ=+&&结点每个时步内解算整套方程如果存在非线性，在时步内迭方法比较显式,时间-进程隐式,静态1.可以遵循非线性定律而无需内部迭代,因为本构计算中位移被“冻结”。2.对同样问题，计算时间增加N3/2。3.物理不稳定不会引起数值不稳定。4.因为无需存储矩阵，用少量内存可以模拟大型问题。5.大应变、大位移和转动模拟无需额外机时。1.整个迭代过程需要遵循非线性定律。2.解算时间增加N2甚至N3。3.模拟物理不稳定性困难。4.需要大内存，或大容量硬盘存储。5.模拟大变形明显需要更多的机时。ƒ如图，一维杆件用数个等尺寸的有限差分网格划分，杆件的密度为ρ，杨氏模量为E。2.6FLAC计算流程举例①②③④i-1ii+112345i-1ii+1结点单元速度、位移应力xxuExxx∂∂=σƒ对于固体材料，微分形式的本构方程为：（1）ƒ运动方程(或平衡方程)为：xtuxxx∂∂=∂∂σρ22ƒ假设杆件无侧向约束。对于单元i，与(1)式对应的中间有限差分公式为：（2）xtutuEtixixixxΔ−=+)()()(1σ（3）式中，(t)表示其变量是在时刻t确定的，上标i表示单元或结点编号。)}()({1)}2()2({1..ttxttuttutixxixxixix−−Δ=Δ−−Δ+Δσσρ（4）ƒ对结点i的有限差分形式运动方程为：)}()({)2()2(1..ttxtttuttuixxixxixix−−ΔΔ+Δ−=Δ+σσρ或写成（5）ƒ在显式算法中，所有有限差分方程右端的值均是已知的。因此，必须先用(3)式算出所有单元的应力，然后再由(5)式和(6)式计算所有结点的速度和位移。ƒ在概念上，这个过程等价于对变量值“同时”更新，而不是像其它方法那样，方程右端混存有“新”、“旧”值，对变量值“依次”更新。积分后得出位移：tttututtuixixixΔΔ++=Δ+)2()()2(.（6）3.FLAC应用实例（1）老虎台矿开采诱发矿震的力学机理和震源机制的模拟计算分析$•地质条件F25F26F1F7-1F18•E5200剖面图•老虎台矿开采历史老虎台矿自1907年开始开采，至今已有近百年的开采历史。•矿震事件统计1988年1月至2000年5月，随着老虎台矿开采深度增大和向断裂构造逼近，矿震频率和震级都呈上升趋势，平均每月发生矿震52.2次，远远超出了抚顺地区天然地震的数量，昀大震级达到3.6ML。矿震震级(ML)1∼1.51.6∼2.02.1∼2.52.6∼3.03.1∼3.6合计矿震数目(次)5721167529863197776•抚顺矿震震级、频度与时间特征ƒ抚顺市矿震有感范围ƒ抚顺市矿震有感范围•计算模型F7-1F18F1F16F1-AF25F263800m3000m1100mxyz1999-2000E5200E6000~1992~20001999-2000E5200E6000不同开采时期断层附近岩体的动力学响应(1992-2000)00.10.20.30.40.50.60.70.80.97.88.08.28.48.68.89.09.29.49.6StepShearDisp.,mm78001-1F25断层露头(1992-2000)00.511.522.537.88.08.28.48.68.89.09.29.49.6StepShearDisp.,mm78001-1F25断层-200m（井田西部）(1992-2000)051015207.88.08.28.48.68.89.09.29.49.6StepShearDisp.,mm78001-1(1992-2000)0510152025307.88.08.28.48.68.89.09.29.49.6StepShearDisp.,mm78001-1F25断层-400mF25断层—煤层（-595m）（井田西部）3g包络域不同开采深度的3g包络域（2）程潮铁矿主溜井特大塌方治理N24m措施井1#主溜井2#主溜井西风井1#、2#主溜井的布置图托斗锚固法治理主溜井塌方技术方案主视图“托斗法”治理主溜井塌方局部剖视图加固前、后主溜井围岩位移矢量场（卸矿水平剖面）加固前加固后加固前、后主溜井围岩位移矢量场（卸矿井筒剖面）加固前加固后（3）鹤壁四矿村庄下厚煤层特殊开采2602N2604F726144F206杨吕寨26064F1054F204260826102612剖面图杨吕寨750m2000m280m700m420mF7430m706m300m4F2044F2064F1055163.55154中央条带两翼长壁开采方案80m150m150m三维计算模型图650m1800m1500m4F1054F2064F204F7中央条带两翼长壁工作面开采后地表下沉0200400600800100012001400x(m)020040060080010001200140016001800y(m)-7000-6500-6000-5500-5000-4500-4000-3500-3000-2500-2000-1500-1000-5000W(max)=-6507.68(mm)mm中央条带两翼长壁工作面开采后地表水平变形200400600800100012001400x(m)2004006008001000120014001600y(m)-12-10-8-6-4-20246810Ex=(-11.8236,11.3979)(mm/m)mm/m(a)X-方向200400600800100012001400x(m)2004006008001000120014001600y(m)-14.0-12.0-10.0-8.0-6.0-4.0-2.00.02.04.06.08.010.012.0Ey=(-12.8346,12.839)(mm/m)mm/m(b)Y-方向(4)综放采场围岩应力分布与稳定性研究600m500m314.17m三维模型走向剖面500m314.17m600m314.17m倾向剖面综放工作面煤层内垂直应力分布三维视图沿煤层走向综放工作面中部围岩主应力场采空区老顶岩层140.8m沿煤层倾斜剖面工作面前方150m围岩主应力场沿煤层倾斜剖面工作面前方15m围岩主应力场沿煤层倾斜剖面工作面处围岩主应力场沿煤层倾斜剖面工作面后方200m围岩主应力场(5)大倾角特厚煤层开采地表塌陷与地表水入渗模拟研究四层煤四层煤四层煤二层煤二层煤F17F16F15F3黄土900m400m+1700m+1800m+1900m+1500m+1600m地表塌陷区地表开裂山体断裂地表抽冒山体错断不同开采水平受采动岩体内的渗流迹线图(a)