FLAC3D,基本知识

FLAC 3D基础知识介绍

一、概述

FLAC（Fast Lagrangian Analysis of Continua）由美国Itasca公司开发的。目前，FLAC有二维和三维计算程序两个版本，二维计算程序V3.0以前的为DOS版本，V2.5版本仅仅能够使用计算机的基本内存64K），所以，程序求解的最大结点数仅限于2000个以内。1995年，FLAC2D已升级为V3.3的版本，其程序能够使用护展内存。因此，大大发护展了计算规模。FLAC3D是一个三维有限差分程序，目前已发展到V3.0版本。

FLAC3D的输入和一般的数值分析程序不同，它可以用交互的方式，从键盘输入各种命令，也可以写成命令（集）文件，类似于批处理，由文件来驱动。因此，采用FLAC程序进行计算，必须了解各种命令关键词的功能，然后，按照计算顺序，将命令按先后，依次排列，形成可以完成一定计算任务的命令文件。

FLAC3D是二维的有限差分程序FLAC2D的护展，能够进行土质、岩石和其它材料的三维结构受力特性模拟和塑性流动分析。调整三维网格中的多面体单元来拟合实际的结构。单元材料可采用线性或非线性本构模型，在外力作用下，当材料发生屈服流动后，网格能够相应发生变形和移动（大变形模式）。FLAC3D 采用的显式拉格朗日算法和混合-离散分区技术，能够非常准确的模拟材料的塑性破坏和流动。由于无须形成刚度矩阵，因此，基于较小内存空间就能够求解大范围的三维问题。

三维快速拉格朗日法是一种基于三维显式有限差分法的数值分析方法，它可以模拟岩土或其他材料的三维力学行为。三维快速拉格朗日分析将计算区域划分为若干四面体单元，每个单元在给定的边界条件下遵循指定的线性或非线性本构关系，如果单元应力使得材料屈服或产生塑性流动，则单元网格可以随着材料的变形而变形，这就是所谓的拉格朗日算法，这种算法非常适合于模拟大变形问题。三维快速拉格朗日分析采用了显式有限差分格式来求解场的控制微分方程，并应用了混合单元离散模型，可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。

FLAC-3D(Three Dimensional Fast Lagrangian Analysis of Continua)是美国Itasca Consulting Goup lnc 开发的三维快速拉格朗日分析程序，该程序能较好地模拟地质材料在达到强度极限或屈服极限时发生的破坏或塑性流动的力学行为，特别适用于分析渐进破坏和失稳以及模拟大变形。它包含10种弹塑性材料本构模型，有静力、动力、蠕变、渗流、温度五种计算模式，各种模式间可以互相藕合，可以模拟多种结构形式，如岩体、土体或其他材料实体，梁、锚元、桩、壳以及人工结构如支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩、界面单元等，可以模拟复杂的岩土工程或力学问题。

FLAC3D采用ANSI C++语言编写的。

二、FLAC3D的优点与不足

FLAC3D有以下几个优点：

1 对模拟塑性破坏和塑性流动采用的是“混合离散法“。这种方法比有限元法中通常采用的“离散集成法“更为准确、合理。

2 即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍。

3 采用了一个“显式解“方案。因此，显式解方案对非线性的应力-应变关系的求解所花费的时间，几互与线性本构关系相同，而隐式求解方案将会花费较长的时间求解非线性问题。面且，它没有必要存储刚度矩阵，这就意味着，采用中等容量的内存可以求解多单元结构；模拟大变形问题几互并不比小变形问题多消耗更多的计算时间，因为没有任何刚度矩阵要被修改。

当然，它也存在以下几个不足之处：

1 对于线性问题的求解，FLAC3D比其他有限元程序运行得要慢；但是，当进行大变形非线性问题

或模拟实际可能出现不稳定问题时，FLAC3D是最有效的工具。

2 用FLAC3D求解时间取决于最长的自然周期和最短的自然周期之比。

三、FLAC3D的特点

1、应用范围广泛

1.1 包含10材料本构模型

Flac3D中为岩土工程问题的求解开发了特有的本构模型，总共包含了10种材料模型：

1. 开挖模型null

2. 3个弹性模型（各向同性，横观各向同性和正交各向同性弹性模型）

3. 6个塑性模型（Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型和修正的cam粘土模型）。

Flac3D网格中的每个区域可以给以不同的材料模型，并且还允许指定材料参数的统计分布和变化梯度。还包含了节理单元，也称为界面单元，能够模拟两种或多种材料界面不同材料性质的间断特性。节理允许发生滑动或分离，因此可以用来模拟岩体中的断层、节理或摩擦边界。

FLAC3D中的网格生成器gen，通过匹配、连接由网格生成器生成局部网格，能够方便地生成所需要的三维结构网格。还可以自动产生交岔结构网格（比如说相交的巷道），三维网格由整体坐标系x,y,z系统所确定，这就提供了比较灵活的产生和定义三维空间参数。

1.2 有五种计算模式

(l)静力模式。这是FLAC-3D默认模式，通过动态松弛方法得静态解。

(2)动力模式。用户可以直接输人加速度、速度或应力波作为系统的边界条件或初始条件，边界可以固定边界和自由边界。动力计算可以与渗流问题相藕合。

(3)蠕变模式。有五种蠕变本构模型可供选择以模拟材料的应力-应变-时间关系:Maxwell模型、双指数模型、参考蠕变模型、粘塑性模型、脆盐模型。

(4)渗流模式。可以模拟地下水流、孔隙压力耗散以及可变形孔隙介质与其间的粘性流体的耦合。渗流服从各向同性达西定律，流体和孔隙介质均被看作可变形体。考虑非稳定流，将稳定流看作是非稳定流的特例。边界条件可以是固定孔隙压力或恒定流，可以模拟水源或深井。渗流计算可以与静力、动力或温度计算耦合，也可以单独计算。

(5)温度模式。可以模拟材料中的瞬态热传导以及温度应力。温度计算可以与静力、动力或渗流计算藕合，也可单独计算。

1.3 可以模拟多种结构形式

(l)对于通常的岩体、土体或其他材料实体，用八节点六面体单元模拟。

(2)FIAC-3D包含有四种结构单元:梁单元、锚单元、桩单元、壳单元。可用来模拟岩土工程中的人工结构如支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩等。

(3)FLAC-3D的网格中可以有界面，这种界面将计算网格分割为若干部分，界面两边的网格可以分离，也可以发生滑动，因此，界面可以模拟节理、断层或虚拟的物理边界。

1.4 可以有多种边界条件

边界方位可以任意变化，边界条件可以是速度边界、应力边界，单元内部可以给定初始应力，节点可以给定初始位移、速度等，还可以给定地下水位以计算有效应力、所有给定量都可以具有空间梯度分布。

2 FLAC-3D内嵌语言FISH

FLAC-3D具有强大内嵌语言FISH，使得用户可以定义新的变量或函数，以适应用户的特殊需要，例如，利用HSH做以下事情:

(l)用户可以自定义材料的空间分布规律，如非线性分布等。

(2)用户可以定义变量，追踪其变化规律并绘图表示或打印输出。

(3)用户可以自己设计FLAC-3D内部没有的单元形态。

(4)在数值试验中可以进行伺服控制。

(5)用户可以指定特殊的边界条件。

(6)自动进行参数分析。

(7)利用FLAC-3D内部定义的Fish变量或函数，用户可以获得计算过程中节点、单元参数，如坐标、位移、速度、材料参数、应力、应变、不平衡力等。

3 FLAC-3D具有强大的前后处理功能

FLAC-3D具有强大的自动三维网格生成器，内部定义了多种单元形态，用户还可以利用FISH自定义单元形态，通过组合基本单元，可以生成非常复杂的三维网格，比如交叉隧洞等。

在计算过程中的任何时刻用户都可以用高分辨率的彩色或灰度图或数据文件输出结果，以对结果进行实时分析，图形可以表示网格、结构以及有关变量的等值线图、矢量图、曲线图等，可以给出计算域的任意截面上的变量图或等直线图，计算域可以旋转以从不同的角度观测计算结果。

四、FLAC3D做计算分析的一般步骤：

与大多数程序采用数据输入方式不同，FLAC采用的是命令驱动方式。命令字控制着程序的运行。在必要时，尤其是绘图，还可以启动FLAc用户交互式图形界面。为了建立FLAC计算模型，必须进行以下三个方面的工作：

1. 有限差分网格

2. 本构特性与材料性质

3. 边界条件与初始条件

完成上述工作后，可以获得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后，进行工程开挖或改变边界条件来进行工程的响应分析，类似于FLAC的显式有限差分程序的问题求解。与传统的隐式求解程序不同，FLAC采用一种显式的时间步来求解代数方程。进行一系列计算步后达到问题的解。

在FLAC中，达到问题所需的计算步能够通过程序或用户加以控制，但是，用户必须确定计算步是否已经达到问题的最终的解

五、FLAC3D分析的使用领域

根据手册中所说，总结如下：

1 承受荷载能力与变形分析：用于边坡稳定和基础设计

2 渐进破坏与坍塌反演：用于硬岩采矿和隧道设计

3 断层构造的影响研究：用于采矿设计

4 施加于地质体锚索支护所提供的支护力研究：岩锚和土钉的设计

5 排水和不排水加载条件下全饱和流体流动和孔隙压力扩散研究：挡土墙结构的地下水流动和土体固结

研究

6 粘性材料的蠕变特性：用于碳酸钾盐矿设计

7 陡滑面地质结构的动态加载：用于地震工程和矿山岩爆研究

8 爆炸荷载和振动的动态响应：用于隧道开挖和采矿活动

9 结构的地震感应：用于土坝设计

10 由于温度诱发荷载所导致的变形和结构的不稳定

12 大变形材料分析：用于研究粮仓谷物流动和放矿的矿石流动

六、后处理

用tecplot绘制曲线

1.第一主应力

2.xdisp、ydisp、zdisp、disp

用excel做曲线

隧道

1做地表沉降槽（zdisp）

2地表横向位移（xdisp）

3隧道中线竖向沉降曲线（zdisp）

4提取位移矢量图，

5显示初期支护结构内力

6显示state（找塑性区）

基坑

1做地表沉降槽（zdisp）

2提取位移矢量图，

3显示初期支护结构内力

4显示state（找塑性区）

边坡

做安全系数和应变图

七、模型最优化

用FLAC3D解决问题时，为了得到最有效的分析使模型最优化是很重要的。这个章节对改进模型的运行提供了一些方法建议。同时，准备计算时需要避免的一些通常出现的缺陷也列了出来。

1.检查模型运行时间

一个FLAC3D例子的运行时间是区域数的4/3倍。这个规则适用于平衡条件下的弹性问题。对于塑性问题，运行时间会有点改变，但是不会很大，但是如果发生塑性流动，这个时间将会大的多。对一个具体模型检查自己机子的计算速度很重要。一个简单的方法就是运行5.1节所给的基准测试。然后基于区域数的改变，用这个速度评估具体模型的计算速度。

2.影响运行时间的因素

FLAC3D有时会需要较长时间才可以收敛主要发生在下列情况下：

(a)材料本身刚度变异或材料与结构及接触面之间的刚度差异很大。

(b)划分的区域尺寸相差很大。

这些尺寸差异越大编码就越无效。在做详细分析前应该研究刚度差异的影响。例如，一个荷载作用下的刚性板，可以用一系列顶点固定的网格代替，并施以等速度。（记住FIX命令确定速度，而不是位移。）地下水的出现将使体积模量发生明显的增加——见理论卷第三章流体-固体相互作用分析。

3.考虑网格划分的密度

FLAC3D使用常应变单元。如果应力/应变曲线倾斜度比较高，那么你将需要许多区域来代表多变的分区。通过运行划分密度不同的同一个问题来检查影响。FLAC3D应用常应变区域，因为当用多的少节点单元与用比较少的多节点单元模拟塑性流动时相比更准确。

应尽可能保持网格，尤其是重要区域网格的统一。避免长细比大于5:1的细长单元，并避免单元尺寸跳跃式变化（即应使用平滑的网格）。应用GENERATE命令中的比率关键词，使细划分区域平滑过渡到粗划分区域。

4.自动发现平衡状态

默认情况下，当执行SOLVE 命令时，系统将自动发现力的平衡。当模型中所有网格顶点中所有力的平均量级与其中最大的不平衡力的量级的比率小于1*10时，认为达到了平衡状态。注意一个网格顶点的力由内力（例如，由于重力）和外力（例如，由于所加的应力边界条件）共同引起。因为比率是没有尺寸的，

所以对于有不同的单元体系的模型，在大多数情况下，不平衡力和所加力比率的限制给静力平衡提供了一个精确的限制。

同时还提供了其他的比率限制；可以用SET ratio 命令施加。如果默认的比率限制不能为静力平衡提供一个足够精确的限制，那么应考虑可供选择的比率限制。

默认的比率限制同样可用于热分析和流体分析的稳定状态求解。对于热分析，是对不平衡热流量和所加的热流量量级进行评估，而不是力。对于流体分析，对不平衡流度和所加流度量级进行评估。

5.考虑选择阻尼

对于静力分析，默认的阻尼是局部阻尼，对于消除大多数网格顶点的速度分量周期性为零时的动能很有效。这是因为质量的调节过程依赖于速度的改变。局部阻尼对于求解静力平衡是一个非常有效的计算法则且不会引入错误的阻尼力（见Cundall 1987）。

如果在求解最后状态，重要区域的网格海域的速度分量不为零，那么说明默认的阻尼对于达到平衡状态是不够的。有另外一种形式的阻尼，叫组合阻尼，相比局部阻尼可以使稳定状态达到更好的收敛，这时网格将发生明显的刚性移动。例如，求解轴向荷载作用下桩的承载力或模拟蠕变时都可能发生。使用SETmechanical damp combined命令来调用组合阻尼。组合阻尼对于减小动能方面不如局部阻尼有效，所以应注意使系统的动力激发最小化（见例3.14）。可以用SET mechanical damp local命令转换到默认阻尼。

6.检查模型反应

FLAC3D 显示了一个相试的物理系统是怎样变化的。做一个简单的试验证明你在做你认为你在做的事情。例如，如果荷载和实体在几何尺寸上都是对称的，当然反应也是对称的。改变了模型以后，执行几个时步（假如，5或10步），证明初始反应是正确的，并且发生的位置是正确的。对应力或位移的期望值做一个估计，与FLAC3D 的输出结果作比较。

如果你对模型施加了一个猛烈的冲击，你将会得到猛烈的反应。如果你对模型作了一些看起来不合理的事情，你一定要等待奇怪的结果。如果在分析的一个给定阶段，得到了意外值，那么回顾到这个阶段所用的时步。

在进行模拟前很关键的是检查输出结果。例如，除了一个角点速度很大外，一切都很合理，那么在你理解原因前不要继续下去。这种情况小，你可能没有给定适当的网格边界。

7.初始化变量

在模拟基坑开挖过程时，在达到目的前通常要初始化网格顶点位移。因为计算次序法则不要求位移，所以可以初始化位移，这只是由网格顶点的速度决定，并有益于用户初始化速度却是一件难事。如果设定网格顶点的速度为一常数，那么这些点在设置否则前保持不变。所以，不要不要为了清除这些网格的速度而简单的初始化它们为零——这将影响模拟结果。然而，有时设定速度为零是有用的（例如，消除所有的动能）。

8.最小化静力分析的瞬时效应

对于连续性静力分析，经过许多阶段逐步接近结果是很重要的——即，当问题条件突然改变时，通过最小化瞬时波的影响，使结果更加“静力”。使FLAC3D解决办法更加静态的方法有两种。

（1）.当突然发生一个变化时（例如，通过使区域值为零模拟开挖），设定强度性能为很高的值以得到静力平衡。然后为了确保不平衡力很低，设定性能为真实值，再计算，这样，由瞬时现象引起的失败就不会发生了。

（2） .当移动材料时，用FISH 函数或表格记录来逐步减少荷载（见1.2节中的例子）。

9.改变模型材料

FLAC3D 对一个模拟中所用的材料数没有限制。这个准则已经尺寸化，允许用户在自己所用版本的FLAC3D中最大尺寸网格的每个区域（假如设定的）使用不同的材料。

10.运行在现场原位应力和重力作用下的问题

有很多问题在建模时需要考虑现场原位应力和重力的作用。这种问题的一个例子是深层矿业开挖_——

回填，此时大多数岩石受很高的原位应力区的影响（即，自重应力由于网孔尺寸的限制可以忽略不计），但是回填桩的放置使自重应力发展导致岩石在荷载作用下可能坍塌。在这些模拟中要注意的重点（因为任何一种模拟都有重力的作用）是网格的至少三个点在空间上应固定——否则，整个网格在重力作用下将转动。如果你曾经注意到整个网格在重力加速度矢量方向发生转动，那么你可能忘记在空间上固定网格了（见例

3.16）。

八、算例

FLAC3D3.0在某隧道工程开挖支护中的应用-命令流.dat

new

set log on

set logfile wang.log

gen zon radcyl p0 0 0 0 p1 9.0 0 0 p2 0 50 0 p3 0 0 8 &

size 4 20 6 4 dim 6 5 6 5 rat 1 1 1 1 group 围岩

gen zon cshell p0 0 0 0 p1 6.0 0 0 p2 0 50 0 p3 0 0 5.0 &

size 4 20 6 4 dim 5.6 4.6 5.6 4.6 rat 1 1 1 1 group 初期支护

gen zon cshell p0 0 0 0 p1 5.6 0 0 p2 0 50 0 p3 0 0 4.6 &

size 4 20 6 4 dim 5.0 4.0 5.0 4.0 rat 1 1 1 1 group 二次衬砌 fill group 原岩

gen zon radcyl p0 0 0 0 p1 0 0 -8.0 p2 0 50 0 p3 9.0 0 0 &

size 4 20 6 4 dim 3 6 3 6 rat 1 1 1 1 group 围岩2

gen zon cshell p0 0 0 0 p1 0 0 -3.0 p2 0 50 0 p3 6.0 0 0 &

size 4 20 6 4 dim 2.6 5.6 2.6 5.6 rat 1 1 1 1 group 仰拱初期支护

gen zon cshell p0 0 0 0 p1 0 0 -2.6 p2 0 50 0 p3 5.6 0 0 &

size 4 20 6 4 dim 2 5 2 5 rat 1 1 1 1 group 仰拱二次衬砌 fill group 仰拱原岩

gen zone reflect normal -1 0 0

gen zone radtun p0 0 0 0 p1 45 0 0 p2 0 50 0 p3 0 0 20 &

size 3 20 3 12 dim 9 8 9 8 rat 1 1 1 1.1 group 围岩3

gen zon reflect dip 0 ori 0 0 0 range x 0 9 y 0 50 z 8 20

gen zon reflect dip 0 ori 0 0 0 range x 9 45 y 0 50 z 0 20

gen zon reflect dip 90 dd 270 ori 0 0 0 range x 0 9 y 0 50 z 8 20

gen zon reflect dip 90 dd 270 ori 0 0 0 range x 0 9 y 0 50 z -8 -20

gen zon reflect dip 90 dd 270 ori 0 0 0 range x 9 45 y 0 50 z -20 20

gen zon brick p0 -45 0 -20 p1 -45 0 -40 p2 -45 50 -20 p3 45 0 -20 &

size 5 20 6 rat 1.1 1 1 group 围岩4

save tun_model.sav

;假设围岩岩体符合mohr-coulomb本构模型，给围岩赋参数命令流如下，

; mohr-coulomb model

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))

b_mod2=E_mod2/(3.0*(1.0-2.0*p_ratio2))

end

set E_mod1=0.6e9 p_ratio1=0.27 E_mod2=0.8e9 p_ratio2=0.26

derive

prop bulk b_mod1 shear s_mod1 cohe 1.8e6 tens 0.8e6 fric 30 range z 4.5 20

prop bulk b_mod2 shear s_mod2 cohe 2.8e6 tens 1.0e6 fric 35 range z -40 4.5

ini dens=2300

set grav 0 0 -10

; boundary and initial conditions

apply szz -1.4e6 range z 19.9 20.1

fix z range z -40.1 -39.1

fix x range x -45.1 -44.9

fix x range x 44.9 45.1

fix y range y 49.9 50.1

hist unbal

hist gp xdis 6.0,0,0

hist gp zdis 0,0,5

hist gp xdis 6.0,50,0

hist gp zdis 0,50,5

plot hist 3

solve

save tun_nature.sav

;对后面计算而言，模型建立时岩体在开挖前认为位移已经终了，因此需要对位移进行“清零”，而应力可以保留。假设隧道先开挖上断面，中间不设支护，直到进尺50m，那么此时位移和应力的分布情况可用如下命令流

ini xdis=0 ydis=0 zdis=0

plot re

model null range group 原岩

model null range group 二次衬砌

model null range group 初期支护

set large

hist unbal

plot hist 3

solve

save tun_ext1.sav

;如果在开挖后适时对隧道进行锚喷混凝土初期支护，该命令流如下：

restore tun_nature.sav

ini xdis=0 ydis=0 zdis=0

plot re

model null range group 初期支护

model null range group 二次衬砌

model null range group 原岩

step 100

hist unbal

plot hist 3

mo el range group 初期支护

prop bulk 1.33e9 shear 0.8e9 cohe 2.2e6 tens 1.2e6 fric 32 range z 4.5 7.0 x -7.0 7.0

prop bulk 10.9e9 shear 8.9e9 dens 2500 range group 初期支护set large

solve

save tun_ext2.sav