FLAC3D原理..

《FLAC3D原理..》由会员分享，可在线阅读，更多相关《FLAC3D原理..（16页珍藏版）》请在装配图网上搜索。

1、2.2 三维数值模拟方法及其原理2.2.1FLAC3D工程分析软件特点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-Hard

4、ening/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.2 FLAC3D分析计算原理计算所采用的数学模型是根据弹塑性理论的基本原理（应变定义、运动定律、能量守衡定律、平衡方程及理想材料的连续性方程等）而建立的。2.2.2.1基本约定在数学及数值模型的表达式中，符号有一定的约定含义，一般A表示张量，A.表示张量A的（ij）分量，a表示矢量，a.表示矢量a的i分量，ai表示a对x.的偏导数。ix.，u.，v.和dv/dt，（i

7、=l,3）分别表示一点的位置矢量分量、位移矢量分iiii量、速度矢量分量和加速度矢量分量。2.2.2.2 数学模型（一）柯西（Cauchy）应力张量与柯西公式对于一个具有体积V的封闭曲面s的物体，在其上取一表面元素As,这个表面元素的单位外法向矢量为n，在某一时刻t，在表面元素对于连续介质中一点，作用着对称的应力张量6.,根据As上作用有力AP，.j则极限TyAPdPT=lim=AtoAsds称为表面力。若用t.表示T的分量，则在三维直角坐标系中可有关系式.t=nb（1）.j这个关系式称为柯西公式，其中，b.j称为柯西应力张量。.j二）应变速率和旋转速率如果介质质点具有运动速度矢量v，则在一个

8、无限小的时间dt内，介质会产生一个由v.dt决定的无限小应变，对应的应变速率分量g.为.j2)3)1dvQvg=（+j）.2QxQx.而其旋转速率分量为.1QvQvw=（4j）.2QxQx.三）运动及平衡方程b与体力b，且.根据牛顿运动定律与柯西应力原理，如果质点作用着应力具有速度v.，则在无限小时间段dt内，它们之间的关系为.Qbdv.+pb=p.Qx.dt.式中，p为质点密度。（4）式称为柯西运动方程。当质点的加速度为零时，上式变为静力平衡方程j+pb二0(5)dxij(四)本构方程上述(4)式与(5)式组成的方程组中含有9个方程，15个未知量，其中12个是应力与应变速率分量，3个是速度分

9、量。其余6个关系式则由本构方程提供，本构方程一般具有如下形式6=H9,E,K)(6)ijijijij式中，才为应力变化速率，H表示一个特定的函数关系，K为与荷载历史有关的参数。2.2.2.3 数值模型FLAC3D的数值剖分网格在计算中是按照四面体进行的，四面体的节点也既是网格剖分的节点，因此，每个计算单元有4个面和4个节点(见图2-1)。(4)35一) 空间微分的有限差分逼近对于一个计算单元，若内部各质点速度为一连续的矢量场v，则根据高斯(Gauss)积分原理有JvndsSij7)dv1idv=dxj式中，n为外法向单位矢量场。由于单元体的应变速率是连续的，因此可以近似认为速度是线性变化的，则

10、(7)式可用下面的求和公式近似逼近兖州矿区建筑物下厚煤层安全开采方法关键问题研究QvV、V i-二乙v（f）n（f）S（f）（8）Qxijjf=i式中，V为单元体体积，S（f）为单元体某一面的面积，f为单元体面数,f=1,4,V（f）为面平均速度的i分量。由于速度场是线性的，则有iv（f）=工vii3ii=1,1十f式中，1为单元体节点数，1=1,4；vi1为1节点的i速度分量。将（9）式代入（8）式可得10)V 乞=-为v1工n（f）S（f）Qx3ijji=1f=1,f知根据正交原理有11)工n（f）S（f）=0jf=1dvidxj则将（10）式两边除以V，并将（11）式代入可得12)1工v

11、in（i）S（i）3Viji=1因此有g=ij13)f（vin（i）+vin（i）S（i）6Vijjii=1二）运动方程的节点公式根据前面对质点运动方程的讨论，对于连续介质，当处于平衡状态时，其14)15)平衡方程为j+pB=0Qxijdvidt介质可由若干个作用着体力B，各自产生一定变形的四面体组成，设节点力为fn，n=1,4，利用虚功原理，假定单元体节点具有速度8vn，内部具有变形速率6g，则由节点力fn和体力B所做的外功率与6所做的内功率应当相ij外功率E可用下式表示E=工Vnfn+f8vBdViiViin=1而内功率I为I=f8cdVVijij由（13）式，对于恒定应变速率单元体可有I

12、=-工（8vqn（i）+8vqnQ）S（i）6iijjjijii=由于应力张量是对称的，定义矢量TlTi=cn（/）S（i）iijj可得I=丄8viTi3iii=将(5)式代入(6)式得E=为8Vnfn+Eb+EIiin=6）7）8）9）（20）（2）式中，Eb和Ei分别是体力pb.和惯性力产生的外功率，对于单元体内恒定的体力pb，Eb可写成22)Eb=pbf8vdV.V.EI-fp8VFV23)而EI可以写成37如前所述，单元体内部速度场以线性变化，为方便描述，选取单元体质心为原点，x、x、2X为坐标轴，可有38v=8vnNn.n=l式中，N(n=l,4)是具有如下形式的线性函数Nn=Cn+

13、CnX+CnX+CnX02233(n=l,4则由下面方程决定NnWj,Xj,Xj)=8l23njcncncncn0、2、324)25)26)8nj为克罗内克尔（KrOneCker）记号根据质心定义可得27)Eb=pbXSvnCnVii0n=1用克莱默1(Cramer)法则解得cn=，则0428)Eb=工5vn竺Vi4n=1El=工VJpNn=liV巴dVdt29)Edn=1pbiV-JpNn巴dV4Vdt30)#在稳定状态下，内功率与外功率必定相等，则有TnpbVdvfn=1+iJpNnidVi34vdt上式最后一项，如果单元体内p恒定，则根据质心定义有dvpVdvJpNndV=(i)nV d

14、t4dt若将pV/4看作是假定的节点质量mn,则上式变为dvdvJpNndV=mn(i-)nV dtdt31)32)33)TnpbVdvfn=+imn()ni34dt根据牛顿定律有34)dvF=M(卜)1=1,nidtn35)（匕）为节点dt1在F作用下产生的加速度，n为包含在全部连续介质中单元节点总数。in设质点质量M=m，不平衡力F为式中，F为节点1所受力的i分量，M为节点1的质量,iF二i+fi36)三）时间微分的显式有限差分对于上面（35）式，可写成dvi)dtFMi1=1,nn37)(55)4138）（39）（40）表示为A（TH*Q,gAt）ijijijij记gAtAsijij则A

15、s.为At时间段内应变增量。ij应力增量由下式确定AoA+Aoc.式中，Aoc由下式定义.ocoo)At.kk.kk.而.fv1nv1nS6.11（41）（42）43）44）45）FLAC3D计算中，对于上式的逼近是在一个足够小At时间段内，采用中心差分格式，即上式变为v(t+)v(t)+Fi2i2Mi同理，节点的坐标差分公式与位移差分公式分别为x(t+At)x(t)+Atv(t+A)iii2u(t+At)u(t)+Atv(t+At)iii2四）本构方程的增量形式FLAC3D中，假定在At时间段内，速度是不变的，本构方程的增量形式可2.2.2.4 本构模型一）弹塑性理论的增量关系一般来讲，弹塑

16、性体破坏的判断准则具有如下形式f（o.）0（46）式中，f为一特定的函数关系，巳为n维应力矢量i分量。i弹塑性体的变形As为弹性变形Ase与塑性变形Asp的总和，即iiiAs=Ase+Asp（47）iii弹性应变与应力之间的关系为Ag=S应ei=1,n（48）iin其中，Si为线性函数。塑性变形可用下式表示AspiQgdg_i式中，九为一常数，g为（巳）的某一函数关系。n将（47）、（49）式代入（48）式得Ac=S(As)-九S(Qgii_niQc对于新应力c+Ac，仍有nnf(c+Ac)=0nn当f(c)为线性函数时，(51)式可写成nfG)+f*(Ac)=0nn其中，f*0=/()f(o

17、)。n由于f（c）=o，再将（50）式代入（52）式可得n=049)(50)(51)52)53)若记则有因此有cn=c+Aciiici=c+S(As)iiinfC)=f*ts(As)nnn，a,、fs(Qg/Qcn)*1-f(0)nnn54)则新应力为cN=QI九S(亜)(56)一i一iidCn二)莫尔库伦(Mohr-Coulomb)模型本构关系本次计算的本构模型选择的是莫尔库伦弹塑性材料模型，弹塑性体可产生弹性及塑性两部分变形，根据虎克定律(Hookeslaw)，应力应变的关系为AQ二EAe(57)其中，E为刚度矩阵。对于弹性变形，应力增量可由下式确定Ac=aAse+a(ase+Ase1 1

18、12(23Ac=aAse+aAse+Ase2 122(13Ac=aAse+aAse+Ase58)313212式中，円、a2可由剪切模量与体积模量得出58)式还可写成S(ase,Ase,Ase)=aAse+aCse+Ase1(123)112(23SAse,Ase,Ase丿=aAse+aAse+Ase2(123)122(13SAse,Ase,Ase丿=aAse+aAse+Ase31231321259)三)破坏准则与流动法则若ccc，贝y123莫尔库伦破坏准贝具有形式fs=0，其中fs虫1-c3Ne+2cJNe对于拉张破坏，有ft=0的形式，其中f二cctt3上两式中，为内摩擦角，c为内聚力，ct为

19、抗拉强度,N=1+sin()1sin()(60)61)(62)对于材料的抗拉强度ct,其值不可能超过c3(见图2-2)，ct的最大值为63)64)65)(66)图2-2莫尔一库伦破坏准则示意图Ctmax对于塑性流动状态，设gs与gt分别为剪切破坏和拉张破坏所对应的函数关系，有gt=C3式中，屮为膨胀角，1+sinG)%=1sinG)#流动准则可定义为函数关系h(c,C)=0，在C1、C3平面内将fs=0和1313ft=0以上部分分为两个区域(见图2-3)，其函数表达式如下兖州矿区建筑物下厚煤层安全开采方法关键问题研究h二-Qt+apC-qp)(67)31其中，ap与qp由下式定义68)ap=.

20、1+N2+Neeqp二QtN2c-Neve从图2-3中可见，若某点的应力状态在1区(fs0,h0,h0)时，则发生拉张破坏。(四)塑性阶段应力修正对于剪切破坏，由(64)式可得do4312do69)将(59)式中的Ase1替，得do屮3Ase和Ae分别用dgs/do、dgs/do和dgs/do代231Sdgsdgsdgs1、do1do2do丿3Sdgsdgsdg2do1do2do)3Sdgsdgsdg31do1do2do丿3=aaN1 2屮=a(N2 屮=aN+a1屮2再由(55)、(56)式，利用f=fs可得on=oi兀CaN)11丫2片on=oi兀aiN丿对于拉张破坏，由(66)式可得(7

21、0)71)(72)=0SQg1Qg2Qg丿3（QcQgsQgsQgs、Qg1，Qg2，Qg3、rQgsQgsQgsQg2=X2S2Qg丿33（Qg1=（X1=X1再利用f=ft可得（73）（74）（75）（76）（77）（78）、k分别由1dgt=0Qg2埜二1Qg3同上，利用（59）式可得Gn=Gikx1 11GN=GIktX2 22GN=GIksX3 31GIGtkt=3X1将（76）式代入（75）式可得GN=GIC?IGt21 13XGN=GICIGt2 23X1GN=Gt3五）塑性应变增量利用（69）式与（73）式，由（49）式可得AsPs=ks1AsPs=ksN3屮AsPs=kt3式中，APs与Ast分别为塑性剪应变增量与塑性张应变增量;ii72）、（76）式给出。#