有限差分法(1)FLAC2D

有限差分法（1）FLAC2D1.FLAC简介1.1FLAC特点FLAC——FastLagrangianAnalysisofContinuaFLAC建立在拉格朗日算法基础上，采用有限差分显式算法来获得模型全部运动方程(包括内变量)的时间步长解，从而可以追踪材料的渐进破坏和垮落，这对研究采矿工程设计是非常重要的。FLAC适用模拟计算岩土材料力学行为，特别适合模拟大变形和扭曲，包括材料的高度非线性(应变硬化/软化)、不可逆剪切破坏和压密、粘弹(蠕变)、孔隙介质的应力—渗流耦合、热—力耦合以及动力学问题等。1.2岩体力学本构模型(1)各向同性弹性材料模型；(2)横观各向同性弹性材料模型；(3)莫尔-库仑弹塑材料模型；(4)应变软化/硬化塑性材料模型；(5)双屈服塑性材料模型；(6)遍布节理材料模型；(7)空单元模型，可用来模拟岩体开挖和开采。1-3莫尔-库仑准则与应变软化加载路径卸载路径压力P体积应变eV岩石充填体岩石应变软化/硬化模型LNNTkkM单元SPA侧B侧sn断层断层模型1.3支护构件与本构模型(1)杆单元模型；(2)梁模型；(3)桩模型；(4)群支架模型。预应力砂浆锚杆的模拟体系cs,mmm锚杆砂浆体锚固结点砂浆剪切刚度砂浆粘结刚度砂浆粘聚力,摩擦力KbSbKbKbSbcs,自由段锚固段FLAC举例(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({:..式中，δij——Kronecker记号；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.模拟大变形明显需要更多的机时。动态衰减•在动态衰减过程中，结点按照牛顿定律运动。结点的加速度与非平衡力成正比。该算法确定这组位移将导致系统趋于平衡，或表征失稳状态。•关于动态衰减的两个重点：①选择时步②阻尼作用时步为满足数值稳定，时步必须满足：：pCxtmin式中，Cp与1/mgp成正比。对于静态分析，结点质量按比例度量，致使局部临界时步等于（t=1），这样提供了收敛的最优时步。然后，对结点惯性质量进行调整来实现稳定性条件：注意：重力质量不受影响。AxGKmpg3)3/4(2max阻尼•速度-比例阻尼导致影响计算的体积力。•FLAC中采用了局部阻尼——结点处的阻尼力与不平衡力大小成正比，正负号的规定确保振动状态被阻滞。局部阻尼•阻尼力Fd为：mtuFFuiiii)(sgn||•在FLAC中，通过监测不平衡力比率（不平衡力Fi与施加力Fm之比）来确定静力状态•当Fi/Fm0.001时，则（计算机缺省）认为模型处于平衡状态。)sgn(iiduFF式中，Fi是不平衡力•在运动方程中引进阻尼力如图，一维杆件用数个等尺寸的有限差分网格划分，杆件的密度为，杨氏模量为E。1.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）设计计算模型的尺寸（2）规划计算网格数目和分布（3）安排工程对象（开挖、支护等）（4）给出材料的力学参数（5）确定边界条件（6）计算模拟3.1建模FLAC基本称谓单元编号结点编号3.2网格生成Gridi,j例如：grid30,20xyi=30j=203.3网格规划：Genx1,y1x2,y2x3,y3x4,y4例如：Gen0,00,1010,2020,01(x1,y1)2(x2,y2)3(x3,y3)4(x4,y4)3.4分区规划网格。例如：GenxI1,yI1xI2,yI2xI3,yI3xI4,yI4i=1,10j=1,21（I区）GenxII1,yII1xII2,yII2xII3,yII3xII4,yII4i=10,20j=1,21(II区）IIIi=10i=1,j=1i=21,j=1i=1,j=21i=10,j=21i=21,j=213.5特殊形状的网格（1）圆形gencirclexc,ycradxc,ycrad3.5特殊形状的网格（2）弧线genarcxc,ycxb,ybthetaxc,ycxb,yb3.5特殊形状的网格（3）直线genlinex1,y1x2,y2x1,y1x2,y2（4）任意形状tab1x1,y1,x2,y2,,xn,yn,x1,y1gentab1(x1,y1)(x2,y2)(x3,y3)(x4,y4)(x5,y5)(x6,y6)3.6赋给单元材料性质mode(弹性)propd1800e-6bu12.5sh5.77i=1,20j=1,10propd2400e-6bu1250sh577i=1,20j=11,20modm（弹塑性Mohr-Coulumb准则）propd1800e-6bu12.5sh5.77c0fri20ten0.015regi,j*特别提示)mod()1(2)mod()21(3ulusshearEGulusbulkEK（1）量纲统一（2）参数换算SIImperialLengthmmmcmftInDensitykg/m3103kg/m3106kg/m3106g/cm3slugs/ft3snails/in3ForceNkNMNMdynesIbfIbfStressPakPaMPabarIbf/ft2psiGravitym/sec2m/sec2m/sec2cm/s2ft/sec2in/sec2StiffnessPa/mkPa/mMPa/mBar/cmIbf/ft3Ib/in33.7赋给模型边界条件（1）固定边界（结点）Fixxi=1,j=1,21Fixyi=1,21j=1Fixxi=1j=1Fixy3.7赋给模型边界条件（2）施加边界力(结点)applyyf-10i=1,21j=21或applysyy-10i=1,21j=21applyxf-5i=21,j=1,21或applysxx-5i=21,j=1,21i=21j=213.7赋给模型边界条件（3）赋单元内应力（单元）inisxx-10i=1,20j=1,20inisyy-5var04i=1,21j=1,21i=21j=213.8计算Setgrav9.81SetlargeStep1000Savetest.sav3.9结果显示Plotgrid显示网格Plotbo显示边界Plotplas显示塑性区Plotsig1fi显示最大主应力1Plotsig2fi显示最小主应力2Plotsdiffi显示主应力差(1-2)Plotstr显示主应力矢量场Plotxdisfi显示X方向位移Plotydisfi显示Y方向位移Plotdisp显示位移矢量场3.10保存与调用结果命令Calltest.txt(或catest.dat)调用数据Savetest.sav保存结果New重新开始Resttest.sav调用结果Quit退出程序4.FLAC解题技巧4.1模型尺寸4.2模拟开挖Modnui=6,15j=5,12(或regioni,j)i=6,15j=5,124.FLAC解题技巧4.3模