有限元分析方法之FLAC程序说明1

FLAC程序说明书

2004-8-26

1

FLAC(version 2.00)

Fast Lagrangian Analysis of Continua

目录

1.0引论 5.6建议及忠告37 1.1 FLAC的技术要求及装6.0 FLAC中的结构模拟39

机„„„„„„3 6.1命令结构39

1.2 绘图机故障分析4 6.2定义结构单元的几何条件及其支2.0立即满意－应用FLAC的一个简单承介质的联动装置40

指导性示例4 6.3实例应用40 2.1建造在非线性之中的壕沟4 7.0例题45 3.0基础知识—显式有限差分法8 7.1例1无摩擦粘土上的毛石基脚45 3.1引论8 7.2例2 粘性摩擦土的边坡稳定47 3.2显式／计算循环8 7.3例3 端部有剪力的弹性悬臂梁3.3有限差分格式9 51 3.3.1导数的表示9 7.4例4弹性，弹塑性及横向各向同性3.3.2运动方程式9 岩石介质中，受初应力作用的圆形隧3.4速度／应变增量方程式9 洞 3.5应力／应变规律10 7.4.1弹性岩石介质51 3.6确定网点处的不平衡力10 7.4.2弹塑性岩石介质53 3.7应力转动修正项11 7.4.3横向各向同性岩石介质55 参考文献12 参考文献56 4.0输入指令 8.0运行FLAC时值得注意的重点及注4.1定义12 意事项57 4.2输入命令12 8.1初始化各变量57 4.3设置你自己的默认条件 8.2改变材料模型57 参考文献 8.3运行含现场应力和重力的问题57 5.0用FLAC解答的问题 附录A本构模型描述59 5.1引论28 A1引论59 5.2运行FLAC28 A2弹性各向同性模型59 5.2.1网格的形成及材料特性的定义A3Mohr－Coulom模型 29 A4空模型 5.2.2应用边界条件33 A5各向异性弹性60 5.2.3应用荷载／变化条件33 A6多处存在结理的模型61 5.2.4数据的打印及绘图33 A7应力软化／强化模型61 5.3特殊问题的考虑36 参考文献63 5.3.1大应变36 附录B利用FLAC时确定平衡条件 5.3.2平面应力36 附录C错误及警告信息66 5.3.3重力36 附录D FLAC中的界面逻辑68 5.3.4图形形状36 FLAC快速查阅命令清单70 5.3.5存入36 FLAC 2.01版补遗72 5.4错误处理37 5.5储存／复原运行37

2

FLAC

Fast Lagrangian Analysis of Continua

(Version 2.00)

@1987

ITASCA Consulting Group, INC.

P.O.Box :14806

Minneapolis, MinneSota 514 ITASCA Consulting Group, INC. 持有执照的FLAC的条款及规则 （Terms and for licensing FLAC）

使用FLAC程序之前，你应当仔细阅读以下各条款及规则．把FLAC插入你的计算机，意味着你已承认这些条款及规则．如果你不这个程序是由Ttasca咨询小组有限公司提供的．给记录的数个体起的名称和给支持所起的名称是使用转递用的，但给程序起的名称．达到你的预想效果，并负责安装，使用及从程序中获得结果．

许可证（License）

在任何一个时间，你只又能在一台计算机使用这个程序． 仅仅为了延伸利用，你可以复制备份程序．

除了本文件提到外，你不能利用，复制，修改或传送本程序或任何复制大部分或一部分．

你不能再执照，出租本程序． 有效期限（Term）

本许可证终止前一直有效．任何时候，你可以通过程序运用备份拷贝来终止它．如果你不遵守（fail to comply with）本的任何条款或条件，它也会终止．你同意这样的终止，使坏程序连同备份拷贝，以任何方式的修改和／或的各部分．

保证书(Warrany)

在本代码，在12个月内将免费改正代码中的任何错误，完整的表列输入及输出文件，并错误的书面形式，给出通知．如果经判断，代码有错误，将视情况免费修改或交换拷贝，或偿还．

责任界限(Limitation of Liability)

概不负责：关于FLAC或任何部分的使用；关于使用而造成的任何损坏或掉失，包括由于使用FLAC而造成的时间，金钱或信誉损失（包括各种修改或改正而造成的）．决不负责因使用FLAC而造成的间接的，错误的，偶然的或随时而发的各种损坏．

1.0引论

FALC是一种显式有限差分代码（explicit finite difference code），它模拟由岩土或其它材料建造的结构物的性能；这些材料达到屈服极限时，可能经历塑流．这些材料是通过构成一个网格的域或者单元来表示的；这个网格由用户调整，以拟合模拟对象的外形。各单元对外力边界约束的反应，是服从某一规定的线性或非线性应力－应变定律．如果应力大得足以使材料发生屈服和流动，实际上，网格将随其所用材料变形和移动．本计算方案被称为＂Lagrangian＂（拉

3

格朗日算法）并完全适宜于模拟大变形．解答的显式特性，指的是任何模型的非线性应力－应变规律都能服从；而不招致显著的过大（overload）（用线性规律相加）

诸如与因岩或土相互作用的隧道衬砌，岩石锚杆或板桩等结构，也可以模拟．因此，就考察支护坑道的稳定性，或考察岩土的不稳定性对地面结构的影响． 1.1 FLAC的技术要求及装机（FLAC Specification and Installation）

FLAC 2.00应具有以下设备： ⑴IBM－PC 或兼容CPU；

⑵0Kb RAM(Random Access Memory 随机贮存器)； ⑶一个360Kb或2Mb软盘； ⑷一个10M或更大的硬盘；

⑸IBM或兼容彩色图像转换器或增强图像转换器； ⑹8007或8029数字处理器； ⑺DOS3.1或更高版本；

⑻（可选择的）Hewlett-Pockerd或兼容的2或6草绘图机及只读器的点阵打印机．

对于上述技术要求，能模拟约2000个域或单元．

FLAC靠一个360K的软盘转递，而必须储存于计算机的硬盘上．提供给你的FLAC点上的ARC程序，用于重新组装你的硬盘上的可执行文件．从软盘将所有文件拷贝到你的磁盘上．然后发出下列命令：

C> ARC x flac.arc flac.exe 这时，计算机将从Flac.arc文件中建立一个：非浓缩的（de-compressed）＂Flac.exe文件．

包含在软盘中的几个例子的文件为：⑴Flac.CGA;⑵Flac.EGA;⑶几个数据文件．文件Flac.CGA及Flac.EGA具有一些Flac命令，这些命令将你的系统置于CGA及EGA图形模式．

如后所述，起动时Flac将在已写入的Flac目录中，寻找称之Flac.INI文件．这个文件包含一系列命令，这些命令设置多次利用本代码正常运行该系统或程序的属性．例如，具有EGA系统的用户，将需要在Flac上建立命令SET.EGA以设置EGA卡的属性．Flac.CGA及Flac.EGA文件具有一些Flac.CGA更名为Flac.INI

C>REName flac.cga flac.ini(设置显示模式)

如果你具有EGA系统，把Flac.CGA更名为flac.INI C>Rename flac.EGA flac.INI

你可以附加任何命令于你想利用任何文本编译器的那些文件．

FLAC能把线图传给Hewlett-packard（或兼容的）式绘图机或带图形ROM的点阵打印机．Flac式绘图机的命令传给用户选择上．多次式绘图机与串联接口COM连接．因此，式绘图机的输在存入Flac之前指向COM．串联接口也必须初始化，以便，缺省值为以（band）把线图传给COMI．

如果你的绘图机与COM2连接，或者如果你想改变波特（baud）速率，直接利用Fｌａｃ命令。ｓｅｔ output=Com 2及ｓｅｔ baud=b 式中b=1200,2400,4800或9600，用以改变确省设置。设置这些值的最简便的方法就是在文件Flac.in.set output及set baud命令。.须记住，绘图机的波特速率，必须与串联接口初始化值一至检查你的绘图机及所有关这些设置的说明．FLAC

4

采用美国的纸张时，应购置相应的绘图机．利用命令

Sｅｔ．ｏｕｔｐｕｔ＝Lｐｔ１

也能将输出传给平行接口．作为代替式绘图机图形，变换图形可以传给点阵式绘图机打印机，必自装入DOS图形程序在装入FLAC前，键入

C>graphics 当使用FLAC时，你可以通过同时按下ｓｈｉｆｔ与ｐｒｉｎｔｓｃｒｅｅｎ键的方式，利用点阵式打印机打出屏幕图高分辨率点阵或EGA屏幕打印可以利用市场上供应的程序包． 1.2绘图机故障分

1.在使用式绘图机中遇到的最常见错误，是不能配置匹配于串联接口与草式绘图机的波特速率，确保草式绘图机上的开关词是正确地调整到的波特速率． 2.应确保绘图机是联接于文件plotter.BAT中所指定的未有串联接口上； 3.如果在草式绘图机的开关上选择了非美国纸张尺寸，刻写在线图上的文字可能变形．应确保选美国纸张尺寸．

2.0立即满意—应用FLAC的一个简单指导性示例

本节是入并开始试验FLAC的用户提供的．选择了一个简单，快速运行的示例题目．它FLAC代码的的状况． 2.1建造在非粘性土中的壕沟

一个既能交互式考察，又能说明FLAC代码某些机能的简单问题，是在否非粘性土介质中开挖的一个壕沟．

键入 C>flac

后，执行flac．计算机将装入程序，并在原有标题下立即显示flac> 为建立最初的有限差分网格，运用GRID命令： flac＞grid 5,5

该命令将创建一个宽5个域（或单元），高5个域（或单元）的网格，现在，应各个域的材料模量及特性，对于本例题，我们采用摩尔－库仑弹塑性模型．键入命令：

flac>model mohr

flac>prop bulk=1e8 Shear=.3e8 fric=35 flac>prop dens=1000 coh=1e10

这里，我们规定Mohr-conlomb模型（如同后面将看到的，不需要键入各命令的全部字母，计算机是能够网格中的每个域可以想像地会具有不同的材料模型和特性．然而，在Model命令后，没有规定域的flac认为所有域均为Mohr-Coulomb．接着给出特性值－包括体积模量bulk Modulus（单位pa，仅写b=或直接写b，后限要求的值，是完全一样的；为了明白起见，这里给出了全名），剪切模量shear modulus；内摩擦角及粘聚力cohesion．你会看出，给出了一个很高的粘聚力值．这又是在体内形成重力应力过程中，采用的初始值．实际上，在最初形成重力应力，我们加载于具有弹性物体．这能避免在模型最初状态期间的任何塑性屈服．只要你懂得明确模拟过程经验，理由会变得很显然．既然已确定了网格及模型特性，便能绘或打印出与模拟有关的数据．发布下列命令：

flac>plot x y

x,y坐标将在网点实际位置以表格形式出现．你会注意到，该表有i行，j列（沿表的左上边从1到6）．每个网点及域与i行，j列的一个对应．在本例中，

5

网点的范围为i=1~6,j=1~6，而域的范围有I=1~5,j=1~5．如果你要更进一步了解，参用第4章的土4－1．为了看清网格图形，发下列命令：

flac>plot grid

依据你的Flac.INI文件而定，你将在屏幕上以低分辩率或EGA模式看到一幅网格图线．看完后，按回车键即可返回flac状态．注意，如果你不规定网格具有坐标（利用GEN或INITIAL命令），那么x,y的规定为坐于网点数减1．例如，在上述网格土中，左下角网点取为层点，并给出坐标（0，0），右下角（网点（6，1））的坐标为（5m，0）．用户可以完全地用GEN及INITIAL命令去规定任意点的坐标．为使本例题简单，我们分网格为5m×5m． 接着，便设定问题的边界条件．在本例中，我们要在底部即两侧边布置轴边界，重力作用于各域并允许它们处于自然状态，产生现场应力．为了使用这些边界（即，在规定的右边无位移或速度），利用下列命令： flac>fix y j=1 flac >fix x I=1 flac>fix x I=6 上述命令具有以下功能： 1.y方向被固定．当flac见到j=1时，便自动取I的范围1～6（即，全范围）．你规定j=1,I=1,6能完成同样功能．

2.左边网点（I=1）及右边网点（I=6）沿x方向均被固定．FLAC再次取j方向的全范围．然后，键入：flac>set grav=9.81

我们设定重力（gravity）式中9.81m/sec2为重力加速度．重力向下为正，向上为负（如果重力认为负的，则物体将会浮起） 如果你想看到模型的某一点的位移，以便观察随时变化的平衡或塌陷，应键入： flac>his hstep=5

flac>his ydis I=2 j=6

这里，我们对表的第一点选取5个时间步长的y向位移跟踪监视．

这时，我们准备好了使初始模型处于平衡．因为flac是一个显式动力代码，我们通过时间跨入模型（允许网格动能逐渐降低的时间（因此，提供了我们寻找的静力解））．允许在物体内形成重力，我们时间跨入模型以------．这里，用SLOVE命令自动检验平衡⑴

flac>Solve

flac花一定时间（几分钟）去＂判断＂．在每个时间生长处，步长数及最大不平衡力将显示在屏幕上，当完成时，Flac将返回一个信息．以示已达到极限．立即再显示flac．

现在，我们能看见在模型内已发生什么．早先要求的y向位移史： flac>plot his 1

一个屏幕图形将以单色图或EGA表观运输．它表示模型约在108个时间步长达到平衡．

我们来考察在物体内形成的重力应力．窗口是自动定义的，如果你想放大或缩小图形，你可以在重新设定窗口时通过键入：

flac>title

title>a simple trench excavation example:the initial gravity stresses

给图形一个标题⑵

6

（如果你用单色模式，脱离彩色键盘），应键入 flac>Set pal=0

（以便设置调色板选择器），然后，键入

flac>plot Syy yei bou gre

你将看见一个边界为绿色（单色屏幕全为绿色），只有黄灰（yellow-brown）色的yy应力图（图２－１）．同样，绘yy应力图应键入：

flac>plot Sxx yel bou gre

你将注意到，重力应力随深度线性增加．这些是能发现的，须键入：

flac>print Sxx Syy

将这个初始状态储存起来是明智的，以便在任何时间可以重新启动来完

成参数研究．为了存入这个，键入：

flac>Save trench.sav

在默认驱动上将建立一个存入文件．然后返回一个flac提示．现在我们

可以在土中开挖一个壕沟．键入：

flac>prop coh=0

由各无粘聚力且带无支护的坚壁，一定会崩塌．因为我们要确切地考察这个过程，在代码中必须设置大应变．为此，键入： flac>Set large

最后，为了绘图，我们只想观察从壕沟开挖起的位移变化，而不是上述的重力调整量――因此，我们可以把x,y的位移调为零⑶：

flac>INIT Xdis=0 Ydis=0 为了开挖壕沟，键入：

flac>model null I=3 j=3,5

因为我们故意取粘聚力低得足以发生崩塌．我们不想用是有一个不平衡

力极限的Solve命令（不平衡力变平衡）⑷，因为我们的模拟决不会收敛于平衡状态．而我们能一个时间步长一次跨过模拟过程，关闭和打印出发生崩塌时的结果．这是显示方法的真正本领－不要求模型在每一个计算后收敛于平衡，因为我们从末解过线性代数方程组，这正如大多数工程师们所熟悉的显示代码一样．我

⑸

们利用step命令：

flac>Step 100

3.0基础知识－显式有限差分法 3.1引论

Flac是一个一般地质力学数值模型，它利用显式有限差分法去解答运动基本方程。一般讲，差分法涉及到被模拟的物体被划分为若干在网点（节点）处相互连接的若干个二维域（单元）。在每个断点处，把被未解的运动方程的格式取为时间步长方式。因此，有可能体系的性能看作是随时间逐渐改变的。

某些著者（Wilkins是x中之一，1963）已经证明，由有限差分法和有限单元法导得的方程式，对于个别问题而言是一致的。使用中，有限差分法如同有限单元法一样，适应性很强，可运用不规则网格，变化材料类型或特性，也可规定不同的边界条件。有限差分法详细用于地质力学问题为Cundall(1976)所讨论。

⑴见第四章Solve命令的设置范围。

⑵标题及图例（Legend）出现在硬拷贝草图上，但不出现在屏幕上。

7

⑶这不会影响计算，因为在计算工序中模型不要求这些位移。保持它们

又是为了方便用户。

⑷更多的细节，见第四章的Solve命令。

⑸Slolve命令也能用于时间步长数量的范围。

3.2显示／计算循环

一般讲，在解答一个特定问题的运动方程时，有两种采用数值码的方法。（求解静力问题所应用的数值方法）隐示法同时求解各网格的未知量值。换言之，要建立一个把未知量与已知量联系起来的方程组，比如，在有限元法某模型中，通过整体刚度矩阵，建立节点力与节点位移间的关系。必须储存并求解这个方程组，从而造成大的计算机内存设备。

显式法利用的概念是：对于一个很小时间步长，某一已知网点处的干扰只能靠它的直接领域来承受。比如，设想在时间t=0时，某个网点的温度升高到了某个值。在一个短时间内，只有各邻近网点到温度已升高了，随时间推进，这个联响将遍布网格，从而在周围各网点处导致较高的温度。

必须合理的选取时间步长，以避免解的数字不稳定性－即时间步长必须小于两相邻网点间传递现象的时间。在flac的力学版本中，用物体内的声速来控制这个时间步长；在代码的超力版本中由热扩散率及对流换算系数来检测时间步长。Flac自动确定保证数字稳定的时间步长。

目前，Flac是通过适当地衰减动力解的办法来求解静力问题。这时“时间步长”不是指的实际时间步长，而是问题的时间步长，其中速度以每一时间步长的长度来度量的。这个解法对于在个人计算机模拟地质材料具有很多优点。

因为从不形成矩阵，代码的内存设备是最小的，而且每个时间步长的计算浪费也是很小的。对于地质材料来说，显式方法的一个很重要方面是能轻而易举地处理非线性本构规律。对于服从非线性应力－应变规律的材料，不要求（可能引起解的明显误差的）迭代法；相应于也给应变改变的应力变化，如同实际发生的一样，可以在已知域处实施。用这种方法时，用适当的物理方式服从非线性规律，而不依赖于迭代法的路径灵敏度。此外，由于不形成矩阵，可以附加少量的计算耗费便能处理诸如大位移这类现象。然而，也有消极的一面，一般来说，显式码比其隐式的静力、弹性问题运行要慢些。

模拟地质材料的含义是什么呢？对于土壤和岩石，这类材料常常遭受破裂或屈服（即，非弹性），且可能出现大位移。因此，一个物体在一个孤立区可能经受崩塌（如，土的滑移），而部分仍稳定的。对于这类问题，在各网点耦合的隐式中，在矩阵求答的过程中，会成为数字上不稳定的。然而，是代码不遭受时间或稳定性的损失。更微妙的优点是，用户能考

运动规律 察屈服或崩塌的形成过程，而不仅能观察最终的

对各个网点系统 平衡状态。

⑴根据域的应力确定所点处的不Flac采用的显式计算循环说明于图3-1，对

平衡力。 于每个时间步长，对网格的每个网点求解运动方

⑵由于平衡确定所点应变分量。 程式。对于一个非平衡状态，出现在每个网点处

的力不是平衡的。根据不平衡力分量及集中于网⑶通过积分，求网点位移 点处的域质量，这格产生网点加速度。积分这些加速度，便得到网点速度，随后，利用这些速度来

应力－应变规律 求应变变化⑴。

对每个域系统 把这些应变用于本构规律，以求该域的相应

应力增量。一旦求运动增量后，把它们在周围各网⑴由网点速度求应变率。 点上引起的力加起来，便求得合成的不平衡力。 ⑵由应力－应变规律计算应力增对每个时间步长，重复这个计算循环。用户量 可以考察该过程任一阶段时，问题的当前状态。如

8

上所述，在Flac内部，上衰减体系的动力响应，便提供一个静力解答。用于该码的衰减方法后面将予以讨论。当随着增加时间步长问题趋于静平衡时，衰减使不平衡力蜕变为零。平衡时，需要的时间步长数目与很多因素有关，其中包括用户所希望的解的精度，材料的屈服程度，以及问题的大小。

图3-1 码运行概述

⑴

3.3有限差分格式(Finite Difference Scheme)

3.3.1 导数的表示－Wilkins（1963）根据偏导数的积分定义，提出了一个差分格式： F1limFnids (3-1) A0Asxi

F某一纯量、矢量或张量；式中

xi位置矢量分量；

A－积分域；

ds－弧长增量；

ni垂直于ds的单位法线分量。

(3-1)式的面积分是连续的，但沿一个有限多边形进行积分时，可以写为一个等效的（虽然近似的）表达式：

n F1NnFikxk (3-2) xiAn1 n式中 N 边数；xi n边n F边n上的F平均值；

A n

xi－边n的矢量长度的分量； 01图3-3 数值积分时边界的离散 ik二维置换张量（permutationtensor)10



用(3-2)式来推导Flac中的所有空间差分方程式。应注意，这种表示没有外形及边数，不同于以矩形网格为依据的许多有限差分表示法。 3.3.2运动方程式——用熟悉的表达式给出的运动方程式为：

iju(3-3) gi txy

式中 —密度；　　　　　　　　　ij—应力张量；m

F→→—速度；　　　　　　　u　　　g—体力分量；

t—时间。 ,u,uu 可以简单例子，试考察受随时变化作用的某质量的运动

方程式(图3-4)： 图3-4受随时间变化

Fu　　　（即Fma)tm 9

的质量

(3-4) 可以用包含半时间步长的速度的中心差分格式⑵表求解(3-4)式，((3-4)式左边的))加速度可以写为

ttuttuu (3-5) 22 tt代入(3-4)式后，得 (3-6) t t t  t  uuFtt22m 

这是半个时间步长时的网点速度。现在，可以用一个附加积分由速度求得位

(tt)t (3-7) 移为：u(tt)u(t)uttt uttutu2

构成图3-5示说明的计算顺序后，使可以求同一时间增量改间的力。

3.4 速度／应变增量方程式(Velocity／Strain Zncrement Equations)

Flac由每个网点处的速度，确定每个时间步长时的一个初应变增量．把这个应变增量用于所选的本构规律，该确定相应的应力增量．在增量形式中，应变张量为：

jiu1u(3-8) ijt2 xjxi式中

ij应变增量张量，i,j1,2；

ii向速度分量，ui,j1,2；

xii向坐标分量，i,j1,2；

t时间步长。

应用运动方程（式（3-3））后，在一个典型四边形单元的每个角隅的网点处有一组x，y向速度。Flac把这个四边形再划分为两对重迭的常应变三角形域。命令为a，b，c，d，（图3-6）。 回想（式（3-8））当计算应变增量张量时，需要速度梯度。根据式（3-2），可以用他们的差分来代替偏导数： 1uiikxku （3-9）

xjA边域

这个求和是沿域的边缘取的，但只知道角隅处的速度。如果沿每个边缘取平均速度分量，当沿边缘的速度为线形变化时，该求和内的项是与由精确积分式（3-1）求得的相同。Flac先三角形a及b求速度梯度，然后对c及d求速度梯度。对于三角形a，式（3-9）的展开式为：

NWSE u1122313uuxuuxuuxiijkkiijkkiijkk xj2A 为此说明导数的一个分量的上述展开式，该考察应变率：

1u12132113 11u12x212u13x21i3x2ux2ux2uux2x12A

 10

（3-11） 展开并消项后，有

1u1232311 211x2113x2ux2ux2x2ux2x12A（3-12） 用同样的方法求所有的其于速度梯度分量，并用于求这些应变增量张量的三个位移。

这里消除了有时在应用有限差分时经历的“水漏”变形状态问题，因为在计算中看到了三角形域，而不是四边形域。 为了防止塑流过程中“锁住”，在三角形a，b间及c，d间采用了函数格式。在网点处取两组三角形的平均力。

3.5应力-应变规律

一旦确定了每个三角形内的应变张量，必须根据该域所选用的本构规律计算相应的应力。在正式考察本构规律之前，需要对原有的应力作两个修正。第一，如果采用大应变逻辑，须就网点位移引起的转动应力作修正。这将在3.6节中讨论。第二，如果须有效应力，必须从直接应力中扣除孔隙压力。应变与应力间的本构关系给成增量形式：

ijfij,ij, （3-12）

2； 式中ij应力增量张量i，j1，

f包含有增量应变，现存总应力状态，材料常数等的本构方程。

本构规律及应变计算应力分量的细节给在附录A中。

最后，由加权平均三角形应力分量的方法来确定四边形的平均应力分量。Flac允许打印或者三角形或均应力分量。

3.6确定网点处的不平衡力

一旦四边形域内的应力张量已经求得，便可以用它来计算各网点处的不平衡力。重写运动方程式：

（3-13） iijup xpgij

注意，它的右边可以用式（3-2）代替，并通除以p。其结果为 ui1p （3-14） Aijjkxkpgi或

iuF

igi m式中 FxmpAiijjkk

取产生的各网点处的所有矢量，因此，上式中的质量m取为由周围各域凝聚成的质量（图3-7）。式（3-14）的求和项取为出现在该图说明的闭和回路。

如果取闭合回路为图示路径，并假设凝聚质量等于由路径所包围的质量，那么（1）由绕各网点的全部这样的路径是精确地覆盖了整个网格区；（2）凝聚质量满足一定的条件，即重心。

虽然，域内的回路路径对质量凝聚是重要的，但它对式（3-14）中的求和项是不重要的。因为，域上的域应力为常项，只要终点不变（在网格线中心），对于任何路径，求和给出相同的结果。因此，如图3-7所说明，网点力变为：

 11

14213243aFi ijjkxkxkbijjkxkxkcijjkxkxkdijjkxkxk （3-15）

这个力在所有应力和坐标，在时刻t，t+△t等时都是已知的。在下一个半时间步长时的网点速度为：

t （3-16） Fiittuittutgi 22m应该指出，如果域中的任何一个搞错了（如，在边界处或如果这个域是空域），相应这个搞错了域的项由式（3-15）中直接省略掉。

如果选择FLAC的大应变模型，把修正的网格坐标用于式（3-16）的新速度：

tt2tttixixiut （3-17）

显然，FLAC不使用位移，但为了方便用户，也用同样的方式对位移矢量进行修正。

3．7应力转动修正项

当利用代码中的大应变逻辑时，由于物体的转动，必须修正应力。应力的正规转动方程式为：

（3-18） ijJiJj式中

cossin (二维)Jij sincos  ij大应变转动时修正了的应力张量.。对于很小的角度，因为cos1,sin,则

1 (二维)Jij 1 因此， 1111J11J1121J12J1112J11J1222J12J12 （3-19） 112112222

11212(因为2112，很小时，20）

这就是参数于与旧轴成角的新轴的应力变化—但我们需要的是参数于旧轴（原始轴）的新应力。因此，11的修正项为

11212 （3-20） 同理，有

2211J21J2121J22J2112J21J2222J22J22

2112211222

（3-21） 22212112

反号并略去2项后，有

22212 （3-22） 还有 JJJJJJJJ12111121211221121122221222 （3-23）

112121222

121122 （3-24）

 12

当采用大应变逻辑时，把这些修正项加到各应力中。

参考文献

13

4.0输入指令(Input Instructions) 4.1定义(Definitions)

在详细讨论输入命令之前，在这里先重温几个定义。关于术语及问题构成的细节见第5.0节。

Zone（域）—有限差分网格的一个单元—类似一个有限单元。

Gridpoint（网点）—有限差分域的角隅节点，由四边形的四角网点确定该四边形域。相邻域在其网点相连接。

Grid（网格）—有限差分域的集合—类似于有限单元网。

在正常运行中，有限差分网格的原点取在网格的左下处。ix及jy轴从这个点开始（图4-1）。因此，每个域及网点具有ix及jy的编号和它们相联系。

因为FLAC是一个显式码，以逐次时间步长来求解运动基本方程式。代码是动态的，它的意思是每个网点具有一个必须衰减以便给出静力解的逐变分量。当代码变到平衡时，这个衰减在FLAC中自动进行。当这些速度达到某些足够小的值时，便得到平衡解。用以确定平衡的手段将在附录B中讨论。 4.2输入命令(Input Commands)

输入Flac命令不同于常见的数值模拟计算机程序，为了操作简单和易于使用，规定了专门的命令。FLAC的操作可以是“交互”式的（即通过键盘发命令），或“文件驱动(file-driven)”式的（即数据存在一个数据文件里，并以软或硬盘的形式读取）。不论那种方式，运行的命令是相同的，由用户的爱好选用。

输入命令是按调序排列的，视其重要，它们包括了一个基本命令字，后跟一些总键及数值输入。下面给出的各命令是通字地键入在输入行上的。应注意，只有前几个字母是大写的。程序只要求键入这些大写字母，以便命令“识别”。许多关键词后跟了一系列关键词要求输入的数值。以小写字母开头的词代表。以I,j,m或n起头的词为整型，分别为实型词（10进制小数）。实数的小数点可以省略，但在任何整数中不能出现小数点。这些值可以用任何个空格或下列定界符中的任一个隔开：（），/ =

你将看到，带其输入参考的附加符号。它们是： 〈 〉—表示必选参数（括号不键入）； „ —表示可以给出任意个这样的参数。 在输入中，“*”或“；”后面的内容为注释内容，不予考虑。如果以分批（即文件驱动）方式运行时，在输入文件中使用这种注释是有用的，因为可以把这些注释掉的内容可显示输出中。

在FLAC中使用以下符号约束，在进入输入时必须记住。

正的值：向下及向右运动时；拉应力及压力时；伸长应变时；剪应变时如右图示

由力矢方向确定力边界条件的正负号（即当指向正坐标轴方向时为正）。虽然我们推荐使用国际（SI）单位，它可以使用任何性能协调的工程单位。在表4-1中，给出了各种单位体系的述评。在程序中，又有摩擦角及膨胀角设作约定，它们按度（°）键入。

表示4-1 单位系统

m 3 密度 Kg/m力 N 应力 Pa 2 重力 m/sec长度

m 33 10kg/mKN Kpa 2 m/secm cm 636310kg/m 10kg/m MN Mdynes Mpa Bar 22m/sec cm/sec 14

ft in 33Slugs/ft Snails/in Lbf Lbf 2Lba/ft psi 22Ft/ sec In/ sec 以下各例具体说明各输入命令的使用。为了固定I=1~10及j=1~10范围内各网点沿x向的位移（或速度），下列各命令会产生相同的结果：

fix x I=1,10 j=1,1 或 fix x j=1 I=1,10 或 fix x I=1,10 j=1

如果I或j之一的范围被略去，这时则取网格的整个范围。

另一个例子，说明利用输入行中的任选参数。为了建立一幅最大长度尺度为0.001的位移矢量图，给下列命令：

plot disp/max=0.001

当程序“认识”“displacement”这个字时，便寻找补充说明或“开关”字。这时，字“max”设置画图的最大位移矢量值。此外，应注意，在一些情形中，在同一输入行中可以给出许多关键词。例如：

plot disp/max=0.001 Sxx/int=10E6/red

将产生一幅如上所述的位移矢量图，并用红色去附加的应力等值线，其间格为10。应注意，上例中的/和=是为了方便用户严格地用作分界符，如愿意的话，可以删去它们。

输入命令一览表

Apply keyword=value

以下的keywords可以利用:

Pressure—压力 XForce—x向力 YForce—y向力

可以把压力作用在各网点行或列的各处，它们或者沿一个边界作用，或者作出在物体内部。当压力作用在i或j的转递增序列时，左起的压力为正。

Xforce及Yforce可以作用于任意列或行处。网点的任选范围是与力作用处的网点的。形式为ii1,i2,jj1,j2的范围必须维持在输入行上，但不管i或j的先后次序。

Call filename

遥控输入一个带文件名的文件可以用Call命令来执行。任何一个组输入指令可以代替该文件以便以遥控或批处理方式运行。命令ReTurn必须为遥控文件中的最末命令，以便把输入回到局部或交互方式。目前，不能Call命令本身作为这个文件。 CUtoff

对于指定为Mohr-Coulomb的所有域，设置Mohr-Coulomb塑性及多处同时存在的拉力截止的节理模型(Ubiquitous joint models to sion cut-off) Fix x

Y

X Y

用这个命令，防止由定义域内的网点改变X或 Y向的速度。具有形式ii1,i2,jj1,j2，不管i或j的先后次序。如果想固定位移，其相应的后初始化为零（这是启动时的默认值）。可以把Fix与INI XV,YV命令一起使用，以便提供一个刚性移动边界条件。注意，如果关键词MARK已给，只有在内标注的网点被固定。

FRee x

Y X Y

15

它与FIX命令相反。它释放网点约束。

GEn x1,y1 x2,y2 x3,y3 x4,y4

对于整个网格或子域来说，可以生成各坐标。如果给出上述GEN命令形式，把（由行、列范围确定的）四边形域再划分区域。范围中的最小行、列号的网点位于坐标（x1,y1）处；其余范围的角点，按顺时针向进程的方式，其坐标（x2,y2），（x3,y3），（x4,y4）。可以用字SAME代替任何一对或对坐标。这时，特殊网点将维持其现行坐标。除非给出了任选关键词RATIO，在四边形域内部，将均匀地配置四边形域内的其余网点。这时，按照参数ri或rj，可以域的间距增大或减少，其中ri或rj分别右i或j方向相继域尺寸间（between succesive zone sizes）的几何比值。例如比例1.1将产生比原有域大1.1倍的相继域(succesive zone)。 GEn Circle xc yc rad

Arc xc yc xb yb theta Line x1,y1 x2,y2

GEn命令形式迫使网格的规定区域符合已知外形，如圆、直线等。这些GEn命令用于定义网格内部的内面形状，这些形状可用来定义开挖或结构物，诸如层状沉降物或物体（bedded deposits or are bodies）。以下形状都是有效的。

ARC—绘制出符合一个圆弧的网格，圆弧的中心为（Xc,Yc）边界起点为（Xb,Yb），并反时针转theta0。

CIRCLE—绘制出半径为rad而中心为（Xc,Yc）的圆形网格。 LINE—绘制出端坐标为（X1,Y1）和（X2,Y2）的直线网格。 注意，构成所给外形部分的各网点都是“有记忆的”（见MARK及UNMARK命令）。于是，可以通过带INITIAL，MODKL及PROPERTY命令的关键词REGION来寻找自带有记号网点所包围的整个区域的地址。例如，可以用命令mod null reg=i,j来取消一个圆内部的区域。其中i,j在圆内的一个域。运用GEN的例子给在

第5.2.1节。

应小心地运用GEN命令。某些几何图形用四边形单元是不可能构成的——例如，由两根相交的线定义的细长三角形区域。此外，用户应知道，一旦网格MARKED（有记号），就不能再用另外的GEN命令来移动它们。如果想移动已MARKED（有记号）的网点，应先返消（UMARKED）该点的记号。 GEn ADJust 网格本身调整，试图获得一个较光滑（smoother）的离散。那些有记号（MARKED）的网点及边界是都不移动的。逐次ADJUST命令可以给出更绘制的网格。其例子可以在第5.2.1节中找到。 Grid icol_jrow

规定了计算网格各域的行列数目，对flac的0K版本，约可采用Mohr-coulomb模型的2000个域。

His

对于网点或域i1,j1,每NSTEP时间步长储存其变量〈Keyword〉历史。在任何时刻有达25个历史，每个历史可储备达1000个。每个历史有顺序地从1起编号。用户记着各历史的次序，因为在绘图或打印时，必须要求提供一个指定的历史编号。对于所有的历史NSTEP必须是相同的，且仅需要在跨入时间步长之前只一次给出；如果没有给出，其默认为10。

历史值是储存硬盘上的直接读取文件中。当终止Flac时（Stop），这个文件被抹去。因此，如果用户想储存（Save）这个文件，必须建立一个储存文件或把这个历史写入一个格式化了的软盘文件（disk file）（见HIS WRITE N命令）。可以跟踪网点或域i,j的历史的keywords为：

Ang—对于i,j域，从正x轴反时针向转至最小主应力的角度； PP—i,j域的孔隙压力； SIG1—i,j域的最大主应力； SIG2—i,j域的最小主应力；

16

SXX—i,j域XX应力； SYY—i,j域YY应力； SXY—i,j域XY应力； X—网点i,j的X坐标； Y—网点i,j的Y坐标；

XDis—网点i,j的X向位移； YDis—网点i,j的Y向位移； Unbal—最大不平衡力；

Xvel—网点i,j的X向速度； Yvel—网点i,j的Y向速度；

XXA—域 i,j的三角形a的xx应力； XYA—域 i,j的三角形a的xy应力； YYA—域 i,j的三角形a的yy应力； XXB—域 i,j的三角形b的xx应力； XYB—域 i,j的三角形b的xy应力； YYB—域 i,j的三角形b的yy应力； XXC—域 i,j的三角形c的xx应力； XYC—域 i,j的三角形c的xy应力； YYC—域 i,j的三角形c的yy应力； XXD—域 i,j的三角形d的xx应力； XYD—域 i,j的三角形d的xy应力； YYD—域 i,j的三角形d的yy应力；

His keyword

以下各HISTORY命令形式允许用户绘、写或重新设置历史。写入屏幕或硬盘上的硬盘文件共nhis个历史数目（nhis=1~总历史数）。可以选择以下keyword。

Dump_nhis—历史数目nhis的历史（时间步长，历史值）被写入屏幕

（Screen）。如果该历史大于页长（page length），将滚过屏幕。利用 ctrl_numlock键来终止滚动。

Write_nhis—把历史数目nhis的历史（时间步长，历史值）写入硬盘上的

Flac.HIS文件。终止Flac后，可以打印或操作该文件。逐次的各HIS Write命令将顺序地相加Flac.His文件。不过，所写的前面文件将冲掉改写现存的flac.His文件。 Reset—清除命令。

Initial keyword_value_<„>_

为某些网格变量赋初值。可以从以下一览表中选择一个或一个以上keyword：

PP—孔隙压力； X—x坐标； Y—y坐标；

SXX—xx应力分量； SYY—yy应力分量； SXY—xy应力分量； XDisp—x向位移； YDisp—y向位移； XVel—x向速度； YVel—y向速度；

为了协助执行INITIAL命令，提供了几个任选keyword。它们为：

Mark—仅对标有记号（MARK）的网点初始化网点变量（其它网点保持原样）。

17

Region i,j—对标有记号（MARK）的网点的连续序列所包围的各个域进行初 始化。域i,j是有记号（MARK）区域内的一个域。

Var xv,yv—在整个指定的范围（或区域）上，给定参数的变量，这里xv为x向变量，yv为y向变量。在具有最小i,j量的域（或网点）与具有最大i,j量的域（或网点）之间，应用这些变量。

用作标志时，可以给出各域或各网点的任选范围。形式为ii1,i2,jj1,j2的范围必须维持在输入行上，但不管i或j的先后次序。

例如，如果想初始化水平应力随深度的变化率（gradient），可以发布下列命令：

flac>ini_sxx=-30e6_(Var=0,_10e6)

这里，网点的xx—应力为-30MPa，而网项为-20MPa，在这些点之间按y线性变化。如果规定i或j的范围，则在整个范围取变量；变量可以为正或负。

可以利用INITIAL命令，把网点移到网格的另一位置。在执行这一措施时，务必小心flac内在地把各四边形域划分为两组重迭的成对三角形域，以便进行计算（见图3-6）。可以用任何方式使四边形变形，但应服从下面两个准则： （1） 四边形的面积必须为正的；

（2） 至少一对三角形域的每个三角形必须具有大于四边形总面积的20%

的面积。

如果这些准则之一不满足，在进行时间步长中，flac将给出错误信息“BAD GENME”以下图例说明可能的域的变形。

在运行前，带用赋零值来初始化位移。这不会影响计算，并在输出时，会给出位移的变化。如果不初始化它们，则输入的位移为总位移。

注意：通过运行来部分地改变坐标，是没有意义的。用户可以规定任意的初应力，但它们可能是不平衡的——这时，应先使模型平衡。 INTerface—n_keyword__ n_keyword=value

用INTerface命令来激活预释放的界面模型（The pr-release

interface model is actirated）。籍助于一个“Interface”可以使一个flac网格的任何部分与任意其它部分相互作用，这“interface”由每个相互作用节点处的法向及剪切弹簧（spings）所组成。实际上，界面相当于一个断层或节理（fault or joint）——即，两个表面间的“软衬层”（soft cushion）。必须事先说明每个界面，其方法是：规定两组可以相互作用的边界网点，并规界面的刚度，粘聚力及摩擦性质。可用于这个命令的keyword为：

Aside i1,j1_ Bside i1,j1 Cohesion Value Friction Value Glued Value KN Value KS Value TBond Value Unglued

Aside及Bside相应于界面的两个侧面，而毗连的各对i,j定义了位于相互作用边界上的各网点。必须安排好各对i,j的顺序，以便代码知道那一边界面包含坚固材料（Solid material）。采用的约定如下：如果人沿边界面“步行”，坚固的网格（Solid grid ）必须在右边。按人步行的方向，从界面上的第一个网点到最后一个网点，对各网点编号。

n是所给界面的唯一标志，是顺序的。

为了构成一个相互作用边界，应先在界面所在网格的两个区域内，表明一

18

行空域。如果用户想使这两个区域处于初始接触，那么须用GEN或INITIAL命令，使网码的这些对立区域移动，以便在跨入时间步长之前，使它们接触。

必须定义每个界面的材料特性（粘聚力、摩擦、法向刚度KN、剪切刚

度KS和拉伸强度TBOND）。粘聚力及拉伸强度的单位按应力表示（即Pa、Psi等）。摩擦用度（°）量测，而刚度为每单位长度上的应力（pa/m,psi/in）。采用的这一组单位必须与问题的其它单位一致。

这两平界面，在任何时候可以“胶合”（glued）在一起，以便构成一个

粘着的界面，不再发生滑移或分离。这种作用的结果是一个标准的Flac网格。GLUE命令的对应命令为UNGLUE。关于利用INTERFACE命令的更详尽细节给在附录D中。 MArk

对所有给range内的各网点作记号。一个有记号的网点在任何情况下不影响计算，但对于INITIAL，PROP及MODEL命令来定界一个区域是起作用的。注意，GEN命令能自动地对各网点作记号。

具有形式i=i1,i2 y=y1,y2。 i或y何者在前不限。 Model keyword__

这个命令把一个本构规律和一个相应于由规定的域范围的网格区或i,j域所在区域联系起来。关于一个Flac网格内部定义的区域的更详细见5.2.1节。在计算循环中，各个域将查阅相应于下列keyword之一本构子程序： ANisotropic—横向各向同性弹性模型； Elastic—弹性各向同性本构模型；

Mohr-Coulomb—Mohr-Coulomb塑性模型；

Null—零域（零域被用于模拟开挖了的材料。零域内的应力都自动赋零）； SS—应变软化；

Ubiquitous—普通存在节理的模型。

具有形式ii1,i2,jj1,j2,i或j何者在前面不限。另外，模型的区域可以定义为REGION=i,j。这时，由一个MARKED网格的连续序列所包围的所有域将赋于该已给模型。

NEW 这个命令允许不脱离Flac的情况下开始一个新问题。 PLot keyword_„>_

这个命令要求在屏幕或Hewlett_Packard笔式绘图机上作一幅画。注意，除非默认窗口已满意，在先用WINDOW命令来定义画面大小。通过在一行上给出几个keyword方法立即可以画出几个变量。在打印或绘图之前，必须（利用MODEL命令）给网格一个材料模型。各keyword及其意义如下： BEam—绘各结构单元的几何形状；

Boundary—绘网格的边界，包括零区域的内边界； Cable—绘钢索锚栓（Cable bolt）几何形状； Disp—绘位移矢量；

E_P—塑性应变（若仅用应变软化材料模型）； Grid—绘有限差分网格；

His nhis—由记录在历史量nhis内的变量组成的一幅时间终历图； Pen—把图形转向Helett_Packard 7470A或7475A笔式绘图机（默认

值为屏幕），如果要笔式绘图，必须作为第一个keyword出现；

PP—孔隙压力等值线；

RE—画成矢量的固定网点处的反力；

STAte—绘在各单元形心处的材料现行屈服状态； STress—主应力矢量；

SXX—xx应力（总应力）等值线； SYY—yy应力（总应力）等值线；

19

SXY—xy应力等值线；

Velocity—绘速度矢量(as arrows)； WAter—绘地下水位线*； XDisp—x位移等值线； YDisp—y位移等值线； XVel—x速度等值线； YVel—y速度等值线；

利用材料特性名称作为一个keyword可以绘任何材料特性的等值线图。 如果不给出更多的参数，将自动选择尺度因子。不过，可以在任何keyword后，跟上一系列“Switch”，它们是为keyword本身赋予某些绘图特征的关键词。每个Switch对它前面的那个keyword起作用。以下的各Switch是可用的： Color—根据前面用SET PAL命令给出的调色板的选择和图形监视器类型，给出线条颜色的选择**。对于CGA调色板非0，线条颜色选择为Green,Red和Yellow-Brown；对于调色板非1；线条颜色选择为Cyan,Magenta和White。对于，颜色选择列表示下：

Palette 0(调色板) keyword 绿色（默认值） GREen 红色 Red 黄色 Yellow 深红色 MAGenta 褐（咖啡）色 Brown 白色 White 淡绿色 LGreen 淡红色 LRed 淡兰色 LBlue 淡深红色 LMagenta 灰白色 GRAy 强白色 IWhite

黑色 Black Polette1 keyword

兰绿色（默认值） Cyan 红色 Red 黄色 Yellow 深红色 MAGenta 褐色 Brown 白色 White 淡绿色 LGreen 淡红色 LRed 淡兰色 LBlue 淡深红色 LMagenta 灰白色 GRAy 强白色 IWhite

黑色 Black

Interval=c—等值线间隔为C。在发布PLOT命令之前，通过直接发布SCLIN

命令可以把等值线的值输给屏幕或笔试图形（见SCLIN命令）；

Max=v—对于所有矢量图，设置相应于最大箭头长度的值为v。而所有其它矢

量将以这个最大长度为尺度；

Noh—消除图形的标题；

Zero—产生被略去的0等值线。

20

在给出PLOT命令之前，利用WINDOW命令可以作一幅放大或缩小的图形。这时，该窗口将保留，而在重新启动时将被记住（The window will then remain set and will be remembered on restart.）。

所提出的建立一个笔试图形的方法应该先建立一幅你已满意的屏幕图形。一旦完成这一步应键入：Plot pen（后跟用于建立屏幕图形相同，一些keyword及Switch）。用keyword.pen来初始化Hewlett-packard 7470A绘图机*，且须先于控制图像的keyword。关于绘图的最后提示——在用Flac作图之前，应定义用作输出装置的绘图机。对于更多的信息，见Set Output及Set Baud命令。

Print keyword__„__ 打印输出用于外荷载，极限尺寸（limits），主要网格变量或扩张变量（Extension Variables）。对于用确定的各网点或域的范围，或位于域i,j内的区域可以作输出。具有形式i＝i1,i2,jj1,j2；i或j的先后次序不限。另外，可以把模型区域定义为REG=ON=i,j；这时，将打印出由标有MARK的各网点的延续系序列所包的域的输出。域i,j是MARKED区域内的一个域。与绘图情形一样，在规定材料模型和特性之前，不会打印网格个各变量。各keyword为：

Apply—作用于模型的力或压力的值和范围； Fix—打印x或y向的各固定网点； Limits—对SOLVE命令的（limits）； MArk—打印有记号的各网点； MEm—占用的内存；

SIrct—打印各结构单元节点处的力弯距和位移；

Interface—打出界面数据，包括节点力及单位法向矢量。 主要网格(Main Grid)keyword

Bulk—体积模量； COhension—粘聚力； Den—物质密度； Dilation—膨胀角；

E_P—总塑性应变（应变软化模型）； Friction—摩擦系数； SHear—剪切模量； X—x_坐标； Y—y_坐标； PP—空隙压力；

JFric—节理（joint）摩擦； JCoh—节理粘聚力；

JAngle—节理角度(Joint Angle)； NUYx—y-x泊松比； NUZx—z-x泊松比； SIG1—最大主应力； SIG2—最小主应力； SXX—xx-应力； SYY—y-y应力； SXY—xy-应力；

Theta—从x_轴反时针向转至最小主应力的角度； XDisp—x-位移；

21

*可以用任何Hewlett-Packard或兼容的笔式绘图机。 YDisp—y-位移；

Xmod—x-模量(Modulus)； Ymod—y-模量； XVel—x-速度； YVel—y-速度；

ASXX—xx-应力（三角形a）*； BSXX—xx-应力（三角形b）； CSXX—xx-应力（三角形c）； DSXX—xx-应力（三角形d）； ASYY—yy-应力（三角形a）； BSYY—yy-应力（三角形b）； CSYY—yy-应力（三角形c）； DSYY—yy-应力（三角形d）； ASXY—xy-应力（三角形a）； BSXY—xy-应力（三角形b）； CSXY—xy-应力（三角形c）； DSXY—xy-应力（三角形d）； State—塑性状态；

0—弹性；

1—屈服处（at yield）；

2—弹性但以前在屈服处（elastic but previously at yield）； 3—在单轴拉伸时，已超过拉力截止（has surpassed tension cut-off）；

4—在屈服处，亦超过拉力截止； 5—在一般拉伸中，已超过拉力截止；

6—目前在屈服处的多处同时出现的节理（ubiqutous joints currently at yield）；

7—在过去已屈服，但目前为弹性状态的多处同时出现的节理（ubiquitous joints have yielded in the past but are currently in the elastic state）；

Iables—将打印储存表格（Stored Tables）； Xreaction—x－反力； Yreaction—y－反力。

PROp keyword=value_<„>_

这个命令为由MODEL命令确定的本构模型赋特性值。为了协助PROPERTY命令的执行，提供了一些任选keyword。它们是：

Region i,j——由MARKED网点的连续序列所包围的各个域将具有规定的

特性。域i,j是MARKED区域内部的一个域；

Var xv yv ——可以给出整个规定范围（或区域）上的特性的变化

（variation）。其中，xv为x向变化；而yv为y向变化。这些变化被应用在具有最小i,j值的域与具有最大i,j值的域之间。

可以给出一个由确定的各域的任选范围。形式 i＝i1,i2,jj1,j2；i或j的范围必须维持在输入行，但i或j先后次序不限。

各特性的keyword，它们的含义及使用它们的本构模型给于下面:

22

*应记得，Flac把每个四边形分为四个三角形域：a,b,c,d.

keyword Description Model

Angle 各向异性角，从x轴逆时针 横向各向同性 向取值。

BUlk 弹性体积模量。 弹性

Mohr-Conlomb 应变软化/强化 普遍存在节理（ubiqutous

joints）

Cohesion 粘聚力 Mohr-Conlomb CTab,DTab,FTab表格数 应变软化/强化 Density 物质密度 应变软化/强化 DIlation 膨胀角（度） 所有模型 Friction 内摩擦角（度） Mohr-Conlomb

应变软化/强化 JAngle 角度（度），对于普遍存在节理， 普遍存在节理 从x轴起逆时针向量取。

JCohesion 普遍存在节理的粘聚力 普遍存在节理 JFriction 多节理面的摩擦角 普遍存在节理 Shear-mod 弹性剪切模量

Xmod 各向同性平面内的弹性杨氏模量 横向各向同性 Ymod 垂直于各向同性面内的弹性杨氏模量 横向各向同性 NUYx 由垂直于各向同性面内的单轴应力

引起的各向同性面内的正应力泊松比 横向各向同性 NUZx 由各向同性面内的单轴应力 横向各向同性 引起的各向同性面内的正应力泊松比

模型类型及其特性，必须在打印或绘制任何数据之前，给以定义。如果所给特性与所选模型不符合，将给出一个警告信息，以通知用户不接受非必需的特性。各模型需要以下特性：

弹性：(1)剪切模型 (2)体积模量 (3)密度

Mohr-Coulomb：(1)剪切模量S (2)体积模量E (3)密度D (4)摩擦角f (5)粘聚力C

(6)膨胀角（任选项） 横向各向同性:(1)剪切模量* (2)x-模量 (3)y-模量 (4)密度 (5)NUYx (6)NUYz *必须确定各向异性弹性的横向剪切模量G。荃于岩石的实验室试验，提示了以下方程式:

假设x-2平面各向同性平面（见附录A） ExEyGXYEx(12xy)Ey

23

普遍存在节理：(1)剪切模量 (2)体积模量 (3)密度

(4)粘聚力(末终损坏的材料intact material) (5)摩擦角(末终损坏的材料intact material) (6)节理粘聚力(节理) (7)节理摩擦(节理) (8)节理角度(节理)

(9)膨胀（末终损坏的材料，任选项） 应变强化/软化：(1)剪切模量 (2)体积模量 (3)密度

(4)摩擦角(初始值－任选项) (5)粘聚力(初始值－任选项) (6)膨胀(初始值－任选项)

(7)ftable(表非与塑性应变间的关系) (8)ctable(表非粘聚力与塑性应变间的关系) (9)dtable(表非膨胀角与塑性应变间的关系)

应变强化及软化的性能由摩擦、粘聚力及膨胀的变化控制，它们是由表格值规定的塑性应变的函数。规定这些值的细节，见附录A。注意，如果将值均为0，这些特性将取所给的值（即，用COHESION,DILATION或FRICTION各Keyword）。

REstore filename

从名为filename的文件中，恢复一个以前储存的问题的状态。 RETurn

这个命令应维持在一个遥控输入文件内（见Call命令）。 SAve filename

在目前状态时，把内存中的现行问题储存在名为filename的文件上。如果已存在这个文件，将被覆盖。 SClin n_x1,y1_x2,y2

这个命令在屏幕或笔试绘图机构成的任何等高值线图上，标注等值线。用户定义一个他想知道等值线的值的“扫描线（scan line）”。用户用命名SCLIN来定义该线的边界点。在扫描线与等值线间的每个交点处，写上字母A~Z。给出一个标题，以表示最小和最大等值线的值，以及等值线间隔。这允许计算两端值间的任意等值线的值。缺省条件是无扫描线。在每个新图形之前，必须给一个新的扫描线。对于一幅图，可以规定多达5根扫描线。其参数如下：

n—扫描线数目（必须1，2，3，4或5）； x1,y1—扫描线始点坐标； x2,y2—扫描线终点坐标。 例如，命令

flac>Sclin _ 1_0,1_0,10

将产生从0，1到0,10贯穿扫描线的各等直线的值。 SEt

这个命令用来设置在一个运行过程中，始终保持不变的全程条件变量。该命令的各keyword为：

Aspect_a——其中a=x:y,x、y是从一个应为方形的非方形图像量得的。根

据图板（graphics beard）和屏幕垂直尺寸调节，a一般为0.5～1.5。 Back_background_color

24

0----黑色，1—蓝色，2—绿色，3—兰绿色，4—红色，5—深红色，6—黄褐色，7—白色。 BAud b—为输出装置设置波特率。b选择为1200，2400，4800及9600bpm。

波特率只能为COM1及COM2设置。

CGA或CSC—设置320×199图像分辨率（pixel resolution）的彩色屏幕

图像模式。这是默认值。

Col_n—n为最大输出列（column）数。注意，在flac运行之前，用户必

须设置适当的打印列模式。有关设置打印机模式的方法，试查阅你的打印机手册。

EGA—把图形模式设置为0×350图形分辨率的增加图形转换器。 Force f—设置SOLVE命令的不平衡力极限。

Gravity_g—用大小为g，方向角为来规定重力，从负的y轴

量起。

Large—大应变（各坐标均为适时修正）。

Log_ON—在默认硬盘驱动区中，打开一个名为flac.log的文件。如果已

存在flac,log文件中，将覆盖它。但打印于屏幕的文件也写入log文件。这对保存交互式对话一段时间的记录，特别有用的。可以编辑这个文件，以建一批数据文件。

Log_OFF—关闭记录功能（Logging function）。它不封闭记录文件

flac.log。如果在这一段时间之后的每一阶段，给log on进行设置，后来的输出将附加到这个文件里。 Mono—设置单色分辨度图形模式。

Output =p—把画出的输出转送给与通道p相连的设备，这里p可以为

COM1,COM2,LPT1或任何其它通道（port）。

Pal 0—0号调色板。给出线条颜色选择，在CGA模式中为红色和黄褐色。 Pal 1—1号调色板。给出线条颜色选择，在CGA模式中为深红色及白色。

对于EGA模式的颜色，是用PLOT命令给出的。默认值=1号调色板。 P_Stress—平面应力（仅对弹性及横向各向同性模型）。 Small—小应变（不适时修正坐标）。

STep n—为SOLVE命令设置步长数的极限值。

Time t—为SOLVE命令设置运行时间的极限值，单位为分。 SOLve_

该命令允许自动检测问题的稳定解，直到达到时间步长，运行时间或不平衡的力的限定条件为止，都执行计算任务。用以下任选的keyword来改变这些：

Step—时间步长； Time—运行时间（分）； Force—不平衡力。

这些keyword的默认值为：

S=500步，T=5分，F=100

应该注意，F取决于所采用的极（system of limits）。如果没作明确的规定，那么将采用旧的极限（old limits），但当设置或重新起动时（when set or on restart），会记得这些的。不过，当给一个New命令时，他们被重新设置为默认值。一个大于1440分（24小时）的时间极限是不会接受的；对于更大的运行时间（如，超过周末），可以依次给出几个SOLVE命令。如果在执行过程中，触动了任何键，在完成现行步长后，Flac将返回问题由用户控制。 STEp_n

执行n个时间步长。如在执行过程中，动了任何键，在完成现行步长后，

25

Flac将返回问题由用户控制。 STop_n

Flac停止，且运行时间结束（the session ends）

STRuct keyword„

STRCUT命令用于定义结构单元（structural elements）的几何条件特性及端部条件。STRUCT命令描述了单元类型及几何条件，单元与flac网格的连系，以及与单位有关的特性类型数。用于定义结构单元的keyword如下： （1） 单元类型keyword 有两种单元类型：

Beam---用几何条件及特性keyword来定义一个梁单元，而它们的值紧

跟在Beam的keyword后面。

Cable---用几何条件及特性keyword来定义一个索单元，而它们的值紧

跟在Cable的keyword的后面。

每个梁或索是由节点起点或终点来定义，或者由主网格联动装置（main grid linkages）来定义。可以把梁或索划分为若干段，在每个段在端部用节点联系起来。梁单元结点用三个自由度来描述：x－及y－向的平移和一个转动。索单元只允许x,y向平移。最后，各个特性通过一个特性类型数来规定。图4-2说明了与梁单元有关的条件。 （2） 单元的几何条件及节点联动装置

为了定义单元几何条件，需要以下两个keyword

Begin keyword－它确定了梁或索的一个端点； End keyword－它确定了梁或索的另一个端点。

单元的各端点，用直接跟在关键词BEGIN及END后的三个可能的keyword及/或数值之一来定义：

Grid_i,j－这里，在单元始或终节点与网点i,j处的土壤或岩石物

质之间，建立了一个联动装置；

Node_n－在建立这个结构时，Flac将自动地给各结构节点一个从1~

节点数的连续号码，其次序按输入值而定。通过发布该命令，可以把当前单元的始或末端与节点n耦合起来；

x,y－在关键词BEGIN后面，可用给定的一对x,y来连接，给各个单

元节点一个空间坐标。可以用该命令来定义任意空间方位处的一切单元。如果x,y都非常接近于一个已存在节点或网点，那么，Flac将迫使这个新节点与已存节点或网点一致。

（3） 单元特性类型keyword

可以为每个梁或索规定特性。这可以用把一个特性代码n与每个梁或索相连系的方法来做到。Keyword的形式为：

Prop_n－把特性代码n分配给单元类型。利用命令struct prop＝n

把各特性与特性代码联系起来。

（4） 单元段keyword（Element Segments）

利用以下keyword可以把每个单元或索划分为若干段：

Seg_k－把单元分为k段。例如，如果始端坐标为（0,0），末端坐标

为（10,0），而k＝10，那么，将把该单元分为长为1.0的10个段。

（5） 索的预拉力keyword

一根索锚栓（Cable bolt）可以预张拉到一个荷载t，这时用keyword：

Tension t

为了进一步定义单元条件，提供以下两个STRuct补充命令: STRuct Node＝n keyword

利用以下定义的keyword，可以对每个结构节点，规定各种边界或几何条件。

26

FIx _ _ _－这个任选项允许把节点n沿x向和/或沿y向位移同

或者固定转动（即,可以把一个梁端锁在适当位置或允许转动）。

FRee _ _ _---这与FIX相反。

Pin---它在节点n形成一个铰接（即,没有弯矩）。 Load _fx,fy,m---这允许用户把x向和/或y向的力，或者弯矩作用于节点n。 x，y---它们是节点n的x,y坐标。这用在建立一个新节点时，这时该节点可以包含在结构单元中。这允许用户在单元形成以前，去定义节点号。 STRuct Prop＝n keyword＝Value„

该命令把特性值与一个特定特性代码n联系起来。图4-3与图4-4分别说明了灌注与锁材料特性的关系。以下的各特性keyword可应用：

E＝Value－梁或索的杨氏模量；

I＝Value－梁的惯性矩（第二次面矩）； Area=Value－梁或索的截面积；

Kbond=Value－水泥浆的粘结刚度（bond stiffness of grout〔力/单

位索长/位移〕）；

SBond=Value－水泥浆粘结度〔力/单位索长〕； Yield=Value－索的屈服强度（力）。

Table n x1,y1 „

该命令为内存中的x和y建立一个供flac利用的表格。把这些表格用于定义地下水位（即，孔隙压力计算用）（见WATER命令），粘聚力和应变软化模型的膨胀表格（见PROP命令）。可以存在多达20个表格；数目n（1~20）是表格标志号。可以给出任意多对x,y；它们是不需按次序排列的，因为flac会把它们按增长的x值分类。如果对于相同的表格号给出更多的TABLE命令，则把新的x,y插入到现有表格的正确位置。如果给出一个与现存x值一模一样的x值，则对表格中的相应y值进行修正（在编辑各表格中，该特性是有用的）。可以用PRINT TABLES命令核实一个表格的内容。

Title——把下一个输入的内容作为一系列笔试图形的打印标题及重新起动文件

的记录标题。

Unmak ——取消在内各任何网点上现有的MARK。 WAter

在每个域内，根据位于所列表格给定的地下水位线以下的域的深度，附着一个孔隙压力。这个压力也正比例于水的密度（如上所给Density）及重力垂直加速度（用SET GRAVITY命令给出）。在地下水位线以上的各域具有的孔隙压力为零。

当用Water命令把孔隙压力附着在网格中时，将增加相同数量的直接应力（Sxx，Syy），因为它们全相应于总应力。可以给出或打印出孔隙压力分布；取压力为正。也可以（用Plot WAter命令）显示地下水位线。

在计算过程中，Flac利用有效应力（即，总应力减去孔隙压力）来构成本构模型。在目前的flac版本中，孔隙压力既不受域的体积变化的影响，也没有任何水流。总应力显示在图和打印输出上。

如果给出命令Water Table=0，则从网格中取消所有的孔隙压力。同样，在运行过程中，可以通过开关一个新表格或改变现行表格的内容的方法，来改变压力分布。

WIndow _

这个命令在屏幕上或笔试绘图机上建立一个供绘图用的虚构窗口，xlo~xhi，ylo~yhi的幅面被映在方形屏幕区。因此，如果窗口面不是方形，会以不同的x,y尺度画出变了形图形。

如果在绘图以前没有定义Window，便选择一个能容纳整个网格的窗口。

27

窗口一直保持到改变设置为止。如果希望自动开窗口，试用Auto关键词。如果窗口小于网格尺寸，则会删去窗口边框处的屏幕图像。用户可以利用这个特点去得到感兴趣的细节部分放大视图。

4.3 设置你自己的默认条件 (Setting your own default conditions)

如果你想不论什么时候当你起动程序时，Flac能接受某些参数或模型，试准备一个名为Flac.INI的文件。这个文件可以包含任何有效的各Flac命令。Flac会自动起动来读这个文件，并处理各个命令。例如，Flac.INI可包含以下命令：

Set force=5e2,CGA，back＝1，pal＝0 Water density＝1000 如果flac.INI文件不存在，flac直接恢复而没有错误信息（flac simply coutinues without error）。注意，在一个flac.INI文件中的某些命令可以产生一个错误信息。例如，如果你在定义网格之前，试图给出一个网格的各特性值，则会出现正常的错误信息。 参考文献（Reference）

（1） Lekhnitskii,S.G.Theory of Elasticity of an Anisotropic

Body.Moscow:

Mir Pubishers,1981.

28

5.0 用FLAC解答的问题 （Problem Solving with flac）

5.1 引论（Introduction）

为了用Flac运行一个模拟试验，必须先确定问题的几何条件，边界条件及材料模型。这些过程与所有应力分析码是相似的。必须为模拟的区域划分为若干个四边形单元或域，用域的四个节点或网点来定义每个域。为了后续计算，在内部，Flac再把每个四边形分为四个三角形单元。一旦网格已定义，必须提供边界及原始条件。这些条件包括：

⑴在x,y向被固定的网点位移或速度； ⑵作用于任何面或网点线的压力； ⑶在任意网点处作用于x,y向的力； ⑷物体中的原始应力； ⑸重力。

对于Flac这样的显式模型，问题的实际解答是稍不同于许多平常的隐式模型的解答的。显式码命名为“时间逼近”（time－marching）。在逼近中，以直到达到平衡为止的逼次渐尽的时间来求解运动基本方程式。实际上，对用户系统，这指的是用一个有物理意义的方法来寻求解问题。

例如，试考虑在靠地面的土中开挖的一个隧洞。这时，由重力引起力是主要的。可以最初设置一个大得足以避免边界效应的尺寸（图5-1）侧面及下面边界是用滚轴固定的，而重力作用于材料。于是，通过若干时间步长，直到各网点处的速度或不平衡力趋于零为止，来“步入”（S stepped）这个模型。即直到这些变量衰减到一个可接受的水平。Flac允许两种基本方法来时间步入每个模拟阶段的平衡。用户可以用一个他所要求的具体的时间步长数，

图5-1靠地面的隧洞

内部网格的最大不平衡力的极限，或求解所允许的真实时间（分）。另外，可以允许程序不受用户干扰去求解平衡问题。前一种方法的优点是由求解方法的明显特性所引起的。当初形成的显式程序是允许求解诸如不可能获得的结构崩塌这类高度非线性问题。这时，它便于能考察当崩塌发生时的崩塌过程。为此，在某一时间时（取决于问题的类型）把解的过程循环若干步，建立一个“Save”文件，绘出图形，并考察数值输出。重新储存这个Save文件，并用这种方法，或继续循环以使用户满意（即，或者到平衡或者到某些破坏状态）。

试考虑图5-1示靠地面的隧洞问题。在本例中，一旦开挖外面的网格，已经出现平衡，物体内的重力应力就已形成，且会处于现场应力状态。这时，通过把开挖域的材料模型改为“空（null）”域的办法来进行开挖。计算过程由再步入平衡来进行。如以前一样，用户能直接让程序求解平衡状态，或人工跨步并考察在中间各步时物体的状态。

因此，用flac的问题解答过程包含在图5-2示的逻辑中。这种处理方法是方便的，因为它实际上表现了发生在自然环境中的过程。几个简单问题的模拟会有助于用户熟系这一求解过程。这些问题是：Prandtl楔形体（wedge），深埋圆形隧洞（kirsch解）及重力荷载下的土体边坡稳定。 5.2运行FLAC（running flac）

FLAC可以用交互式（通过键盘进入命令）或文件驱动（成批）式两种方式之一进行运行。首先，用户必须通过软驱动器将代码拷贝到硬盘，并通过给在第一章的以下各指令解压（de-compressit）。然后，通过键入：

c>flac

29

来启动这个程序。计算机将调用flac，并显示以下信息：

FLAC(Version 2.00)

Fast Lagrangian Analysis of Continua （C）Copyright Peter Cundall （1987） Itasca Consulting Group.Inc.（1987）

Strike return to continue

如果存在一个flac.INI文件，当击Return键时，程序将读取并执行各个命令。当flac已接受这些命令时，显示以下提示符： flac> 如果你想交互式地运行某个码，试开始键入前面已介绍，并在以下各节进一步讨论各个命令。当按Return键时，Flac将执行每个命令。如果出现一个错误，给用户的警告将出现于屏幕。

也可以通过一个事先准备好的文件来运行这个代码，该数据文件包含一组与交互式时应进入的相同的一组命令。这个数据文件可以用任何DOS编辑器来建立，并扩展名为.DAT的文件来储存。为了再调用并执行这个数据文件，利用CALL命令: flac> Call file.dat

式中file是你给Data.file取的名字。你会看到，当FLAC每读取一行数据时，在屏幕上数据入口处向上卷动。最后，通过利用任选项的事先储存的状态，可以重新启动这个程序。对于这三种输入方案的典型数据组给于下一节中。

5.2.1网格的形成及材料特性定义（Mesh generation and material properties definition）

必须首先输入形成有限差分网格的各命令。控制这个功能的三个基本命令为：

flac> grid icol，jrow flac> gen〔2〕

flac> initial x or y

GRID命令生成一个宽icol个域，高jrow个域。i-轴平行于问题的x-轴，而j-轴平行于问题的y-轴（图4-1）。

GEN命令用于提供网格坐标，并形成所要求的任何外形的网格。GEN命令的特殊性能包括“分级”（grading）网格，以便表示远的边界。

这个命令也用于建立网格内部（或在网格内部）的各种外形的边界。预先编辑的外形包括圆、弧及直线。直线函数能进一步用于形成任意形状的孔洞。下列是利用GEN命令生成各形状的例题。 例1——方形域

在它的最简单的形式中，GEN命令能为网格提供新的坐标。命令

flac> grid 10，10 flac> m e

将建立一个10域×10个域的方形网格（11个网点×11个网点）并对各个域指定了一个弹性模型。如果打印这阶段的网格各坐标，需发布的命令为： flac> pr x y

你将注意到，x,y都从0.0~10.0运行——即flac规定了一个网点之间相同一个单位间距的方形网格。如果网格的实际坐标为x向0.0～500.0，而y向为0.0～1000.0，那么采用以下GEN命令 flac> gen_0.0,0.0_0.0，1000.0_500.0，1000.0 _500.0，0.0_i=1,11 _j＝1,11 注意，由i=1,11,j=1,11确定的网格区从网格的左下角开始，并以顺时针方式绕各外角围成。位于这些外角点内部的所有网点，将根据各点坐标重新确定其坐标。现在，再打印出这些坐标，以便看出这些坐标，的确已经改变了。注意，正好能给出一个网格区的各个新坐标。这个网格区是由i，j的范围定义的。在上述例题中，如果把范围规定为i=5，11 ; j＝5,11，这只影响网格的右上角。范围

30

必须按顺时针方式规定。

例2——开矩形孔洞

能用GEN命令来形成网格中的畸变。例如，试用以下命令：

flac> new

flac> grid 20,20 flac> m e

flac> gen 0,5 0,20 20,20 5,5 I＝1,11 flac> gen same same 20,0 5,0 I=11,21 flac> plot grid

在本例中，唯一的一个网格还被用每个GEN命令畸变了。第一个GEN命令，由半个网格形成一个畸变的四边形，而第二个GEN命令。“复盖”（wrap）了网格的其余部分构成一个矩形孔洞。注意，逐次的GEN命令是相加的――即只一次改变，直到再用GEN或INI命令来改变，这些网格坐标才保留为新坐标。

例3――按等比级数逐渐变大（小）的网格

为了表示各边界，可以用GEN命令来逐次变化一个网格（to grade a mesh to represent for boundary）例如，在许多情形下，在岩石物质的极深处进行开挖。须确定开挖周围应力及位移的详细信息，这些地方的干扰是大的。但是，在远离开挖处，不必细微。在以下例题中，网格的左下部分（is left tighty discretized），在外面把边界分为x,y向的（the boundaries are graded outward in x，and y,direction）。试发布以下命令：

flac> new

flac> grid 20,20 flac> m e flac>

flac> plot grid

GEN命令将使网格沿x,y向以1.25传于前一个网跨的比率扩展至1000个单位。注意，如果比例介于0～1之间，网格的尺寸将随坐标值的增加而减少。例如，发布以下命令：

Flac> new

Flac> grid 10，10 Flac> m e

Flac> gen –100，0 –100，100 0，100 0，0 rat .80，1.25 Flac> plot grid

你将会看见沿x,y向逐渐变小的网格。

例4――带孔洞的网格

常需要在网格中进行开挖，要构成复杂形状的开挖，是十分麻烦的，特别是用直接移动各单独的网点来构成圆弧时。这里给出了利用GEN命令来建立开挖形状的两个例子。

首先，建立一个圆形开挖。试用以下命令：

Flac> grid 20,20 Flac> m e

Flac> gen circle 10,10 5 Flac> plot grid

这个命令自动构成一个外轮廓在网格内部的圆孔，圆孔中心在（x＝10，y＝10）处，而半径为5.0。注意，网格其余部分仍为方形――即，单元的各角均为900。为了使网格更好的与新开孔一致，键入：

31

Flac> gen adjust Flac> plot grid

逐次的各GEN ADJUST命令将使网格变形（warp），以便日益显著地弄通（to increasingly greater－levels）

当利用Circle，Arc或Line命令，在网格内部构成内部形状时，flac把由有记录的闭合路线形成的网格的不同区域间予以区分。在上例中，Gen Circle命令，在由圆边界构成的网格内部，建立了两个区域：边界内的区域和边界外部的边界。如果你想看见网格的边界在那里，试键入：

Flac> pr mark

那些由flac调整了以便形成边界的各个域，在打印输出中，用一个“M”来表示。 应小心注意――如果由闭合的轮廓线把两个区域隔开，那么才能形成两个区域。换句话说，始末点在网格内部而不构成一个闭合边界的一条线段，只能形成一个区域。对于上例，可以用关键词REGIN去规定二个域x的模型、特性或原始值。例如，为使上例隧洞内部无效，只要健入命令：

Flac> mod null region＝10,10

Flac> plot grid

你会注意到，隧洞中的各个域已被挖去。这时可以用以下命令代替它们。

Flac> m e region＝10,10

注意，可以用任意模型及与模型相协调的特性来代替开挖。更多细节，见Model及Prop命令。

例5――移动网点

可以用Initial命令，把一个点或若干点从当前位置移到一个新位置。以下各命令建立了一个网格，并用Initial命令来使它变形：

Flac> grid 5,5

Flac> model elastic

Flac> gen 0,0 0,10 10,10 10,0（沿x,y向规定网点坐标为0～10） Flac> Ini x＝－12 I＝1 j＝6（把左上角水平移动－2个单元）

Flac> Ini x＝12 I＝6（把各右边界网点各移2个单元，注意，因为没

有给j的范围，取为整个范围）

注意，可以用Initial命令，把任何网点移至任何位置。当然，各单元不能重迭。如果重迭了将给出一个为“BAD.GEOMETRY”的警告信息，直到校正了网格结构中的错误为止，flac不会继速执行的。为了得到合理的解，应把各域的高宽比的实际范围，保持在约1:10或更小。

一旦定义了网格，便必须对物体内部的各个域，规定一个或更多的材料模型及所要求的模型。这可利用3.0节中所叙述的Model及Property命令来达到。 5.2.2 应用边界条件（Apply the Boundary Conditions）

在生成网格以后，便应用边界及/或原始条件。可能的边界条件如表5-1所示。

表5-1边界及初始条件命令小结

Command Effect

Initial

Sxx 初始化域的Sxx应力 Syy 初始化域的Syy应力 Sxy 初始化域的 Sxy应力 Xdis 初始化域的x－位移 Ydis 初始化域的y－位移 Xvel 初始化域的x－速度 Yvel 初始化域的 y－速度

32

Fix

X 固定一个网点处的x－向速度（或位移） Y y X Y x,y Free

X 与Fix命令相反 Y X Y Apply

Pressure 沿整个一行或列各网点作用的压力 Xforce 作用于网点的行或列上的x－力 Yfoce y－力

例如，试考虑受重力荷载的网格的形成。在这个网格内部，开挖了一个靠地表的圆形隧洞。可以键入以下命令序列，供交互式考察GEN命令的效应：

Flac> grid 10,10 Flac> model elastic

Flac> gen circle 5,5 2（在初始化应力之前，形成隧洞的网格） Flac> plot grid

Flac> gen adjust（调整网格） Flac> plot grid

Flac> prop s＝.3e8 b＝1e8 d＝1600 （特征值） Flac> set grav＝9.81 Flac> fix x I＝1

Flac> fix x I＝11（在两侧面及底部让土壤均匀沉降的滚轴边界条件） Flac> fix y j＝1

Flac> * 在物体内形成的应力，这时，开挖隧洞并观察反力 Flac> pr mark（显式在物体内部的各个区域的边界） Flac> mod null region＝5，5（使隧洞内部的区域失效） flac> plot grid

Flac> set 300

Flac> plot str bou . . . 视需要继速模拟。。

上述问题说明了一个初始应力问题――即在开挖以前，存在于物体中的应力。因为在时间步长开始以后，网格不能改变，它必须通过变形，以适应初应力处于平衡前的图形边界。这种典型的问题，与物体受应力以前就存在隧洞的？大不相同。对于无支撑的开挖转形极限应力状态是相同的；不过隧洞处的位移是不同的。这类似于把在一个预应力板中剥开的孔洞与在一个板中先开洞而加应力的情况相比较。在这种情况下，用户必须决定那种类型的问题是适宜。

5.2.3 应力荷载/变化条件（Applying loading/changing condition） flac允许在求解的过程中，改变任意点的模型条件。这些变化可以为以下形式：

材料的开挖；

网点荷载或压力补充或删除； 任何域的材料模型或特性的变化； 固定或解除任意网点的位移或速度。

为了在交互式中完成这些操作，连接在屏幕提示处键新的命令。

33

作为在一个运行过程中，能改变各种条件的例子，试考虑一个随挖随添的简单例子，这里先进行开挖，而后立即回填。试运行以下简化问题：

flac> grid 10,10

flac> m e ﹡ elastic rock mass

flac>prop s=5.7eq b=11.1eq d=2000

flac>fix x y i=1 flac> fix x y y=1

flac> fix x y i=11 固定所有边界 flac> fix x y y=11

flac>int sxx=－30e6 syy=－20e6(场地应力) flac>his nstep=1

flac>hia xdis i=5 y=5(回采工作与壁处的 x位移) flac>salve

flac>﹡ new excavate

flac>mad null i=3.4 y=3.6(开挖非1) flac>plot grid

flac>step zoo(步入平衡) flac>plo his 1(校核平衡)

flac>plot gird str (绘应力，网格)

flac>mod mohr i=3.4 y=3.6(用sangfill特征) flac>mod null i=3.4 y=7(开挖非2)

flac>ini xd=0 yd=0(初始化位移，以致于看出开挖引起的变化) flac>step zoo(再次步入接*衡) flac>plo his 1

flac>plot grid str(绘应力及网格) flac>plot coh(绘不同的材料类型)

在本例中，对于一个精确解来说，边界太近了。但它说明了在运行过程中，flac是能改变模型条件的。应注意——如果用户用其他材料（在本例中，为sandfill）来代替空域后的网格，各取消的单元将被代替。这说明，即使各个域给为空材料，且不绘或打印出来，但它们仍然存在并能在任何时候被改变。材料特征的一幅轮廓图也能通过在所要求的特征后面，键入plot来绘出。

这里利用history plot ,是使用户能判别什么时候问题足以接*衡。这时，开挖壁处网格的x向位移是用作移向平衡的指示器。当逐次时间步长间的值（惟一，应力等）有微小变化时，则模型大体上处于平衡了。

5.2.4 数据的打印及绘图（Printing and ploting of data）

如前面所讨论，直接利用print命令，可以得到一个适当的数据的屏幕打印：

flac> print keyword

打印出的数值矩阵类似于网格本身的逻辑格式。各keyword是显然的，但以下两点例外：

1. 如前面讨论，flac把用户定义的各矩形域在分为名为a,b.c及d的四个三

角形子域。用户可以利用以下的命令请求打印每个子地域的应力 Flac>print asxx,....etc.

34

2. flac保存了一份每个域的屈服史记录。它常常用于了解那些域现在正在

屈服，那些过去已经屈服但现在为弹性，以及那些域现在及过去均为弹性。可以用print sta 来提供每个域的现行弹性状态。打印在矩阵域中的各个数值具有如下意义：

0 弹性 1屈服处

2弹性但以前屈服

3已超过单轴拉伸的拉力cut-off 4屈服，也超过拉力cut-off

5已超过一般拉伸的拉力 cut-off

6现在屈服的多处存在的节理（ubiquitous joints） 7过去屈服但现在为弹性状态的多处存在的节理

用户可以用SET LOG ON 及 SET LOG OFF装置，来得到打印于屏幕的任何文件的打印拷贝。在一个flac数据组中（或者来自批文件，或者一个交互式文件），这命令

flac>set log=on 将打开一个硬盘文件（FLAC.LOG），该文件将记录后续的所以命令及 flac的各个响应。要关闭该记录装置，应键入： flac>set log=off 在退出flac后那么可以关闭这个文件，或者把它拷贝给打印机，以作为永久记录，直接利用DOS COPY或PRINT命令。

网格的几何条件及各种输出变量的图形，可以在PC显示器或一台Hewlett-Packard或兼容绘图机上完成。适用于用户的现有图形类型，编列如下：

1. GRID绘制有限差分网格，当时间步长继续时，GRID命令将绘出

有限差分网格变形后的形状（在大应变模型时）。这对于大应变逻辑是特别重要的。因为，破坏方式常常从边了形的图形中变得明显（since,often,failure modes become evident from the defarmed plots） 2. Displacement及Velocipy给出各个网点处的位移及速度矢量图。 3. Boundary绘制网格边界及开挖。它常常是便于把边界图重叠在其

他不同的图上，诸如矢量及等值线图。 4. sxx,syy及sxy绘制xx, xyy或xy应力。 5. Xdisp及YDISP绘制x位移及x位移等值线。 6. XVEL及 YVEL绘制 x速度及y速度。 7. STRESS绘制主应力矢量。 8. PP绘制孔隙压力等值线 9. HISTORY绘制由HIS命令所选择的变量的时间经历（time history）. 10. STATE绘制线的屈服状态。 11. BEAM绘制结构单元。

12. CABLE绘象铺（cable bolt）的几何条件。

必须首先用SET命令初始化图形屏幕。2.00版要求一个IBM（或兼容机）的CGA或EGA图形卡（graphics card）或hercules单色卡。SET命令定义为： 1.屏幕类型﹡——低分辨度彩色（CGA），高分辨度彩色（EGA），或高分辨度单色的。默值为低分辨度单色屏幕。屏幕的选择点，以下命令之一设置：SET CGA

35

（低分辨度彩色图形装置时），SET MONO（单色时），SET EGA（增强图形转换装置）

﹡ 单色图要求彩色图形转换装置，但是，将以任选的高分辨度单色图形方式，绘在彩色屏幕上。一个的flac版本(a separate version of flac) 是适用于 Hercuies Graphics board绘图的情况。

2.背景颜色——用SET back=iback来设置背景颜色，其中iback的范围为0~7：

0=黑色（blake） 1=兰色(blue) 2=绿色(dreen) 3=兰绿色(cyan) 4=红色(red) 5=深红色(magenta) 6=黄褐色(yellow-brown) 7=白色(white) 默认值为单色方式的绿色，因为背景为线条颜色。例如，如果你想在单色方式的黑色屏幕上的红色图形，那么Set back=4。 CGA及EGA转换器（board）要求一个黑色背景（或back=0）。如果你有这些任选项之一，最好用set back =0建立一个flac.ini文件

3.调色板的选择——利用SET-PAL命令时，有两种可能的调色板（palette）选择设置(ipal=0或1)。每种调色板有三种可能的前景颜色。ipal=0时，CEA的颜色选择为褐色、红色及黄褐色。ipal=1时，CGA的颜色选择为兰褐色、深红色及白色。EGA时，14种颜色选择是可能的。关于颜色选择的细节，参见PLOT（第四节）。通过PLOT命令的颜色keyword，来定义前景颜色。

设置一个彩色，有黑色背景底分辨度屏幕，而调色板为非0的示例性命令为： flac> set pol=0 back=0 CGA

一旦图形设备已经说明，便可以设置绘图窗口。窗口应该是方形的，以便清除图形的畸变。命令形式为：

flac>wind xlow,xhigh ,ylow ,yhigh

如果没有给出WINDOW命令，将自动选择一个尺度。注意，通过改变窗口尺寸，可以把图形放大以便细致的考察网格。窗口应该始终保持为方形，否则，图形会畸变。即使为方形，图形也可能稍微有点畸变。这是由于每个屏幕的调节器的变化所引起的。为了清除任何畸变，试调节显示器背面的调节器钮，或利用SET ASPECT。

现在可以直接键入后给各重叠图形的keyword的PLOT命令，来把这些图形绘到屏幕上。一个示例性图形命令为：

flac>plot Sxx disp boun

它将重迭Sxx应力等值线，位移矢量及网格和开挖边界，在图形的左上角，生成一个描述图形的标题。图形更一步增强（enhancements to the plot）可以用“sub-keywords”或开关来获得。这些开关可以用于：

（1） 定义最大矢量长度（max==value）； （2） 定义等值线的间隔（Intenval==value）； （3） 定义线条颜色；

（4） 取消标题（remave headings）。

例如，上述命令能给成：

flac>plot Sxx Int=1e4 Cyan disp max=0.001 mag-boun 这里，Sxx应力等值线绘成为蓝绿色（取调色板非1），具有1E04的等值线间隔；深红色的位移量是以最大矢量长度为0.001来度量的； 网格及开挖的边界均以

36

背景颜色来表示。

在绘图以前，可以用SCLN命令把各等值线标号。

一旦一幅图已在屏幕上构成，并感到满意时，可以在一个Hewlett-Packard笔式绘图机上生成一幅图形的拷贝。通过在PLOT的keywords前面，键入命令PEN，得到一幅这样的图形。例如命令：

flac>plot pen Sxx disp boun 将把这幅图形的内存信息转移于笔式绘图机。

也可以把这个图形传给图形点矩打印机（graphics dot matrix printer）。在把flac装入之前，通过发布DOS命令GRAPHICS，既可做到这一点。当一幅图绘在屏幕上时，可以用print screen键把它转移到打印机上。一般说来，图形在打印机上的高宽比是相同的。为了建立一幅在打印机上不致畸变，而在屏幕上的“换算”（corrected）畸变图形，用一个非方形窗口作某些试验是必要的。一些商业软件包（commercial package）是适用建立单色或彩色CGA或EGA转换器的点矩图形的。

5.3 特殊问题的考虑（special considerations）

flac有一些注释特性(feature of note)，在此予以回顾。 5.3.1 大应变（large strain）

当物体遭受应变时，正常的假设是：位移的大小小得足以使物体的坐标保持不被改变。这是flac中默认的小应变逻辑情况。然而，如果应变都很大，当应变发生时，各网点的坐标必须修正，而各域的应力也应校正。当处理塑性本构模型及可能的结果崩塌模型时，大应变逻辑是不可少的。用命令：

flac>set large

来调用flac中的大应变任选项，就是这种情形。如果你想返回到小应变模型，其命令：

flac>set small

那么，在进一步的计算中，将对网点坐标不作修正。 5.3.2 平面应力 （plane stress）

flac的默认解是平面应变。不过，平面应力假设可用以下命令来设置： flac>set p-stress

在运行中间从平面应变转向平面应力是不合理的，因此，这个命令应仅用于网格初始下时。

5.3.3 重力 （gravity）

通过设置重力，可以把重力力包含在物体内： flac>set grav 9.8 0

值9.8是重力引起的加速度，而0指的是重力竖向作用。（默认值为0，且是不必要的）。重力的符号约定如下；

+== 向下 －== 向上 5.3.4 图形形状 （Aspect）

一个方形网格应以一个方形呈现在显示屏幕上。如果办不到，有两种改正方法。

第一个方法是用竖向尺寸控制或你的PC监视口的调节器。如果这个方法不能解决问题，那么在屏幕上绘一个方形网格（ 即，grid 2,2）并用直尺手工测量其水平及竖向尺寸。跨高比为水平/竖向用以下命令设置flac的跨高比（aspect ratio ）

37

flac>set asp=a 5.3.5 存入（logging）

利用“logging on”和“logging off”可以得到？的flac对话期间的一个硬盘文件。其命令为：

flac>set log =on

及 flac>set log =off

它们分别表示打开存入（log on）及关闭存入文件（log off）。在这两个命令之间打入屏幕的任何信息，将被存入名为flac.log硬盘文件中。别忘了，当你开始一个新的flac运行时，一个存在的flac.log 文件已被覆盖。如果你想保留这个log文件，在你开始新的 flac 运行以前，将它更名. 5.4 错误处理（Error handling）

flac具有一个错误矫正逻辑，它能识别很多常见错误，并在交互对话方式时，为用户给出一个不中止程序的信息。控制系统返回到提示程序。各种错误的信息一览表给在附录C中。

5.5保存/复原运行(Saving/Restoring)

在模拟过程的任意处，flac 允许用存入一个状态，然后又复原它。通过键入以下命令，来建立一个保存文件

flac>Save Z:file.Sav

其中“file”是你想给的保存文件名；“Z”是驱动器说明。为了复原一个已存在的文件，直接键入命令：

flac>rest Z:File.Save

其中“file”是以前存入状态的名称。自由的把存入文件用于参数研究，绘图及打印总是最好的主意。

5.6 建议及忠告 （Suggestions and Advice）

1. flac采用常应变。如果应变梯度是强烈的，你需要很多域来表示非均匀分布。为了校核，试用更多的域来运行相同的问题。对于塑流，采用常应变是因为使用许多低价单元比几个高价单元的好。

2. 试使域的划分尽可能的均布。避免细长域或很畸变的域。 3. 对于一个新问题，总是以少量的几个域先作试运行，以便对于各种响应及可能的困难，获得一个快速感受。当你了解了试运行结果时，才增加域的数目，以便得到更准确的结果。

4. flac将花一个很长的时间收敛，如果（a）有巨大的刚度变化；（b）有巨大的域尺寸对比。

5. 一个很刚性的荷载板，常常可以用一系列给出常速度的固定网点来代替（重新调FIX命令表固定速度，而不是位移）。 6. 为了确定一个崩塌荷载，经常最好是在“应变控制”条件下，而不是在“应力控制”条件下来进行——即，应用一个常速度并测量各反力，而不是应用力并量测各位移。当外加荷载趋于崩塌荷载时，一个崩塌的系统变成难以控制（这对一个真实体系以及一个模型体系是正确的）。

7. 只要有可能，试利用对称条件，以节省计算机内存及运行时间，例如，如果一个体系是对称于一竖轴的，你可以把对称线表示 为X向具有固定网点的竖向边界（但在Y向自由的）。

8. 频繁的利用save文件。例如当进行参数研究时，保存中间状态。如果由于某种原因运行终止，你将有这些中间状态。

38

9. 处理一个flac模型正如一个物理模型一样。在一个flac运行中，试图重现自然中实际上会发生的各个阶段。别忘记，对于一个非弹性体系，没有唯一的平衡状态。可能有许多满足平衡的可能阶段；你所得到的阶段取决于历史。

10. flac怎样表现一个体系的性能。常用一些简单试验去校核你正在作你想做的事。例如，如果加载条件和几何条件是一致的，试校核反应是对称的，或者，在使荷载变化或其他变化`，试先执行少数几个步长（如5）；以便检验其初试反应具有正确的符号及位于正确的位置。

Might alse do back-of-the-envelope estimates of expected order of magnitude stress or displacements and compare them to flac out put 11. 如果把一个强烈的冲击（violent shock）作用于一个体系，你会获得一个强反应。如果你把非物理的合理的东西加于这个体系，你必须预料一些奇怪的结果。

12. 在继续进行模拟之前，试以鉴定的眼光考察输出。例如，如果除在一个角域处大速度外，一切正常，则直到你了解原因为止，便不要继续下去了。这时，你可以让一个“固定”的网点自由。 13. FLAC不直接给出一个“安全系数”。如果你需要一个安全系数，可以通过取真实值与引起破坏的值的比的办法，来对你认为重要的安全参数来定义安全系数。例如：

FW==

引起崩塌的地下水位

实际地下水位tan(实际摩擦角）

tan（引起破坏的摩擦角）F=

Fl=

破坏荷载

设计荷载注意，总是用较小的值去除较大的值（假设在实际条件下，体系不会破坏） 14．如果你想flac通宵运行，但次日上午仍有计算机可以利用，试直接设置一个适合大步长数。当你需要计算时，可以键入任何键，以便返回到命令(command level)。万一电源中断，在包含一个调用(called)频繁save file的控制下，作长运行是比较保险的。

15． 试用历史图形(history plots)去确定什么时候一个模型处于平衡。

6．0 FLAC中的结构模型（Structural Modeling in Flac）

任意几何条件及特性结构，以及它们同土壤或岩石介质的相互作用，是可以用FLAC来模拟的。提供了两类结构单元：

（1） 梁单元（它可以同拧一个梁单元和/或土壤或岩石介质连接，

以模拟各种问题，比如，挡土墙，缺口开挖的支撑（Support struts in an open-cut excavation），以及隧洞的混凝土或喷射混凝土（Shotcrete）衬砌，或诸如荃脚（footing）这样的地面结构。）

（2） 索或岩栓（Rackbolt）单元（能用它来表示全灌浆

39

（fully-grouted）或点锚（point-anchor）岩栓，索栓（Cavlebolt），或牵索栓（Tie-back anchor）

所有情形，定义这些结构符必须的各命令是相当简单的，只不过调用（invoke）一个十分强有力且适应性强的结构逻辑（Structural lagic）。因为这个结构逻辑是用如同其余代码（与矩阵结构法不同）一样的相同有限差分逻辑，其结构可以遭受大位移。

6．1 命令结构 （Cmmand Structure）

为了模拟结构与土壤及岩石介质的相互作用，用户必须定义支撑物的几何条件（除了物体外），支撑物的特性，以及各结构单元间和它们与支撑介质间的连动装置的类型。

为了完成这些功能，提供三个基本命令。命令 STRuct

BeamCable Begin

Gridi,Mode,xbjnybend

Gridi,Mode,xejn ye

定义了单元类型，它的始点联动装置，结构的节点或X，Y坐标，参数终点，单元分为若干段的数目，以及单元的特殊类型。一根索栓可以预张到某一荷载T，用keywordTENSION。 命令

STRuct Prop=k keywords......

定义了梁或索栓的特征类型。可以利用以下特征keyword: Aera 截面积

E 扬氏模量

I 梁的惯性矩（二次面矩） Kbond 灌浆的粘结刚度 Sbond 灌浆的粘结强度 Yield 索的屈服强度

最后，命令 STRuct Node=n keyword

定义了各结构单元节点处的联动装置。用户选择各节点处的固定或自由位移或转动，作用于X，Y向的俩或弯矩，或压个节点处规定一个铰链。默认值时，把各节点考虑为自由的。以下keyword可用于定义联动装置：

Fix <2> 节点n处，固定x或y位移或转动 FRee <2> 与FIX相反

Load fx fy m 用户提供的取决于方向的力或弯矩 Pin 在节点n处建立铰链

x，y 用户提供的节点n的x及y坐标 用这三个命令，能规定任意几何条件、边界条件及各个结构。

6.2 定义结构单元的几何条件及其直承介质的联动装置（Defining the Structural Geometry and Its Linkage to the Support Medium）

对于每给梁或索必须定义其几何条件及联动装置。可以把一根单梁或索分为若干段，已改进计算精度，采用了上述基本命令/位于括号内的keyword都是任选项，虽然必须规定每一个keyword。

40

各梁单元都是用于任何弯曲刚度重要的结构构件。应用梁单元的例子包括： (1) 诸如基础或荃脚这类地面结构（使用的逻辑是通用的，因此可以把地面

结构表示梁系（assemblage of beams） (2) 各种支持结构（retaining structures）

(3) 诸如混凝土或喷射混凝土衬砌，或钢木支撑这一类内部连续隧洞支撑。 应用索单元的例子包括：

（1）支撑结构的牵索锚（点或全灌浆锚固（anchorage）） （2）索锚（点或全灌浆锚固） （3）岩栓（rackbolt）（点或全灌浆锚固）

用关键词Begging或End 来定义单元的始末端联动装置。对于begging 及 End的几个适合选项为：（1）grid=i，j (2)node=n 及(3) Xb,Yb 或Xe,Ye

关键词grid=i，j表示梁或索的始（或末）端，与岩石或土壤介质的网点i，j相连接。令node=n把梁或索的始（或末）端与结构的另一个节点（即，另一个梁或索单元）相连接。通过一对表示梁或索始（或末）端坐标的x,y，用户可以把一个结构单元安置在网格内部或外部的任意位置。

每个不同类型的老梁或索，用命令STRuct Prop=n来规定其特征。例如，如果某根梁 一种结构截面（one structual cross-section）取n=1，对于W6×25梁（SI单位），其命令为：

flac>strac prop=1 E=200e9 I=2.3e5 A=4.8e-3 应指出，所有量的必须用等价的单位系绘出（见表4-1），在计算荷载时，代码不考虑结构自重，因此，如果结构自重在问题是重要性，应把等价于自重的竖向力作用于结构。

补充命令struct node=n keyword提供了描述结构单元间连接的任选项。任选项包括：

（1） 自由或固定x,y位移或转动。 （2） 铰节点

（3） 外荷载或弯矩 （4） 节点坐标。

这些联动装置任选项由以下适宜的选项keyword给出。 Fix x,y,r 固定分给节点处的x,y位移和/或转动。 Free x,y,r 释放

Pin 在一个节点处提供一个铰链（即，释放弯矩）

Load fx,fy,m 在一个节点处作用一个x,y向荷载，或弯矩。荷载沿卡笛尔

坐标为正。反时针向弯矩为正。

x,y 规定一个新节点号的坐标。注意，如果你用一个任选项而不

是node=n命令已建立了各节点，FLAC将自动把节点号分配给它？建立的那些节点。为了了解现行各节点号，试发布命令PRINT STRUCT

6.3 实例应用 （Exmaple Applications）

为了具体说明FLAC中各结构单元的使用，给出几个简单例子。第2.0节中的例题具体说明了一个崩塌的壕沟。这里，我们用两根支撑开挖壁的结构表示这个壕沟。以各命令序列是对把这些支撑土即放置在开挖上的情形而言的。

flac>gird 5,5

41

flac>m mohr

flac>prop s=0.3e8 b=1e8 d=1600 fric=20 coh=0 flac> fix y j=1 flac> fix x i=1 flac> fix x i=6 flac>set large flac>hist nstep=1

flac>hist xdis i=3 j=6 flac>set grav=9.81 flac>solve

flac>*excavate trench and install braces flac>model null i=3 j=3,5 flac>*properties for w6×25 beam in SI units flac>struc prop=1 E=200e9 I=2.3e-5 area=4.8e-3 flac>struc beam beg gr=3,6 end gr=4,6 s=3 pr=1 flac>struc beam beg gr=3,4 end gr=4,4 s=3 pr=1 flac>step 50

flac>plot gird beam 反复考察平衡移动，绘结构 flac>plot rf bound 图，反力图，位移图 flac>plot bound dis beam

图6-1说明了两根支撑的效应（与图2-2和2-3相比）。壕沟的崩塌仍会发生（对该模型），但破坏范围减少了。补充的支撑和/或板桩（由竖梁表示的）可能试图稳定壕沟。

第二个例子说明了地面结构基础荷载。这里，一个简单的交叉撑平台（）建造在放于土壤介质上的混凝土板上。结构受有作用于各支撑柱上的竖向集中荷载。本项目是考察结构中的荷载及弯矩，以及土壤介质中的应力及位移。

以下命令结构用于准备和运行该问题：

flac>*a Simple cross-braced stracture on a sail flac>grid 10,10 *the sail mass flac>m e

flac>prop s=0.3e8 b=1e8 d=1600 flac>fix x i=1 flac>fix x i=11 flac>fix y j=1 flac>srt grav=9.81 flac>salve

flac>*let sail equilibrate under gravity flac>*build structure flac>*concrete slab

flac>strac prop=1 E=17.58e9 I=0.0104 a=0.5 flac>strac prop=2 E=200e9 I=2.3e-5 a=4.8e-3

flac>struc beam beg gr 5,11 end gr 7,11 s=2 pr=2 flac>struc beam beg node=1 end 4,13 s=2 pr=2 flac>struc beam beg 4,13 end 6,13 s=2 pr=2

42

flac>struc beam beg 6,13 end 4,10 s=2 pr=2 flac>struc node =8 5.0 ,11.5

flac>struc beam beg node=8 end node=1 s=1 pr=2 flac>struc beam beg node=1 end node=4 s=1 pr=2 flac>struc beam beg node=1 end node=6 s=1 pr=2 flac>struc beam beg node=1 end node=2 s=1 pr=2 flac>struc beam node=1 fix r

flac>struc beam node=2 fix r 4 flac>struc beam node=4 load 0 -1e6 0——似应为 node=3

flac>struc beam node=6 load 0 -1e6 0 2(译者著)

flac>plot beam *check strcture flac>pr struc *check linkage flac>step 100

flac>pr struc

flac>plot beam boun ydis 反复考察朝平衡的移动 flac>plot rf beam boun

本问题的原始图形如图6-2示。荷载对带结构物的土壤的影响说明于图6-3。产生在结构物中的力及力矩，是各外荷载和与弹性网格相互的结果，并可用PRINT STUCT来。

例三：

第三个例题中，考察产生在喷射混凝土及混凝土中的荷载，该衬砌位于双向应力场中，是作用？？。因此用梁单元表示与岩石介质连接触的衬砌。各界面单元也可用于模拟衬砌与岩石间的划移效应（见附录D）。

对于这种情形，网格及衬砌的形成如下： flac>grid 15,15

flac>m mahr (mohr-coulomb model) flac>gen 0,0 0,30 30,30 30,0

flac>gen 14,14 14,60 60,60 60,14 rat 1.2 1.2 i=8,16 j=8,16 flac>gen 14,-30 14,14 60,14 60,-30 rat 1.2 0.833 i=8,16 j=1,8 flac>gen -30,30 -30,14 14,14 14,-30 rat 0.833 0.833 i=1,8 j=1,8 flac>gen -30,14 -30,60 14,60 14,14 rat0.833 1.2 i=1,8 j=8,16 flac>gen circle 14,14 4 flac> gen adjust

flac>*rock properties

flac>prop s=5.75e9 b=6.6e9 d=2000 coh=1e7 fr=35 flac>*boundary conditions flac>fix x y i=1 flac>fix x y i=16 flac>fix x y j=1 flac>fix x y j=16

flac>ini sxx=-60e6 syy=-30e6(初应力)

flac>solve（处于平衡，坑道隧洞及设置支撑；然后让它达到平衡） flac>mod null region=8.8

43

flac>pr mark

flac>struc beam beg gr 7 9 end gr 8 9 flac> struc beam beg gr 8 9 end gr 9 9 flac> struc beam beg gr 9 9 end gr 9 8 flac> struc beam beg gr 9 7 end gr 9 7 flac> struc beam beg gr 8 7 end gr 7 7 flac> struc beam beg gr 7 7 end gr 7 8 flac> struc beam beg gr 7 8 end gr 7 9 flac> plot gird beam flac> step 100

flac>plo bou disp red beam

flac>pr struc 反复考察朝平衡的移动 flac>pr state

各衬砌单元中的荷载及弯矩（从PRINT STRUCT）列表6-1。衬砌通风？周围的最后应力状态如图6-4所示

表6-1 隧洞衬砌例题PRINT STRUCT结果 structural node data.... ID X Y X-disp Y-disp ang-dis xfix yfix rfix 8 1.000E+01 1.400E+01 1.955E-02 9.416E-05 -7.235E-05 no no no 7 1.117E+01 1.117E+01 1.279E-02 3.815E-03 -1.863E-03 no no no 6 1.400E+01 1.00E+01 3.403E-04 5.737E-03 -1.306E-04 no no no 5 1.6E+01 1.091E+01 -1.188E-02 3.793E-03 2.095E-03 no no no 4 1.800E+01 1.400E+01 -2.128E-02 2.984E-0.4 2.018E-05 no no no 3 1.683E+01 1.683E+01 -1.334E-02 -2.5E-03 -2.178E-03 no no no 2 1.400E+01 1.800E+01 5.852E-04 -6.271E-03 1.286E-04 no no no 1 1.091E+01 1.6E+01 1.434E-02 -3.678E-03 1 .809E-03 no no no

ID X-load Y-load moment i j pin i ang 8 0.000E-01 0.000E-01 0.000E-01 7 8 no 7 0.000E-01 0.000E-01 0.000E-01 7 7 no 6 0.000E-01 0.000E-01 0.000E-01 8 7 no 5 0.000E-01 0.000E-01 0.000E-01 9 7 no 4 0.000E-01 0.000E-01 0.000E-01 9 8 no 3 0.000E-01 0.000E-01 0.000E-01 9 9 no 2 0.000E-01 0.000E-01 0.000E-01 8 9 no 1 0.000E-01 0.000E-01 0.000E-01 7 9 no Structural element data... ID node-1 node-2 prop F-shear F-axial mom-1 mom-2 8 8 1 1 beam -3.087E+03 4.079E+06 -6.855E+03 -1.457E+03 7 7 8 1 beam 3.002E+03 4.073E+06 2.336E+03 6.855E+03 6 6 7 1 beam –9.753E+01 8.277E+06 2.037E+03 -2.336E+03 5 5 6 1 beam 6.595E+02 9.351E+06 4.351E+03 -2.037E+03 4 4 5 1 beam -3.912E+03 4.343E+06 -9.037E+03 -6.351E+03 3 3 4 1 beam 4.091E+03 4.520E+06 3.437E+03 9.037E+03 2 2 3 1 beam -3.761E+02 9.382E+06 2.336E+03 -3.487E+03

44

1 1 2 1 beam -2.568E+02 8.185E+06 1.457E+03 -2.336E+03 Structural properties...

Prop NO. E I Area

1 1.380E+10 2.800E-04 1.500E-01

最后的例题是考察受跨中荷载的轻质钢筋混凝土梁（lightly reinforced cohcrete beam ）的性能。这里，有索单元表示钢筋。该问题的输入命令为：

flac>*a simple reinforced bean flac> gird 10,2 flac>m e

flac>prop s=9e9 b=10e9 d=2400 flac>set large p-stress flac>*boundary conditions flac>fix y j=1 i=1 flac>fix y j=1 i=11

flac>apply yforce -1e7 j=3 i=6

flac>struc cable beg -1 -1 end 9.9 -1 s=12 prop=1

flac>struc prop=1 yi=1e6 kb=1e9 sb=1e7 e=200e9 a=2e-4 flac>his ydisp i=6 j=3

flac>salve f=1e5 1000 t=15 flac>save beam.sav flac>ret

这时，跨中竖向位移约为12mm。

7.0 例题 （Example Problems）

下面提供了一系列岩土力学领域的问题，已作为说明FLAC码的各种特征。可能的话，把这些结果与解析相比较，已验证FLAC码的效果。

7.1 例1：无摩擦粘土上的毛面荃脚（Raugh footing an a friction cohesive

sail）

能为数值模型进行颇为严格检验的一个标准土力学问题，是预测支承在粘土上的一个荃脚的崩塌荷载（SLon及Randolph）。从“Prandtl’s Wedge”解得到的土的承载能力为：

q=(2+π)c=5.14c

式中q为破坏荷载，而C为材料的粘结力。可知道发生破坏的模型。

一个颇粗略的网格被用于该问题的FLAC模拟最初规定为10×10方形网格，但把右边及底部延伸可一点，以便提供合理的逼近于一个半无限边界。在网格的左边利用了对称条件，而底部及右边边沿，沿x及y方向均是固定的，可以用三种可能的方法，来模拟荃脚荷载：

（1） 把一个压力边界条件应用在跨越所要的各单元的负方向。

（2） 把等价于的荃脚应力的网点力，作用在跨越荃脚宽度的各表面网点的

负y方向。

（3） 把一个常速度边界，应用于跨越荃脚宽度的负y方向。

因为我们的目标是把预测解与精确解相比较，选择了第三种替换方法，通过监视破坏边界下方格网中的力的变化，可以计算出承载能力。这就看出了对

45

外边界条件进行到发生崩塌为止的试验的必要性。

把具有以下材料特性的mahr-coulomb本构定律用于本模拟： 密度 1000kg/m3 粘结力 1.0×105pa

剪切模量 0.3×108pa 摩擦角 0 体积模量 1.0×188pa 此外，为了与精确协调一致，采用了小应变逻辑。 用以下命令序列给出输入：

set log=on tit

prandtl’s wedge gr 10,10

*mohr-coulomb model m m

prop s=0.3e8 b=1e8 fri=0 coh=1e5 *Extend grid boundaries a bit

gen 0,0 0,15 20,15 20,0 rat 1.2 0.8 (第n域长：0.8×1.2) *apply the boudary conditions fix x i=1 fix x y j=1 fix x y i=11

* rigid footing moving at constant velocity fix x y i=1,4 j=11 int yv=-1e-3 i=1,4 j=11

* displacement history under footing his nste=1

his ydis i=6 j=11

*solve with new limits solve f=2e2 t=40 s=1000 set log =off save b:foot.sav

应注意，表示荃脚的速度条件，是 利用fix x y (i=1,j=11)命令的后面，紧跟初始化这些网点，y速度为-1×103m/步长的方法得到的。对solve命令的为，甚至了不平衡力的极限，时间极限为40分，或步长极限1000时间。保存了荃脚下面的y位移时间史，并将在行结束时绘制成图形。

1000步以后，运行停止，发表命令： flac>print state

将给出以下输出：

State (译注：1—屈服，0，2—弹性)

J i 1 2 3 4 5 6 7 8 9 10 10 1.000 2.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 9 2.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 8 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000

46

7 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 6 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 5 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.000 4 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.000 3 2.000 2.000 2.000 2.000 2.000 2.000 1.000 2.000 0.000 0.000 2 2.000 2.000 2.000 2.000 2.000 2.000 2.000 0.000 0.000 0.000 1 2.000 2.000 2.000 2.000 2.000 2.000 2.000 0.000 0.000 0.000 数值1表示该域的材料现为屈服。数值2和0表示材料现为弹性，通过绘位移矢量及速度等值线图，可以得到一幅更好的已屈服材的几何图形。首先，设置屏幕为彩色方式。如果你仅有单色屏幕，为了建立着幅图应不理睬以下两个命令： flac>set CGA(或Set EGA) flac>set pal=0

而键入

flac>plot xv z int=1e-4 yel dis red bo

开关字Z消除了零等值线区间，而区分类int=1e-4把该等值线区间设置为

1×104m/sec。

为了获得一幅笔式图（假使你已经恰当的连接了笔式绘图机——见POLT命令），应在上一个命令行的POLT后面，了解插入PEN字。

崩塌时，在荃脚下面的y向力。由以下命令给出： flac>print yreaction 也可以用以下命令绘出图形： flac>plot rf 屈服的平均应力，可取为y力除以截面积，与精确解5.14×105pa相比较，共值为q=5.17×105pa。因此，即使相当粗略的离散，取误差离真实解小于1%。 时间关系如图随着时间而增加的位移，表示模型已没有达到平衡，而是现行崩塌。

有关这个问题输入方面的几点最后说明是必要的。Solve命令的解的极限，已对默认极限作重新定义。最后，应注意，这个数据文件能否建立在硬盘上，先在运用程序时调用，而不是以交互式方法使用。利用编辑器把相同的数据组键入一个文件（如Foot.DAT），并通过以下命令来调用： flac>Call foot.DAT 当flac调用它时，用户会看到依次执行的每一行。如果你希望将出现在哦上的信息，记录在一个硬或软盘中，直接利用命令 flac>Set log=on

所有打印到屏幕上的数据，将存入硬盘文件flac.log中。

最后，命令save foot.sav将建立一个再启动现行模型状态的文件。为了再启动一个保存的状态，在默认驱动器上直接利用Resto命令(如，flac>Rest foot.sav)。

7.2 例2: 粘性摩擦土的边坡稳定（Slope Stability for a Cohesive and

frictional Soil）

在工程土力学中遇到的另一个常见问题，上有摩擦土体的边坡稳定。在本例

47

中，我们运行两个实例研究。第一，用陡于砂的休止角（angle of repose）的初试边坡，当模拟无粘接力的砂。当然，这个边坡应崩塌，而我们希望给出当它发生崩塌时的进程。第二，我们给这个材料附加一个晓得粘接力，为了确定它是否稳定再进行考虑。

我们将利用一个简单的例子，并认为土为均匀的，这里，可以利用解析及图解法，来确定边坡的稳定及安全系数。flac 的本域在于它能考察更复杂的边坡几何条件，比如，可以存在几层不同材料的土限。通过对不同的域规定不同的材料模型，考察这类问题并不费力。

本例说明了一些网格生成功能，以及域的“开挖”方法的使用。以下命令序列是用于初始化边坡及建立一个重新启动文件。（ 边界条件及/或材性可能不同于原有的）：

tit

*Granular Slape under gravitational load sr 20,10

*Mohr-Coulomb Model m m

*Warp grid to form a slope—first use line command to form slope face gen line 5,3 9,10

*Mark area of slope in front of face so that the regians are formed mark i=1,6 j=4

*null rigion in front of slope mad null reg=1,10

*sail properties—note large cohesion to force initial elastic *behavior for determining inmtial state. This will prevent *Slope failure when initializing the gravity stresses prop s=0.3e8 b=1e8 d=1500 fri=20 coh=1e10 *displacement history of slope his nste=1

his ydis i=10 j=10

*displacement boundary cindition x fix x i=1 fix x i=21 fix x y j=1 *apply gravity Set grav=9.81

*sovle for initial gravity stresses Slove

*reset displacement companmets to zero ini xdis=0 ydis=0 *Save initial state Save hill1.sav

*set cohesion to 0 Prop coh=0

48

*Use large strain logic Set large Step zoo

Save hill2.sav Step zoo

Save hill3.sav Return

设置了20×10个域的初始方形网格。注意，除非用户另外定义，默认值指定该初始的尺寸为x向20个单位，y向为10个单位，规定各域为一个Mohr-Coulomb本构模型（因为没有给定范围），并是有以下特征：

密度density 1500kg/m3 剪切模量shear modulus 3×108pa 体积模量 bulk modulus 1×108pa 摩擦角friction angle 20° 粘结力 cohesion 1×1010pa

你会注意到，最初我们规定了一个很大的粘结力。通过复查显示模型工作方式，能看出这样做的理由。这时，我们必须形成一个初始网格，并把重力作用于各网点，且让重力应力达到平衡。因为这个模型事实是动力的，我们通过时间来步入，使材料团结（consolidating the material），直到各网点的速度很少为止。对于大多数问题，我们希望尽可能快地发生这个团结过程，通过在团结过程中，要求具有弹性的材料，就可以做到这一点。一旦各应力已平衡，我们就可以规定选用任何材性、开挖、施加的荷载等，并继续其模拟过程。在这里具体说明的例题中，最初，我们指定了一个塑性本构模型，但设置了很高的粘结力，以便迫使材料具有弹性行为。一旦这些应力已达到平衡，我们可以重新定义要求的粘结力值。这就免得去重新定义平衡后的本构模型及性质，以节省一些输入步骤。如果先把一个初始的弹性模型用于网格的团结，而后改变为Mohr-Coulomb模型、并定义材料的摩擦及粘结力，也能起到相同的作用。用作初始化一个重力荷载网格的上述两个方法，都是可疑的，由用户自己选用。不过，重要是用户用遵循图所说明的一般求解过程。

我们已经定义了基本网格，现在必须将它变形为边坡形状，且让土壤在边坡的下面。这可用QEN命令办到。

首先，我们定义一根我们将用于定义可能发生的边坡边界的内部直线。其命令为： flac> gen line 5,3 9,10

它构成一根始点（x,y）=(5,3)而终点为(x,y)=(9,10)的直线。你可以绘出网格以验证所发生的那根内部直线。应注意，如果生成实线，该线的两个端点应与已有的网点重合如果它们不在正确位置，你可以在利用GEN命令以前，（用INI命令）来改变个网点的位置。其次，我们连接使没有切断边坡脚（slope toe）的边破左边的区域失效。这由建立一个“region”来实现（即，我们需要把网格划分为由一个边界分开的两个区域）。你将注意到，定义边坡面的直线,始于网格内部，且必须延伸至左边界。有两种方法能延伸这个边界。可以用另一个GEN LINE命令，网格的左侧面向边坡脚延伸一条水平线，或可以对这根水平线作记号（MAKE），以通知FLAC是一个存在的边界。这里，我们将利用MAKE命令。输入（the entry）

49

flac>mark i=1,6 j=4

将把沿j=4 从i=1,6的一条边界线作上记号。现在我们成功的定义了两个区域：边坡及边坡上面已取消区域。这可用发布Pr MAKE命令来验证。用以下命令来拆除位于边坡以上的区域：

flac>mad null region=1,10

它将使域1，10所在的区域失效。现在，试绘出网格以便看一下最后的结果。 下面，把重力加速度设置为9.81m/sec2(向下为正)。所应用的边界条件，包括模型前面和后面的滚轴边界，以及一个固定的底边。现在我们要得到一个平衡重力应力的初始弹性状态。用默认极限值的Solve命令可以办到。当达到100N的不平衡时，得到平衡。作为考察运行进展的一种方式，用户可以要求边坡顶点某一网点的y位移史。用命令HIS来实现它。当命令Solve已达到它的极限值时，试绘出这个位移史，以验证网格的确处于平衡状态，（即，plot his 1）通过以下命令，可以得到弹性应力的一个故值或图形表示法： flac>print Sxx Syy

flac>plot Sxx yellow bound* flac>plot Syy yellow bound*

这时，建立一个保存当前文件平衡状态的再起文件是聪明的（smart）。以防万一在进一步运动时，改变材料参数或本构模型，这样做便可节约时间。完成这些研究仅须重新起动弹性状态，因此，免除了重新形成平衡状态的必要性**。以下命令将在默认驱动器上建立一个名为Hill1.sav的重新起动文件

flac>Save hill1.sav

这时，我们应仃止（quit）flac，并通过直接键入以下命令，重新把保存的状态存在后一个时期处：

flac>restore hill1.sav

在完成时间步长处，该模拟仍在RAM中，因此，我们可以从这里继续运行。

对于模拟的下一步阶段，当考虑可能的破坏过程时，我们要设置实际土壤的材性及时间步长。在这一过程期间，将绘出边坡递增位移图形。为了避免分析数据中的任何混淆，我们不能只了解位移变化，而不顾知道模拟开始以来的累计位移（comulative displacements）（即，时间步长=1以来的位移）。在计算过程中，显示法不要求位移，但作为方便用户，保留了每个网点的累计总位移。因此，在计算过程中，我们可把任何点处的位移调整为零（即，初始化），而不会影响最后得到的解答。为此，利用命令： flac>ini xdis=0 ydis=0

这样，绘出的位移，或打印出的位移，将仅表示从前状态以来的位移变化。

其次，利用PROP命令来设置各个域的材性。对于目前由土壤组成的各个域来说，粘结力设置为零。再次用Region命令来实现。最后，我们把计算模型设置为大应变，以便提供随着边坡破坏呈现出的精确的几何表示。因为边坡崩塌将由于土的休止角（Angle of Repose）出现，我们不要利用Solve命令（不会达到平衡的）***。利用Step命令，并同时直接时间一步入一个小的时间步长数的模拟来停止打印或绘出这些中间阶段，这是很容易的。这里，我们看到了显示法在能处理在整个时间中决不会收敛于一个平衡状态的高度非线性问题方面的能力。

试以200个时间步长间隔进入步长，接着便打印及绘制其结果，其命令为： flac>Step 200

50

flac>print xv yv yd state

flac>plot xv z yellow displ yed bound

以下一组图形，表示出了在75.475及1275个时间步长时，以上仍给出的Plot命令的各个结果。这些图形清楚地说明了边坡的逐渐崩塌。

* 假设在绘图前，屏幕已初始为彩色屏幕，调色板为0。 ** 这里值得提到，特别是当模拟复杂开挖或材性变化时，建设足够的利用重新起动任选两页。

*** 当然，如果采用了较大的不平蘅力，或一个较小的时间步长数，或总运行时间，可以利用Slove命令。

图7-5表示由以下命令产生的初始网格图：

flac>plot pen gird(用自动窗口绘网格) 然后以zoo 补充时间步长问题。这时键入

flac>step zoo

图7-6表示了时间步长为475时，在位移矢量及边坡边界上的x速度矢量。这个图是由以下命令产生的：

flac>plot pen xv z red disp max=0.4 gre boun gre

你能清楚的看出边坡崩塌的各起点。X速度等值线用图形说明了屈服面的发展。图7-7是；另一个800时间步长以后的图形。这里，企图达到土的休止角（angle of repose）时，边坡正在崩塌。显然，随着继续延伸时间步长，各网点的位移变成与事实不符的（这是一个连续介质码）。FLAC能自动地校核过大网格？，如果发出并显示出不损伤运行错误信息，它会停止运算过程。

用户必须了解到，各大畸变是有意的是有限度的。当各域的形状变的畸形到太大的程度时，所得到的计算可能变得不准确。

完成以上模拟以后，应键 flac>new

而一个新的运行 将被初始化。通过规定土的一个小的粘结力来进行另一个模拟，可以得到一个有趣的比较。因为我已保存了初始状态，我们只需从这个时间步长开始重新启动：

flac>rest hill1.sav

随后紧接早先用过 的方法，应键入： flac>ini xdis=0 ydis=0

flac>prop coh=7e4 regin =10,1 flac>set large flac>step zoo

如前所述那样的问题，以便考察变化各特征的后果。

7.3例三：端部有剪力的弹性悬臂梁（Elastic Catilever Beam with End Sheer）

以下一个例题，分析了一个悬臂梁的弯曲，在它的端面作用一个剪力。虽然象flac这样连续介质代码不最适宜于（be not best suited）分析包括弯矩的结果问题，但这个问题提出了一个颇难以证实的试验（this problem presents a fairly difficult）

所给出的输入数据给于以下命令序列中：

51

set log on tit

beam with end shear gr 30,6 *

*5:1 beam with 6 zones across width *

m e

set p-stress col=132

prop d=1000 balck=1e8 shear=0.5e8

*clamp end of beam fix x y i=1 *apple end shear

apple yforce=-0.142857e5 i=31 j=1.7 *(load of 1e5 units) step 3500 set log off save colan.sav return

一个30×6个域的网格代表了平面应力梁。输出也按132列配置。梁的一端x,y向固定，而其余的为无约束。20×105N的负y-力作用于梁的另一端。该荷载分布于7个端点中，每个网点为-0.142857×105N。虽然这个例题进行颇慢，约3500时间步长后得到了收敛。图7-8说明了沿渠的最后y-位移 最后位移为0.41m，与解析解0.39相比，误差为5%。

7.4 例4：弹性、弹塑性及横向各向冈性岩石介质中，受初应力作用的圆形隧洞 （Cirular Tunnel In an Elastic, Elasto-Plastic, and

Transversely-Isotropic Rock Mass Under Initial Stresses）

本例说明了几种本构模型以及GEN命令的应用，这里，用GEN命令去生成 一个网格坐标系 ，应在平面应变的无限岩石介质中，构成一个隧洞的原型边界。 通过仅考察四分之一洞孔，显示了把对称条件引入模型的方法。在本例中，也说 明了利用初始银及应用边界应力，来提供现场应力（in-situ field stresses） 7.4.1 弹性岩石介质 （Elastic Rack Mass）

本例的输入数据给予了以下清单中。该文件可以建立在硬盘中。并可以利 用CALL命令读出，或可以用交互式将它键入。 tit

Hole in elastic medium

* generate the grid

gr 20,20 * elastic model m e

*give initial grid coordinates and grade mesh to boundaries

gen 0,0 0,6 6,6 6,0 rat 1.2 1.2

52

*generate a quarter circle in lower left hand corner of grid *using the arc function gen arc 0,0 1,0 90 * elastic properties

prop s=2.9e9 b=3.9e9 dens=2500 *initial stresses

ini Sxx=-30e6 Syy=-30e6

*boundary conditions for a quarter symmetry problem fix x i=1 fix y j=1

fix x y j=21 fix x y j=21

*monitor history of displacements of tunnel periphery his nste=1

his ydis i=1 j=12 his xdis i=12 j=1

*solve for initial stress state solve

*now excavate hole mod null region=1,1

*time step to equilibrium step 500

*create a save file save ehole.sav

*return to interactive mode return

对于弹性本构规律规定了一个20×20域的网格（400个单元）。我们希望有一个半径为1米的隧洞，且将该隧洞周边（periphery）精细地离散，这样，其结果是合理的。为此，我们想要在包围隧洞区域内，具有网距约为1/5半径的x-及y-网格，但我们也需要相距约为5倍半径的边界（boundaries about 5 rad away），以便避免孔洞对位移的影响。为了建立一个网距0.2米的网格，用GEN命令来在定义x-及y-向0~6米的网格7坐标。如果你想证实GEN的结果，应利用print x及print y命令。

用GEN ARC命令来构成隧洞周边。该命令的效果可以通过设置窗口 flac>wind -0.5,2 -0.5,2

来观察。并可以通过绘出网格来观察。在模型中的初应力已平衡后，我们将开挖隧洞周边内部的各个域。

其次，各边界条件必须设置成体现出沿左边及下边界的对称条件，及位于“无限远”处的固定位移边界（顶边及右边）。为此，沿i=1,j=1,21（左边）及j=1,i=1,21（右边）边界，布置滚轴，而沿i=21,j=1,21(右边)j=21,i=1,21（顶边）边界布置固定点。有了已固定在空间中的网格后，可以利用命令initial Sxx= Syy=来施加？地应力(in-situ stresses)。

如同上一个例题一样，在开挖以前这些应力必须先平衡（因为是客观事实）。这只需要一个时间步长，因为FLAC直接把这些应力装入其数组的适当位置。这

53

可以用以下命令之一来完成

flac>step 1 及

flac.>solve

你可以用prin Sxx 或print Syy命令来核实所得到的原有应力状态*。此时，作好建一个save.file，这样，可以从这个初始状态来重新起动该文件。

再其次，通过使隧洞各域失效的方法来建立开挖。这可以用一组model null命令来办到。为了检验你已失效的那些原有域，试绘出网格*。整个网格及开挖区的放大图，如图7-9示。这一图形由以下命令生成： flac>plot pen gird * plots full gird flac>wind -0.5,2 -0.5,2(建立小窗口) flac> plot pen gird(绘开挖周围的网格)

现在，把问题步入平衡。对于小应变逻辑的20×20网格，在带慢速8087（4.7MHz）芯片（chip）的标准PC机上，约需30秒/时间步长，才能达到平衡解（取决于你定义的平衡）。

最好的做法是，绘制遭受最大位移变化各处（即，开挖边界处）的位移史图，以确定模型朝平衡的进展。你将能迅速的看出问题的收敛特征，你可以用PLOT命令在屏幕上绘出应力及位移。记住，如果你想用彩色屏幕，你必须键入SET CGA（或SET EGA），并使用你所酷爱的调色板选择。

把这些结果与平面应变无限大板中有孔洞的解析解（kirsch）相比较，其结果如图7-11所示。如你所能看到的，这些解比较起来是相当好的。引起的误差颇趋向于隧洞的域的划分。用一个30×30域的模型作同样的运行，其吻合度在5%以内。

7.4.2 弹塑性岩石介质（Elasto-plastic Rock Mass）

下面，对同一问题进行比较，但利用弹塑性本构规律。唯一需要调整的输入数据是本构定律反其特征定义。为此，应把读取的输入文件的第二、第三行改为

flac>mohr

flac>prop s=2.8e9 b=3.9e9 d=2500 fric=30

用solve命令来使各应力达到平衡，并需要一个时间步长。用以前的mode null region命令来挖隧洞。最后，用inital xd=0 yd=0命令可以重新设置各位移。对于平衡的模拟，需要500~600时间步长——或在标准PC机上约为20~3小时求解时间。在AT机上，或在带序协处理器（6-12MHz）的等数机上，或当利用有提高了速度的处理板的商业用计算机时，其时间将显然地减少。

对于塑性情形，孔口周围的应力分布图7-12示。该图明确的说明了隧洞周围形成的切向应力数值，已超过它的极限。理论的和预测的分布图与理论解（BRAY 1967）相当匹配。破坏区半径由下式给出：

2pqu[1tg2(45)ccot2rbr=a [1tg2(45)]ccot2式中 a——隧洞半径 c——粘结力

p ——远端场的流体压力

1/ （7-1）

tg(45/2)1

tg(45/2) qu——单轴压力强度

——摩擦角

塑性域中的切向及径向应力为：

ccot[(r/a)1] rccot[(r/a)1] （7-2） 而弹性区的应力为：

p(/r2) rp(/r2) （7-3）

[tg2(45/2)1]pqu2式中 rbr 2tg(45/2)1这时，用（7-1）式的破坏区半径为1.85m。你可以考察用flac来预测的屈

服区半径范围，这时，键入

flac>print state

或

flac>polt state

各状态指示器表明，屈服区约为1.6米，且与解析是合理的。再后，误差是由相对正常网个引起的。一个30×30域的模型吻合在7%以内。

该问题的主应力矢量，位移矢量及剪应力等值线绘在图7-13，7-14、7-15中。这时图形是由以下命令产生的：

flac>set cga pal=0(如果你有一台彩色监视器)

flac>wind -0.5,2 -0.5,2 flac>plot pen str bou

flac>plot pen disp max=0.3 red bou green flac>sclin 1 0,0 5,5

flac>plot pen Sxy z boun green

如果你想只有一幅屏幕图形，便去掉关键字PEN。 7.4.3 横向各向同性岩石介质 最后一个算例是对横向冈性岩石介质中的隧洞进行的，它受有双轴初始应力。隧洞的位移解是由EISSA给出的。在本例中，用于使初值达到平衡的边界压力如图示。如前两个隧洞例题一样，我们用同一网格从弹性岩石介质开始。键入的数据与前两个例题一样，不同处为：材料模型及性质必须改变，而且不是用固定边界条件，应施加边界压力，以维持现场应力。应注意，这里的对称性只在各向异性角等于0时才有效；如果它等于某个角度，其对称性不切合实际的。

通过固定沿i=1线的x-速度及沿j=1线的y-速度，再次取1/4对称部分。为了提供边界荷载，我们在各对称边界上利用Apply命令：

flac>apply press=1.0e6 i=1,21 j=21

flac>apply press=0.5e6 i=21 j=21,1

这些命令把一个双轴压力作用在边界上。应力的方向由受应力作用的各节点的编

55

号来控制。（祥见APP-命令）。我们也用命令Init Sxx=--0.5e6 Syy=-1e6来建立初应力状态，以致于在1~3个时间步长内平衡。通过键入print Sxx Syy可以随时考察这些应力。

一旦这些应力已达平衡，应利用以下命令来 材料模型及给出新的特征： flac>model anis(把各域设置为各向异心房模型)

flac>prop s=0.1103e6 xm=0.976e6 ym=0.953e6 flac>prop nuy=0.167 nuz=0.165 dens=2000

如以前一样，现在应使各位移为零，并开挖隧洞。利用step或solve命令。这时可以让问题？到平衡。约需要1000个时间步长，在标准IBMPC机上，每个步长约花30秒。

所得到的主压力及位移矢量图形，重叠在边界轮廓线上，如图7-17及7-18所示。这些图形由以下命令组产生的：

flac>wind -0.5,2 -0.5,2(设置一个小宽口)

flac>set cga pal=0(设置彩色屏幕。0号调色板) flac>plot pen str boun green(绘主应力)

flac>plot pen max=4 red bo gre(绘位移矢量)

如果只要求一幅屏幕图形，取消关键字PEN。如果你有一抬双笔绘图机（Two-pen ploter），命令行中的各种颜色，以便改变每一个新图形所遇到的笔颜色（pen color）。主压力总是以1#颜色绘制，而主拉应力将用2#颜色绘制。

图7-19绘出了计算与解析位移，它们约为于x轴所成角度的函数。本例的最大误差接近1%，它主要是由闭合的邻近应力边界及粗糙的网格引起的。

8.0运行FLAC时值得注意的重点或注意事项（Important points or acaautions

note when running FLAC）

本节祥述在运行FLAC中的某些难点，我们已发现这些难点，对某些用户来说，可能构成错误，随着用户反馈的继续，我们将扩充及本节内容。 8.1 初始化各变量（Initializing Variables）

通常的做法是，在两次运行之间初始化各网点的位移，有助于一个模型的互相贯通，该模型采用了许多不同的开挖阶段。这是可以办到的，因为在计算中，代码不要求各种位移——方便用户，可以由网点的速度求出位移。

然而，初始化速度是一件难事，如果把网格的速度设置为常数值，直到重新设置为止，我们将仍为这个常数值。因此，不准把网点速度初始化为零，这会直接清除它们——影响模拟结果。

8.2 改变材料模型（Changing Material Models）

在一个模拟过程中，FLAC对你可以利用的不同材料模型的个数，是有一个的。代码已选定的允许用户对每个域（如果你愿意）具有一个不同的材料（对你的FLAC版本的最大号码网格而言）。不过，如果对于最大型问题（larger-Sized problem），你要改变整个网格的材料模型，可能回出现内存的问题。

8.3 运行含现场应力和重力的问题（Running Problems with in-situ field

stresses and gravity）

有许多问题必须把现场应力及重力作用于模型。深埋的随挖随填的采矿，就是这类问题的一个例子。这时，岩石的质量遭受到强现场应力作用（即，可略去限定网格范围的重力应力），但其就位的回填支柱，将在深井开采中产生重力应力， 在荷载下矿井可能崩塌。在这些模拟中应注意的重点（如同作用上重力的

56

任何模拟一样），至少是网格上的两个点必须固定于空间中——否则，整个网格会因重力而平移。如果你曾注意到整个网格沿负y-方向平移，那么，可以判定，你忘记了把网格固定在空间中。

如果你想只有随深度变化的场地应，应利用含深度变化任选项的INITIAL命令：

flac>ini Sxx=-30e6 (Var=0,10e6)

这将产生从网格底部处（i=1,j=1）为-30e6，变化到网格顶部处（i=n,j=n）为-20e6的xx-应力，因为i,j范围没有给出，FLAC取为整个网格。在本例中，x-变化（x-variance）设置为零。

57