FLAC3D原理

....2.2三维数值模拟方法及其原理2.2.1FLAC3D工程分析软件特点FLAC3D是由美国ItascaConsultingGroup,Inc.为地质工程应用而开发的连续介质显式有限差分计算机软件。FLAC即FastLagrangianAnalysisofContinua的缩写。该软件主要适用于模拟计算岩土体材料的力学行为及岩土材料达到屈服极限后产生的塑性流动，对大变形情况应用效果更好。FLAC3D程序在数学上采用的是快速拉格朗日方法，基于显式差分来获得模型全部运动方程和本构方程的步长解，其本构方程由基本应力应变定义及虎克定律导出，运动平衡方程则直接应用了柯西运动方程，该方程由牛顿运动定律导出。计算模型一般是由若干不同形状的三维单元体组成，也即剖分的空间单元网络区，计算中又将每个单元体进一步划分成由四个节点构成的四面体，四面体的应力应变只通过四个节点向其它四面体传递，进而传递到其它单元体。当对某一节点施加荷载后，在某一个微小的时间段内，作用于该点的荷载只对周围的若干节点（相邻节点）有影响。利用运动方程，根据单元节点的速度变化和时间，可计算出单元之间的相对位移，进而求出单元应变，再利用单元模型的本构方程，可求出单元应力。在计算应变过程中，利用高斯积分理论，将三维问题转化为二维问题而使其简单化。在运动方程中，还充分考虑了岩土体所具有的粘滞性，将其视作阻尼附加于方程中。FLAC3D具有一个功能强大的网格生成器，有12种基本形状的单元体可供选择，利用这12种基本单元体，几乎可以构成任何形状的空间立体模型。FLAC3D主要是为地质工程应用而开发的岩土体力学数值评价计算程序，自身设计有九种材料本构模型：（1）空模型（NullModel）（2）弹性各向同性材料模型（Elastic,IsotropicModel）（3）弹性各向异性材料模型（Elastic,anisotropicModel）（4）德拉克-普拉格弹塑性材料模型（Drucker-PragerModel）（5）莫尔-库伦弹塑性材料模型（Mohr-CoulombModel）....（6）应变硬化、软化弹塑性材料模型（Strain-Hardening/SofteningMohr-CoulombModel）（7）多节理裂隙材料模型（Ubiquitous-JointModel）（8）双曲型应变硬化、软化多节理裂隙材料模型（BilinearStrain-Hardening/SofteningUbiquitous-JointModel）（9）修正的Cam粘土材料模型（ModifiedCam-clayModel）除上述本构模型之外，FLAC3D还可进行动力学问题、水力学问题、热力学问题等的数值模拟。在边界条件及初始条件的考虑上，FLAC3D软件十分灵活方便，可在数值计算过程中随时调整边界条件和初始条件。FLAC3D具有强大的后处理功能，用户可以直接在屏幕上绘制或以文件形式创建或输出打印多种形式的图形、文字，用户还可根据各自的需要，将若干个变量合并在同一幅图形中进行研究分析。FLAC3D软件还可对各种开挖工程或施加支护工程等进行数值仿真模拟，软件自身设计有锚杆、锚索、衬砌、支架等结构元素，可以直接模拟这些支护于围岩（土）体的相互作用。FLAC3D拥有可以自行设计的FISH语言，用户可根据自身需求，自己设计材料的本构模型、屈服准则、支护方案、复杂形状的开挖方式等工作。特别注意的是，岩石是一种脆性材料，当外荷载达到岩石强度后，材料发生断裂破坏，产生弱化现象，应属于弹塑性体。在FLAC3D中，一般对于弹塑性材料，判断其破坏与否的基本准则有两个，即Drucker-Prager准则和Mohr-Coulomb准则。根据室内岩石力学性质试验结果，其典型应力应变曲线反映出岩体破坏包络线符合莫尔—库伦屈服准则，故本次建立的本构力学模型选择莫尔—库伦弹塑性材料模型为宜。2.2.2FLAC3D分析计算原理计算所采用的数学模型是根据弹塑性理论的基本原理（应变定义、运动定律、能量守衡定律、平衡方程及理想材料的连续性方程等）而建立的。2.2.2.1基本约定在数学及数值模型的表达式中，符号有一定的约定含义，一般[A]表示张....量，Aij表示张量[A]的（i,j）分量，[a]表示矢量，ai表示矢量[a]的i分量，,i表示对xi的偏导数。xi，ui，vi和dvi/dt，（i=1,3）分别表示一点的位置矢量分量、位移矢量分量、速度矢量分量和加速度矢量分量。2.2.2.2数学模型（一）柯西（Cauchy）应力张量与柯西公式对于一个具有体积V的封闭曲面s的物体，在其上取一表面元素s，这个表面元素的单位外法向矢量为n，在某一时刻t，在表面元素对于连续介质中一点，作用着对称的应力张量ij，根据s上作用有力P，则极限dsdPsPTs0lim称为表面力。若用ti表示T的分量，则在三维直角坐标系中可有关系式ijiint（1）这个关系式称为柯西公式，其中，ij称为柯西应力张量。（二）应变速率和旋转速率如果介质质点具有运动速度矢量[v]，则在一个无限小的时间dt内，介质会产生一个由vidt决定的无限小应变，对应的应变速率分量ij为)(21ijjiijxvxv（2）而其旋转速率分量ij为)(21ijjiijxvxv（3）（三）运动及平衡方程根据牛顿运动定律与柯西应力原理，如果质点作用着应力ij与体力bi，且具有速度vi，则在无限小时间段dt内，它们之间的关系为dtdvbxiijij（4）式中，为质点密度。（4）式称为柯西运动方程。....当质点的加速度为零时，上式变为静力平衡方程0ijijbx（5）（四）本构方程上述（4）式与（5）式组成的方程组中含有9个方程，15个未知量，其中12个是应力与应变速率分量，3个是速度分量。其余6个关系式则由本构方程提供，本构方程一般具有如下形式),,(][ijijijijH（6）式中，][为应力变化速率，H表示一个特定的函数关系，为与荷载历史有关的参数。2.2.2.3数值模型FLAC3D的数值剖分网格在计算中是按照四面体进行的，四面体的节点也既是网格剖分的节点，因此，每个计算单元有4个面和4个节点（见图2-1）。（一）空间微分的有限差分逼近对于一个计算单元，若内部各质点速度为一连续的矢量场[v]，则根据高斯（Gauss）积分原理有SjiVjidsnvdvxv（7）式中，[n]为外法向单位矢量场。由于单元体的应变速率是连续的，因此可以近似认为速度是线性变化的，则（7）式可用下面的求和公式近似逼近41)()()(fffjfijiSnvxvV（8）单元面单元节点图2-1四面体单元示意图....式中，V为单元体体积，)(fS为单元体某一面的面积，f为单元体面数，f=1,4，)f(iv为面平均速度的i分量。由于速度场是线性的，则有4,1)(31flllifivv（9）式中，l为单元体节点数，l=1,4；vil为l节点的i速度分量。将（9）式代入（8）式可得4,1)()(4131lffffjllijiSnvxvV（10）根据正交原理有041)()(fffjSn（11）则将（10）式两边除以V，并将（11）式代入可得)()(4131lljllijiSnvVxv（12）因此有)(41)()()(61llliljljliijSnvnvV（13）（二）运动方程的节点公式根据前面对质点运动方程的讨论，对于连续介质，当处于平衡状态时，其平衡方程为0ijijBx（14）)(dtdvbBiii（15）介质可由若干个作用着体力[B]，各自产生一定变形的四面体组成，设节点力为[f]n，n=1,4，利用虚功原理，假定单元体节点具有速度[v]n，内部具有变形速率[]，则由节点力[f]n和体力[B]所做的外功率与ij所做的内功率应当相等。外功率E可用下式表示dVBvfvEiVininni41（16）....而内功率I为VijijdVI（17）由（13）式，对于恒定应变速率单元体可有41)(61illiijljljijliSnvnvI（18）由于应力张量是对称的，定义矢量TllljijliSnT（19）可得4131lliliTvI（20）将（15）式代入（16）式得41nIbniniEEfvE（21）式中，Eb和EI分别是体力bi和惯性力产生的外功率，对于单元体内恒定的体力bi，Eb可写成ViibdVvbE（22）而EI可以写成ViiIdVdtdvvE（23）如前所述，单元体内部速度场以线性变化，为方便描述，选取单元体质心为原点，1x、2x、3x为坐标轴，可有41nnniiNvv（24）式中，Nn(n=1,4)是具有如下形式的线性函数3322110xcxcxccNnnnnn（25）nc0、nc1、nc2、nc3（n=1,4）则由下面方程决定njjjjnxxxN321,,（26）nj为克罗内克尔（Kronecker）记号，根据质心定义可得410nnniibVcvbE（27）....用克莱默（Cramer’s）法则解得410nc，则414ninibVbvE（28）而41nVinniIdVdtdvNvE（29）则414nVinininidVdtdvNVbfvE（30）在稳定状态下，内功率与外功率必定相等，则有VinininidVdtdvNVbTf43（31）上式最后一项，如果单元体内恒定，则根据质心定义有niVindtdvVdVdtdvN)(4（32）若将V/4看作是假定的节点质量mn，则上式变为ninVindtdvmdVdtdvN)(（33）故ninininidtdvmVbTf)(43（34）根据牛顿定律有nlillinldtdvMF,1)(（35）式中，Fil为节点l所受力的i分量，Ml为节点l的质量，lidtdv)(为节点l在Fil作用下产生的加速度，nn为包含在全部连续介质中单元节点总数。设质点质量llmM]][[，不平衡力[F]l为liliilifVbTF]43[（36）（三）时间微分的显式有限差分对于上面（35）式，可写成....nlillinlFMdtdv,11)(（37）FLAC3D计算中，对于上式的逼近是在一个足够小t时间段内，采用中心差分格式，即上式变为lilliliFMtttvttv)2()2(（38）同理，节点的坐标差分公式与位移差分公式分别为)2()()(tttvtxttxlilili（39）)2()()(tttvtuttulilili（40）（四）本构方程的增量形式FLAC3D中，假定在t时间段内，速度是不变的，本构方程的增量形式可表示为),(*tHijijijij（41）记ijijt（42）则ij为t时间段内应变增量。应力增量由下式确定Cijijij（43）式中，Cij由下式定义tkjikkjikCij)(（44）而llliljljliijSnvnv4161（45）2.2.2.4本构模型（一）弹塑性理论的增量关系一般来讲，弹塑性体破坏的判断准则具有如下形式0)(if（46）式中，f为一特定的函数关系，i为n维应力矢量i分量。弹塑性体的变形i为弹性变形ei与塑性变形pi的总和，即....pieii（47）弹性应变与应力之间的关系为niSen