AXES.h 声明一个与坐标轴变换有关的数据类型Axes和相关函数
CONTABLE.h 声明了一个ConTableList类，该类为本构模型提供了一张表，该表用来存储模型单元或节点ID号 USERMODEL.h 用户自定义本构模型派生类的声明 USERMODEL.cpp 用户自定义本构模型的由应变增量获得应力增量的实现
在编写自定义本构模型的C++文件时，需要在FLAC3D自带的莫尔-库仑模型基础上增加新的材料参数变量。如中尉编写的节理岩体损伤莫尔-库仑模型增加如图所示的材料参数。
FLAC3D本构模型的主要功能就是根据应变增量得到新的应力增量，应力应变关系表达式很可能为矩阵表达式，在计算时，可能需要对矩阵求逆。对于逆矩阵的求取，可根据逆矩阵定义或矩阵的初等行变换来进行，建议采用一种在数值计算方法中常用的全主元高斯－约当（Gauss-Jordan）消去法。该方法具有数值稳定而且精度较高的特点。由于岩体的柔度张量矩阵元素都非常小（10的-7、-8次方级），在进行乘法计算时，将产生极大的舍入误差，而高斯－约当消去法的消去过程是在全矩阵中选主元（绝对值最大的元素）来进行，故可使舍入误差对结果的影响减小到最小。将编写好的头文件jdmohr.h和源文件jdmohr.cpp导入到工程文件（.dsw）中，通过编译（compile）和链接（link），形成了节理岩体损伤莫尔-库仑本构模型的动态链接库文件jdmohr.dll，将其复制到FLAC3D的安装目录下，如图
程序的调试是编程过程中不可缺少的重要步骤，尤其是对于数值计算，需要实时监控计算流程和各种变量，以防出现死循环和结果严重失真。在VC++中可设置将FLAC3D中的
10
EXE文件路径加入到程序的调试范围中，并将自定义的FLAC3D的.dll文件加入到附加动态链接库(Additional DLLs)中，然后在.cpp文件里的Initialize()或Run()函数中设置断点，进行调试。使用自定义模型时，首先在命令流中使用命令CONFIG cppudm对FLAC3D进行配置，使其能接收动态链接库文件，然后通过MODEL load命令将自定义本构模型（即动态链接库文件）加载到FLAC3D中，这样，FLAC3D就可以识别出新的模型名和属性名。同样，在恢复（命令：RESTORE）一个使用自定义模型的文件时，也必须使用CONFIG cppudm命令和MODEL load命令。 如图
6．FLAC3D自定义本构模型
对应manual之Optional Features中最后关于自定义本构模型的部分，给出以下翻译。 FLAC3D中可以用C++编写动态链接库文件（.dll文件）来自定义本构模型，其主要功能是根据应变增量返回新的应力张量。
C++语言的特点是采用面向对象方法，使用类来代表对象进行编程。与对象有关的信息被封装在类中，这些信息在类的外部是不可见的。与对象的通信通过成员函数操作封装数据来完成。一个新的对象类型可以从基类派生而来，基类的成员函数也可被派生类的同名函数覆盖。这种方式便于程序的模块化设计。
用C++自定义FLAC3D本构模型的过程主要包括：基类、成员函数的定义，模型注册，模型与FLAC3D间的数据传递以及模型状态指示。 本构模型的基类
基类提供了实际本构模型的框架，其它类都来源于该类。基类命名为ConstitutiveModel，由于它定义了一系列“虚”（语法上以0为标志）的成员函数，因此又被成为“抽象”类。这也是说不能创建该类的任何对象，而其它继承自该类的类的对象必须给出实际的成员函数来代替基类的“虚”函数。基类的头文件（.h）中定义了这些“虚”的成员函数。 成员函数
ConstitutiveModel类的成员函数及功能如下：
const char *Keywork() 返回一个指向包括模型名称的字符数组的指针，以便用户在FLAC3D中使用MODEL命令时，C++能够识别。
const char *Name() 返回一个指向包含模型名称的字符数组的指针，以识别用户使用如PRINT Zone等命令。
const char **Properties() 返回一个指向包含模型力学指标名称的字符串数组的指针，并用一个空指针代表字符串数组的结束。这样程序能识别用户输入的PROPERTY命令。
11
const char *States() 返回一个包含单元状态名称的字符串数组的指针，并用一个空指针代表字符串数组的结束。这样程序就能输出或显示用户自定义的模型内部状态（如：塑性流动、屈服、受拉）。
SetProperty(unsigned n, const double &dVal) dVal的值由FLAC3D命令PROP > Double GetProperty(unsigned n) 返回力学指标顺序为n的值，供模型使用。
const char *Copy(const ConstitutiveModel *cm) 该成员函数首先调用基类的Copy函数，然后从指向cm的模型对象中拷贝所有的数据。如果出现错误，返回一个指向错误字符串的指针；否则，返回0。在调用Initialize()函数时并不需要拷贝被计算的数据。
const char *Initialize(unsigned uDim, State *ps) 该函数在使用命名CYCLE命令和当大应变更新坐标时对每一个单元使用。模型对象可能对属性变量和状态变量进行初始化，也可能不初始化。uDim为计算维数（在FLAC3D中为3）。ps是结构子，包含了当前被计算的子单元的当前应力分量和前一步的应变增量分量。如果出现错误，返回一个指向错误字符串的指针；否则，返回0。注意，在调用Initialize()时应变没被定义。单元的平均应力分量在状态结构子中可以得到；不能通过Initialize()来改变它们。
const char *Run(unsigned uDim, State *ps) 该函数在FLAC3D计算循环的每一步对每一个子单元使用。模型必须根据应变增量对应力张量进行更新。ps结构子包含当前被计算的子单元的当前应力分量和前一步的应变增量分量。当调用Run()函数时，应力分量已经包含了应力修正。如果出现错误，返回一个指向错误字符串的指针；否则，返回0。
double ConfinedModulus(void) 模型对象必须返回一个对最大受限模量做出最好估计的值，这被用来计算稳定的计算时步。
double ShearModulus(void) 返回一个对当前切线剪切模量做出最好估计的值，FLAC3D据此来确定动态模式下静态边界系数。
double BulkModulus(void) 返回一个对当前切线体积模量做出最好估计的值，FLAC3D在当前步不调用该函数。
double SafetyFactor(void) 返回某个安全系数值，FLAC3D在当前步不调用该函数。
ConstitutiveModel *Clone(void) 创建同一个类的新对象作为当前对象，返回一个指向ConstitutiveModel类的指针，该函数在单元使用该模型时被调用。
const char *SaveRestore(ModelSaveObject *mso) 该函数在使用命名SAVE或RESTORE时被调用。模型对象首先调用基类的SaveRestore()函数对文件进行保存和恢复。其数据类型只允许整数和实数变量，其它数据类型需要被转换成这两种类型。派生类函数必须调用mso->Initialize(nd, ni)，这里nd为双精度型数据个数，ni为整数型数据个数。这些变量然后通过调用mso->Save(ns, var)，这里ns为变量顺序（从0到nd-1或0到ni-1，取决于保存整数或实数）。var是要保存的变量。对于双精度或整数型变量，分别有Save()函数来处理。 模型类定义也应该包含一个能调用基类结构子的结构子。如果基类结构子由bRegister调用，其值设为true，则派生的模型就被FLAC3D注册了。跟模型唯一对应的类型号（uTypeIn）也必须被传递。这样就能够在从保存文件中恢复时重新调用正确的模型。建议使用较大的类型号（100或更大）以避免和内在模型（从1开始编号）发生冲突。数据成员的初始化可以通过结构子来完成。 模型的注册
每一个用户自定义本构模型有模型名称、力学参数名称和状态指示器。FLAC3D通过调用合适的成员函数来确定这些信息。FLAC3D通过调用模型对象的静态全局实例获得用户定义模型的信息。true值使得基类的结构子将新模型注册并将其加入到模型列表中。一个特定模型只有一个静态注册产生，方便将其放到模型的C++源码中，这样当相应的.dll文件被装载时，模型被注册。
12
模型与FLAC3D之间的信息交换
在FLAC3D和用户自定义模型间的最重要的连接就是成员函数Run(unsigned nDim, State *ps)，它在模型的计算循环时计算其力学响应。State结构被用来传输信息和生成模型。State结构的成员主要如下：
Run()的主要功能就是根据应变增量计算出新的应力。在非线性模型中，它也用来传递模型的内部状态。如补充模型显示当前处于屈服状态或者以前已经发生屈服。每一个子单元设定变量mState以一串二进制代码记录模型状态。 模型状态指示器
FLAC3D中的单元是由四面体子单元所组成，每一个四面体有记录其当前状态的成员变量。该成员变量共有16位，能够代表最多15种不同的状态。对于用户定义本构模型，用户可以命名一种状态并为其分配特定的位。FLAC3D对组成单元的每个四面体单元调用Run()函数来更新其应力值。在这个过程中，状态指示器也通过本构模型更新。对于内置模型，状态指示器显示了四面体单元的破坏状态。本构模型对四面体状态指示器变量和当前破坏状态使用逻辑“或”运算符来进行更新。
7．数值计算中初始应力场的模拟
中尉认为，在做数值模拟时，岩体初始应力场就是在本次开挖前岩体的应力状态，这种应力状态应该是稳定的（数值模拟中以最大不平衡比率来衡量）。但在实际中，岩体的初始应力状态由于受到多种因素的影响而非常复杂。以一个露天煤矿边坡岩体为例，首先，重力作为外加载荷将产生重力场，如果该岩体在过去的地质历史中经历了各种动力运动，由于温度不均、水压梯度、地表剥蚀或其它物理化学变化等也可引起相应的应力场。除了自重应力场和构造应力场外，由于露天煤矿采用爆破开采，这种动载荷对边坡岩体应力场的扰动也不可忽略。现在要对在露天煤矿边坡岩体下的井工开采进行数值模拟，考虑上述各种因素，则初始应力场的确定就十分困难了。如果有条件做地应力测量的话，当然更好。如果没有条件，就只能做简化。中尉的论文只考虑井工开采对露天边坡的扰动，因此忽略了动载荷的影响，也没有考虑构造应力，只加入自重应力场。开挖前的应力状态首先是一种平衡状态，因此自重应力作为外力应与岩体中存在的应力场产生平衡。由于考虑本次开挖的影响，因此开挖前的应力状态一般也应不产生破坏（塑性区）而且位移场也应为零。因此可认为：在进行数值模拟时，岩体在自重应力场下达到平衡状态时的应力状态就是初始应力场。但同时还要对位移和塑性状态清零。以FLAC3D为例，可有如下获取初始应力场并开挖计算的方法： (1) 对模型施加重力，采用弹性模型计算至平衡状态，保存； (2) 恢复保存文件，位移清零，采用弹塑性模型，施加重力； (3) 开挖，计算，保存结果。
这种方法由于第一步采用弹性模型可避免出现塑性区，但第二步位移场仍然需要清零，值得注意的是，第二步仍然需要施加重力，因为这样才能保证此时的平衡状态（重力与第一步计算出的应力场平衡）。
8．FLAC3D应变分析
按塑性力学理论，一点的应变张量可以分为与体积变化有关的球形应变张量和与物体形状变化有关的应变偏量。FLAC3D里称球形应变张量为体积应变张量，应变偏量为剪应变张量。在进行理想塑性分析时，常假设体积不变，此时球形应变张量为零，应变偏量等于应变张量
13
（尤其对于金属材料）；或者考虑体积力只引起体积变形，剪切力只引起剪切变形。但在进行岩土塑性分析时，尤其是对于土体，具有压硬性和剪胀性的特性。压硬性指的是静水压力（体积力）与剪切变形的耦合作用，即静水压力也会引起剪切变形。剪胀性是指切应力与体积应变的耦合作用，即切应力也会引起体积改变（扩容或体积收缩）。对于岩土材料来说，静水压力不仅产生弹性和塑性的体应变，而且由于静水压力的存在，还会引起剪变形刚度的增大而使切应变变化；而切应力不仅会产生弹性或塑性的切应变，而且会引起剪胀或剪缩。因此与弹性和理想弹塑性材料相比，岩土的本构关系显然要复杂得多。但FLAC3D中常用的莫尔－库仑模型还是采用的理想塑性分析。 FLAC3D的塑性本构模型采用的是增量理论，因此其计算结果里并没有应变值，而出现的是体积应变增量(vsi)、剪应变增量（ssi）和体积应变率(vsr)、剪应变率(ssr)。由于应变增量张量是一个三维张量，上述各项应该出现的是由多个值组成的三维状态量，但FLAC3D里给出的都只有一个值，实际上对体积应变给出的是平均应变增量；对剪应变给出的是最大剪应变增量值。
尽管FLAC3D的计算结果不能直接得到全应变，但后处理程序语言FISH中还提供了全应变增量fsi和全应变率fsr的变量z_fsi(p_z,arr)和z_fsr(p_z,arr)，其中，arr表示六个分量。可根据这两个变量写FISH代码得到全应变增量和全应变率。 好像实际分析中用到体积应变的很少，对剪应变的分析多一些，如在边坡分析中通过剪应变率形成的条带判断滑动面的生成。其原理在于：滑动面各点的剪应变改变速率要快于沿该点滑移面法线方向上的各点。
9．FLAC3D的调参
以下转自中尉发到SIMWE论坛FLAC\\FLAC3D版面的帖子 数值模拟与经验方法的不同之处在于前者具有普适性，后者则局限在相似的条件下使用。一般经验方法需要的参数较少并且容易获得和比较准确，如可通过有效的观测手段直接得出。但对于数值模拟的参数就不一样了，因为对于岩土工程来说，限于测试手段，一般很少进行大型的岩体原位三维试验而在实验室中对取自现场的岩样进行试验。这样得到的岩样参数往往并不能代表岩体参数，人们已经对岩体参数进行了一些研究，但这些研究很多仍然靠一些模糊的分类来界定参数，还没达到定量化或者实际使用效果并不是很好。数值模拟说到底只是一种数学方法，任何参数的变化均可能得出跟实际完全相悖的结果。因此对于数值模拟，调参就是一项很重要的工作，它决定了模拟结果的可信性；同时由于数值模拟中的参数众多，为得到可信结果，参数怎么调就成为一个比较棘手的问题。针对不同的要求需要调哪些参数？调大还是调小？调一个还是调多个？这些问题是不能用排列组合和穷举法一一试验的。这些问题对我们使用数值模拟方法就提出了一个前提条件，即：我们必须对所模拟的问题有一个相当的了解，对其模拟结果有个大概的估计，这种了解和估计源于其他已被验证过的理论或经验知识。我们通过调参后得到的数值模拟结果不能违反这些经验知识。
数值模拟的调参应该根据计算原理来分析。以FLAC3D中常用的莫尔库仑模型为例，这是一个假设单元在弹性阶段为线弹性材料，在塑性阶段为理想塑性材料的弹塑性准则。在线弹性阶段，显然按虎克定律，应力和应变之间只通过变形模量（以弹模和泊松比表示）发生关系。变形模量的大小反映物体抗变形能力的大小。变形模量越大，发生的变形（位移）就越小，反之则反。变形模量在FLAC3D中又分为剪切模量和体积模量。前者反映抗剪切变形能力，后者反映抗体积变形能力。因此，如果基于经验知识，发现位移过大或过小，可通改变模量值来调整。而在塑性阶段，应力不再简单地根据应变按弹模和泊松比来调整，还要根据流动法则来确定（拉破坏服从关联流动，剪破坏服从非关联流动），这时，粘聚力C、内
14
FLAC3D中一些问题的讨论
