FLAC3D中一些问题的讨论
目录
1．FLAC3D的固流耦合计算模式--------------------------------------------------------------------------------- 1 2．FLAC3D固流耦合学习小结 ------------------------------------------------------------------------------------ 5 3．关于流固耦合的问题-------------------------------------------------------------------------------------------- 6 4．也谈采用FLAC3D对地下采矿的模拟 ---------------------------------------------------------------------- 8 5．FLAC3D本构模型开发 ------------------------------------------------------------------------------------------ 8 6．FLAC3D自定义本构模型 ------------------------------------------------------------------------------------- 11 7．数值计算中初始应力场的模拟----------------------------------------------------------------------------- 13 8．FLAC3D应变分析 ----------------------------------------------------------------------------------------------- 13 9．FLAC3D的调参 ---------------------------------------------------------------------------------------------------- 2 10．开采沉陷垂直剖面等值线的生成-------------------------------------------------------------------------- 3 11．FLAC3D的应变硬化软化模型 ------------------------------------------------------------------------------- 4 12．FLAC3D的塑性流动格式-------------------------------------------------------------------------------------- 5 13．FLAC3D的动画制作 -------------------------------------------------------------------------------------------- 5 14．地下连续墙基坑开挖支护------------------------------------------------------------------------------------ 6 15．一个汇的小例子------------------------------------------------------------------------------------------------- 9 16．用3DEC生成岩体随机节理网络 ------------------------------------------------------------------------- 11 17．固结小算例 ----------------------------------------------------------------------------------------------------- 12 FLAC3D的固流耦合计算模式英文原文----------------------------------------------------------------------- 14
1．FLAC3D的固流耦合计算模式
https://www.wendangku.net/doc/8514126524.html,/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．关于流固耦合的问题
https://www.wendangku.net/doc/8514126524.html,/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 &dV al) 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) 该函数在使用命名SA VE或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、内摩擦角、剪涨角和抗拉强度都要参与进来，分析起来就很复杂了。但在塑性破坏阶段，我们更多地不是关心应力值大小，而是关心破坏区域的大小。显然，破坏区域的大小只与屈服函
数中的粘聚力、内摩擦角和抗拉强度有关。通过对最大主应力和最小主应力图上表示的莫尔库仑模型分析，增大粘聚力、内摩擦角和抗拉强度都能使破坏区域减小。但这时还有一个参数敏感性的问题，即：哪个参数对结果影响最大。这可通过调整某一个值，同时固定另外两个值来进行试验。
最后要说明的一点是：上述参数的调整并不是任意的！应根据试验结果进行合理的调整。比如对于变形模量值，它是由弹模和泊松比求出来的。反之，弹模和泊松比也能够由变形模量表示出来。注意到，泊松比的取值是有范围的，一般不超过0.5！因此，变形模量的赋值应该满足泊松比的取值范围。其他如抗拉强度也不能任意加大，岩体本身就是一种不抗拉材料！
10．开采沉陷垂直剖面等值线的生成
以下转自中尉发到SIMWE论坛FLAC\FLAC3D版面的帖子
在应用FLAC3D对煤矿开采沉陷的研究中，经常需要得到某一垂直剖面上的位移、应力的等值线分布。等值线可以由其他软件如tecplot或surfer软件作出，论坛里也有很多相关帖子，这里不再赘述。问题的关键在于要把垂直剖面上的单元或节点数据导出到数据文件中以便由绘图软件画出等值线。由于煤系地层均为层状，建模时一般采用逐层建模，各层单元的ID号就是连续的，如果要得到某一层或地表的位移等值线则可很方便地采用FISH遍历该层各单元，读出各单元数据就行了；但某一垂直剖面上的各单元的ID号并不是连续的，这时可以通过FISH中的内存变量mem()将该剖面上的所有单元连接起来形成一个链表，然后对该链表操作读出剖面上各单元的数据。以下分别是找出剖面上各单元的findsection函数和读出单元数据的readdata函数。
def findsection
head1 = null
p_z = zone_head
loop while p_z # null
if z_ycen(p_z)=310 ;该垂直剖面为y=310
new1 = get_mem(2)
mem(new1) = head1
mem(new1+1) = p_z
head1 = new1
endif
p_z =z_next(p_z)
endloop
end
findsection
def readdata
ad=head1
loop while ad#null
p_z=mem(ad+1)
p_smax=z_sig3(p_z) ;将剖面上所有单元的最大主应力放到p_smax中
ad=mem(ad)
endloop
end
然后就可用FISH编写读取p_smax变量的函数了。
下图为该例子的剖面云图和等值线。
11．FLAC3D的应变硬化软化模型
对于岩土体的常规静力计算，FLAC3D的莫尔库仑模型应该是用得最多的了，但这个模型在塑性阶段采用的是理想塑性模型，简单说，就是材料达到屈服后，其应力不再随应变而变化，塑性阶段的应变只有塑性应变。这跟一般的岩石力学试验得出的结果是不同的，岩石在很多情况下都表现出了应变软化的性质，即在屈服点以后，应力随着应变的增大而减小，在塑性阶段发生的应变包括弹性应变和塑性应变两部分。FLAC3D将这一现象体现在了一种本构模型中，称之为应变软化模型。同理，屈服点以后，应力随应变增大而增大的模型称为应变硬化模型。（其实在塑性力学中，材料只分为理想塑性材料和硬化材料。其定义为：随着塑性应变等的出现和发展，屈服面的大小和形状不发生变化的材料称为理想塑性材料，反之叫硬化材料。这里的硬化材料就包括上述软化性质）。
对于应变软化或硬化，FLAC3D的处理方法仍然是在莫尔库仑模型的理论框架下进行计算，
只不过需要定义一个随应变而变化的材料属性（粘聚力、内摩擦角和抗拉强度）函数，在材料达到屈服后的计算的每一时步，程序都根据这些函数弱化或强化这些材料属性，在每一步计算时，材料的本构关系仍然遵从莫尔库仑模型。这种模型的好处还在于可以直接从单元的材料属性中得出剪切塑性应变和拉伸塑性应变，这样就可以很方便地根据全应变张量得到单元的弹性应变（全应变减去塑性应变为弹性应变）。当然，该模型的难点也显而易见，那就是如何给出这个材料属性软化或硬化的函数。
12．FLAC3D的塑性流动格式
FLAC3D的二次开发大有可为！中尉的毕业论文主要进行了FLAC3D的二次开发工作，其中本构模型的开发是其中最重要的一项工作。仔细研读manual之theory and background,里面有各个模型的原理介绍，从塑性模型最简单的莫尔库仑模型入手，可以对其计算流程有个深入了解，为自定义本构模型奠定基础。
中尉根据manual翻译了FLAC3D中塑性流理论增量方程一节。
FLAC3D中的塑性流动遵守如下5条准则。(1) 破坏准则。即为屈服函数，即当塑性流发生时的应力组合关系。该函数为应力空间中的一个曲面，所有在该面之内的应力点保持弹性行为。(2) 应变分解准则。该准则将应变增量分解为弹性和塑性部分。(3) 弹性增量关系：应力增量为弹性应变增量的线性函数。(4) 流动法则该法则指定了塑性应变增量向量与塑性势函数的正交性。(5) 新的应力矢量分量满足屈服函数。
FLAC3D中，t+1时刻的弹性假设应力状态由t时刻的应力状态加上由总应变增量采用弹性应力－应变关系求出的应力分量增量获得。如果该弹性假设应力状态满足屈服函数，则用塑性流动法则得到位于屈服曲线上的应力。否则，新的应力状态由t时刻的弹性假设应力给出。莫尔－库仑模型
FLAC3D中的莫尔－库仑模型遵从莫尔－库仑准则（压剪破坏）和拉伸破坏准则。强度包络线上的应力点在剪破坏时服从非关联流动法则，拉破坏时服从关联流动法则。
FLAC3D在运行莫尔－库仑模型时，首先通过总应变增量由虎克定律计算出应力增量，再将应力增量加到应力分量上计算出弹性假设应力。主应力及其方向根据一点应力状态和主应力关系即可得到。
如果主应力满足组合屈服准则，则必发生剪切破坏或者拉伸破坏。对于剪切破坏，根据剪切流动法则得到新的应力增量并得到新的应力；对于拉伸破坏，根据拉伸流动法则得到新应力增量并得到新的应力。
如果应力点位于组合屈服包络线下方，则当前时步没有塑性流动发生，新的主应力由弹性关系给出。
说明：由于博客里的HTML文件粘不上公式，具体文件在中尉的电脑里共享。
13．FLAC3D的动画制作
FLAC3D的3.0版本有制作动画功能，其原理很简单，就是根据设置的特定时步记录一次作为关键帧，计算结束后，对关键帧之间进行插值，生成动画。
主要命令如下：
命令流开头：
set movie avi size 700 300 step 20 file dispvector3-2.avi；设置象素、记录步数和保存文件名movie star；记录动画开始
命令流末尾：
movie finish；记录动画结束
中尉做了一个边坡体下开采的岩层位移矢量动画。
14．地下连续墙基坑开挖支护
https://www.wendangku.net/doc/8514126524.html,/old/dispbbs.asp?boardid=23&Id=388
geotechman
我用壳结构单元模拟地下连续墙,梁单元模拟钢支撑,分步开挖,分步支护,怎么算出来的钢支撑轴力老是拉力,并且壳单元上的弯矩很小,为几吨,本基坑开挖深度为17米,计算结果与实际不符,请各位大侠帮我看一下!壳单元到底能不能模拟地下连续墙,该怎样模拟?局部坐标如何定义?请高人指点!! config fluid
gen zone brick p0 -1 0 -40 p1 15 0 -40 p2 -1 20 -40 p3 -1 0 0 size 16 20 40 gen zone brick p0 16 0 -40 p1 40 0 -40 p2 16 20 -40 p3 16 0 0 size 24 20 40
sel shell id 20 ran x 15.9 16.1 y 0 20 z -26 0 sel shell prop iso(28.5e9,0.17) thick 1 dens 2500 sel group wall ran id 20 sel node local xdir 0 0 1 ydir 0 1 0 ran sel group wall sel link attach zd=rigid xd=free yd=free xrdir=free yrdir=free zrdir=free ran sel group wall
group soil1 ran z -2 0 group soil2 ran z -20 -2 group soil3 ran z -22 -20 group soil4 ran z -25 -22 group soil5 ran z -40 -25
model mohr ran group soil1 def derive1 s_mod1=y_mod1/(2.0*(1.0+p_ratio1)) b_mod1=y_mod1/(3.0*(1.0-2.0*p_ratio1)) end set y_mod1=5.9e6 p_ratio1=0.4 derive1 prop bul b_mod1 shear s_mod1 cohe 4e4 fric 12
model mohr ran group soil2 def derive2 s_mod2=y_mod2/(2.0*(1.0+p_ratio2)) b_mod2=y_mod2/(3.0*(1.0-2.0*p_ratio2)) end set y_mod2=5.3e6 p_ratio2=0.36 derive2 prop bul b_mod2 shear s_mod2 cohe 2.8e4 fric 19
model mohr ran group soil3 def derive3 s_mod3=y_mod3/(2.0*(1.0+p_ratio3)) b_mod3=y_mod3/(3.0*(1.0-2.0*p_ratio3)) end set y_mod3=8.5e6 p_ratio3=0.3 derive3 prop bul b_mod3 shear s_mod3 cohe 1.4e4 fric 29
model mohr ran group soil4 def derive4 s_mod4=y_mod4/(2.0*(1.0+p_ratio4)) b_mod4=y_mod4/(3.0*(1.0-2.0*p_ratio4)) end set y_mod4=9.7e6 p_ratio4=0.25 derive4 prop bul b_mod4 shear s_mod4 cohe 1.3e4 fric 26
model mohr ran group soil5 def derive5 s_mod5=y_mod5/(2.0*(1.0+p_ratio5)) b_mod5=y_mod5/(3.0*(1.0-2.0*p_ratio5)) end set y_mod5=14.9e6 p_ratio5=0.24 derive5 prop bul b_mod5 shear s_mod5 cohe 1.0e4 fric 34
model fl_iso ran group soil1 prop perm 1e-12 poro 0.47 model fl_iso ran group soil2 prop perm 2.8e-11 poro 0.467 model fl_iso ran group soil3 prop perm 1.2e-10 poro 0.467 model fl_iso ran
group soil4 prop perm 1.9e-8 poro 0.475 model fl_iso ran group soil5 prop perm 2.1e-9 poro 0.44 ini fdensity 1e3 ini fmod 2e9 model fl_null ran sel group wall model fl_null ran group soil1
ini pp 2e4 grad 0 0 1e4 range x 0 40 y 0 20 z -40 -2 fix pp range z -2.1 -1.9
ini dens 1.88e3 ran group soil1 ini dens 1.95e3 ran group soil2 ini dens 1.96e3 ran group soil3 ini dens 1.81e3 ran group soil4 ini dens 1.97e3 ran group soil5
int 1 face ran plane norm 1 0 0 origin 15 20 20 dist 0.01 z -26 0 int 1 maxedge 1 int 1 prop ks 1e8 kn 1e8 fric 20 cohe 100 ini x add 1 ran x -1 15 y 0 20 Z -40 0 attach face ran x 15.9 16.1 y 0 20 z -40 -26
fix x ran x -0.1 0.1 fix x ran x 39.9 40.1 fix y ran y -0.1 0.1 fix y ran y 19.9 20.1 fix z ran z -40.1 -39.9
ini szz 0 grad 0 0 1.88e4 ran z -2 0 ini sxx 0 grad 0 0 1.22e4 ran z -2 0 ini syy 0 grad 0 0 1.22e4 ran z -2 0
ini szz 1.4e3 grad 0 0 1.95e4 ran z -20 -2 ini sxx -1e3 grad 0 0 1.17e4 ran z -20 -2 ini syy -1e3 grad 0 0 1.17e4 ran z -20 -2
ini szz 3.4e3 grad 0 0 1.96e4 ran z -22 -20 ini sxx 2e2 grad 0 0 1.176e4 ran z -22 -20 ini syy 2e2 grad 0 0 1.176e4 ran z -22 -20
ini szz -2.96e4 grad 0 0 1.81e4 ran z -25 -22 ini sxx -5.942e4 grad 0 0 0.905e4 ran z -25 -22 ini syy -5.942e4 grad 0 0 0.905e4 ran z -25 -22
ini szz 1.04e4 grad 0 0 1.97e4 ran z -40 -25 ini sxx -3.942e4 grad 0 0 0.985e4 ran z -40 -25 ini syy -3.942e4 grad 0 0 0.985e4 ran z -40 -25
apply sxx 0 grad 0 0 1.22e4 ran x 39.9 40.1 z -2 0 apply sxx -1e3 grad 0 0 1.17e4 ran x 39.9 40.1 z -20 -2 apply sxx 2e2 grad 0 0 1.176e4 ran x 39.9 40.1 z -22 -20 apply sxx -5.942e4 grad 0 0 0.905e4 ran x 39.9 40.1 z -25 -22 apply sxx -3.942e4 grad 0 0 0.985e4 ran x 39.9 40.1 z -40 -25 apply syy 0 grad 0 0 1.22e4 ran x 39.9 40.1 z -2 0 apply syy -1e3 grad 0 0 1.17e4 ran x 39.9 40.1 z -20 -2 apply syy 2e2 grad 0 0 1.176e4 ran x 39.9 40.1 z -22 -20 apply syy -5.942e4 grad 0 0 0.905e4 ran x 39.9 40.1 z -25 -22 apply syy -3.942e4 grad 0 0 0.985e4 ran x 39.9 40.1 z -40 -25 apply sxx 0 grad 0 0 1.22e4 ran y -0.1 0.1 z -2 0 apply sxx -1e3 grad 0 0 1.17e4 ran y -0.1 0.1 z -20 -2 apply sxx 2e2 grad 0 0 1.176e4 ran y -0.1 0.1 z -22 -20 apply sxx -5.942e4 grad 0 0 0.905e4 ran y -0.1 0.1 z -25 -22 apply sxx -3.942e4 grad 0 0 0.985e4 ran y -0.1 0.1 z -40 -25
apply syy 0 grad 0 0 1.22e4 ran y -0.1 0.1 z -2 0 apply syy -1e3 grad 0 0 1.17e4 ran y -0.1 0.1 z -20 -2 apply syy 2e2 grad 0 0 1.176e4 ran y -0.1 0.1 z -22 -20 apply syy -5.942e4 grad 0 0 0.905e4 ran y -0.1 0.1 z -25 -22 apply syy -3.942e4 grad 0 0 0.985e4 ran y -0.1 0.1 z -40 -25
apply sxx 0 grad 0 0 1.22e4 ran y 39.9 40.1 z -2 0 apply sxx -1e3 grad 0 0 1.17e4 ran y 39.9 40.1 z -20 -2 apply sxx 2e2 grad 0 0 1.176e4 ran y 39.9 40.1 z -22 -20 apply sxx -5.942e4 grad 0 0 0.905e4 ran y 39.9 40.1 z -25 -22 apply sxx -3.942e4 grad 0 0 0.985e4 ran y 39.9 40.1 z -40 -25 apply syy 0 grad 0 0 1.22e4 ran y 39.9 40.1 z -2 0 apply syy -1e3 grad 0 0 1.17e4 ran y 39.9 40.1 z -20 -2 apply syy 2e2 grad 0 0 1.176e4 ran y 39.9 40.1 z -22 -20 apply syy -5.942e4 grad 0 0 0.905e4 ran y 39.9 40.1 z -25 -22 apply syy -3.942e4 grad 0 0 0.985e4 ran y 39.9 40.1 z -40 -25 set gravity 0 0 -10 set large plot cont zdisp set fluid on hist unbal cyc 500
ini xdisp=0 ydisp=0 zdisp=0 fix pp ran x 16 40 y 0 20 z -2.1 -1.9 fix pp 0 ran x -0.1 16.1 y -0.1 20.1 z -3.1 -2.9 group excavate ran x 0 16 y 0 20 z -3 0 model fl_null ran group excavate model fl_null ran sel group wall model fl_null ran group soil1 model null ran group excavate set fluid on cyc 500
fix pp ran x 16 40 y 0 20 z -2.1 -1.9 fix pp 0 ran x -0.1 16.1 y -0.1 20.1 z -7.1 -6.9 group excavate
Flac3D命令--完整经典版
实例分析命令： 1. X ，Y ，Z 旋转 Shift+ X ，Y ，Z 反向旋转 Gen zone ……；model ……；prop ……（材料参数）；set grav 0,0,-9.81（重力加速度） plot add block group red yellow 把在group 中的部分染成红色和黄色 plot add axes black 坐标轴线为黑色；print zone stress% K 单元应力结果输出 ini dens 2000 ran z a b （设置初始密度，有时不同层密度不同）；ini ……（设置初始条件）；fix ……（固定界面） set plot jpg ；set plot quality 100 ；plot hard file 1.jpg 图像输出（格式、像素、名称） plot set magf 1.0视图的放大倍数为1.0；plo con szz z 方向应力云图 2. ini z add -1 range group one 群one 的所有单元，在z 方向上向下移动1m ；然后合并 命令 gen merge 1e-5 range z 0此命令是接触面单元合并成一个整体，1e-5是容差 3. (基坑开挖步骤)：Step 1: create initial model state （建立初始模型）Step 2: excavate trench （开挖隧道） 4. group Top range group Base not 定义（群组Base 以外的为）群组Top 5. plot blo gro 使得各个群组不同颜色显示 6. （两个部分间设置界面；切割法）：gen separate Top 使两部分的接触网格分离 为两部分；interface 1 wrap Base Top 在（Base 和Top ）这两部分之间添加接触单元；plot create view_int 显示，并创建标题view_int ；plot add surface 显示表面；plot add interface red 界面颜色红色 7. （简单的定义函数及运行函数）new ；def setup 定义函数setup ；numy = 8定义常 量numy 为8；depth = 10.0 定义depth 为10；end 结束对函数的定义；setup 运行函数setup 8. （隧道生成）上部圆形放射性圆柱及下部块体单元体的建立，然后镜像。 9. 模拟模型的材料问题时为什么要去定义某个方向上的初始速度？— 10. 渐变应力施加：apply nstress -1e6 gradient 0,0,1e5 range z 3.464,0 plane dip 60 dd 270 origin .1 0 0；施加法向应力：apply nstress -1e6 range plane dip 60 dd 270 origin .1 0 0 11. d ip dd 确定平面位置使用：（纠结） 12. p rint gp position range id=14647 输出节点坐标 13. a pply sxx -10e6 gradient 0 , 0, 1e5 range z -100 , 0在这个求解方程中，z 为变量，所以xx σ为：65=-1010+10xx z σ?? ；原点（0,0,0） 14. f ree x range x -.1 .1 z 6.9 10.1放松x=0 平面上，z=7，10 这一部分在x 方向的约 束（可以在此处产生破坏） 15. 体积模量K 和剪切模量G 与杨氏模量及泊松比v 之间的转换关系如下： =3(1-2v)E K G=2(1+v) E 16. 一般而言，大多数问题可以采用FLAC 3D 默认的收敛标准（或称相对收敛标准），即当体 系最大不平衡力与典型内力的比率R 小于定值10-5；（也可由用户自定义该值，命令：
FLAC3D 实例命令流1
第1部分命令流按照顺序进行 2-1定义一个FISH函数 new def abc abc = 25 * 3 + 5 End print abc 2-2使用一个变量 new def abc hh = 25 abc = hh * 3 + 5 End Print hh Print abc 2-3对变量和函数的理解 new def abc hh = 25 abc = hh * 3 + 5 End set abc=0 hh=0 print hh print abc print hh new def abc abc = hh * 3 + 5 end set hh=25 print abc set abc=0 hh=0 print hh print abc print hh 2-4获取变量的历史记录 new gen zone brick size 1 2 1 model mohr prop shear=1e8 bulk=2e8 cohes=1e5 tens=1e10 fix x y z range y -0.1 0.1 apply yvel -1e-5 range y 1.9 2.1 plot set rotation 0 0 45
plot block group def get_ad ad1 = gp_near(0,2,0) ad2 = gp_near(1,2,0) ad3 = gp_near(0,2,1) ad4 = gp_near(1,2,1) end get_ad def load load=gp_yfunbal(ad1)+gp_yfunbal(ad2)+gp_yfunbal(ad3)+gp_yfunbal(ad4) end hist load hist gp ydis 0,2,0 step 1000 plot his 1 vs -2 2-5用FISH函数计算体积模量和剪砌模量 new def derive s_mod = y_mod / (2.0 * (1.0 + p_ratio)) b_mod = y_mod / (3.0 * (1.0 - 2.0 * p_ratio)) end set y_mod = 5e8 p_ratio = 0.25 derive print b_mod print s_mod 2-6 在FLAC输入中使用符号变量 New def derive s_mod = y_mod / (2.0 * (1.0 + p_ratio)) b_mod = y_mod / (3.0 * (1.0 - 2.0 * p_ratio)) end set y_mod = 5e8 p_ratio = 0.25 derive gen zone brick size 2,2,2 model elastic prop bulk=b_mod shear=s_mod print zone prop bulk print zone prop shear 2-7 控制循环 New def xxx sum = 0 prod = 1 loop n (1,10)
FLAC3D 实例命令流1
第1部分命令流按照顺序进行2-1定义一个FISH函数 new def abc abc = 25 * 3 + 5 End print abc 2-2使用一个变量 new def abc hh = 25 abc = hh * 3 + 5 End Print hh Print abc 2-3对变量和函数的理解 new def abc hh = 25 abc = hh * 3 + 5 End set abc=0 hh=0 print hh print abc print hh new def abc abc = hh * 3 + 5 end set hh=25 print abc set abc=0 hh=0 print hh print abc print hh 2-4获取变量的历史记录 new gen zone brick size 1 2 1 model mohr prop shear=1e8 bulk=2e8 cohes=1e5 tens=1e10
fix x y z range y -0.1 0.1 apply yvel -1e-5 range y 1.9 2.1 plot set rotation 0 0 45 plot block group def get_ad ad1 = gp_near(0,2,0) ad2 = gp_near(1,2,0) ad3 = gp_near(0,2,1) ad4 = gp_near(1,2,1) end get_ad def load load=gp_yfunbal(ad1)+gp_yfunbal(ad2)+gp_yfunbal(ad3)+gp_yfunbal(ad4) end hist load hist gp ydis 0,2,0 step 1000 plot his 1 vs -2 2-5用FISH函数计算体积模量和剪砌模量 new def derive s_mod = y_mod / (2.0 * (1.0 + p_ratio)) b_mod = y_mod / (3.0 * (1.0 - 2.0 * p_ratio)) end set y_mod = 5e8 p_ratio = 0.25 derive print b_mod print s_mod 2-6 在FLAC输入中使用符号变量 New def derive s_mod = y_mod / (2.0 * (1.0 + p_ratio)) b_mod = y_mod / (3.0 * (1.0 - 2.0 * p_ratio)) end set y_mod = 5e8 p_ratio = 0.25 derive gen zone brick size 2,2,2 model elastic prop bulk=b_mod shear=s_mod print zone prop bulk print zone prop shear
FLAC3D快速入门及简单实例
FLAC3D快速入门 及简单实例 李佳宇编 LJY指南针教程
前言 FLAC及FLAC3D是由国际著名学者、英国皇家工程院院士、离散元的 发明人Peter Cundall博士在70年代中期开始研究的，主要面对岩土工程的通 用软件系统，目前已经在全球70多个国家得到广泛应用，在岩土工程学术界 和工业界赢得了广泛的赞誉。前国际岩石力学会主席 C.Fairhurst(1994)对 FLAC程序的评价是：“现在它是国际上广泛应用的可靠程序。” 我从研二(2010年)开始接触FLAC3D，最初的原因是导师要求每一个人至 少学会一个数值计算软件，而他嘴里每天念叨最多的就是FLAC，自己当时对数 值计算一无所知，便答应老师要学会FLAC3D。第一次打开软件界面，我心里 就凉了大半截，面对着一个操作界面跟记事本无异的所谓“功能强大”的岩土工 程专业软件，半点兴趣也提不起来。年底，从项目工地回到学校准备论文开题， 老师对我的开题报告非常不满意，当着全教研室师生的面，劈头盖脸大批一顿， 第二天又找谈话。在巨大的压力和强烈的自尊心驱使下，我硬着头皮开始啃 FLAC3D，一个半月之后，终于有了初步的计算结果，对老师有个交代，我也能 回家过年了。 前面这一段过程可能是大多数FLAC3D初学者的必经阶段，或者是即将 开始软件学习的人惧怕的事情。毫无疑问，FLAC3D极其不友好的界面是阻碍 初学者前进的很大障碍，当然还包括它是一个全英文的软件。但是当你费尽周折 的走进FLAC3D的世界，你就会发现它独特的魅力，比如简洁的界面，快捷的 命令流操作，高效的计算方法，不易报错等等。另外一个拿不上台面的优点就是 它非常小巧，包括Manual在内一共才几十兆大小，而且已经被破解成绿色版， 只要把它和命令流装进U盘，你就可以随便找一个身边功能最强大的电脑开始 计算了，如果你有过ANSYS、ABAQUS等大型软件痛苦的安装经历，你便能 毕业之后，本以为不用再接触数值计算，但工作需要使得我又一次开始与理解“绿色版”的含义，当然还请大家尊重知识产权，支持正版。 FLAC3D进行亲密接触。我的领导给了我很多新思路和很大的支持，如今我的 水平比研究生时有了不少提高。于是，我想把我的经验总结成文，希望对初学 者起
FLAC3D命令流(挺实用)
1、怎样查看模型？ 答：plot grid 可以查看网格，plot grid num 可以查看节点号。 2、请问在圆柱体四周如何施加约束条件？ 答：可以用fix ... ran cylinder end1 end2 radius r1 cylinder end1 end2 radius r2 not，其中r2 《FLAC 原理实例与应用指南》FLAC3D 实例分析教程 1 FLAC3D 实例分析教程 刘波韩彦辉（美国）编著《FLAC 原理实例与应用指南》 北京：人民交通出版社，2005.9 Appendix(附录) 版权所有：人民交通出版社，Itasca Consulting Group, USA 说明： 1．本实例分析教程是为方便读者学习、应用FLAC 和 FLAC3D 而编写的，作为《FLAC 原理、实例与应用指南》一书的 附录。 2．计算算例参考了Itasca Consulting Group 的培训算例，命 令流的解析旨在方便读者理解FLAC 和FLAC3D 建模及求解问题的一般原则与步骤。 3．实例分析的算例中，FLAC 算例是基于FLAC 5.00 版本、FLAC 3D 算例是基于FLAC 3D 3.00 版本实现计算分析的。读者在学习和研究相关算例时，请务必采用Itasca 授权的合法版本进行分析计算。 4．本实例分析教程仅供读者参考，读者在参考本教程算例 进行工程分析时，编者对可能产生的任何问题概不负责。 编者 《FLAC 原理实例与应用指南》FLAC3D 实例分析教程 2 实例分析1：基坑开挖 图1：基坑开挖的位移等值线图 例1： ; simple tutorial – trench excavation ;简单的隧道开挖方法 ; Step 1: create initial model state ;第一步：初始模型的建立 ; create grid ; 建立网格 gen zone brick size 6 8 8 ;建立矩形的网格区域，其大小为6×8×8 ; assign material ;指定材料的性质 model mohr ;采用摩尔-库仑模型 prop bulk 1e8 shear 0.3e8 ;模型的材料性质：体积模量1×108 剪切模量0.3×108 prop fric 35 coh 1e3 tens 1e3 ;摩擦角35°粘聚力1×103 抗拉强度1×103 ; set global conditions ;设置全局参数 set grav 0,0,-9.81 《FLAC 原理实例与应用指南》FLAC3D 实例分析教程 3 ;设置重力加速度，z 坐标正方向为正，故为-9.81 ini dens 2000 ;初始密度为2000 ; set boundary conditions ;设置初始边界情况 fix x range x -0.1 0.1 ;在x 方向上固定x = 0 边界，为滚动支撑 fix x range x 5.9 6.1 ;在x 方向上固定x = 6 边界，为滚动支撑 fix y range y -0.1 0.1 ;在y 方向上固定y = 0 边界，为滚动支撑 fix y range y 7.9 8.1 ;在y 方向上固定y = 8 边界，为滚动支撑 fix z range z -0.1 0.1 FLAC3D模拟实例循环开挖与支护n res ini.sav set geometry=0.001 ini ydis0 ini xdis0 ini zdis0 ini yvel0 ini xvel0 ini zvel0 m m prop bulk 4.0e9shear 2.5e9fri32coh 2.0e6& range grou diban-shayan ; prop bulk 1.8e9shear 1.2e9fri25coh 1.0e6& range grou diban-niyan any grou hangdao any ; prop bulk 1.2e9shear0.8e9fri22coh0.8e6& range grou diban-gentuyan ; prop bulk 1.9e9shear 1.3e9fri24coh 1.0e6& range grou diban-tniyan ; prop bulk0.7e9shear0.8e9fri21coh0.7e6& range grou mc any grou gzm any ; prop bulk 3.0e9shear 2.5e9fri30coh 1.8e6& range grou dingban-fenshayan ; prop bulk 1.5e9shear 1.2e9fri25coh 1.1e6& range grou dingban-niyan ; prop bulk 3.5e9shear 2.5e9fri34coh 1.4e6& range grou dingban-shayan ;添加接触面 gen separate gzm interface1wrap mc gzm interface2wrap dingban-fenshayan gzm interface1prop kn20e9ks10e9tens1e9 interface2prop kn20e9ks10e9tens1e9 set mech ratio=5e-4 def excav_mc loop n(excav_p,excav_p_z+cut_liang) 有接触面模型 gen zone brick p0 0 0 0 p1 5 0 0 p2 0 5 0 p3 0 0 2 size 10 10 4 group 1 gen zone brick p0 0 0 2 p1 5 0 2 p2 0 5 2 p3 0 0 4 size 10 10 4 group 2 gen zone brick p0 0 0 4 p1 5 0 4 p2 0 5 4 p3 0 0 6 size 10 10 4 group 3 gen separate 2 int 1 wrap 1 2 gen separate 3 int 2 wrap 2 3 plot block group int red fix z range z -0.1 0.1 fix x range x -0.1 0.1 fix x range x 4.9 5.1 fix y range y -0.1 0.1 fix y range y 4.9 5.1 model mohr prop bulk 2.4e9 shear 1.5e9 coh 2.5e6 fri 32 tension 1.2e6 ran group 1 prop bulk 2.4e9 shear 1.5e9 coh 2.5e6 fri 32 tension 1.2e6 ran group 2 prop bulk 2.4e9 shear 1.5e9 coh 2.5e6 fri 32 tension 1.2e6 ran group 3 ini dens 2660 ran group 1 ini dens 2660 ran group 2 ini dens 2660 ran group 3 inter 1 prop kn 0.88e11 ks 0.88e11 fri 15 coh 2e6 ten 1e6 inter 2 prop kn 0.88e11 ks 0.88e11 fri 15 coh 2e6 ten 1e6 set grav 0 0 -10 ini szz -7e6 range z 0 6 ini sxx -8.4e6 range z 0 6 ini syy -8.4e6 range z 0 6 app nstress -7e6 range z 5.9 6.1 x 0 5 y 0 5 plot con zd plot add int black solve 无接触面对比模型 gen zone brick p0 0 0 0 p1 5 0 0 p2 0 5 0 p3 0 0 2 size 10 10 4 group 1 gen zone brick p0 0 0 2 p1 5 0 2 p2 0 5 2 p3 0 0 4 size 10 10 4 group 2 gen zone brick p0 0 0 4 p1 5 0 4 p2 0 5 4 p3 0 0 6 size 10 10 4 group 3 plot block group int red fix z range z -0.1 0.1 fix x range x -0.1 0.1 fix x range x 4.9 5.1 fix y range y -0.1 0.1 fix y range y 4.9 5.1 model mohr n gen zone brick p0 0 0 0 p1 2500 0 0 p2 0 2000 0 p3 0 0 56 size 100 80 2 group 煤层底板; gen zone brick p0 0 0 56 p1 2500 0 56 p2 0 2000 56 p3 0 0 60 size 100 80 1 group 煤层; gen zone brick p0 0 0 60 p1 2500 0 60 p2 0 2000 60 p3 0 0 76 size 100 80 1 group 直接顶;顶板gen zone brick p0 0 0 76 p1 2500 0 76 p2 0 2000 76 p3 0 0 100 size 100 80 1 group 老顶;顶板gen zone brick p0 0 0 100 p1 2500 0 100 p2 0 2000 100 p3 0 0 280 size 50 40 4 group 中粒、粉砂岩互层;顶板 gen zone brick p0 0 0 280 p1 2500 0 280 p2 0 2000 280 p3 0 0 480 size 50 40 4 group 泥岩、砂泥岩互层;顶板 gen zone brick p0 0 0 480 p1 2500 0 480 p2 0 2000 480 p3 0 0 640 size 50 40 4 group 表土层;顶板 modelmohr prop bulk 6.94e9 shear 5.21e9 coh 24.6e7 fric 38 tens 3e7 range group 煤层底板 prop bulk 6.25e9 shear 0.35e9 coh 1e7 fric 25 tens 1e7 range group 煤层 prop bulk 1.33e9 shear 0.44e9 coh 1.5e7 fric 32 tens 1e7 range group 直接顶 prop bulk 7.44e9 shear 5.12e9 coh 8.56e7 fric 37 tens 3e7 range group 老顶 prop bulk 8.26e9 shear 4.72e9 coh 6.25e7 fric 33 tens 2.6e7 range group 中粒、粉砂岩互层prop bulk 8.83e9 shear 1.89e9 coh 1.5e7 fric 33 tens 1e7 range group 泥岩、砂泥岩互层prop bulk 0.017e9 shear 0.0036e9 coh 0.1e7 fric 20 tens 1e7 range group 表土层 initial density 2.6e3 range group 煤层底板 initial density 1.4e3 range group 煤层 initial density 2.4e3 range group 直接顶 initial density 2.5e3 range group 老顶 initial density 2.64e3 range group 中粒、粉砂岩互层 initial density 2.6e3 range group 泥岩、砂泥岩互层 initial density 1.8e3 range group 表土层 fix z range z -0.1 0.1 fix x range x -0.1 0.1 fix x range x 2499.9 2500.1 fix y range y -0.1 0.1 fix y range y 1999.9 2000.1 setgrav 0 0 -9.8 histunbal plotcontzdis 第5部分命令流按照顺序进行 new gen zon radcyl p0 9 0 0 p1 18 0 0 p2 9 36 0 p3 9 0 9 & size 4 24 8 4 dim 3 3 3 3 rat 1 1 1 1.2 group outr gen zone cshell p0 9 0 0 p1 12 0 0 p2 9 36 0 p3 9 0 3 & size 1 24 8 4 dim 2.7 2.7 2.7 2.7 rat 1 1 1 1 group liny fill group inr gen zon reflect dip 90 dd 90 orig 9 0 0 gen zon reflect dip 0 dd 0 ori 9 0 0 ;1/2隧道地层部分 gen zon brick p0 18 0 -9 p1 36 0 -9 p2 18 36 -9 p3 18 0 9 size 10 12 8 group outr1 gen zon brick p0 18 0 9 p1 36 0 9 p2 18 36 9 p3 18 0 13 size 10 12 2 group outr2 gen zon brick p0 0 0 9 p1 18 0 9 p2 0 36 9 p3 0 0 13 size 8 12 2 group outr2 gen zon brick p0 18 0 -23 p1 36 0 -23 p2 18 36 -23 p3 18 0 -9 size 10 12 8 group outr3 gen zon brick p0 0 0 -23 p1 18 0 -23 p2 0 36 -23 p3 0 0 -9 size 8 12 8 group outr3 gen zon reflect dip 90 dd 270 orig 0 0 0 group linz range group liny x -12.1 -5.99 ;绘制模型图 plot add axes red plot set rotation 20 0 45 plot block group ;0-自重应力场求解 set gravity 0 0 -10 ;给定边界条件 fix z range z -23.01,-22.99 fix x range x -36.01,-35.99 fix x range x 35.99,36.01 fix y range y -0.01 0.01 fix y range y 35.99,36.01 model mohr ;围岩的密度 ini density 1800 ;体积、剪切、摩擦角、凝聚力、抗拉强度 prop bulk=1.47e8 shear=5.6e7 fric=20 coh=5.0e4 tension=1.0e4 step 3000 save Gsol.sav ;第1步 rest Gsol.sav initial xdisp=0 ydisp=0 zdisp=0 model null range group inrz y 0 6 作一个盾构隧道开挖，地下水位设为z=12.6m，（隧道中心线设为z=0m）。分了7层土，给定初始条件的代码如下： config fluid （建模过程省） model fl_isotropic prop perm 1e-12 poro 0.48 range z 9.1 12.6 prop perm 1e-12 poro 0.55 range z -0.9 9.1 prop perm 1e-12 poro 0.5 range z -7.4 -0.9 prop perm 1e-12 poro 0.5 range z -11.4 -7.4 prop perm 1e-12 poro 0.45 range z -18.4 -11.4 prop perm 1e-12 poro 0.45 range z -30.4 -18.4 ini fdens 1e3 ini fmod 8.5e7 ini sat 0 range z 12.6 15.1 ini sat 1 range z -30.4 12.6 ini pp 0 grad 0 0 -1e4 range z -30.4 12.6 fix pp range x -.1 .1 fix pp range x 39.9 40.1 fix pp range y -.1 .1 fix pp range y 119.9 120.1 fix pp range z 12.5 12.7 fix pp range z -30.5 -30.3 ;material mechanic parameters model mohr def derive s_mod1=E_mod1/(2.0*(1.0+p_ratio1)) b_mod1=E_mod1/(3.0*(1.0-2.0*p_ratio1)) s_mod2=E_mod2/(2.0*(1.0+p_ratio2))FLAC原理实例与应用指南FLAC3D实例分析教程
FLAC3D模拟实例循环开挖与支护
FLac3D接触面简单实例
FLAC3D条带开采案例(代码)
FLAC3D 实例命令流5
FLAC3D隧道耦合实例