FLAC3D中一些问题的讨论

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