FLAC3D原理..

1、2.2三维数值模拟方法及其原理2.2.1 FLAC3D工程分析软件特点FLAC3D是由美国ItascaConsultingGroup,Inc.为地质工程应用而开发的连续介质显式有限差分计算机软件。FLAC即FastLagrangianAnalysisofContinua的缩写。该软件主要适用于模拟计算岩土体材料的力学行为及岩土材料达到屈服极限后产生的塑性流动，对大变形情况应用效果更好。FLAC3D程序在数学上采用的是快速拉格朗日方法,基于显式差分来获得模型全部运动方程和本构方程的步长解，其本构方程由基本应力应变定义及虎克定律导出，运动平衡方程则直接应用了柯西运动方程，该方程由牛顿运动定律导出。

2、计算模型一般是由若干不同形状的三维单元体组成，也即剖分的空间单元网络区，计算中又将每个单元体进一步划分成由四个节点构成的四面体，四面体的应力应变只通过四个节点向其它四面体传递，进而传递到其它单元体。当对某一节点施加荷载后，在某一个微小的时间段内，作用于该点的荷载只对周围的若干节点(相邻节点)有影响。利用运动方程，根据单元节点的速度变化和时间，可计算出单元之间的相对位移，进而求出单元应变，再利用单元模型的本构方程，可求出单元应力。在计算应变过程中，利用高斯积分理论，将三维问题转化为二维问题而使其简单化。在运动方程中，还充分考虑了岩土体所具有的粘滞性，将其视作阻尼附加于方程中。FLAC3D具有一个

3、功能强大的网格生成器，有12种基本形状的单元体可供选择，利用这12种基本单元体，几乎可以构成任何形状的空间立体模型。FLAC3D主要是为地质工程应用而开发的岩土体力学数值评价计算程序，自身设计有九种材料本构模型：(1)空模型(NullModel)(2)弹性各向同性材料模型(Elastic,IsotropicModel)(3)弹性各向异性材料模型(Elastic,anisotropicModel)(4)彳惠拉克-普拉格弹塑性材料模型(Drucker-PragerModel)(5)莫尔-库伦弹塑性材料模型(Mohr-CoulombModel)(6)应变硬化、软化弹塑性材料模型(Strain-Har

4、dening/SofteningMohr-CoulombModel)(7)多节理裂隙材料模型(Ubiquitous-JointModel)(8)双曲型应变硬化、软化多节理裂隙材料模型(BilinearStrain-Hardening/SofteningUbiquitous-JointModel)(9)修正的Cam粘土材料模型(ModifiedCam-clayModel)除上述本构模型之外，FLAC3D还可进行动力学问题、水力学问题、热力学问题等的数值模拟。在边界条件及初始条件的考虑上，FLAC3D软件十分灵活方便，可在数值计算过程中随时调整边界条件和初始条件。FLAC3D具有强大的后处理功能，

5、用户可以直接在屏幕上绘制或以文件形式创建或输出打印多种形式的图形、文字，用户还可根据各自的需要，将若干个变量合并在同一幅图形中进行研究分析。FLAC3D软件还可对各种开挖工程或施加支护工程等进行数值仿真模拟，软件自身设计有锚杆、锚索、衬砌、支架等结构元素，可以直接模拟这些支护于围岩(土)体的相互作用。FLAC3D拥有可以自行设计的FISH语言，用户可根据自身需求，自己设计材料的本构模型、屈服准则、支护方案、复杂形状的开挖方式等工作。特别注意的是，岩石是一种脆性材料，当外荷载达到岩石强度后，材料发生断裂破坏，产生弱化现象，应属于弹塑性体。在FLAC3D中，一般对于弹塑性材料，判断其破坏与否的基本

6、准则有两个，即Drucker-Prager准则和Mohr-Coulomb准则。根据室内岩石力学性质试验结果，其典型应力应变曲线反映出岩体破坏包络线符合莫尔一库伦屈服准则，故本次建立的本构力学模型选择莫尔一库伦弹塑性材料模型为宜。2.2.2FLAC3D分析计算原理计算所采用的数学模型是根据弹塑性理论的基本原理(应变定义、运动定律、能量守衡定律、平衡方程及理想材料的连续性方程等)而建立的。2.2.2.1 基本约定在数学及数值模型的表达式中，符号有一定的约定含义，一般A表示张量，Aij表小张量A的（i,j）分量，a表小矢量，ai表小矢量a的i分量,表小对Xi的偏导数。Xi,Ui,vi和dvi/dt,

7、（i=1,3）分别表示一点的位置矢量分量、位移矢量分量、速度矢量分量和加速度矢量分量。2.2,2.2数学模型（一）柯西（Cauchy）应力张量与柯西公式对于一个具有体积V的封闭曲面s的物体，在其上取一表面元素s,这个表面元素的单位外法向矢量为n,在某一时刻t,在表面元素对于连续介质中一点，作用着对称的应力张量ij,根据s上作用有力巳则极限1sm0dPds称为表面力。若用ti表示T的分量，则在三维直角坐标系中可有关系式tiniij(1)这个关系式称为柯西公式，其中,（二）应变速率和旋转速率ij称为柯西应力张量。如果介质质点具有运动速度矢量v,则在一个无限小的时间dt内，介质会产生一个由vidt决

8、定的无限小应变，对应的应变速率分量ij为ij1产X）2XjXi而其旋转速率分量ij为ijJ/上）2Xjx(2)(3)（三）运动及平衡方程根据牛顿运动定律与柯西应力原理，如果质点作用着应力具有速度vi,则在无限小时间段dt内，它们之间的关系为ij,dvibiXjdtij与体力bi,且（4）式中，为质点密度。（4）式称为柯西运动方程。当质点的加速度为零时，上式变为静力平衡方程bi0（5）Xj（四）本构方程上述（4）式与（5）式组成的方程组中含有9个方程，15个未知量，其中12个是应力与应变速率分量，3个是速度分量。其余6个关系式则由本构方程提供，本构方程一般具有如下形式ijHj（0,0,）（6）式

9、中，为应力变化速率，H表示一个特定的函数关系，为与荷载历史有关的参数。2.2.2.3数值模型flac3D的数值剖分网格在计算中是按照四面体进行的，四面体的节点也既是网格剖分的节点，因此，每个计算单元有4个面和4个节点（见图2-1）图2-1四面体单元示意图（一）空间微分的有限差分逼近对于一个计算单元，若内部各质点速度为一连续的矢量场（GausS）积分原理有VidvvinjdsVXjSjv,则根据高斯式中，n为外法向单位矢量场。由于单元体的应变速率是连续的，因此可以近似认为速度是线性变化的,则（7）式可用下面的求和公式近似逼近V*4Vixjf1(f)n(f)q(f)njS(8)式中，v为单元体体积

10、,f=1,4,Vi为面平均速度的(9)Vixj13i4lVi1f4(f)nj1,flS(f)(10)S为单元体某一面的面积，f为单元体面数,分量。由于速度场是线性的，则有-(f)1iVi-Vi311,if式中，l为单元体节点数，l=1,4;vi1为l节点的i速度分量将（9）式代入（8）式可得(11)则将（10）式两边除以Vixj13ViV,并将(4vln(l)S111）式代入可得(12)根据正交原理有4（f）Q（f）njSf1因此有146Vi1/i(l)i(VinjVjZ)ni)S(13)（二）运动方程的节点公式根据前面对质点运动方程的讨论，对于连续介质，当处于平衡状态时，其平衡方程为一LBi

11、0（14）xjBi（bi*（15）dt介质可由若干个作用着体力B,各自产生一定变形的四面体组成，设节点力为fn,n=1,4,利用虚功原理，假定单元体节点具有速度vn,内部具有变形速率,则由节点力fn和体力B所做的外功率与ij所做的内功率应当相外功率E可用下式表小EvnfinVvBdV(16)n1而内功率I为IVjjdV(17)由(13)式，对于恒定应变速率单元体可有1 4I1(vl"vjn)Sl(18)6 i1由于应力张量是对称的，定义矢量TlTi1ns1(19)可得14I-vlTil(20)3 l1将(15)式代入(16)式得4EvfinEbEI(21)n1式中，Eb和E1分别是体

12、力bi和惯性力产生的外功率，对于单元体内恒定的体力bi,Eb可写成EbhVvdV(22)而E1可以写成EIvi皿dV(23)Vdt如前所述，单元体内部速度场以线性变化，为方便描述，选取单元体质心为原点，x、x2、*3为坐标轴，可有4vv：Nn(24)n1(25)(26)式中，Nn(n=1,4)是具有如下形式的线性函数nnnnnNC0CiXiC2X2C3X3Cnnncnc0、c1、c2、c3(n=1,4)则由下面方程决止NnX1j,X2j,X3jnjnj为克罗内克尔(Kronecker)记号，根据质心定义可得用克莱默Ebbivinc0Vn1(Cramer's)法则解得EbEI在稳定状态下

13、,nC04,则nbiVVi，"4nVin1nVifinNn>biVTVn"Mv内功率与外功率必定相等，则有fnTinbiVi34上式最后一项，如果单元体内ndViNnLdV-VdtndViNnidVVdt恒定，则根据质心定义有V(dVi)n4dt若将V/4看作是假定的节点质量mn,则上式变为ndVinLdVdVinfinbiV4n,dvi、nmE根据牛顿定律有l1,nn式中，Fi<l>为节点l所受力的分量，M<l>为节点l的质量,(28)(29)(30)(31)(32)(33)(34)(35)(也)l为节点dtl在Fi<l>作用下产

14、生白加速度，nn为包含在全部连续介质中单元节点总数。设质点质量Mlml,不平衡力F<>为lFil6等fil(36)34(三)时间微分的显式有限差分对于上面(35)式，可写成(37)dvil1l()Fil1,nndtMlFLAC3D计算中，对于上式的逼近是在一个足够小t时间段内，采用中心(38)差分格式，即上式变为Vi(t)Vi(t)-Fi22Ml同理，节点的坐标差分公式与位移差分公式分别为Xil(tt)Uil(tt)iltXi(t)tvi(t)2uJ(t)tvj(t')(39)(40)(四)本构方程的增量形式FLAC3D中，假定在t表小为时间段内，速度是不变的，本构方程的增

15、量形式可*ijHij(41)jtj则ij为t时间段内应变增量。应力增量由下式确定(42)式中,ijijiC由下式定义(43)(ikkjikkj)t(44)14l6nVinjllolvjniS(45)2.2.2.4本构模型(一)弹塑性理论的增量关系(46)一般来讲，弹塑性体破坏的判断准则具有如下形式f(-i)0式中，f为一特定的函数关系，一为n维应力矢量i分量。弹塑性体的变形为弹性变形与塑性变形p的总和，即ep_i_i_i弹性应变与应力之间的关系为_e_iSi_ni1,n其中，S为线性函数。塑性变形可用下式表示p_g(47)(48)(49)式中，为一常数，g为（_n）的某一函数关系。将（47）、

16、（49）式代入（48）式得g_iSi_nSi（-）（50）一n当f（_n）为线性函数时,（51）式可写成_*_f_nf_n0其中，f.f.f0no由于f（_n）0,再将（50）式代入（52）式可得_*_*_gfSn_nf?()0-n若记则有因此有N_i_i_i1s_i_iSi_n_If_nfSn(_n)(51)(52)(53)(54)InfSng/_nf0n对于新应力_n_n,仍有则新应力为(二)莫尔库伦(Mohr-Coulomb)模型本构关系Hooke'slaw),应力应变的关系为(57)本次计算的本构模型选择的是莫尔一库伦弹塑性材料模型，弹塑性体可产生弹性及塑性两部分变形，根据虎克

17、定律(E其中，E为刚度矩阵。对于弹性变形，应力增量可由下式确定式中,e1e2e3e2e1e1e3e3e2(58)2可由剪切模量与体积模量得出4G32G3(58)式还可写成S2S3e1,e1,e1,e2,e2,e2,e3e3e3e1e2e3e2e1e1e3e3e2(59)(三)破坏准则与流动法则若123,则莫尔一库伦破坏准则具有形式fs13Nfs=0,2cN其中(60)对于拉张破坏，有ftft=0的形式，其中t(61)上两式中，为内摩擦角，C为内聚力,N1sin()t为抗拉强度,1sin()对于材料的抗拉强度:其值不可能超过3(见图2-2),(62)t的最大值为tmaxctg()对于塑性流动状态

18、，设g,与gt分别为剪切破坏和拉张破坏所对应的函数关系,有gsi3Ntg3式中，为膨胀角，(64)(65)(66)N流动准则可定义为函数关系h（1，3）0,在1、3平面内将fs0和图2-3莫尔一库伦模型塑性流分区ft0以上部分分为两个区域（见图2-3）,其函数表达式如下.tPPh3a1其中，ap与p由下式定义N2NPtN2cN(68)从图2-3中可见，若某点的应力状态在1区（fs0,h0）时，则发生剪切破坏，在2区（fs0,h0）时，则发生拉张破坏。（四）塑性阶段应力修正对于剪切破坏，由sg(64)式可得(69)3将（59）式中的gs/1、gs/2和gs/3代Si2NS2(70)1N再由(55)、(56)N1N2N