FLAC3D原理.docx

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

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

3、附加于方程中。FLAC3D具有一个功能强大的网格生成器，有12种基本形状的单元体可供选择，利用这12种基本单元体，几乎可以构成任何形状的空间立体模型。FLAC3D主要是为地质工程应用而开发的岩土体力学数值评价计算程序，自身设计有九种材料本构模型：（1）空模型（Null Model）（2）弹性各向同性材料模型（Elastic, Isotropic Model）（3）弹性各向异性材料模型（Elastic, anisotropic Model）（4）德拉克-普拉格弹塑性材料模型（Drucker-Prager Model）（5）莫尔-库伦弹塑性材料模型（Mohr-Coulomb Model）（6）应变

4、硬化、软化弹塑性材料模型（Strain-Hardening/Softening Mohr-Coulomb Model）（7）多节理裂隙材料模型（Ubiquitous-Joint Model）（8）双曲型应变硬化、软化多节理裂隙材料模型（Bilinear Strain-Hardening/Softening Ubiquitous-Joint Model）（9）修正的Cam粘土材料模型（Modified Cam-clay Model）除上述本构模型之外，FLAC3D还可进行动力学问题、水力学问题、热力学问题等的数值模拟。在边界条件及初始条件的考虑上，FLAC3D软件十分灵活方便，可在数值计算过程中

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

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

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

8、小的时间dt，介质会产生一个由vidt决定的无限小应变，对应的应变速率分量 ij为 （2）而其旋转速率分量 ij为 （3）（三）运动及平衡方程根据牛顿运动定律与柯西应力原理，如果质点作用着应力 ij与体力bi，且具有速度vi，则在无限小时间段dt，它们之间的关系为 （4）式中， 为质点密度。（4）式称为柯西运动方程。当质点的加速度为零时，上式变为静力平衡方程 （5）（四）本构方程上述（4）式与（5）式组成的方程组中含有9个方程，15个未知量，其中12个是应力与应变速率分量，3个是速度分量。其余6个关系式则由本构方程提供，本构方程一般具有如下形式 （6）式中，为应力变化速率，H表示一个特定的函数

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

10、数，l=1,4；vil为l节点的i速度分量。将（9）式代入（8）式可得 （10）根据正交原理有 （11）则将（10）式两边除以V，并将（11）式代入可得 （12）因此有 （13）（二）运动方程的节点公式根据前面对质点运动方程的讨论，对于连续介质，当处于平衡状态时，其平衡方程为 （14） （15）介质可由若干个作用着体力B，各自产生一定变形的四面体组成，设节点力为fn，n=1,4，利用虚功原理，假定单元体节点具有速度 vn，部具有变形速率 ，则由节点力fn和体力B所做的外功率与 ij所做的功率应当相等。外功率E可用下式表示 （16）而功率I为 （17）由（13）式，对于恒定应变速率单元体可有 （

11、18）由于应力量是对称的，定义矢量Tl （19）可得 （20）将（15）式代入（16）式得 （21）式中，Eb和EI分别是体力 bi和惯性力产生的外功率，对于单元体恒定的体力 bi，Eb可写成 （22）而EI可以写成 （23）如前所述，单元体部速度场以线性变化，为方便描述，选取单元体质心为原点，、为坐标轴，可有 （24）式中，Nn(n=1,4)是具有如下形式的线性函数 （25）、（n=1,4）则由下面方程决定 （26） nj为克罗克尔（Kronecker）记号，根据质心定义可得（27）用克莱默（Cramers）法则解得，则 （28）而 （29）则 （30）在稳定状态下，功率与外功率必定相等，则

12、有 （31）上式最后一项，如果单元体 恒定，则根据质心定义有 （32）若将 V/4看作是假定的节点质量mn，则上式变为 （33）故 （34）根据牛顿定律有 （35）式中，Fi为节点l所受力的i分量，M为节点l的质量，为节点l在Fi作用下产生的加速度，nn为包含在全部连续介质中单元节点总数。设质点质量，不平衡力F为 （36）（三）时间微分的显式有限差分对于上面（35）式，可写成 （37）FLAC3D计算中，对于上式的逼近是在一个足够小 t时间段，采用中心差分格式，即上式变为 （38）同理，节点的坐标差分公式与位移差分公式分别为 （39） （40）（四）本构方程的增量形式FLAC3D中，假定在 t

13、时间段，速度是不变的，本构方程的增量形式可表示为 （41）记 （42）则 ij为 t时间段应变增量。应力增量由下式确定 （43）式中，由下式定义（44）而 （45）2.2.2.4 本构模型（一）弹塑性理论的增量关系一般来讲，弹塑性体破坏的判断准则具有如下形式 （46）式中，f为一特定的函数关系，为n维应力矢量i分量。弹塑性体的变形为弹性变形与塑性变形的总和，即 （47）弹性应变与应力之间的关系为 （48）其中，Si为线性函数。塑性变形可用下式表示 （49）式中， 为一常数，g为的某一函数关系。将（47）、（49）式代入（48）式得 （50）对于新应力，仍有 （51）当为线性函数时，（51）式可

14、写成 （52）其中，。由于，再将（50）式代入（52）式可得 （53）若记则有 （54）因此有 （55）则新应力为 （56）（二）莫尔库伦（Mohr-Coulomb）模型本构关系本次计算的本构模型选择的是莫尔库伦弹塑性材料模型，弹塑性体可产生弹性及塑性两部分变形，根据虎克定律（Hookes law），应力应变的关系为 （57）其中，E为刚度矩阵。对于弹性变形，应力增量可由下式确定 （58）式中， 1、 2可由剪切模量与体积模量得出（58）式还可写成 （59）（三）破坏准则与流动法则若，则莫尔库伦破坏准则具有形式fs= 0，其中 （60）对于拉破坏，有ft = 0的形式，其中 （61）上两式中， 为摩擦角，c为聚力， t为抗拉强度，（62）对于材料的抗拉强度 t，其值不可能超过 3（见图2-2）， t的最大值为 （63）对于塑性流动状态，设g s与g t分别为剪切破坏和拉破坏所对应的函数关系，有 （64） （65）式中， 为膨胀角， （66）流动准则可定义为函数关系，在 1、 3平面将和以上部分分为两个区域（见图2-3）