1、目 录1 FLAC3D 的固流耦合计算模式 -12 FLAC3D 固流耦合学习小结 -53关于流固耦合的问题 -64也谈采用 FLAC3D 对地下采矿的模拟 -85 FLAC3D 本构模型开发 -86 FLAC3D 自定义本构模型 -117数值计算中初始应力场的模拟 -138 FLAC3D 应变分析 -139 FLAC3D 的调参 -1410开采沉陷垂直剖面等值线的生成 -1511 FLAC3D 的应变硬化软化模型 -1612 FLAC3D 的塑性流动格式 -1713 FLAC3D 的动画制作 -1714地下连续墙基坑开挖支护 -1815一个汇的小例子 -2116用 3DEC 生成岩体随机节理
2、网络 -2317固结小算例 -24FLAC3D 的固流耦合计算模式英文原文 -2611FLAC3D 的固流耦合计算模式http:/ 的计算模式中是否需要做孔压分析取决于是否采用 config fluid 命令。1 无渗流模式（不使用 config fluid）即使不使用命令 config fluid，仍然可以在节点上施加孔压。这种模式下，孔压将保持为常量。如果采用塑性本构模型的话，材料的破坏将由有效应力状态来控制。节点上的孔压分布可由 initial pp 命令或 water table 命令来设定。如果采用 water table 命令，由程序自动计算水位线以下的静水孔压分布。此时，必须施加
3、流体密度（water density）和重力（set gravity） 。流体密度值和水位位置可以用命令print water 显示。如果水位线是由 face 关键字来定义的，则可用命令 plot water 命令显示水位。这两种情况，单元的孔压都由节点孔压值平均求出，并在本构模型计算中用作有效应力。这种计算模式下，体积力中不反映流体的出现：用户必须根据水位线以上或以下相应地指定干密度和湿密度。使用命令 print gp pp 和 print zone pp 可分别得到节点或单元孔压。plot contour pp 命令可绘出节点孔压云图。2 渗流模式（使用 config fluid）如果使用
4、命令 config fluid，则可进行瞬时渗流分析，孔压改变和潜水面的改变都可能出现。在 config fluid 模式下，有效应力计算（静态孔压分布）和非排水计算均被执行。除此之外，还可进行全耦合分析，这种情况下，孔压改变将使固体产生变形，同时体积应变反过来影响孔压的变化。如果采用渗流模式，单元孔压仍由节点孔压平均求出。但这种模式，用户只能指定干密度（不论是水位以上还是以下） ，因为 FLAC3D 将流体的影响考虑到了体积力的计算中。采用渗流模式时，渗流模型必须施加到单元上，使用命令 model fl isotropic 模拟各向同性渗流，model fl anisotropic 模拟各向
5、异性渗流， model fl null 模拟非渗透物质。注意，力学模型为空的单元并不代表渗流模型为空。流体性质（参数）可施加到单元或节点上。各向同性渗透率、孔隙率、比奥系数和非排水热系数等单元流体性质由命令 property 施加。对于各向同性渗流，渗透率通过 perm 关键字赋予。对各向异性渗流，渗透率的 3 个主值采用关键字 k1,k2,k3 赋予，主方向由关键字 fdip,fdd,frot 确定。渗透率的主方向服从右手系统。fdip 和 fdd 分别为 k1 和 k2 确定的平面的倾向和倾角。frot 为 k1 轴和倾角矢量的旋转角。如果不特别指定，比奥系数默认为 1，孔隙率默认为 0.
6、5。节点的渗流性质由命令initial 指定。这些性质包括流体重度、流体体积模量、比奥模量、流体抗拉强度和饱和度。每种性质在空间上都可以变化。流体重度也可以用 water 命令给出。在渗流模式里，有必要知道可压缩性被定义在以下两种参数中：（1）比奥系数和比奥模量；（2）流体体积模量和孔隙率。第一种参数表征的是固体颗粒的可压缩性（对不可压缩颗粒，比奥系数设为 1） 。对第二种参数，固体颗粒被认为是不可压缩的。单元属性可由命令 print zone property 显示，节点属性由 print gp 命令显示。流体重度，如果随着水位位置被确定，则可由 print water 命令显示。渗流性质可
7、由命令 plot bcontour property 显示。对于各向异性渗流，渗透率的各球形分量可通过使用单元的属性关键字kxx,kyy,kzz,kxy,kxz,kyz 来显示（注意，这些球形分量不可被直接初始化） 。初始节点孔压分布的施加对于渗流模式和非渗流模式都是一样的（如，要么用 initial 2pp 命令或用 water table 命令） 。在指定节点可用命令 fix pp 或 free pp 对孔压固定或释放。流体涌入或渗漏或可由命令 apply 施加。渗流计算由命令 set fluid 和 solve 控制。如，set fluid on 或 off 命令开启或关闭渗流计算模式。
8、具体使用开启或关闭模式取决于渗流分析的耦合程度。渗流分析结果以下面这些命令给出。命令 print gp pp 和 print zone pp 分别给出节点和单元孔压。节点和单元孔压历史可由命令 histroy gp pp 和命令 history zone pp 进行监测。对于瞬时计算，孔压与时间的关系可由命令 history fltime 监测。命令 plot contour pp 绘出节点孔压云图。命令 plot contour saturation 绘出饱和度云图。命令 plot fluid 绘出流量矢量图。渗流模式的所有信息由命令 print fluid 命令给出。FISH 还提供了一些
9、渗流变量。其中一个与节点有关的变量 gp_flow，只能通过 FISH 函数使用。该变量描述了通过节点的净流入或流出量。因为可以提供一个系统总的流入或流出量，这些流量的统计在孔压固定的边界是很有用的。渗流边界条件，初始条件FLAC3D 默认为不透水边界，即认为所有节点上的孔压随着从邻近单元流入或流出的量发生自由变化。可以使用命令 fix pp 将节点上的孔压设为 “自由” ，也可使用 free pp 使节点上的孔压“固定” 。如果孔压固定，流体可以在外边界上流入或流出节点。下面总结这两种边界条件的影响：1，孔压自由这是默认的不透水边界条件。节点与外界之间不发生流量交换。系统根据当前饱和度值和流
10、体是否形成涡凹现象来计算压力和饱和度变化。2，孔压固定这是一种流体通过外界流入或流出的边界条件。如果设定孔压为 0，饱和度才可能变化。否则，饱和度被设为 1（FLAC3D 假设孔压只在完全饱和材料中存在） 。孔压不能被固定在低于拉力极限的值，如果出现这种情况，FLAC3D 会将其设定到拉力极限值。如前所述，边界条件不是任意的。FLAC3D 在进行计算前会“检查”并“修正”这些条件。可使用 fix pp 命令将孔压固定在某个值，也可在外边界或内边界上使用命令 apply pp。如果边界条件被用于一个非表面节点，则必须加关键字 interior。apply 命令具有可以用“历史”命令进行监测的优点
11、。渗流边界条件可以通过 apply 命令用在单个或部分节点、单元面或单元上。命令 apply pwell 为边界节点指定了一个流入或流出井。如果加上 interior 关键字，则该条件用于内部节点。命令 apply discharge 和 apply leakage 为边界单元的表面分别指定了涌出和渗漏边界条件。命令 apply vwell 为指定区域内的单元提供一个流速。这些边界条件除了 apply leakage 外，均可使用 history 监测命令。具有固定孔压节点就好像是流入源或流出源。没有直接的命令显示这些节点的流入或流出量。但可通过 FISH 变量 gp_flow 来记录。孔压的
12、初始分布，孔隙率，饱和度和流体属性可通过命令 initial 或 property 施加。如果还加了重力，则孔压初始分布应与重力梯度，水的重度和节点饱和度和孔隙率相容。如果这些初始分布不相容，则计算开始时所有单元中将出现流体流动。因此，应在模拟开始时设一定的计算步来检验初始条件是否相容。如果模型中含有接触面，有效应力将沿着这些接触面进行初始化（即：在节点应力初始化时，认为接触面应力包含孔压） 。water table 命令将包含沿着接触面的孔压，这是因为定义在单元节点上的孔压也在接触面节点上。如果接触面的上下两面连在一起，在没有阻力时，将发生穿越接触面的流体流动。但程序不对沿着接触面的流体流动
13、（裂隙流）进行计算。3单渗流与渗流耦合问题FLAC3D 既能进行单渗流分析，也能进行固流耦合分析。耦合分析可由 FLAC3D 内置力学模型完成。但要注意，渗流模型中的空单元并不是力学空单元。必须用命令 model fl_null 给单元赋予流体空属性。对于耦合过程，FLAC3D 提供了几种计算模式。其中之一是假设孔压一旦被赋予便不再改变。该方法并不要求任何额外空间存储计算过程。除此之外涉及到渗流的计算模式都要求使用命令 config fluid。命令 model fl_iso 使所有单元中都能发生渗流。不同的耦合计算模式在下面讨论。一般情况下，在能跟所模拟问题的物理过程相似的情况下，应使用尽可
14、能简单的模式。计算模式的选择根据以下几个方面确定。时间比例对所需模拟的渗流或耦合问题用 FLAC3D 估计与涉及的不同进程相关的时间比例是非常有用的。对有关研究问题的时间度量和扩散性的认识有助于估计最大网格宽度、最小区域尺寸、时步大小和计算可行性。如果不同进程的时间比例相差太大，则很可能采用一种简单的（非耦合）方法。时间比例可用特征时间给出。以下这些由量纲分析得出的定义，都是基于解析的连续源理论表达式。它们可用于得出 FLAC3D 分析的大致时间比例。力学过程特征时间、流体扩散过程特征时间流体扩散率FLAC3D 中使用了取决于控制过程的储水系数的几种形式：流体存储系数、地下潜水相存储系数、弹性
15、存储系数以上定义，有几点特性值得注意：（1）因为 FLAC3D 中显式的时步对应于最小区域中信息从一个节点传到下一节点所需要的时间，时步的大小可用计算特征时间公式中特征长度的最小区域来估计。重要的是注意FLAC3D 中在用流体扩散率（即使是在耦合模拟中）计算显式流体时步。因此，时步的大小可用特征长度的最小区域尺寸来估计。（2）在饱和流体问题中，简化的体积模量不但导致时步的增加，同样导致到达稳定状态时间的增加，所以总步数增加，该总步数可用模型和最小区域的特征长度来估计。（3）在部分饱和流体流动问题中，可通过调整流体体积模量加速收敛以趋于稳定状态，但要注意不可将体积模量减小太多以至产生数值不稳定。
16、数值稳定条件能由流体储量在一个特征长度区域的高度上必须保持低于地下潜水储量的要求推导而出。（4）为避免扩散问题中的边界效应，模型的特征长度必须大于某个尺度。同样，最小模拟时间由某个关系式控制。（5）在耦合流体问题中，实际扩散率由流体刚度与岩土介质的刚度比来控制。完全耦合模拟方法的选择用 FLAC3D 进行完全耦合的准静态固流耦合分析通常要耗费大量时间，且有时候并不必要。很多情况下，可使用不同程度的非耦合方法简化分析并加快计算速度。下面的例子给出了对应于流固耦合的不同水平的 FLAC3D 模拟方法。选择计算方法时有 3 个主要的因素需要考虑：（1）模拟时间比例和扩散过程的特征时间；（2）耦合过程
17、中强制扰动特性；（3）流固刚度比。时间比例首先通过从扰动的开始阶段计算时间来考虑时间比例因素。定义分析所需要的时间（模拟时间） ，对应于耦合扩散过程的特征时间。4短期行为（不排水）如果对应于耦合扩散特征时间，分析所需时间非常短，在模拟结果中流体流动的影响几乎可以忽略不计，则可采用不排水模拟（config fluid, set fluid off） 。数值模拟中不涉及真实的时间，但如果给流体体积模量一个实际值，则体积应变将导致孔压的变化。长期行为（排水）如果分析所需时间大于耦合扩散特征时间并在模拟时间到达时排水，则孔压场可不耦合到立场中。稳定状态的孔压场可用单纯流动模拟确定（set fluid
18、on, set mech off） （不代表扩散率） ，然后力学场可通过在设置流体模量为 0 的力学模式中（set mech on, set fluid off）计算到平衡状态获得。 （严格说，这种方法仅对弹性材料有效，因为塑性材料力学行为是与路径有关的） 。另外一种描述时间比例的方法是不排水行为和排水行为。严格地说，不排水表示与模型外界无流体交换。排水则是与模型外界有完全的流体交换，这就意味着流体压力能在各处达到平衡。由于不排水试验通常所需时间很短，而排水试验则需要很长的时间以使多余的流体压力消散，因此， “排水” 和“ 不排水”这两个词常分别和“短期行为”和“长期行为”联系在一起。在现场，
19、 “短期行为” 通常意味着流体流动可以忽略，而“长期行为”则意味着几乎所有压力降都变为 0（这需要一个很长的过程） 。注意在无渗流计算模式（不使用 config fluid）和短期行为模式（使用 config fluid,set fluid off）的模拟过程中，由于施加的孔压发生变化而产生的总应力的修正并不是通过程序在内部执行。但是，孔压增量可用 FISH 函数监测，并用于减小循环到力学平衡前的总法向应力。如果地下水位已在网格内部移动，同时需要调整饱和及非饱和的质量密度。强制扰动到耦合进程的特性将扰动强加到完全耦合的固流系统可能导致流体流动边界条件和力学边界条件的改变。如，流向位于层间含水层
20、内的井的瞬时流体流动是由井内孔压变化引起的。作为公路路堤建设成果的饱和地基的固结是由路堤高度确定的力学载荷控制。如果扰动是由于孔压的变化，很可能流体流动进程可不与力学过程耦合。如果是固体产生的扰动，非耦合的程度取决于流固刚度比。刚度比相对刚度比对用于解决固流耦合问题的模拟方法有重要影响：相对刚性岩土介质（相对刚度比远小于 1）如果岩土介质骨架刚度很大（或者流体是高压缩性的）且相对刚度比很小，孔压的扩散方程可以不耦合，因此扩散率有流体控制。建模方法取决于流体或固体扰动的力学机制：（1）在固体控制的模拟中，孔压可假设保持不变。在弹性模拟中，固体表现的力学行为好像流体不存在；但在塑性分析中，孔压压力
21、的出现可能导致破坏。这种模拟方法在边坡稳定性分析中使用。（2）在孔压控制的弹性模拟中（如由于流体被挤出导致的沉降） ，体积应变不显著影响孔压场，且流体的计算可独立进行（set fluid on, set mech off） （这种情况下，扩散率是精确的，因为对于相对刚度比小于 1，总压缩系数等于流体扩散率） 。一般地，孔压变化会影响应变，且这种影响可以通过随后在力学模式中将模型循环到平衡状态来加以研究（set mech on, set fluid off） 。相对柔性岩土介质（相对刚度比远大于 1）如果岩土介质骨架刚度很小（或流体不可压缩） ，且相对刚度比很大，则岩土骨架控制系统扩散率的耦合。
22、模拟方法也取决于控制的力学机制。（1）在力学控制的模拟中，计算可能比较耗时。FLAC3D 的显式时步可能会很小，为增加时步，可减小流体模量。注意，流体模量不应该设置得大于流体的实际值。5（2）在多数孔压控制系统的实际例子中，经验表明，孔压场和力学场的耦合是微弱的。如果介质是弹性的，可用单纯流动模式（set mech off, set fluid on）计算，然后在单纯力学模式（set mech on, set fluid off）中计算到平衡。必须注意为保持系统的扩散率（以及特征时间比例） ，流体模量必须在流体计算阶段调整到某一值，且在力学计算阶段为 0，以防通过体积应变进一步调整。对于建模方
23、法的选择，建议根据以下步骤进行。首先，对于特定的问题条件和特性，确定扩散进程的特征时间，且将此时间同所关注的实际时间进行对比。其次，考虑对于系统的扰动是由孔压控制还是由固体控制。第三，确定流体刚度对固体基质骨架的刚度比。最后，基于这 3 个因素综合考虑选择合适的建模方法。在建模时需注意：（1）为建立无地下水流动的有效应力分析的初始条件，用 water table 或 initial pp 命令，或者用 FISH 函数建立稳定状态的孔压分布。指定正确的位于地下水位以下区域的湿密度和地下水位以上区域的干密度。（2）为建立地下水流动的有效应力分析的初始条件，如果地下水位位置未知，用 initial
24、命令或 FISH 函数建立稳定孔压分布，或者指定 set fluid on 和 set mech off 并逐步计算到稳定状态。将流体模量设为一个较小的值以加快部分饱和系统的收敛速度。注意设定的流体模量值应满足数值稳定。（3）为建立孔压驱动分析的初始条件，如果地下水位位置未知，用 initial 命令或 FISH 函数建立稳定孔压分布，或者指定 set fluid on 和 set mech off 并逐步计算到稳定状态。将流体模量设为一个较小的值以加快部分饱和系统的收敛速度。注意设定的流体模量值应满足数值稳定。（4）非流固耦合方法推荐用于孔压驱动系统，且在相对刚度比远大于 1 时谨慎使用。注
25、意在单纯流动分析阶段调整流体模量的值以满足耦合扩散率是正确的。（5）完全耦合分析中，注意对于相对刚度比远大于 1 时的情况，如果流体模量调整得低于某个值，时间响应将会接近于无限大的流体模量的时间。2FLAC3D 固流耦合学习小结通过几天来对 FLAC3D 手册固流耦合部分的翻译和学习，对 FLAC3D 在流体方面的应用有了大致的了解。其计算模式主要分为几种情况：无渗流、非耦合方法以及各种不同程度的固流耦合方法。这些方法都有各自的适应条件，应根据需要模拟问题的具体情况加以选择。总体来说，固流耦合涉及到了两种不同形态物质的力学行为研究甚至二者的相互作用，这比单一介质的分析复杂多了，其参数也比单一介
26、质分析要多得多。众所周知，数值模拟的一个最大难题便是参数的确定，参数越多，各种不可预测因素就越多，直接影响到数值计算的有效性和准确性。尤其是对于岩土工程，各种不确定因素太多，很多人由此对数值模拟结果产生了怀疑，认为那是花架子，放在项目报告里只有唬人的作用。中尉对此不敢苟同。任何技术的发展均遵从实践理论实践的过程，经验以及由经验得出的方法确实在很多工程中起到了不可抹杀的作用，但经验如果上升不到理论，就永远是经验，不能是“放之四海皆准” 的真理，真理虽然没有绝对，但在一定时期内它是推动整个技术进步的源动力。没有理论，就不能把实践上升到更高程度的实践。数值模拟目前在所有领域内均取得了广泛的应用便是其
27、具有强大生命力的最好证明。诚然，岩土工程领域是有很多不确定因素，数值模拟结果的不尽人意应该是多种因素造成的，我们为何不检讨获取数值模拟参数的现场观测手段、实验室手段，而去质疑这种方法本身？在这几天的学习过程中，还有一点对翻译的额外体会。以前也看过很多中国人直接翻6译自老外的教材，以计算机教材最多，看的过程中总感觉翻译的痕迹太重，可当自己翻译出来时，才发现无论怎么修改，润色，仍然晦涩不堪。中尉的一位朋友对这个问题曾有评论，中尉很以为然，援引之：中文的表达方式惯用短句，层层推进，像一条河，后浪推前浪，引人入胜。而英文的句子像一棵棵枝繁叶茂的大树，句子套句子，得透过无数的从句、定语去寻找主干。翻译者
28、或许是在两种思维方式中穿梭得迷了路，翻译过来的句子披着中文的外衣却完全是英文的表达方式，看到主语之后就是一大堆莫名其妙的短句、分号、连接符，最后找到谓语、宾语时已经不知所云。3关于流固耦合的问题http:/ rad_tunnel=(-1)*rad_tun width_m=1 dep_r1=20 dep_r11=(-1)*dep_r1 dep_r2=15 dep_r22=-dep_r2 depth_r=-(dep_r1+dep_r2) wid_add=wid_x-dep_r1 h_add=-dep_r2-dep_r1-dep_r1-h_model end para_model gen zon c
29、yl p0 0 0 depth_r p1 add rad_tun 0 0 p2 add 0 width_m 0 p3 add 0 0 rad_tun size 6 1 8 rat 1.2 1 1 group tunnel gen zon radcyl p0 0 0 depth_r p1 add dep_r1 0 0 p2 add 0 width_m 0 p3 add 0 0 dep_r1 assign names to groups of zones group rock1 range group tunnel not gen zone brick p0 0 0 dep_r22 p1 add
30、dep_r1 0 0 p2 add 0 width_m 0 p3 add 0 0 dep_r2 -fluid flow model- model fl_iso range group rock1 prop perm 1.0e-9 poro 0.3 biot_c 1 model fl_iso range group tunnel prop perm 1.0e-9 poro 0.3 biot_c 1 model fl_iso range group rock2 prop perm 1.0e-7 poro 0.3 biot_c 1 ini fmod 2e9 ini fdensity 1e3 ; -
31、mechanical model - model elas ini dens 1500 prop bulk 5e10 shear 3.0e9 range group rock1 prop bulk 5e10 shear 3.0e9 range group tunnel prop bulk 5e10 shear 3.0e9 range group rock2 ; assign boundary conditions fix x range x 69.9 70.1 fix x range x -70.1 -69.9 fix z range z -100.1 -99.9 fix y ; - init
32、ial stress state - ini pp 2.0e5 grad 0 0 -10000 fix pp range z -0.1 0.1 fix pp range z -100.1 -99.9 ; set fluid biot on ;该处取消，则产生超空隙水压力 set fluid off set grav 0 0 -10 set large 8solve plot sketch plot add cont ppqinjianshe个人愚见：水应力不需要乘以 Porosity4也谈采用 FLAC3D 对地下采矿的模拟以下转自中尉发到 SIMWE 论坛 FLACFLAC3D 版面的帖子地
33、下采矿，以层状的煤层开采为例，目前很多都是采用全部垮落法，煤层开采后，采空区岩层从弯曲下沉到出现裂隙并最终垮落，按开采沉陷“上三带”的说法，就是采空区上覆岩层形成了垮落带、裂隙带和弯曲下沉带。对于这样一个会出现完整岩层垮落成松散矸石的过程，采用基于连续介质的 FLAC 来模拟，本身就并不十分恰当！因为采空区上覆岩层垮落后，已经不再是连续介质。如果采用 FLAC 模拟，在大范围的开采（长壁开采一般范围都很大）后，绝大多数情况下是算不到平衡状态的，但有可能出现稳定的塑性流动状态，从计算原理上来说，这是非常正常的，而不是什么参数选取的问题。采矿工程不同于其他地下岩土工程之处就在于，它不对采空区进行任
34、何支护处理，而其他工程不论是地下厂房，隧道开挖等都会进行及时支护，保持了围岩的完整性和稳定性。虽然从开采沉陷角度来说，岩层移动最终会有个稳定过程，但由于 FLAC 本身不能模拟岩层垮落后形成的散体对上覆岩层起到一定的支撑作用，因此如果不对采空区进行另外的处理，一般是不会计算到系统默认的平衡状态（最大不平衡力比率达到 1105）的，除非开采范围很小。于是，论坛里出现了很多关于用 FLAC 模拟地下开采后对采空区进行处理方法的讨论，如对采空区进行弹性充填等。对于这些方法，本人也曾做过一些试验，效果并不好，毕竟采空区垮落不是充填，这样做有一种“亡羊补牢”的感觉。如果非要用 FLAC 来模拟地下开采，
35、本人建议研究重点放在裂隙带及以上的弯曲下沉带到地表这段范围，因为这部分岩层还比较完整、连续，对于涉及到采场覆岩破坏结构方面的研究，建议还是采用离散元进行。以上观点，是。己应用 FLAC 做开采沉陷两年来的体会，如有不妥，欢迎各位指正并交流！5FLAC3D 本构模型开发FLAC3D 基于显式有限差分格式，计算过程首先调用运动方程，由初始应力和边界条件计算出新的速度和位移，然后，由速度计算出应变率，进而根据本构方程获得新的应力或力。显式有限差分法计算流程如图。9FLAC3D 本构模型的主要功能就是根据应变增量（应变率）返回新的应力张量。 FLAC3D 采用面向对象的语言 C+编写。C+语言的特点是
36、采用面向对象方法，使用类来代表对象进行编程。与对象有关的信息被封装在类中，这些信息在类的外部是不可见的。与对象的通信通过成员函数操作封装数据来完成。一个新的对象类型可以从基类派生而来，基类的成员函数也可被派生类的同名函数覆盖。这种方式便于程序的模块化设计。 FLAC3D 中的所有本构模型都是以动态链接库文件（.dll 文件）的形式提供，自定义本构模型也不例外。动态链接库文件必须采用 VC+6.0(SP4)或更高版本编译得到。用 VC+编写自定义 FLAC3D 本构模型的过程主要包括：基类、成员函数的定义，模型注册，模型与FLAC3D 间的数据传递以及模型状态指示。基类的描述基类提供了实际本构模
37、型的框架，用户自定义模型类都继承于该类。基类命名为ConstitutiveModel，由于它定义了一系列“虚”（语法上以 0 为标志）的成员函数，因此又被称为“抽象”类。这也是说不能创建该类的任何对象，而其它继承自该类的类的对象必须给出实际的成员函数来代替基类的“虚”函数。基类的头文件（.h 文件）中定义了这些“虚”的成员函数。模型成员函数 ConstitutiveModel 类有很多成员函数，主要有： const char *Keyword() 返回一个指向包括模型名称的字符数组的指针，以便用户在 FLAC3D 中使用 MODEL 命令时，C+ 能够识别。 const char *Name(
38、) 返回一个指向包含模型名称的字符数组的指针，以识别用户使用如 PRINT Zone 等命令。 const char *Properties() 返回一个指向包含模型力学参数名称的字符串数组的指针，并用一个空指针代表字符串数组的结束。这样程序能识别用户输入的 PROPERTY 命令。 const char *States() 返回一个包含单元状态名称的字符串数组的指针，并用一个空指针代表字符串数组的结束。这样程序就能输出或显示用户自定义的模型内部状态（如：塑性流动、屈服、受拉） 。 const char *Run(unsigned uDim, State *ps) 该函数在FLAC3D 计算循
39、环的每一步对每一个子单元使用。模型根据应变增量对应力张量进行更新。这个函数是本构模型的核心，不同的本构模型通过不同的 Run()函数得到不同的应力张量。所有本构模型的 Run()函数都继承自 ConstitutiveModel 类的 Run()函数。模型注册每一个用户自定义本构模型有其自己的模型名称、力学参数名称和状态指示器。FLAC3D 通过调用模型对象的静态全局实例获得用户定义模型的信息并使用结构子（constructor）将新模型注册到模型列表中。一个特定模型只有一个静态注册产生，方便将其放到模型的 C+源码中，这样当相应的.dll 文件被装载时，模型被注册。数据传递在 FLAC3D 和
40、用户自定义模型间的最重要的连接就是成员函数 Run(unsigned nDim, State *ps)，它在模型的计算循环时计算其力学响应。State 结构被用来传输信息和生成模型。同时，在非线性模型中，Run() 函数也用来传递模型的内部状态。状态指示 FLAC3D 中的单元是由四面体子单元所组成，每一个四面体有记录其当前状态的成员变量。该成员变量共有 16 位，能够代表最多 15 种不同的状态。对于用户定义本构模型，用户可以命名一种状态并为其分配特定的位。头文件和源文件在 C+中，主要有头文件（.h）和 C+源文件（.cpp）两种文件类型，头文件作为一种包含功能函数、数据接口声明的载体文件
41、，用于保存程序的声明(declaration)，是用户应用程序和函数库之间的桥梁和纽带。而源文件（.cpp）用于保存程序的实现(implementation)，10是具体实现代码。先在 VC+6.0 环境中建立一个工程（project）文件 jdmohr.dsw，工程文件中包含了空的头文件和源文件。 FLAC3D 的本构模型需要用到如表所示的一些头文件和源文件： FLAC3D 本构模型的 C+头文件和源文件名称 功能 CONMODEL.h 声明与本构模型交换数据的变量数据类型 State，保存结果的数据类型ModelSaveObject 和本构模型基类 ConstitutiveModel ST
42、ENSOR.h 声明存储对称张量（应力或应变张量）STensor 类 AXES.h 声明一个与坐标轴变换有关的数据类型 Axes 和相关函数 CONTABLE.h 声明了一个 ConTableList 类，该类为本构模型提供了一张表，该表用来存储模型单元或节点 ID 号 USERMODEL.h 用户自定义本构模型派生类的声明 USERMODEL.cpp 用户自定义本构模型的由应变增量获得应力增量的实现在编写自定义本构模型的 C+文件时，需要在 FLAC3D 自带的莫尔 -库仑模型基础上增加新的材料参数变量。如中尉编写的节理岩体损伤莫尔-库仑模型增加如图所示的材料参数。 FLAC3D 本构模型的
43、主要功能就是根据应变增量得到新的应力增量，应力应变关系表达式很可能为矩阵表达式，在计算时，可能需要对矩阵求逆。对于逆矩阵的求取，可根据逆矩阵定义或矩阵的初等行变换来进行，建议采用一种在数值计算方法中常用的全主元高斯约当（Gauss-Jordan）消去法。该方法具有数值稳定而且精度较高的特点。由于岩体的柔度张量矩阵元素都非常小（10 的-7、-8 次方级） ，在进行乘法计算时，将产生极大的舍入误差，而高斯约当消去法的消去过程是在全矩阵中选主元（绝对值最大的元素）来进行，故可使舍入误差对结果的影响减小到最小。将编写好的头文件 jdmohr.h 和源文件jdmohr.cpp 导入到工程文件（.dsw
44、）中，通过编译（compile）和链接（link ） ，形成了节理岩体损伤莫尔-库仑本构模型的动态链接库文件 jdmohr.dll，将其复制到 FLAC3D 的安装目录下，如图11程序的调试是编程过程中不可缺少的重要步骤，尤其是对于数值计算，需要实时监控计算流程和各种变量，以防出现死循环和结果严重失真。在 VC+中可设置将 FLAC3D 中的 EXE 文件路径加入到程序的调试范围中，并将自定义的 FLAC3D 的.dll 文件加入到附加动态链接库(Additional DLLs)中，然后在.cpp 文件里的 Initialize()或 Run()函数中设置断点，进行调试。使用自定义模型时，首先在命令流中使用命令 CONFIG cppudm 对 FLAC3D 进行配置，使其能接收动态链接库文件，然后通过 M
