2.2 三维数值模拟方法及其原理
2.2.1 FLAC3D工程分析软件特点
FLAC3D是由美国Itasca Consulting Group, Inc. 为地质工程应用而开发的连续介质显式有限差分计算机软件。FLAC即Fast Lagrangian Analysis of Continua 的缩写。该软件主要适用于模拟计算岩土体材料的力学行为及岩土材料达到屈服极限后产生的塑性流动，对大变形情况应用效果更好。
FLAC3D程序在数学上采用的是快速拉格朗日方法，基于显式差分来获得模型全部运动方程和本构方程的步长解，其本构方程由基本应力应变定义及虎克定律导出，运动平衡方程则直接应用了柯西运动方程，该方程由牛顿运动定律导出。
计算模型一般是由若干不同形状的三维单元体组成，也即剖分的空间单元网络区，计算中又将每个单元体进一步划分成由四个节点构成的四面体，四面体的应力应变只通过四个节点向其它四面体传递，进而传递到其它单元体。当对某一节点施加荷载后，在某一个微小的时间段内，作用于该点的荷载只对周围的若干节点（相邻节点）有影响。利用运动方程，根据单元节点的速度变化和时间，可计算出单元之间的相对位移，进而求出单元应变，再利用单元模型的本构方程，可求出单元应力。在计算应变过程中，利用高斯积分理论，将三维问题转化为二维问题而使其简单化。在运动方程中，还充分考虑了岩土体所具有的粘滞性，将其视作阻尼附加于方程中。
FLAC3D具有一个功能强大的网格生成器，有12种基本形状的单元体可供选择，利用这12种基本单元体，几乎可以构成任何形状的空间立体模型。
FLAC3D主要是为地质工程应用而开发的岩土体力学数值评价计算程序，自身设计有九种材料本构模型：
（1）空模型（Null Model）
（2）弹性各向同性材料模型（Elastic, Isotropic Model） （3）弹性各向异性材料模型（Elastic, anisotropic Model） （4）德拉克-普拉格弹塑性材料模型（Drucker-Prager Model） （5）莫尔-库伦弹塑性材料模型（Mohr-Coulomb Model）
— 29 —
兖州矿区建筑物下厚煤层安全开采方法关键问题研究 （6）应变硬化、软化弹塑性材料模型（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软件十分灵活方便，可在数值计算过程中随时调整边界条件和初始条件。
FLAC3D具有强大的后处理功能，用户可以直接在屏幕上绘制或以文件形式创建或输出打印多种形式的图形、文字，用户还可根据各自的需要，将若干个变量合并在同一幅图形中进行研究分析。
FLAC3D软件还可对各种开挖工程或施加支护工程等进行数值仿真模拟，软件自身设计有锚杆、锚索、衬砌、支架等结构元素，可以直接模拟这些支护于围岩（土）体的相互作用。
FLAC3D拥有可以自行设计的FISH语言，用户可根据自身需求，自己设计材料的本构模型、屈服准则、支护方案、复杂形状的开挖方式等工作。
特别注意的是，岩石是一种脆性材料，当外荷载达到岩石强度后，材料发生断裂破坏，产生弱化现象，应属于弹塑性体。在FLAC3D中，一般对于弹塑性材料，判断其破坏与否的基本准则有两个，即Drucker-Prager准则和Mohr-Coulomb准则。根据室内岩石力学性质试验结果，其典型应力应变曲线反映出岩体破坏包络线符合莫尔—库伦屈服准则，故本次建立的本构力学模型选择莫尔—库伦弹塑性材料模型为宜。 2.2.2 FLAC3D分析计算原理
计算所采用的数学模型是根据弹塑性理论的基本原理（应变定义、运动定律、能量守衡定律、平衡方程及理想材料的连续性方程等）而建立的。 2.2.2.1 基本约定
在数学及数值模型的表达式中，符号有一定的约定含义，一般[A]表示张
— 30 —
量，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，则极限
T?lim?PdP?
?s?0?sds称为表面力。
若用ti表示T的分量，则在三维直角坐标系中可有关系式
ti?ni?ij （1） 这个关系式称为柯西公式，其中，?ij称为柯西应力张量。 （二）应变速率和旋转速率
如果介质质点具有运动速度矢量[v]，则在一个无限小的时间dt内，介质会产生一个由vidt决定的无限小应变，对应的应变速率分量?ij为
1?v?vj?ij?(i?) （2）
2?xj?xi而其旋转速率分量?ij为
1?vi?vj?ij?(?) （3）
2?xj?xi（三）运动及平衡方程
根据牛顿运动定律与柯西应力原理，如果质点作用着应力?ij与体力bi，且具有速度vi，则在无限小时间段dt内，它们之间的关系为
??ijdv??bi??i （4） ?xjdt式中，?为质点密度。（4）式称为柯西运动方程。
— 31 —