FLA3D-杂

FLA3D命令集锦(新手进阶必备)

GROUP命令 .......................................................................................................................................................................... 2 FLAC3D的背景颜色 .............................................................................................................................................................. 2 RANGE命令 .......................................................................................................................................................................... 2 MACRO命令 ......................................................................................................................................................................... 2 GROUP—RANGE—MACRO 不同命令实现同一功能 ............................................................................................................ 2 视图操作 ............................................................................................................................................................................... 3 ;初始化条件 .......................................................................................................................................................................... 3 ;计算平衡条件 ...................................................................................................................................................................... 3 HIST命令 .............................................................................................................................................................................. 3 LOG命令 ............................................................................................................................................................................... 3 切片 ...................................................................................................................................................................................... 3 ANY-并集运用 ....................................................................................................................................................................... 3 FISH 主要语句 .................................................................................................................................................................... 3 运算 ...................................................................................................................................................................................... 4 单元/节点遍历 ...................................................................................................................................................................... 4 STRING连接（保存文件） ................................................................................................................................................... 4 TABLE操作 ............................................................................................................................................................................ 4 全部节点的应力输出 ............................................................................................................................................................ 4 OPEN、CLOSE 用法 .............................................................................................................................................................. 4 READ(AR, N) 将N个记录写到数组AR ................................................................................................................................ 5 INTERFACE 的变量命令 ...................................................................................................................................................... 5 GP ZONE命令 ........................................................................................................................................................................ 5 用户自定义的 节点、单元额外变量 ................................................................................................................................... 6 TITLE命令 ............................................................................................................................................................................. 6 应力梯度 ............................................................................................................................................................................... 6 基本形状网格 ....................................................................................................................................................................... 6 后处理 ................................................................................................................................................................................... 7 剖面应力SZZ和位移ZDISP数据都导出来 .......................................................................................................................... 8 MOVIE命令 ........................................................................................................................................................................... 8 结构单元 ............................................................................................................................................................................... 9 DIM 内部区域的尺寸EDGE FILL ........................................................................................................................................ 9 FLAC3D自定义本构模型翻译 ............................................................................................................................................... 9 自定义新的本构模型 .......................................................................................................................................................... 10

1

三维显示有限差分基本原理 .............................................................................................................................................. 12 工程实例 ............................................................................................................................................................................. 21 1、FLAC3D 锚杆建模中出现的问题 ..................................................................................................................................... 21 2、锚杆受拉超限 ................................................................................................................................................................... 21 3、FLAC3D计算完毕关机 ...................................................................................................................................................... 23 4、隧道模拟 ........................................................................................................................................................................... 23 5、显示最大位移点的坐标参考 ........................................................................................................................................... 24 6、固结问题输出某一节点随时间的孔压变化 ................................................................................................................... 24 7、初期支护中的“钢拱架”建立 ............................................................................................................................................ 25 8、最大、最小主应力的分布说明了什么 ........................................................................................................................... 26 9、一维固结 ........................................................................................................................................................................... 26 10、基坑桩锚支护-预应力问题 ............................................................................................................................................ 27 11、FLAC3D模拟金属矿采矿与充填应用注意事项 ............................................................................................................ 30 12、FLAC3D塑性区的讨论 .................................................................................................................................................... 30

GROUP命令

group soil range z 1 2 ;定义group

prop bulk 7.8e6 shear 3.0e6 coh 10e3 fric 15 ran group soil ;给group赋值

model null range group soil y y1 y2 ;开挖group soil y方向范围y1 y2（坐标） model elastic range group soil

plot block group range group dam;只显示group dam ;**************************************** FLAC3D的背景颜色 pl set back white RANGE命令

（1）range name trench x 0 1 y 0 4 z 0 2 ;重新命名组trench，范围为x 0 1 y 0 4 z 0 2

model null range trench ;开挖新定义组trench （2）model null range x=2,4 y=2,6 z=5,10 range name Big_Brick x -3 3 y -2 2 z -1 1 model elastic range Big_Brick

prop bulk 1e8 shear 1e8 range Big_Brick

（3）range name Layer1 plane dip 0 dd 0 ori 0 0 0 above

range name Layer2 plane dip 0 dd 0 ori 0 0 0 below ;以面dip 0 dd 0 ori 0 0 0分成Layer1、Layer2

（4）range cylinder end1 x1 y1 z1 end2 x2 y2 z2 radius r

cylindrical range with one end of the cylinder axis (end1) at location (x1, y1, z1), the other end (end2) at location (x2, y2, z2), and with a cylinder radius of r '由(x1, y1, z1) 、(x2, y2, z2)两点确定 旋转轴

;****************************************

MACRO命令

macro Sand 'bulk 1e8 shear 0.5e8 coh 0 tens 0 fric 35' macro Clay 'bulk 1e7 shear 0.3e7 coh 1e7 tens 0 fric 0' prop sand range Layer1 prop clay range layer2 macro Pt0 'p0 0 0 0'

macro Pt1 'p1 add 10 0 0' macro Pt2 'p2 add 0 10 0' macro Pt3 'p3 add 0 0 10' macro Model_Size 'size 4 5 6'

macro Big_Brick 'gen zone brick Pt0 Pt1 Pt2 Pt3 Model_Size' Big_Brick

macro 'Pt0' 'p0 15 15 15'

gen Big_Brick ; this will cause an error

;**************************************** GROUP—RANGE—MACRO 不同命令实现同一功能 （1）using a RANGE object ...

range name Big_Brick x = (-3,3) y = (-2,2) z = (-1,1) model null range Big_Brick （2）using a GROUP object ...

group Big_Brick range x = (-3,3) y = (-2,2) z = (-1,1) model null range group Big_Brick （3）using a macro object ...

macro Big_Brick 'x = (-3,3) y = (-2,2) z = (-1,1)' model null range Big_Brick

（4）gen zone brick size 10,10,10

macro SiltySand 'bulk 1.5e8 shear 0.3e8'

macro ClayeyGravel 'bulk 1.5e8 shear 0.6e8 fric 30 coh 5e6 ten 8.66e6'

model elas range z 5 10 prop SiltySand range z 5 10

2

;**************************************** 视图操作

旋转操作：x y z m ctr+r 还原

ctr+p 打印图片 ctr+z 鼠标选择

ctr+g 彩图变为灰色图

shift + （x y z m） 旋转 缩放 edit---copy to clipboard----粘贴到 word

;**************************************** ;初始化条件

initial xdis 0 ydis 0 zdis 0 ;xyz方向位移为零

initial xvel 0 yvel 0 zvel 0 ;xyz方向位移速率为零 initial state 0 ;模型全部状态归复弹性状态 set grav 0 0 -9.81 ;设置模型加速度

;**************************************** ;计算平衡条件

set mech force=50 ;最大不平衡力小于50n，停止计算

solve force 1e3 ;最大不平衡力小于1e3n，停止计算 solve ratio 1e-6 ;计算位移速率小于1e-6，停止计算 ;**************************************** HIST命令

Hist reset ;清除已有历史信息

set hist_rep 1 ;每一个时步记录一次 hist n=5 ;每5个时步记录一次

hist write 1 3 vs 2 begin 150 end 375 file xx.txt ;记录从150到375单元的1、3项与2项的关系数据，并保存为文件xx.txt hist gp zdisp 4,4,8

hist id=1 gp zdis 1 1 3 ;ID为监测变量的编号，默 hist id=2 gp zdis 1 2 3 认的ID是从1开始编号。建议对监测变量进行编号，以便后处理、调用

plot his 2 3 vs 4 ;plots histories 2 and 3 versus history 4; history 4 plots along the abscissa（横坐标）. plot his 2 3 ;显示监测ID编号2 3的曲线 pl his xla 'string' ylab 'string' ;设置x、y轴的名称 pl his xmaximum var1 xminimum var2 ;设置x轴的最大为var1 最小刻度为var2

pl his ymaximum var1 yminimum var2 ;设置y轴的最大为var1 最小刻度为var2

;**************************************** LOG命令 set logfile xx.log set log on

print …… ……

set log off

;**************************************** 切片

Pl create S ;定义新视图S

pl set plane ori 0 1.5 0 norm 0 1 0 ;定义剖面位置 pl con zd plane ;显示变量 pl add ske

pl add dis plane pl add axe

pl con sxx outline on ;查看水平应力云图， outline on表示显示模型的边界，默认为 off

;**************************************** ANY-并集运用

（1）print gp dis range id 517 any id 533 any '输出2个节点的变形值 any表示‘并集’ （2）指定边界条件

fix x range x -0.1 0.1 any x 5.9 6.1 any '这个表示在x=0和x=6.0方向上固定x方向的位移。any的意思是表示\且\

的意思，也就是说本来要写两行语句的：fix x range x -0.1 0.1，fix x range x 5.9 6.1，现在用any就可以只写一条语句。

fix y range y -0.1 0.1 any y 5.9 6.1 any

;**************************************** FISH 主要语句 （1） if….then ………… else

………… endif command ………… endcommand

loop n（1,21） ………… endloop

loop while ………… ………… endloop

caseof ………… ………… case n1

3

………… case n2 ………… endcase

;**************************************** 运算

= # > < >= <=

负数在运算时， 要加（），以免和减号 - 混淆

degrad 'π/180 第2行： x1 y1 第3行： x2 y2 ……

空行 '

注：在表的文本文件最后，需要有一个回车换行符，否则会出现\的错误，

完成的文件需要读入操作才可以供flac3d调用。采用 table read 命令进行读入： table 1 read xxx.dat

读入后，可以使用plot table 和print table命令查看生成的表文件

apply xvel 1.0 hist table 1 range x 5.9,6.1 y 0,6 z 0,2 pi 'π

ngp '节点总数 nzone '单元体总数

;**************************************** 单元/节点遍历

p_z = zone_head

loop while p_z # null …………

p_z = z_next(p_z) endloop

p_gp = gp_head

loop while p_gp # null …………

p_gp = gp_next(p_gp) endloop

;**************************************** STRING连接（保存文件）

save_file=string(n)+'_step.sav' D3DECFile = 'D3dec_Model.dat' FlacFile = 'Flac3D_Model.dat'

file_name='7-6_add'+string(n)+'.sav' save file_name cal flacfile

;**************************************** TABLE操作

table 1 name load_settlement '创建新表格 xtable(n,s)=…… 对编号为n的表 的第s行、x列进行赋值

ytable(n,s)=…… 对编号为n的表 的第s行、y列进行赋值

plot table 1 line

plot table 1 both '点、线同时显示

常用编辑文本文件的方法进行表的读入与调用 第1行： 表的名称

'1.0表示 表格中的 y向数据的因子

TABLE n x1 y1 . . . erase erases all entries in table n. insert

name 'string' 'changes the name of table number n to 'string'. The table ID number is not changed. read filename 'reads file, filename (in the format described below), and places it in table n. ;**************************************** 全部节点的应力输出

将FLAC中节点以及节点对应的坐标全部输出到txt文件中？以及将每个节点的应力输出 set log on print gp pos set log off

;**************************************** OPEN、CLOSE 用法

;close 当前文件关闭，返回值为0表示成功关闭

open(filename,wr,mode) '先打开文件，然后 可以读，可以往里写内容

;===filename可以是带双引号的字符串，也可以是变量名。 ;===wr必须是个整数：0-只读打开，文件要求存在；1-文件为写入打开，已经存在的文件将被覆盖 （overwritten）。2-文件为写入打开，存在的文件被改写(appended to) ;===mode必须为一个整数值，0-读写fish变量，只有整型，浮点型和字符串的数据类型被传输，不传递变量名。1-读写ASCII数据。可读入一行数据，行之间用 CR/LF分开。每行最大80字符。

;===2-设定为binary读模式，任何文件都以binary的读取模式打开.

;===返回值： 0表示文件打开成功；1-表示文件名不是字符串；2-文件名是一个空字符串；3-wr或mode不是整型；4-mode参数错误（不是0或1）；5-wr参数 错误（不是0或1）；6-不能打开所要读取的文件，比如文件不存在；7-文件已经打开；8-不是fish模式的文件

4

;read(ar,n) ---放入数组。读取n个记录放入数组ar，每个记录要求是一行ASCII数据，或者是一个FISH变量。数组至少有n个元素大小。

;====返回值0表示无误；-1：表示读取错误，比如到了文件的末尾；n-表示读入n行后，恰好到了文件的末尾。

;****************************************

READ(ar, n) 将n个记录写到数组ar

read(ar, n) reads n records into the array ar. The

array ar must be an array of at least n elements. The

returned value is:

0 requested number of lines were input

without error

-1 error on read (except end-of-file)

n positive value indicates that end-of-file

was encountered after

reading n lines

write(ar, n) 将数组ar的前n个记录 写到文件，

详细：

write(ar, n) writes n records from the first n

elements of the array ar. The array ar must be an

array of at least n elements. The returned value is:

0 requested number of lines were output

without error

-1 error on write

n positive value (in ASCII mode) indicates

that the nth element was not a string

parse(s, i) This function scans the string s and

decodes the ith item

pre_parse(s, i) This function scans the string s and

returns an integer value according to the type of the

ith item, as follows. 0 missing item 1 integer 2 float 3 string missing (unable to interpret as int or float) type(e) 数据类型函数，1=整型，2=浮点，3=

字符串，4=指针，5=数组

;**************************************** Interface 的变量命令

interface 1 face range x x1 x2 y y1 y2 z z1 z2

Interface 变量

Interfaces may be identified with the following

functions.

i_find(id)————address of the interface with ID id; returns null if not found (pointer)

i_id(p_i) ————ID of the interface at address p_i (integer) The interfaces may also be scanned starting at i_head and stepping through the list using i_next.

i_head————address of the first interface in the list of interfaces (pointer) i_next(p_i)————address of the pointer to the next

interface in the list (pointer)

wrap

For example, the following command would find the

twinned faces between group \and group \and

put interface elements on these \faces. Only faces

with centroid within the range x 50.0 75.0 would be

considered.

interface 1 wrap rock soil range x 50.0 75.0

Structural Elements 的变量命令参见ftd128.pdf（64/114）

;****************************************

nd_pos( np, p, dof ) position p (p ∈{1,2} denotes current

or reference position, respectively; dof-component, dof ∈

{1,2,3}). The reference position which stiffness is the configuration for

matrices

have been formed and does not

change during a small-strain analysis. The current position is

updated after each timestep. During a large-strain analysis,

the reference position is set equal to the current position

during each large-strain update.

例子：xt0 = nd_pos(_nd, 2, 1)

yt0 = nd_pos(_nd, 2, 2)

zt0 = nd_pos(_nd, 2, 3)

nd_id(np) ID number of node np. Each node has a

unique ID number.

;****************************************

GP ZONE命令 p_gp =gp_near(x,y,z) 'address of gridpoint closest to (x, y, z) p_z =z_near(x,y,z) 'address of zone closest to (x, y, z) gp_head zone_head gp_next(p_gp)

z_next(p_z) p_gp=find_gp(id) 'address of gridpoints with ID number id

p_z=find_zone(id) 'address of zone with ID number id gp_group(p_gp, ind) ????

gp_region(p_gp, ind) ????

gp_zdisp(p_gp):地址为p_gp的节点的z向变形

5

gp_xpos(p_gp) 'x-coordinate of gridpoint ;--- test of stress rotation --- px=gp_xpos(p_zp) gp_xpos(p_zp)=k*px config gpextra 2 gp_xvel(p_gp) 'x-velocity at gridpoint call fishcall.fis

gen zone brick size 1 1 1 gp_near(x,y,z) '得到靠近坐标(x,y,z)的节点地址

mo el gp_id(p_gp) '节点ID号,是整数

maxdisp=gp_id(p_gp) print gp pos ran id 59 '输出59号节点的坐标位

置 (2.0,2.0,3.0)

注意： 节点指针的循环是从ID号最小的地址

开始的。

zone 函数

z_sig3(p_z) '最小主应力 z_sig1(p_z) '最大

主应力

z_prop(p_z,'young') 'z_prop(p_z,'young')=2e6 z_group(p_z) 'zone group name 是单元所在的组

名变量(string)

z_id(p_z) 'zone ID number (integer) z_model(p_z) ' string类型

z_prop(p_z,string) '可以赋值，可以取值 z_xcen(p_z) 'zone的质心的x坐标

z_sxx(p_z) 'zone的xx-stress

stid=2131

z_p=find_zone(stid)

xtable(3,m)=z_id(z_p)

ytable(3,m)=-z_szz(z_p) ;**************************************** def install pnt = zone_head loop while pnt # null z_depth = -z_zcen(pnt) y_mod = y_zero + cc * sqrt(z_depth) z_prop(pnt, 'shear') = y_mod / (2.0*(1.0+P_ratio)) z_prop(pnt, 'bulk') = y_mod / (3.0*(1.0-2.0*P_ratio)) pnt = z_next(pnt) end_loop end set p_ratio=0.25 y_zero=1e7 cc=1e8 install plot block prop bulk ;**************************************** 用户自定义的 节点、单元额外变量

config zextra n '首先进行单元 额外变量的

配置

config gpextra m '首先进行节点 额外变量的

配置

gp_extra(p_gp,n) z_extra(p_z,n)

prop she 300 bu 300

def ini_coord pnt = gp_head loop while pnt # null gp_extra(pnt,1) = sqrt((gp_xpos(pnt)-xc)^ 2+(gp_zpos(pnt)-zc)^ 2) gp_extra(pnt,2) = atan2((gp_xpos(pnt)-xc),(zc-gp_zpos(pnt))) pnt = gp_next(pnt) endloop

end set xc=0 zc=0 ini_coord

;**************************************** TITLE命令 JOB title：title 'Mesh for trench example'

VIEW title：plot set title text 'Mesh for trench example'

;****************************************

应力梯度 ini dens 2000 ini szz -40e3 grad 0 0 20e3 ran z 0 2 ini syy -20e3 grad 0 0 10e3 ran z 0 2 ini sxx -20e3 grad 0 0 10e3 ran z 0 2 set grav 0 0 -10 **** ini xdis=0 ydis=0 zdis=0 hist gp xdisp 1,0,0 hist gp zdisp 0,0,2 **** apply syy = -20e6 grad 0,0,20e5 range y -20.1,-19.9 z 0,10 With this command, a yy-stress component is applied to boundary faces located between y = -20.1 and -19.9. The yy-stress varies linearly with z from σyy = 0 at z = 10 to σyy = -20e6 at z = 0. S=S(0)+gx*x+gy*y+gz*z ' gx gy gz 分别对应grad

0,0,20e5 中的 三个分量 ;****************************************

基本形状网格 gen zone brick size 6,8,8 p0 -10, -10, -20 & p1 10, -10, -20 & p2 -10, 10, -20 & p3 -10, -10, 0

6

plot surf

gen zone brick size 6,8,8 p0 -10, -10, -20 & p1 10, -10, -20 p2 -10, 10, -20 & p3 -10, -10, 0 p4 10, 10, -20 & p5 -10, 10, 10 p6 10, -10, 0 & p7 10, 10, 10 plot surf

gen zone radbrick &

p0 (0,0,0) p1 (10,0,0) p2 (0,10,0) p3 (0,0,10) & size 3,5,5,7 & ratio 1,1,1,1.5 &

dim 1 4 2 fill 'The fill keyword fills the brick region with zones plot surf

gen zone radbrick &

p0 (0,0,0) p1 (10,0,0) p2 (0,10,0) p3 (0,0,10) & size 3,5,5,7 & ratio 1,1,1,1.5 & dim 1 4 2 fill

gen zone reflect dip 0 dd 90 gen zone reflect dip 90 dd 90

gen zone reflect normal (xv yv zv) origin (xv yv zv) ;

; identify the trench

gen zone radc & dim 3 3 3 3 & ratio 1 1 1 1.2 & size 3 8 8 5 & edge 10 &

p0 100 95 100 & fill

gen zone reflect dip 90 dd 90 ori 100 100 100 gen zone reflect dip 0 dd 0 ori 100 100 100

;**************************************** 后处理

plot create GravV

plot set plane dip=90 dd=0 origin=3,4,0 plot set rot 15 0 20

plot set center 2.5 4.2 4.0 plot add bound behind plot add bcont szz plane

plot add axes plot show

****

plot hold hist 1 vs -2 'plot hist m vs n中m代表y轴； n代表x轴；负值是对其值的镜像

pl sk magf 20 'sk 显示模型的外围网格线 pl con szz ou on magf 20

****

实体建模的显示0 E a' B- k, Z, q, plot set rotation 0 0 454 v& n' l# A* a. s. G, T! c5 K

plot block group ;显示块体& F w% |/ ~1 ^9 o plot add axes red

plot add surface yellow8 @6 _) D; {! ` plot show0 P l. M' c m/ c4 N# c ;(1.形成切面的操作) ;以0为中心取切面

plot set plane ori 0 0 0 norm 0 1 0 ;输出szz方向的云图 plot con szz plane 添加网格线 plot add ske

;在切面上添加位移矢量 plot add dis plane

;在切面上增加坐标系$ w5 j& i/ V. J\ plot add axe

;在切面上绘制竖向位移云图: ]\ plot con zd plane% S$ p% e. x1 d2 e5 R, {/ t

;(2.结果输出命令:print(输出关键词)，gp（节点），zone(单元)，history（监测历史）。)# J n4 ?6 n! K4 \\) ;单元应力结果输出

print zone stress% K( T% v. d! x. ?# I

;节点变形的信息$ K% Z% `2 q+ P+ F2 W5 n' Y9 v% V print gp dis n( l) x+ ~4 n+ s

;存储节点信息(可以用记事本打开.log) set log on set 文件名.log print zone stress

print gp dis) W- |6 K: V5 I( { set log off1 \\' ^8 b, r; [8 u6 i

;生产变形矢量图 plot sk dis

;输出变形网格

plot sk magf 20;(变形倍数)( O8 Q ;变形后的网格与应力云图结合起来 plot con szz ou on magef 20

;（4.形成一半的剖面图: e) F- [1 ~% B\ plot set plane ori 0 0 0 norm 0 1 0 plot add surf plane behind yell9 ;在剖面后显示位移云图% H. d& I. c\ plot add cont disp plane behind shade on;: v; 7

;显示单元的id的网格模型 plot block group id on

;显示最大、最小、中间主应力云图/ X Q\@& p5 x( e# y6 plot con smax plot con smin plot con smid

;**************************************** if x_pos = 10.0 then new = get_mem(2) mem(new) = head mem(new+1) = p_gp head = new endif 例子：

def find_add ;定义fish函数find_add head = null ;给head赋值

p_gp = gp_head ;第一个网格结点的指针赋给p_gp loop_while p_gp # null ;当p_gp值不为null时作循环 x_pos = gp_xpos(p_gp) ;将指针为p_gp的结点的x坐标值赋给x_pos

if x_pos = 10.0 then ;如果x_pos = 10.0 则(执行) new = get_mem(2) ;从主内存空间里得到2个fish变量对象并返回第一个对象的开始地址

mem(new) = head ;将head类型和数值置于地址为new的fish变量 mem(new+1) =p_gp

head = new ;将new值赋给head endif

p_gp = gp_next(p_gp) ;将结点指针为p_gp的下一个结点的指针赋给p_gp endloop;结束循环 end ;结束fish函数

实际上这个fish函数为满足条件（x坐标为10的）的zone的地址开辟一定的地址空间，各地址之间存在一定的联系;找到第一个符合条件的zone地址后，用new = get_mem(2) 从主内存空间里得到2个fish变量对象并返回第一个对象的开始地址，并用下面的mem(new) = head ，将第一次开辟的两个变量的第一个变量存储地址head（注意第一次head＝null），第二个变量存储第一个符合条件的zone地址，并将第一个变量的地址赋予head（head＝new），第一次循环结束；下次循环，同样开辟两个变量对象，第一个变量对象记录上次循环开辟的第一个变量的地址，第二个记录第二个符合条件的zone地址，其余循环依次类推，这样子就建立了一个符合条件的zone地址链条，方便以后使用。 在调用时，用ad=head，就将最后一个循环开辟的两个变量的第一个变量的地址赋予ad，进行调用

时，后找到的zone地址将被先调用，最后一直循环到最先开辟的两变量，因最先开辟的两变量的第一个变量的地址为null，因此可以控制循环结束。

;****************************************

ini state 0 ' 是把单元初始化为弹性状态，在施加初应力时可能使得某一部分单元进入了塑性状态。因此把它改回弹性。

;**************************************** id,cid是什么意思？

id是指在整个结构中的编号，而cid是指在某一类比如说cable中的编号。拿cable 中的一个单元来说，它既有自己在整个结构中的cd,又有自己在cable中的cid 剖面应力szz和位移zdisp数据都导出来

restore xx.sav' set log on1 set logfile xxx.txt. print ZONE stress range id xxx print ZONE stress range id xxx

print ZONE stress range id xxx7 i% @. S/ f, w! ===============================% \\+ _4 ?. restore ???.sav set log on

set logfile ???.txt

print gp disp range id ???4 B1 Q4 m\ print gp disp range id ??? print gp disp range id ??? set log off

我用FLAC3D中计算由得到单元应力，怎么计算节点力或节点应力啊？或可以直接输出节点力或节点应力吗？单元应力是平均应力吗？和节点力有什么关系啊？

flac3d可以直接得到zone形心应力，网格点应力是插值出来的，不能直接输出。

;**************************************** MOVIE命令

rest 6-1.sav

ini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0

app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 plot set rot 20 0 30

plot con szz ou on magf 10 plot add hist 1

set movie avi step 1 file 6-5.avi movie start solve

movie finish

;**************************************** CONFIG keyword The following keywords apply.

8

cppudm C++ user-defined models (only available with the C++ user-defined model option) creep creep material analysis (only available with creep model option)

dynamic fully dynamic analysis (only available with dynamic model option)

fluid fluid-flow analysis (see see Section 1 in Fluid-Mechanical Interaction (FLAC3D Manual)) gpextra n

n extra gridpoint variables for FISH use (see FISH REFERENCE)

thermal thermal analysis (only available with thermal model option) zextra n

n extra zone variables for FISH use (see FISH REFERENCE)

;**************************************** ;**************************************** 结构单元

SEL keyword value cable

cable begin . . . end . . . id nseg

pretension property

cable keyword pretension value example：

sel cable beg 1.0,0.4,1.5 end 5.0,0.4,1.5 nseg 4 sel cable beg 1.0,0.4,0.5 end 5.0,0.4,0.5 nseg 4 sel cable beg 1.0,1.2,1.5 end 5.0,1.2,1.5 nseg 4 sel cable beg 1.0,1.2,0.5 end 5.0,1.2,0.5 nseg 4 sel cable prop emod 2e9 ytens 1e8 xcarea 1.0 & gr_coh 1e10 gr_k 2e9 gr_per 1.0

applies given pre-tension force to all cableSELs in the range and with an ID number of id; if id is not given, then all cableSELs in the range are considered. A positive pre-tension force places a cableSEL into tension. The given force is added to the current force being carried by each cableSEL.

property keyword value . . .

assigns the specified property to all cableSELs in the range and with an ID number of id; if id is not given, then all cableSELs in the range are considered. The following properties are available.

densitydensity (needed if dynamic mode or gravity is active)

emod Young's modulus

gr_coh grout cohesive strength (force) per unit length

gr_fric grout friction angle (degrees) gr_k grout stiffness per unit length gr_per grout exposed perimeter slide large-strain sliding flag

slide_tol large-strain sliding tolerance thexp thermal expansion coefficient xcarea cross-sectional area

ycomp compressive yield strength (force) ytens tensile yield strength (force) ;**************************************** dim 内部区域的尺寸EDGE FILL

gen zone radcylinder size 5,10,10,5 &

p0 (0,0,0) p1 add (400,0,0) p2 add (0,200,0) p3 (0,0,200) &

dim (20,20,20,20) fill not

del range group 2 not '模型中只保留 group 2 ，其他的都删除

;****************************************

out(s) 在屏幕上显示s的字符信息，s必须是字符类型。 执行成功，返回0，否则为1.

in(s) 键盘输入函数。等待用户键盘输入 ;**************************************** FLAC3D自定义本构模型翻译

对FLAC3D3.0的udm部分进行了初略的翻译，请注意这里面的编译环境改变。

*.sln包含了一个或多个工程(*.vcproj)，下面演示工作过程：

登陆Visual Studio .NET(Visual C++ 7.0)，选择File/open Solution，选择到udm.sln文件。 1. 编译工程

选择Build/Rebuild solution，将在Release文件夹中产生一个名为userssoft.dll文件。 2. 创建一个新的工程

对用户定义的模型要重新建立一个工程，将现有的工程(*.vcproj)拷贝到一个新的文件夹，具有有三个步骤： （1）拷贝udm.vcproj到mymodel.vcproj（mymodel为自己取得模型名字）& g' （2）在udm solution图表上右键，选择Add/Existing project

（3）选择mymodel.vcproj( 可以看到在solution explorer中有两个相同名字的工程（都成为udm），右键选择第2个工程，将其名字改成mymodel。) G+ t# b( v' H' m

9

3.选择Release/Debug编译选项& B. _- @: j2 n+ e+ S7 N

（1）选择Build/Configuration Manager

（2）选择Release或者Debug: C: I; S1 v2 U3 l( H 注1：vcmodels.lib是一个Release编译文件（不包含任何调试信息），但是通过Debug选项也可以建立DLL文件并读入FLAC3D，虽然它不能全部使用且可能会降低一些优化选项的运行（具体参照Debug和Release的区别）

注2：下面的摄制过程仅对当前编译环境有效。4. 改变定义模型的输出文件名

（1）在mymodel的工程上右键，选择Properties / Configuration Properties / Linker / General

（2）在Output File输入\5. 添加自定义模型的资源和头文件到工程中 . （1）右键usersoft.cpp，选择Remove，同样对usersoft.h进行移出操作

（2）在上述的位置填入新的.cpp和.h文件（选择Add / Add Existing Item）（3）对输出dll（上述的第4步）进行更名，重新对第1步进行编译即可以得到所需的自定义模型文件.DLL。 自定义新的本构模型

FLAC3D自定义本构模型跟FLAC手册中讲到的用FISH来自定义本构模型一样。然而在FLAC3D中不支撑FISH语言来自定义本构模型，自定义本构模型开发必须用C++语言，且编译成DLL文件（动态链接库），动态链接库文件能在需要的时候随时加载上去。本构模型的主函数主要是返回新的应力，给出应变增量。然而自定义本构模型也必须给出一些其他的信息：比如模型名称和读入、写出保存文件等操作。

C++语言是一种面向对象的程序设计语言，使用类（classes）来代表对象（objects）。对象的数据被对象封装起来，在对象的外面数据是不可见的。通过成员函数来访问对象，而成员函数可以对封装的数据进行操作。另外C++语言强烈支持对象的等级结构，新的对象性质可以从一个基本对象产生，基本对象的成员函数可以被派生出来的对象的成员函数替代。这些特点使程序更加模块化。举个例子：主程序需要在程序代码的不同地方建立与派生类的不同变量之间的接口关系，但是这些仅仅只关系到基类，与派生类无关。运行时间系统自动的调用适当的派生类的成员函数。对C++比较好的介绍来自Stevens (1994)；它假象读者有一定得编程语言知识，特别是对C语言的了解。

在4.2部分将介绍怎么用C++语言开发自定义本构模型。这节主要包括基类、成员函数、本构模型编

号、自定义本构模型与FLAC3D之间的传递信息，本构模型状态指示器。在4.3节中将介绍怎样生成DLL本构模型。这一节主要包括自定义本构模型的支持函数，实例本构的源代码，FISH支持的用户自定义本构，和怎样生成和加载一个DLL文件。在这一节中所有的参考文件被包含在“\\ITASCA\\Models\\UDM”文件夹下面的“UDM.ZIP”这个压缩包文件里面。

注意：FLAC3D 3.0版本是用Microsoft Visual C++（VC++）7.1版本编译的。用户自定义的DLL文件最好采用与其相当的编译器来编译，以使用户自定义的DLL文件能与FLAC3D兼容。

1自定义本构的方法 （1）自定义本构的基类

以上介绍的方法是FLAC3D自定义本构支持的方法。基类为从基类派生出来的实际的本构模型提供框架。这个基类叫ConstitutiveModel类，被称为“绝对”的类，因为他声明了许多完全虚有的成员函数（通过=0语法附加到函数原型）。这意味着这个基类不能产生任何对象，以及从这个基类派生出来的任何对象都必须提供真实的成员函数，以替代ConstitutiveModel类中的虚有成员函数。例子4.1，提供了ConstitutiveModel（包含在文件“CONMODEL.H”中）这个类的部分代码，ConstitutiveModel类中的一些成员函数，像公共函数在例子4.1中省略掉了。公有函数的使用（像YoungPoissonFromBulkShear）是不用证明的，有关他们使用的例子可以在提供的本构模型源程序中找到，FLAC3D使用其它的函数来操作和访问本构模型，用户可以毫无理由的使用和重新定义这些。 class ConstitutiveModel { public: EXPORT ConstitutiveModel(unsigned uTypeIn,bool bRegister=false);

EXPORT virtual ?ConstitutiveModel(void);

// ROUTINES THAT MUST BE SPECIFIED BY THE DERIVED TYPE

virtual const char *Keyword(void) const=0; virtual const char *Name(void) const=0;

virtual const char **Properties(void) const=0; virtual const char **States(void) const=0;

virtual double GetProperty(unsigned ul) const=0; virtual ConstitutiveModel *Clone(void) const=0; virtual double ConfinedModulus(void) const=0; virtual double ShearModulus(void) const=0; virtual double BulkModulus(void) const=0; virtual double SafetyFactor(void) const=0; virtual unsigned Version(void) const=0;

virtual void SetProperty(unsigned ul,const double &d)=0; EXPORT virtual const char *Copy(const ConstitutiveModel

10

*cm)=0;

virtual const char *Initialize(unsigned uDim,State *pst)=0;

virtual const char *Run(unsigned uDim,State *pst)=0; EXPORT virtual const char *SaveRestore(ModelSaveObject *mso)=0; }; （2）成员函数

任何派生出来的本构模型必须提供真实的成员函数来代替ConstitutiveModel类中的虚拟函数。这些函数的功能在以下描述。

const char *Keyword()表示返回字符数组的指针，这个字符数组包括本构模型的名字，用户在使用MODEL这个命令时需要用到。举个例子：“elastic”在C++中是个有效地字符串。

const char *Name()表示返回字符数组的指针，这个字符数组包括的本构模型名称是用来在输出时使用的（举个例子：PRINT zone这个命令当中可能用到）。它的名字可能和Keyword()这个成员函数的名字相同，也有可能不同，但是要注意FLAC3D在输出的时候会截去长的字符串。举个例子：“Linear/elastic”是个有效地字符串。 const char **Properties()，表示返回字符串数组的指针，这个字符串数组包括模型的材料参数名称，和一个指向这个字符串末尾的空指针。接下来的例子是一个有效地字符串数组：{“shear”，“bulk”，0}，给出的材料参数名字将会被PROPERTY命令识别。注意到字符串数据就像上面一样必须以0为结尾。 const char *States()表示返回字符串数组的指针，这个字符串数组表示单元状态名称，以及一个指向字符串末尾的空指针。这个名字用在输出或是显示用户自定义本构模型内部状态（比如塑性流动）。接下来的例子是一个有效的字符串数组：{“yielding”，“tension”，0}。变量mState在4.2.4中介绍。注意到字符串数据就像上面一样必须以0为结尾。 SetProperty(unsigned n, const double &dVal)，变量dVal的值在FLAC3D中以“PROP name=dVal”这个命令赋值。变量n是一个从1开始的序列数，通过调用Properties()这个函数

来表示材料参数名称。本构对象需要在他适当的私有内存中存取以上提供的数值。

double GetProperty(unsigned n)返回本构模型材料参数一个序号n（在Properties()函数中定义的，n=1从第一个材料参数开始）

const char *Copy(const ConstitutiveModel *cm)，这个成员函数首先调用基类Copy函数，然后通过cm从对象本构中复制所有必须的数据（假定目前的本

构模型与其是相同的派生类）。假如发生错误，将返回一个描述错误的字符串，反之返回0。当调用Initialize()这个函数时，就没有必要复制数据重新计算。

const char *Initialize(unsigned uDim, State *ps)，当CYCLE这个命令使用时或是大变形模式时，这个函数在每个模型对象中仅仅被调用一次（例如每个单元）。这个对象能够初始化材料参数或是单元状态变量，或是什么都不初始化。求解问题的维数用uDim这个变量来表示。结构体ps（4.2.4）包含本构模型中单元当前的信息。当检测到错误信息时，必须有一个指针指向一个表示错误的字符串，反之就返回0。注意，当调用Initialize（）这个函数时，并不给出应变。对整个单元的平均应力分量可以从状态结构体重得到，他们不能被Initialize（）这个函数的成员函数调用。

const char *Run(unsigned uDim, State *ps)当FLAC3D对每个单元尽心扫描时，这个函数被调用起来主要是针对自结构（每个单元最多可以划分为10个子结构）。这个模型通过应变增量来更新应力张量。Ps这结构体包括目前的应力分量，并且计算每个子结构的应变增量分量。当这个函数调用时，同时也考虑了由于旋转而产生的应力误差。当发现错误时，应当返回一个表示错误的字符串指针，反之就返回0.

double ConfinedModulus(void)这个对象返回一个值，这个值可以更好的估计最大的受压模量。它应用在FLAC3D计算稳定时间步中，对于一个线性的弹性模型来说，这个受压模量是K + 4G/3。

double ShearModulus(void)这个对象返回一值，这个值可以更好的估计当前的正切剪切模量。这个使用在FLAC3D动力本构模型中粘滞静态边界系数。

double BulkModulus(void) FLAC3D目前还没有使用这个对象，但是这个对象可以很好的返回对当前的正切体积模量的估计。

double SafetyFactor(void)这个对象目前也还没有使用，它必须返回一些值，像10.0。

unsigned Version(void)这个对象返回本构模型版本号。这个用来处理早期本构模型版本保存的结果文件，而这早期的本构模型可能忽略了一些变量。

ConstitutiveModel *Clone(void)创建一个新的对象，必须与当前类相同，返回一个ConstitutiveModel类指针。不管FLAC3D什么时候在一个单元中调用本构，它都被调用。 const char *SaveRestore(ModelSaveObject *mso)当SAVE或是RESTORE命令给出时，这个被函数调用。本构模型必须首先调用基类的SaveRestore()这个函数。SaveRestore()允许本构模型保存或是恢复每个对象的数据。仅仅允许保存整数类型或是浮点数类型，其他的数据类型必须要转换成这两种类型。派生类函数必须首先调用mso->Initialize(nd,ni)，这里nd表示双精度存取（或是恢复）数据字节数，ni表示整型存取（或是恢复）数据字节

11

数。变量通过mso->Save(ns,var)这个函数识别，这里ns是变量的序号（从0到nd-1或是从0到ni-1），而这主要取决于实际整型或是浮点型数据存取或恢复的位数，var是需要保存的变量。有单独的Save()去处理整型或是实型变量。这个结构体类ModelSaveObject是不怎么重要的，除了以上提到的函数。它被定义在“CONMODEL.H”这个头文件中。 重新定义的类也必须包含构造函数，而这个构造函数调用基类的构造函数。如果这个bRegister变量的值为true，基类的构造函数就被调用，然后派生出来的本构模型就被FLAC3D注册登记过。一个与模型相关的数（uTypeIn）也必须通过；这使得当从一个文件中恢复出来的单元能够重新设置正确的本构模型。一般取一个比较的值来表示这个（比如100或是更大），以避免与从1开始的模型内置变量产生冲突。在其他的所有情况中，派生类构造函数在被调用时，应当不含任何参数，就像Clone这个成员函数一样。通过构造函数可以初始化成员数据，就像例4.2所示。在这个例子中本构模型的特别号是整型变量mnUserMohrModel（具体看例4.5），这些符号“dBulk”、“dShear”等等是派生类的数据成员。

Example 4.2 Typical model constructor

UserMohrModel::UserMohrModel(bool bRegister) :ConstitutiveModel(mnUserMohrModel,bRegister), dBulk(0.0), dShear(0.0), dCohesion(0.0), dFriction(0.0), dDilation(0.0),

dTension(0.0), dYoung(0.0), dPoisson(0.0), dE1(0.0), dE2(0.0),

dG2(0.0), dNPH(0.0), dCSN(0.0), dSC1(0.0), dSC3(0.0), dBISC(0.0), dE21(0.0) { } （3）本构模型的编号

每个用户自定义本构模型都包含它自己的名称，它的材料参数的名称，以及单元应力状态指示器。FLAC3D通过调用合适的成员函数来实现以上这些，就像在4.2.2节中介绍的一样。FLAC3D通过一个静态的、全局的本构模型例子来调用构造函数，以识别用户自定义的本构模型（参见例子4.3）。当FLAC3D被加载（主要是内置本构模型），或是当一个DLL文件被加载上（外置本构模型），这时这个构造函数被创建。这个参数的值为true时导致了基类构造函数对新本构模型的注册，然后把它加在本构模型列表中。只有一个静态的、注册的、特别的本构模型被定义，它是为了方便插入到C++本构模型源代码中，以至于当这个模型相关的DLL文件被加载时，它能够被注册。当FLAC3D需要本构模型的任何信息，或是需要去实例化一个模型时（使用Clone函

数）时，这个静态的实例化模型被调用。

Example 4.3 Global instantiation of a model object static ElasticModel modelInstance(true);

// ... forces a constructor call to the model registry （4）计算中本构模型和FLAC3D软件之间的交流 三维显示有限差分基本原理

当FLAC3D达到平衡或是稳定的塑性流动时，它通过显示有限差分来模拟三维连续介质的力学行为。监控的力学响应主要是通过特殊的数学模型和数值计算过程得到。接下来介绍这两方面。 1、数学模型描述

介质的力学行为主要来源于一般原理（应变定义、运动规律），和理想材料的本构关系。这个数学结果表达式通常是一些偏微分方程，涉及到力学（应力）和运动学（应变率、速度）变量。这些偏微分方程联合个别的几何关系、材料参数，以及给定的边界条件和初始条件就可以求解。

虽然FLAC3D在平衡状态附近，主要关注介质的应力状态和变形，但是必须要注意到该数学模型中的运动方程。

（1）符号约定

在FLAC3D中采用拉格朗日算法，介质中的一个点，通过矢量xi，ui，vi和dvidt，i?13，来定义一个点的坐标，位移，速度和加速的。记号ai表示矢量?a?的第i个分量，在笛卡尔坐标系中；Aij表示张量?A?的第（i，j）个分量。a，i表示变量对xi的偏导数。（变量a可以使标量，矢量和张量）

默认结构受拉为正，变形伸长为正。爱因斯坦求和记号只针对下标，i，j，k（i，j，k=1，2，3）。应力介质中一已知点的应力状态是通过对称应力张量?ij来表示。任意斜面上的应力矢量?t?可以通过柯西公式得到（拉为正），如下：

ti??ijnj

（2.37）

?n?表示任意斜面上的单位法向矢量

（2）应变率和转动率

假设介质的离子以张量?v?运动。在一个无限短时间

dt内，介质产生一个无限小的应变为vidt，相关的应变

率张量可以写成如下：

12

?ij?12?vi,j?vj,i? （2.38） 第一应变率张量不变量描述了体积单元的的膨胀程度。张量?ij中没有包含变形率，由于速度矢量的平移和角速度的转动，一个体积单元会产生一

个瞬间的刚体位移，如下：

??1i?2eijk?jk （2.39）

eijk表示置换符号，矢量???表示转动率张量，定义如下：

?ij?12?vi,j?vj,i? （2.40） （3）运动平衡方程

采用连续介质的动量原理和柯西公式，平衡方程如下：

?iij,j??bi??dvdt （2.41） ?为介质的密度，?b?表示单位体力，d?v?dt表

示速度矢量对时间的导数。当力开始施加到介质上，平衡方程贯穿在整个的数学模型中，以及单元介质的运动中。在静力计算过程中，d?v?dt为0，公式2.41简化为如下偏微分方程：

?ij,j??bi?0 （2.42）

（4）边界条件和初始条件

应力边界条件主要是通过公式2.37来表示，位移边界主要是通过指定边界的速度分量为0来实现。初始条件中体力也是可以施加的，需要指定初始应力状态。 （5）本构方程

公式2.41和公式2.38中包含有9个方程，15个未知量，这15个未知量是6个应力分量和6个应变分量，以及3个速度分量。通过本构方程可以提供额外的6个方程，这个6个方程描述介质的应力和应变之间的关系。一般定义如下：

????ij?Hij??ij,?ij,?? （2.43） ????ij为共轭应力张量，?H?为一已知函数，?为考虑加载历史变量，

????ij?d?ijdt??ik?kj??ik?kj （2.44）

这里d???dt表示应力矢量???对时间的实导数，???表示转动率张量。 2、数值方程

FLAC3D通过以下三步骤来求解：

（1）有限差分逼近（变量的一阶导数、时间导数用有限差分来逼近，假定变量在很短的空间内和时间间隔内线性变化）

（2）离散逼近（将连续介质离散为与之相当的网格，在这个网格中，所有的力包括外力和内力，都作用在单元节点的三个方向上）

（3）动态求解方法（在平衡方程中引入惯性定律，使得系统慢慢达到平衡）

连续介质的运动定律通过以上三步骤，转化为离散单元节点上的牛顿定律。一般普通的微分方程可以通过时间显示差分求解。

介质的偏微分平衡方程中涉及到的空间导数，出现在以速度来定义的应变率中。为了进一步的定义速度变量和相关的空间间隔，连续介质被划分为常应变率的四面体单元，这个四面体单元的顶点就是网格单元的顶点。如图2.2所示。

图2.2 四面体单元

（1）空间导数的有限差分形式

四面体的应变率张量分量的有限差分方程是通过节点平衡方程得到。四面体共有1-4个节点，节点n所对的面，称之为面n。通过高斯散度定律可得如下方程：

?Vvi,jdV??Svinjds （2.45）

通过高斯散度定律将四面体上的体积分转化为四个面上

的面积分。假设四面体是常应变率，所以速度矢量是线性变化的，每个面的法向方向也是常量，公式2.45简化为：

4Vvi,j??v(f)in(f)jS(f) （2.46）

f?1上标（f）表示面f，vi表示变量vi的平均，假设速度线性变化，可以得出：

13

v(f)13?4i?vli （2.47） l?1,l?f这里的上标l指的是节点。将公式2.47带入2.46可得：

14vl4Vvf)f)i,j?3?i?n(jS( （2.48）

l?1f?1,f?l假设在公式2.45中vi=constant，我们得到散度定理如下：

?4n(f)(f)jS?0 （2.49）

f?1利用时2.49可以将2.48化简为：

v1l)(l)i,j??3V?4vl(injS （2.50）

l?1因此应变率张量就可以表示为：

?14l(l)ll)ij??6V??vinj?vjn(l)i?S( （2.51）

l?1节点平衡运动方程

节点运动平衡方程通过虚功原理来推导，在任一个瞬时，得到一个等价静态问题。通过在节点惯性平衡方程中引入虚功原理，来求解整个结构。

固定时间t，我们通过平衡方程来研究静力等价状态问题

?ij,j??Bi?0 （2.52）

体力的定义如公式2.41

B?dv?i????bi?idt?? （2.53） 在有限差分的框架中，介质由一些连续的承受体力?B?的常应变四面体单元代替。在静态平衡方程中，作用在单个四面体上的节点力?f?n,n?1,4，以及

四面体应力和体力都是通过虚功方程推导出来。假设四面体有一个虚拟的速度??v?n（它在四面体内会产生一个线性变化的速度场??v?和一个常应变率????）

，我们假定外力虚功由节点力?f?n和体力

?B?产生，内力虚功由虚速度下的应力?ij产生。

外力功率为：

4E???vnifni??V?viBidV （2.54）

n?1内力虚功率为：

I??V??ij?ijdV （2.55）

利用公式2.51，公式2.55可以写成常应变率的形式：

??16?4I??vll)(l)i?ijn(j??vlj?ijni?S(l) （2.56）

l?1由于应力张量是对称张量，现在定义一个矢量Tl如下:

Tl(l)i??ijn(l)jS （2.57）

由式2.57，式2.56可写成如下：

??13?4I?vliTli （2.58）

l?1将公式2.53带入公式2.54中，可得：

4E???vnnIifi?Eb?E （2.59）

n?1Eb为体力?bi对外力虚功率的贡献，EI为惯性力对虚

功率的贡献。对于一个四面体常体力：

Eb??bi?V?vidV （2.60）

EI?????dviVvidtdV （2.61） 根据以前的假设，在四面体内速度场是线性变化的。为了描述它，现在定义一个新的局部坐标系，坐标系原点在四

面体形心，坐标轴为x1?,x2?,x3?。四面体内任一点速度可以用四面体4个顶点的速度差值得到，如下：

?4vni???vniN （2.62）

n?1Nn,n?1,4为线性函数

Nn?cn0?cn1x1??cnx?2?cn23x?3 （2.63） cnnnn0,c1,c2,c3,n?1,4,为常数，可以通过求解以下方程得

出：

Nn?x1?j,xj2?,x3?j???nj （2.64） 14

?nj为克罗内克符号，通过定义质心，这项

?Vx?jdV?0，将公式2.62和2.63带入2.60中得：

Eb4??bncni??vi0V （2.65） n?1由式2.64和形心的性质可知cn0?14，带入式2.65中，得：

Eb??4?vn?bViin?14 （2.66）

将式2.62带入2.61中，得：

4EI????vnnin?1??NdviVdtdV （2.67）

最后将式2.66和2.67带入2.59中可得：

4E???vn??bVdv?i

n?1??fni?i4??V?NnidtdV??（2.68）

在四面体的静力平衡方程中，在任何虚速度下，外力虚功率等于内力虚功率：由式2.58和2.68可知：

nTn?fi?bVdvi?3?i4??V?NnidtdV （2.69）

假定在四面体内，加速度是不变的，则式2.69中最

后一项可改写为：

n?ndvi?dvV?NdtdV??i??dt???V?NndV （2.70）

又因为?为常数，?Vx?njdV?0，c0?14，则式2.7为：

dvn??Nni?V?dvi?VdtdV?4??dt?? （2.71） 定义?V4为虚拟的节点质量mn，节点质量可由式2.72确定，节点质量可以保证在平衡的过程中数

值计算的稳定。式2.71改写为：

n?ndvi?dv?V?NdtdV?mn?i?dt?? （2.72） 因此式2.69改写为：

n?fnTni3??bVi4?mn??dvi?i??dt?? （2.73）

至此等价系统的平衡方程已建立，如式2.73，在每个节点，总的静态平衡力??f?，包括所有四面体的贡献、节点贡献荷载?p?以及集中力为0。为了描述以上规律，我们约

定如下记号：如果一变量带有上标，就表示节点变量，表示该节点在全局节点编号中的编号。???.???l??表示与

节点l相连的所有四面体对节点l的贡献之和。采用以上记号，节点牛顿运动定律可以写成如下形式：

?l?F?l??M?l???dvi?i?dt?? l?1,nn （2.74）

n?n为介质的节点总数，节点质量M?l定义如下：

M?l????l???m??? （2.75）

不平衡力?F??l?定义如下：

?l?F?l?i?????l???Ti?bV???3?i4?????Pi （2.76）

当系统达到平衡时，不平衡力渐渐趋于0。 时间导数的显示有限差分

考虑本构方程2.43和变形率和节点速度之间的关系2.51，公式2.74可以写成如下形式：

dv?l?idt?1M?l?F?l?i?t,?v?1??2??3??p??l?i,vi,vi,,vi?,?? l?1,nn （2.77） 记号{}?l?表示计算节点l速度的速度子集。在FLAC3D中

求解2.77一般采用对时间t的显示有限差分。在这个方法

中假定节点速度在时间间隔?t内线性变化。公式2.77中左边的导数采用中心差分来代替。中心差分时，当计算位移、力时只存取一半时间步的速度。节点速度通过以下迭代公式计算：

v?l?(t??t)?v?l??t?t?1??2??3??l?ii(t?2)??l?M?l?Fi?t,?vi,vi,vi,,v?p?2i?,?? （2.78）

同理节点坐标更新也是通过中心差分的形式得到：

x?l?l?i(t??t)?x?i(t)??tv?l?i(t??t2) （2.79） 15

从式2.78和2.79中可知，一阶误差项不存在，因此上述迭代方法是二阶精度的节点位移迭代公式如下：

u?l??l??ti(t??t)?u?l?i(t)??tvi(t?2) （2.80） u?l?i(0)?0

本构方程的增量形式

在FLAC3D里面假定在时间间隔?t，速度为常量，因此本构方程2.43可以写成如下形式：

????H*ijij??ij,?ij?t? （2.81）

???*ij为共轭应力增量，Hij为一已知函数。 假定位移在时间?t内，线性变化：

?ij?t???ij （2.82）

??ij为时间t时的应变改变量。 应力增量通过???ij计算得到，如下： ??ij????Cij???ij （2.83） ??Cij为一应力修正，定义如下： ??Cij???ik?kj??ik?kj??t （2.84）

转动率张量由公式2.40或是2.51计算

大小应变模式

以上描述的微分方程都是描述大变形，涉及到大位移，以及位移的线性变化和转动。在FLAC3D里面称之为大变形模式。

在实际应用中，转动足够小，以至于分量?ij??ij与整体相比很小，可以忽略不计。张量

???可以用?I?来代替，公式2.83中的应力修正项可以忽略。当然在小位移和位移变化中，应变率张量中涉及的空间导数，可以通过初始值估计。节点坐标也不再由2.79式更新。

在FLAC3D里面小变形假定小位移、位移梯度和转动，节点坐标不再更新，以及不考虑应力转动修正。 保持数值稳定的力学时间步

公式2.78能求解出正确解答的前提是数值计算的稳定。通过一个理想的质点-弹簧系统，可以得出一些观点。质点-弹簧系统的平衡方程如下（矩阵符号

表示）：

?P*???K??u???M???dv??dt?? （2.85）

这里大括号表示节点矢量，?P*?表示外力，?K?表示弹

簧系统的刚度矩阵，?M?表示分布质量矩阵。如果我们能解释公式2.74中的不平衡力作为外力施加到弹簧系统上，以及公式2.85中弹簧的反力，那么理想材料的近似就是可取的。

在采用有限差分来分析质点-弹簧振荡系统时，时间步不能超过临界时间步，而这临界时间步与整个系统的最小固有周期有关。因此在有限差分中必须提供时间的一个上限，以此来保证数值计算的稳定。

推导临界时间步的关系，需要一些固有周期方面的知识。然而在实际中，全局的特征值分析是不可行的，也不经常采用这种方法来得到固有周期。在FLAC3D里面采用在稳定系统中的局部扰动方法来实现一个目标。在整个数值计算过程中，整个系统都是采用一不变的单位时间步?t?1。

假定公式2.74中右边节点质量为变量，并且参与整个局部平衡过程。

图3 质点弹簧系统1

首先考虑一维自由度弹簧系统，见图3.质点在给定初始位移情况下的运动方程如下：

?kx?md2xdt2 （2.86）

k为弹簧的劲度，m节点质量。这个方程的二阶有限差分临界时间步如下：

?t?T? （2.87）

这里T为系统周期。

T?2?mk （2.88） 16

图4 质点弹簧系统2

现在假定无穷个质点和弹簧见图4，由于结构对称，这个组合系统的力学行为可以采用图4所示简化方法分析。图4中（a）所示系统可以简化为一刚度矩阵为4k的单一弹簧系统。从式2.87和2.88中可推出极限稳定临界公式如下：

m?k??t?2 （2.89）

假定?t?1，如果节点质量大于或是等于弹簧劲

度，则系统式稳定的。在局部分析中，公式2.89的有效性，通过解释m作为节点l的节点质量贡献

ml，以及k和相关的节点刚度贡献kl扩展到四面

体。节点质量贡献通过无穷的准则来提供一个时间

步的上限来导出。节点刚度通过局部劲度矩阵的简单对角化得到。

四面体在节点l的内力贡献为：Tli3，这个力可以假定弹簧的恢复力?kllijuj（式2.85）。假定在时间间隔dt内，可得如下公式：

dTlill3??kijvjdt （2.90） 将公式2.57带入2.90中可得：

d?ij(l)l)3n??klljS(ijvjdt （2.91）

假定在节点l的q方向上有一单位速度，其他方向上的速度为0，则从2.91式可得：

kd?qj(l)(l)qqdt??3njS （2.92）

式2.92中重复指标q并不累计求和。引入小应变的增量虎克定律来描述应力应变之间的关系，式2.92可得：

k?qq?1dtqqdt??3n(l)S(l)q （2.93）

式2.93中?1?K?43G，K为体积模量，G为剪切模量。由于速度已经确定，则由式2.51可得：

?qq??13Vn(l)(l)qS （2.94） 将式2.94带入2.93中得：

k?1(l)(l2qq?9V??n)qS?? （2.95） 现在定义节点刚度矩阵的上限为：

kl?max(k11,k22,k33) （2.96）

从式2.89可得四面体的节点质量贡献为：

ml??19Vmax??2?n(l)iS(l)??,i?1,3? （2.97） 节点质量是为了保证数值计算的稳定。 力学阻尼

平衡方程需要提供阻尼，以得到静态解或是准静态解（非惯性）。在FLAC3D中局部非粘性阻尼为默认的静力计算阻尼。在动力分析中也是可以应用局部阻尼的。

在FLAC3D中求解均匀的稳定的运动时需要选择一种阻尼。这种情况发生在蠕变计算和桩的极限承载力计算中。这个阻尼称之为组合阻尼。组合阻尼在消减动能方面比局部阻尼效率要高。

在动力计算中一般采用瑞利阻尼和人工粘性阻尼。 局部阻尼：

在平衡运动方程2.74中施加局部阻尼力，则平衡方程变为：

?l?F?l??f?l?ii(i)?M?l???dv??dt?? l?1,nn （2.98）

f ?l ? ?l?(i)为阻尼力，通过普遍的平衡力Fi和速度v?l?i定义如下：

f?l?i)???F?l??l?(isign?vi???sign(y)??1,y?0; （2.99）

??1,y?0; ??0,y?0这个阻尼力通过阻尼常数?来控制，?默认为0.8。 局部阻尼具有如下特点：

只有存在加速度时才会产生阻尼力，稳态运动是不会产生

17

阻尼力。

?阻尼常数为一无量纲量。

阻尼力是与频率无关的，在一个集合体中，不同的固有周期，采用相同的阻尼常数。

局部阻尼与滞后阻尼有点相似，每次循环结束时能量损失都是独立的，或是不同的。

由式2.99可知，阻尼力是与与总合成力成比例的。而不像粘性阻尼是与速度的大小成比例的。当F?l?i和v?l?i同号时，2.98可写为：

dv?l?i?F?l?i(1??)dtM?l? （2.100） 当F?l??l?i和vi异号时，2.98可写为：

dv?l?i?F?l?i(1??)dtM?l? （2.101） 以上两式可以写成表观质量的形式如下：

dv?l??l?iF?l?dt?F?l?idvim?或是dt?im? （2.102） ?式2.102中m??M?l??l?M1??,m?1??

图5 单自由度弹簧系统的运动

单自由度的弹簧系统的运动图3，当给定一个初始位移，在t=0，x=a见图5所示。像这样的运动，在每个循环中，当速度为0时，表观质量增加两次；当速度到达最大值时，其减小两次。阻尼力在每个循环中消减动能两次，在两个速度的峰值，部分质量被删减。注意：加速度在整个过程中是连续的，当质量被移除时，加速度瞬时值为0。

每个循环过程中移除的能量是通过放弃B点的动能或是：

?W?1??2?m??l?2?(i)?2??m(i)(i)??v(i)??? （2.103）

瞬间移除的平均动能为：

W1???l?2(i)?4?m(i)?m(i)??v(i)? （2.104）

每个循环中最大储存能量的损失率称为单位能耗，在FLAC3D中，单位能耗可以采用阻尼常数来表述。对一个单自由度系统，或是单模态的振动，动能的峰值与储存能量的峰值一样，因此单位能耗可由下式表达：

?W(i)m??(i)?m(i)?W?4?(i)?m?(i)?m?? （2.105）

(i)??4

临界阻尼比D可以在小的阻尼中，建立与单位能耗之间的关系：

D??W(i)W(i)?4??? （2.106） ?的默认值为0.8，所以D的默认值为0.25

公式2.105可以通过单个单元或是一列三个单元在突然施加重力荷载的情况下，系统振动到平衡来证实。两个实例都表明?WW值都很接近，说明局部阻尼式与频率无关的。然而对一个多模态同时作用下的系统，局部阻尼对每个模态都采用相同的阻尼力。 混合阻尼

式2.99中考虑的阻尼力，仅仅考虑当速度分量的符号改变时才改变。在实际，有一重要的匀速运动（与需要衰减的振动大小相比），速度可能没有0点，因此能量并不消耗。通过一个单自由度系统承受谐波振动，可以推导出一可选方程（低效率）来考虑上述问题。这个系统的不平衡力是与速度的负值成比例。例如：

v?asin??t? （2.107）

d2vdt2??a?2sin??t? （2.108） 因此dFdt是与速度v成比例的。在一个振荡系统中，

dFdt的平均值为0，即使速度v的平均值不为0。在

FLAC3D里面采用一种交替形式的阻尼，称之为混合阻尼来实现这个目的。阻尼力由速度和dFdt两项来构成，如下：

18

f??F?l???dF?l??lii??l??(i)?2?sign????dt?-sign??vi??? （2.109）

?这种阻尼适合于一振荡耗散运动中，有一个刚体运

动的系统。然而混合阻尼在能量消散方面要比基于速度的局部阻尼效率低。因此局部阻尼依然是首先的阻尼。

2.5.3 网格离散

在三维常应变率单元中，四面体有不会产生沙漏变形的优点（由节点速度产生的变形形式，没有应变率，因此没有节点力增量）。当在塑性变形的框架下，这些单元不能够提供足够的变形模式。在某些特殊情况下，当某些重要的本构方程，要求体积没有变化时，这时他们不能个别的变形。在这些情况，单元响应比理论预期值要显得过渡坚硬。在FLAC3D中，采用一种混合离散技术来处理上述问题。 混合离散技术的原理是通过适当调整四面体第一应变率张量不变量，使得单元体积具有更大的灵活性。在这个方法中，单元中最粗糙的离散是采用四面体单元来离散，单元中某个四面体的第一应变率张量不变量做为整个单元中所有四面体的体积平均值。这个方法见图6。在草图的特殊变形模式中，个体的常应变率将会产生一个体积改变，而这与不可压缩的塑性流动理论不符。在这个例子中，集合四面体（长方体）的体积仍然为常数，采用混合离散技术能使每个四面体都能反映单元这个性质，因此将会与理论预见值比较相符。

图6 混合离散技术的单元变形

在FLAC3D中，一个单元是由nt四面体组成，如图7所示nt=5。假定一个单元，假定四面体l的应变率单元式首先考虑的，然后将他们分解为偏量和静水压力项：

??l??l???l?ij??ij?3?ij （2.110）

式2.110中??l?应变率张量偏量，??l?为第一应变率张量不

变量。

??l????l?ii （2.111）

图7 8节点单元两种离散方法（每种五个四面体）

第一应变率不变量，由每个四面体对体积的平均值得到：

?z?ntk?1??k?V?k???nt? （2.112）

k?1V?k这里V?k?为四面体k的体积。因此得应变率张量的计算式

如下：

??l???zij???lij?3?ij （2.113）

当屈服发生时，膨胀本构模型将会导致平均正应力的变

化。同理，由应变率增量张量导出的第一应力不变量张量也需要在整个单元上进行体积平均。在这过程中，四面体l的应力张量也是分解为偏张量和静水压力项：

??l?s?l??l?ij?ij???ij （2.114）

s?l?为应变率张量，??l?为平均正应力。

??l??1?l?3?ii （2.115）

单元的第一应力不变量，通过在所有四面体中进行体积平均得到：

ntkk?z??k?1???V???ntV?k? （2.116） k?1由上2.116，式2.114可写为：

??l??l?ij?sij??z?ij （2.117）

在FLAC3D中，离散的过程在初始网格形成时：当单元被

19

定义时,在内部就把其离散化为四面体。例如一个8节点单元，可以离散为两种不同的5个四面体的组合体，见图7。节点应力的计算（基于应变率和应力的计算）可以使用其中一种方法，或是两种方法的一个组合。两种离散方法的优点在于，对称单元在对称荷载作用下，产生对应的响应。在这中情况下，就采用这两种方法结合的混合离散方法。节点力由这两种方法计算出来进行平均得到。 2.5.4 数值计算流程 2.5.4.1网格离散

在FLAC3D中一般由用户将模拟区域进行网格离散，每个单元将会由程序自动离散为一些列四面体单元。用户能够选择是采用一种离散方法还是采用混合离散方法。一般在应力和位移变化比较剧烈的地方推荐采用混合离散方法。程序默认采用混合离散的方法。

2.5.4.2初始条件很边界条件

求解问题的边界条件包括面力、集中荷载以及位移边界。另外体力和初始应力条件也是需要施加的。程序开始时，所有的应力和节点速度都为0，然后初始应力开始施加。集中荷载指定在表面的节点上，位移边界条件是由节点的速度来精确控制的。体力和面力在内部被转化成一些列等效的节点荷载。以上构成了数值计算的初始状态。 2.5.4.3主要的计算步骤

FLAC3D采用显示的时间历程有限差分，对每一时间步，主要的计算过程如下： 由节点速度计算新的应变率；

通过应变率，采用本构方程计算新的上一步的应力；

由应力和力，通过运动方程导出新的节点速度和位移。

每个时间步都重复上述三个步骤，在求解的过程中监测最大不平衡力。这个力要不趋于0，表示整个系统达到平衡状态，要不等于一非0常数，表明部分或是整个系统达到一个材料的稳定塑性流动。为了分析计算结果，计算过程可以随时打断。 2.5.4.4应变率计算

从一个已知的速度场，开始计算应变率，对单元中的每个四面体来说，使用以下有限差分方程来计算：

?4)l(l)ij??1l(l6V??vinj?vjni?S(l) （2.118）

l?1采用混合离散技术以后，通过式2.112和2.113计算新的对角应变率张量。对一个四面体l来说，公式如下：

z??l??l??ij??ij?3?ij （2.119）

?z为单元的第一应变率不变量平均值，计算公式如下： ?nt?k??k?z?k?1?V??nt?k? （2.120）

k?1Vn?k?t为单元中四面体的组合数，V为四面体k的体积。

2.5.4.5应力计算

在一个单元中，通过增量本构方程H*ij来计算每个四面体的应力增量。

??ij?????Cij??ij （2.121） ????????H*ijij??ij,??ij? （2.122）

4???t(l)l)ij??6V??vlinj?vl(jni?S(l) （2.123）

l?1在小变形模式中，不考虑应力修正项??Cij，在大变形中，需要修正：

??Cij???ik?kj??ik?kj??t （2.124）

?4??1(l)(l)l)ij6V??vlinj?vljni?S( （2.125）

l?1新的应力通过附加应力增量得到。采用混合离散技术后，通过式2.116和2.117计算新的对角化的应力张量。

??l?s?l?ij?ij??z?ij （2.126）

?nt?k??k??z?k?1?V?nt?k? （2.127）

k?1V2.5.4.6节点质量计算

四面体节点l的质量贡献公式如下：

ml??19Vmax??2?n(l)iS(l)??,i?1,3? （2.128） 注意所有的时间步都是统一的。全局节点质量M?l?由所

有与之共点的四面体采用式2.128求和所得。

M?l????l???m??? （2.129）

在小变形模式中，节点质量值在循环开始前只计算一次，在大变形中，每10步计算一次。

20

(2) 在程序C++文件(usermodel.cpp)中修改模型结构(UserModel::UserModel(bool bRegister): Constit -utive Model)的定义，这是一个空函数，主要功能是给(1)中定义的所有私有成员赋初值，一般均赋值为0.0。

(3) 修改const char **UserModel: roperties()函数，该函数包含了给定模型的参数名称字符串，在FLAC3D的计算命令中需要用到这些字符串进行模型参数赋值。7 A4 H8 (4) const char **UserModel::States()函数是单元在计算过程中的状态指示器，可以按照需要进行修改指示器的内容。

(5) 按照派生类中定义的模型参数变量修改double UserModel::GetProperty()和void UserModel:: SetProperty()函数，这两个函数共同完成模型参数的赋值功能。

(6) const char * UserModel::Initialize()函数在执行CYCLE命令或大应变模式下对于每个模型单元(zone)调用一次，主要执行参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值。值得注意的是，Initialize()函数调用时没有定义应变分量，但可以调用应力分量，但不能对应力进行修改。$ (7) const char * UserModel::Run()是整个模型编制过程中最主要的函数，它对每一个字单元(sub-zone)在每次循环时均进行调用，由应变增量计算得到应力增量，从而获得新的应力。在计算过程中，要根据单元应力情况对单元状态指示器进行赋值。当进行塑性模型编制时，需对达到塑性的应力状态进行修正。6 D. t* U% J7 E* ~1 D& P

(8) 修改const char * UserModel::SaveRestore()中的变量，修改方法同(2)和(5)，该函数的主要功能是对计算结果进行保存。 }# u* q9 R N (9) 程序的调试有两种方法。①在VC++的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(Additional DLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试；②在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。

8、最大、最小主应力的分布说明了什么

清楚地了解最大主应力与最小主应力的分布对于岩土工程的数值模拟十分的重要！比如：边坡开挖稳定性分析或地下厂房的开挖稳定性分析等知道了最大主应力的方向就可以避免使地下厂房的直

立边墙与其垂直，应该使其夹角尽量地小，这样才能保证高边墙和地下厂房的稳定！还有就是边坡的开挖面也尽量要与最大主应力的方向保持小的角度这样也有利于边坡的稳定，通过分析可以为设计提供依据。另一方面，反过来，对于数值分析计算评价已经施工的高边坡和地下厂房的稳定性，就要知道其内最大主应力受施工扰动或顺序影响下的方向变化来辅助判断它们的稳定性。对于最小主应力当然如楼上fengzhijunw 所说还可以判断坡体内是否存在拉应力，因为岩体或土体的抗拉强度一般都比抗压强度低的多，最容易发生受拉破坏，所以了解岩土体内的拉应力区的方向或大小十分的重要。当然判断或指导边坡或地下厂房或各种工程结构物的设计、施工或长期稳定性评价的可依据因素很多：变形速率、塑性应变区等，但是最大主应力与最小主应力的分布情况对于辅助判断仍然十分重要！！ 9、一维固结

; 单级等速大面积竖向均布加载的单层均质地基压缩问题6 S6 w, |+ K& ]2 I. ]# p. b

;第一步 平衡自重应力$ d# o8 b; t6 Q0 q new

Config fluid

gen zone brick size 1 1 10 p0 0 0 -10 p1 1 0 -10 p2 0 1 -10 p3 0 0 0

pl add sur r pl sh

model fl_iso

prop perm 1e-12 poro 0.3 ini density 9.81e2

ini fmod 2e9 ftens -1e-3 set flow off set mech on

;fix pp range z -10.1 -9.9 model elas

prop bulk 2512562.814 shear 1152959.262 dens 1335.7 fix x range x -0.1 0.1- g9 P8 p& N\ fix x range x 0.9 1.0 \ fix y range y -0.1 0.1

fix y range y 0.9 1.0 |% s/ }# i# ]* y) A0 }; F

fix z range z -10.1 -9.9! d4 Q, _. y\ ini szz 0 grad 0 0 15990.3 range z -10 0 ini syy 0 grad 0 0 12604.69 range z -10 0 ini sxx 0 grad 0 0 12604.69 range z -10 0 ini pp 0 grad 0 0 -9810 range z -10 0 ;fix pp range z -0.1 0.1 ;fix pp range z -10.1 -9.9 plot surface

set gravity 0 0 -9.81\26

solve( plot con pp save gujie0.sav ;第二步流固耦合求解 set flow on set mech on model fl_iso

prop perm 1e-12 poro 0.3 ini fdensity 981

ini fmod 2e9 ftens -1e-3

;app nstr -5e3 range z 0 x 0 1 y 0 1 ;ini pp 1e4 fix pp 0 range z 0.1 -0.1 plot surface

;set mech substep 1000 ;solve age 302400 plot con zdis app nstr -1e4 range z 0 x 0 1 y 0 1 solve age 604800 app nstr -2e4 range z 0 x 0 1 y 0 1 solve age 1209600

app nstr -3e4 range z 0 x 0 1 y 0 14 solve age 1814400 app nstr -4e4 range z 0 x 0 1 y 0 1 solve age 2419200 app nstr -5e4 range z 0 x 0 1 y 0 1' solve age 3024000

app nstr -6e4 range z 0 x 0 1 y 0 1 solve age 3628800

app nstr -7e4 range z 0 x 0 1 y 0 1 solve age 4233600\ app nstr -8e4 range z 0 x 0 1 y 0 1 solve age 4838400

app nstr -9e4 range z 0 x 0 1 y 0 10 solve age 5443200

app nstr -10e4 range z 0 x 0 1 y 0 1solve age 6048000$ X solve age 6652800: T6 solve age 7257600 solve age 78624009 solve age 8467200 solve age 90720001 solve age 9504000 solve age 9676800 solve age 10886400 solve age 12096000' solve age 13305600 solve age 14515200

solve age 15724800

solve age 181440006 {5 J4 ]* ~0 {7 T solve age 20563200U8 ?- h solve age 22982400 solve age 25401600, 时间/d 解析解/cm 本人的解/cm 解析解/本人的解4 7 0.294 0.1905 1.5433070877 O9 t5 e- 14 0.835 0.53475 1.561477326* t4 F) 21 1.526 1.022 1.4931506853 n; C* 28 2.35 1.6428 1.430484539

35 3.284 2.3881 1.3751517943 S' ^$ 42 4.316 3.2495 1.328204339

49 5.439 4.2192 1.28910694# s7 d; m7 56 6.644 5.2899 1.255978374

63 7.926 6.4548 1.227923406\ 70 9.28 7.7075 1.204022056 77 10.405 8.867 1.173452126

84 11.352 9.9483 1.141099484& X3 91 12.203 10.957 1.113717258 98 12.978 11.898 1.090771558

110 13.694 12.775 1.0719373789 v% 112 14.361 13.365 1.074523008 126 15.568 15.068 1.033182904 140 16.63 16.35 1.017125382 154 17.566 17.466 1.005725409

168 18.394 18.437 0.997667733: y* p\ 182 19.125 19.281 0.991909133% r, 210 20.342 20.654 0.984893967 238 21.294 21.693 0.98160697 266 22.036 22.479 0.980292718

294 22.617 23.074 0.9801941584 v- e;这是计算结果，开始误差很大，在50%，随着加载和时间的推移到294天是能够比较接近，请高人帮忙看一下，这个是什么原因啊？ 10、基坑桩锚支护-预应力问题

FLAC3D基坑桩锚支护中关于锚索预应力的问题

最近看了jinxiao7758的题为 《3D平面应变问题基坑桩锚支护命令流》的帖子，很受启发，对学习桩锚支护的学习很有帮助。自己尝试做了一个立体开挖的例子，基坑为长20m,宽10m,深10m。在基坑的4个角点布置4根桩，基坑长度一半处布置1根桩，共六根桩，基坑开挖1.5m后，在每根桩上加锚索，并施加预应力40KN。命令流如下： new

set log on

set logfile w.log O# j, g! a6 W' E- m$ X' D\ ;建立模型

gen zone brick size 20 10 5 p0 0 0 0 p1 100 0 0 p2 0 50

27

\0 p3 0 0 30

group soil1 range z 25 30

group soil2 range z 15 257 f7 T, F$ U6 h x0 B+ A) P7 group soil3 range z 0 15

;定义参数 m m: s prop bulk 7.0e6 shear 3.2e6 fric 25 cohesion 2e3 dil 0 ten 1e10 density 1400 range z 25 30 prop bulk 18.6e6 shear 9e6 fric 25 cohesion 1.8e3 dil 0 ten 1e10 density 1900 range z 15 25\ prop bulk 63.e6 shear 38.6eric 29 cohesion 6e3 dil 0 ten 1e10 density 2100 range z 0 15

;定义应力\ def ini_szz6 O1 W* V( H/ N/ h _grad = -1400*(-9.81) _szzl = -_grad*(30-25)

_szz0 = -_grad*302 Y4 b8 f' P% t; f\ command

ini szz add _szz0 grad 0 0 _grad range z 30 25! f* ini szz add _szzl range z 0 30

end_command8 r2 `4 U2 V! c- N( v _grad = -1900*(-9.81)

_szzl = -_grad*(25-15). `8 e7 O+ @8 `4 b% q _szz0 = -_grad*25; z; Command

ini szz add _szz0 grad 0 0 _grad range z 15 25, g/ m2 ini szz add _szzl range z 0 25' ~ end_command

_grad = -2100*(-9.81)

_szz0 = -_grad*15: I! v2 V% i( X+ i% n+ d Command! b, F, ~) E2 E\

ini szz add _szz0 grad 0 0 _grad range z 0 15% a' }- e! y. I2 `) t5 o

end_command End ini_szz

def ini_conf pnt = zone_head

loop while pnt # null val=0.43*z_szz(pnt) z_sxx(pnt) = val z_syy(pnt) = val pnt=z_next(pnt) end_loop end ini_conf

;定义边界条件# D4 B0 |1 v- fix z range z -0.1 0.1

fix x range x -0.1 0.1 fix x range x 99.9 100.1 fix y range y -0.1 0.1: fix y range y 49.9 50.15 [( M( C\ set gravity 0 0 -9.81 set mech ratio 5e-5 sol4 N' [1 k: H2 ?4 |& v, y save w1

prop tens 0.0

solve( O* Y/ \\3 r: O2 v6 t g save w2.sav9 K8 I9 q& t, u\

;定义初始条件

new. m- V( X; n( [$ T# p/ s) K; r res w2.sav ini state 0

ini xdis 0 ydis 0 zdis 0* }3 F, S% f* u; U5 s, ? ini xve 0 yve 0 zve 0* d# {. ?4 w. e. o8 x5 s his reset

;加桩8 X+ w* D8 o* s

sel pile begin 40 20 30 end 40 20 10 nseg 20

sel pile begin 40 30 30 end 40 30 10 nseg 20. U; p% |4 B1 Q& l! S; S

sel pile begin 50 20 30 end 50 20 10 nseg 20\A4 U' i- D. X' m+ J

sel pile begin 50 30 30 end 50 30 10 nseg 20

sel pile begin 60 20 30 end 60 20 10 nseg 20/ @4 k, m5 @: I2 \\\

sel pile begin 60 30 30 end 60 30 10 nseg 202 w0 d/ P' g3 F; m- g

sel pile prop dens 2.5e3 emod 28.0e9 nu 0.2 xcarea 1.76625 sel pile prop xciy 0.248378906 xciz 0.248378906; _9 C( z9 }/ `4 ^4 h! m5 N: Z, ]

sel pile prop per 4.71 cs_ncoh 10e3 cs_nfric 20 cs_scoh 10e3% c2 ^4 K& y# c3 k

sel pile prop cs_sfric 20 cs_nk 50e5 cs_sk 50e5 9 ~, z- A+ d% f- [3 }0 k( X9 Z: r set mech ratio 1e-5

his unbal l0 }6 W5 s) T# d( k9 j5 e! Y sol

save w3.sav ;2 L. r* y# B8 Y new: o

restore w3.sav- k model null range x 40 60 y 20 30 z 28.5 30 ;加第一排锚索 sel cable id=1 beg 40, 20, 29 end 40 ,15, 29 nseg 5! sel cable id=1 beg 40, 15, 29 end 40 ,10, 29 nseg 55 S sel cable id=1 pretension 40e3 range y 15 20

28

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 15 20

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 10 15/ J3 P. C; {\D' p% p$ |

sel del link range id 127) x5 b' K\

sel link id=127 127 target node tgt_num 3% p& c9 E* k) X9 ~2 x' E9 e0 I

sel link attach xdir=rigid ydir=rigid zdir=rigid &

xrdir=rigid yrdir=rigid zrdir=rigid range id 127\X# r7 @( h0 l7 o: ]! V! P3 u' u e: T C, }7 y& a ~7 K% X

sel cable id=3 beg 50, 20, 29 end 50 ,15, 29 nseg 5 4 x# h1 }# M4 H\

sel cable id=3 beg 50, 15, 29 end 50 ,10, 29 nseg 5, n' e3 P% ~3 k

sel cable id=3 pretension 40e3 range y 15 20

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 15 207 X% B$ `$ d6 j# P+ p

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 10 15 sel del link range id 139

sel link id=139 138 target node tgt_num 453 _' n. t, q% q+ o

sel link attach xdir=rigid ydir=rigid zdir=rigid &; V( B6 q. V1 {3 K& \\/ d, v3 H

xrdir=rigid yrdir=rigid zrdir=rigid range id 139

sel cable id=5 beg 60, 20, 29 end 60 ,15, 29 nseg 5 sel cable id=5 beg 60, 15, 29 end 60 ,10, 29 nseg 5

sel cable id=5 pretension 40e3 range y 15 20

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 15 20

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 10 15

sel del link range id 151/ }: l) [7 n1 r- T) [/ v; d\

sel link id=151 149 target node tgt_num 87. g( ?5 R' A# L Q

sel link attach xdir=rigid ydir=rigid zdir=rigid & xrdir=rigid yrdir=rigid zrdir=rigid range id 151 ' X& Z7 d: E$ `, x6 u

sel cable id=2 beg 40, 30, 29 end 40 ,35, 29 nseg 5

sel cable id=2 beg 40, 35, 29 end 40 ,40, 29 nseg 5 sel cable id=2 pretension 40e3 range y 30 35

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 30 35 sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 35 40) n7 x\ sel del link range id 1633 \\- ^$ \\; X$ `8 d. m/ B\ sel link id=163 160 target node tgt_num 24 sel link attach xdir=rigid ydir=rigid zdir=rigid & xrdir=rigid yrdir=rigid zrdir=rigid range id 163

sel cable id=4 beg 50, 30, 29 end 50 ,35, 29 nseg 5 ' Z- g! @. _1 j$ h( u5 |+ u

sel cable id=4 beg 50, 35, 29 end 50 ,40, 29 nseg 5$ z& N. } e* z& u% w

sel cable id=4 pretension 40e3 range y 30 35

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 30 35

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 35 40# K7 H( o/ i/ B% d% K$ Y4 P sel del link range id 175

sel link id=175 171 target node tgt_num 66' G& ?5 F+ q( G- a0 [

sel link attach xdir=rigid ydir=rigid zdir=rigid & xrdir=rigid yrdir=rigid zrdir=rigid range id 175

sel cable id=6 beg 60, 30, 29 end 60 ,35, 29 nseg 57 A5 Q! Q9 y! ~+ N$ C

sel cable id=6 beg 60, 35, 29 end 60 ,40, 29 nseg 5 sel cable id=6 pretension 40e3 range y 30 35

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 0 gr_k 0 gr_fric 0.0 range y 30 35: }& g& x( H; H: `* t S

sel cable prop emod 2.05e11 ycomp 1.0e5 ytens 15.34e5 xcarea 5.54e-3 gr_coh 1.5e10 gr_k 0.56e9 gr_fric 25.0 gr_per 0.264 range y 35 40 sel del link range id 187

sel link id=187 182 target node tgt_num 108# d* i) |: T( r\ sel link attach xdir=rigid ydir=rigid zdir=rigid & xrdir=rigid yrdir=rigid zrdir=rigid range id 187 plot sel pile fo fz% s1 N; a U0 M- g0 d3 L sol( h& |0 `+ O6 C7 H$ W& k5 J print sel cable force

print sel cable stress8 I9 g2 O) ^' Y2 S print history

save w4.sav) A1 b2 l( Q9 ^8 r

上述命令流模拟的是开挖1.5米后，施加锚索并施加预应

29

力的过程，由于支护桩数量很少，继续开挖土体发生破坏，计算无法收敛，故没有继续开挖。针对上述命令流小弟有几个问题请教：

1.上述例子中桩锚之间的连接是单元与单元之间的连接，我只会一根一根的进行连接，那位高手能教教小弟怎样用fish命令，寻找桩锚之间的节点，并进行联结。

2.通过set mech ratio 1e-5 语句对计算进行控制，FLAC中默认的为多少，一般计算设成多少比较合适。$ \\; C0 p1 d9 t: i# ^+ x/ p

3.按上述命令流计算后，print sel cable force, print sel pile force,其中锚索自由段受力最大为6.538KN，但施加的预应力为40KN，锚索预应力施加后进行平衡，预应力会有所降低，但会相差这么大吗？《3D平面应变问题基坑桩锚支护命令流》中也有这个问题，施加预应力后40KN变为3.361KN，不知道是什么原因，希望各位指点。

11、FLAC3D模拟金属矿采矿与充填应用注意事项 FLAC3D在模拟金属矿山采矿中的优势很突出，首先因为金属矿是一种块体单元的组合，且渗流，蠕变及流变等性质可以不考虑，所以操作起来相对比较简单。但是想进行很好的模拟也要注意许多问题。根据本人几年来对金属矿采矿和充填的数值模拟的实践，特写下在计算模拟中应注意的一些问题。 1. 金属矿山建模相对比较简单，但是当金属矿断层比较多时，应根据具体的情况适当的进行简化和忽略，有些人为了模拟的完全符合，结果把所有的东西都堆进去，那样使模型变得不协调，计算的结果可能会发生比较大的偏差。故合理简化模型，使其称为互相协调的且能相对较大程度反映主要问题的模型为最佳模型。- Z* R4 v+ C/ U6 t\ 2. 模型构建完毕之后，进行边界的确定和荷载的加载，此时也要根据你分析问题的区域，首先把模型先在草稿纸上进行简化分析，看哪些地方需要固定，哪些面是可以左右移动的，然后根据你计算的需要进行加载。有些模型可以直接进行重力加载而忽略了地应力的因素，有些模型以地球构造应力的分布为主。所以应对其进行认真的区分，看你所构建的模型属于哪种情况。

3. 大家都知道，加载完毕之后，开挖（采矿）之前，应对模型进行初始平衡，在初始平衡的时候，也要注意，把初始平衡的点设置的的很小，那样会无止境的计算平衡下去，一般根据实际情况平衡到10e-3就可以了，也可以根据经验设置步数来平衡，使其尽快完成。/ x' |, X$ d2 V$ ]

4. 在采矿过程中，应根据实际的力学原理，分析哪种开采顺序是比较合理的，然后有选择的进行比

较几种可行的采矿顺序的计算情况。千万不宜随意的设置很多开采顺序进行，然后看应力变化云图等来分析哪种更好，这种时候往往是分析不出来的，因为有些顺序的变化，给你的分析变得没有标准。所以首先应从理论上确定一个比较优的开采顺序，然后基于此顺序进行适当调整，来改变其变化，看其应力及屈服情况，才能较好的进行比较，得出最佳方案。故基准或标准开采顺序的确定相当必要。. 5. 开采完毕进行充填时应注意，充填分为以下几种，在金属矿山常用的是胶结充填和非胶结充填。胶结充填分为废石胶结充填，尾砂胶结充填，膏体充填等，非胶结充填分为尾砂充填和块石充填。所以在进行模拟充填时，应根据不同的充填类型进行分别设计模拟，采用尾砂胶结充填和废石胶结充填及膏体充填可以考虑采用块体进行模拟，即假设其为单元块体而不是流体，测出其粘结力，摩擦角及内聚力，弹性模量等，对其空区进行充填。若是非胶结充填，最好用PFC来模拟其充填，即采用颗粒流法模拟。 12、FLAC3D塑性区的讨论

FLAC3D塑性区是否应该包括-p部分，看了一些解释，好像不统一，-p的是历程中出现了塑性，那么现在是否还塑性状态？绝大部分人应该认为-p的是现在不处于塑性状态，那么塑性区就不应该包括-p的，还有个问题，塑性区的塑性应变应该是按什么标准？好像按剪切塑性应变出的图与ABAQUS、ANSYS的按塑性应变出的图不太一样。希望大家讨论

在FLAC3D中采用塑性本构模型进行模拟时，通常要显示应力符合屈服准则的区域，即塑性区。塑性区有剪切破坏塑性区（shear failure）和拉伸破坏（tension failure）,其中每一个又有两种状态，即应力正位于屈服面上（-n，new）和曾进入过屈服状态，而现在已经退出（-p，past）。

30

工程实例

1、FLAC3D 锚杆建模中出现的问题

请教各位大虾：

我用FLAC3d建模计算锚杆的应力和位移，发现几个问题不知如何处理？

1. 计算出来的位移值偏大，达到10cm，而施加

荷载还未达到设计值。

2. 锚杆周边的位移和应力云图好像不太正常。 new Title

gen zone brick p0 0 0 0 p1 6 0 0 p2 0 0 -12 p3 0 6 0 size 15 30 15 ratio 1.0 1.0 1.0

gen zone reflect normal 0 -1 0 origin 0 0 0 gen zone reflect normal -1 0 0 origin 0 0 0 model mohr

prop bulk 9.52e6 shear 4.65e6 coh 0 fric 25 ini dens 2000

fix z range z -11.9 -12.1 fix y range y 5.9 6.1 fix y range y -6.1 -5.9 fix x range x 5.9 6.1 fix x range x -6.1 -5.9 set grav 0 0 -10 set large hist unbal

sel cable id 1 begin 0.037 0 0.9 end 0.037 0 -3.6 nseg 15

sel cable id 2 begin -0.037 0 0.9 end -0.037 0 -3.6 nseg 15

sel cable prop xcarea 0.0167 emod 200e9 ytens 335e6 gr_k 16.7e6 gr_per 0.459 gr_coh 200

sel cable id 1 pre 5e4 range x 0.0369 0.0371 y -0.001 0.001 z 0.69 0.71

sel cable id 2 pre 5e4 range x -0.0371 -0.0369 y -0.001 0.001 z 0.69 0.71

plot set rotation 30 0 40 plot set center auto

plot add contour disp outline on plot add sketch plot add disp

plot add sel cable force red yellow plot add sel cable grout stress blue plot add sel cable grout slip green plot show step 2000

2、锚杆受拉超限

锚杆能够承受拉力，而不能承受超过规定的剪应力，岩体的变形超过了锚杆抗剪范围，会出现上述警告。 锚杆的端头是否经过预应力锚杆处理的，关于预应力锚杆的模拟，你可以到仿真论坛 lakerwater的帖子看看。 ;计算中锚杆的长度发生了变化!!!!!!!!!!!!!

;如果不施加初锚力，锚杆的长度不会变，但是一旦施加了初锚力，锚杆的长度

;就会发生变化，能直接从图中看出，并且越靠近隧道口（y=0附近），锚杆长度

;变化越明显。隧道尾（y=20附近）处锚杆的长度变化不明显。

;同时，当初锚力达到6t时，锚杆变化异常，出现了很多警告。如下：

;+++ WARNING: axial dir. of cables and piles or normal dir. of liners

;and geogrids at shared node out of alignment by more than 30 degrees.

;这是怎么回事？请大侠赐教。。。。。。。

rest Model.sav ;施加2t的初锚力 set log on

set logfile 锚杆2t.log ;材料模式：摩尔-库伦 m m

;========围岩参数========= ;加载层(细砂岩)

pro b 10.9e8 s 8.7e8 f 40 c .93e5 t .93e4 ran z 47.12 70 ini d 2500 ran z 47.12 70 ;-----------------------

;老顶(中粒砂岩)

pro b 12.1e9 s 10.8e9 f 44 c 1.18e6 t 1.18e5 ran z 36.91 47.12

ini d 2500 ran z 36.91 47.12 ;----------------------- ;直接顶(细砂岩)

pro b 10.9e8 s 8.7e8 f 40 c .93e5 t .93e4 ran z 33.32 36.91 ini d 2500 ran z 33.32 36.91 ;----------------------- ;伪顶(碳质泥岩)

pro b 7.1e8 s 4.7e8 f 33 c .53e5 t .53e4 ran z 32.8 33.32 ini d 2500 ran z 32.8 33.32 ;----------------------- ;煤层(二号煤层)

pro b 7.1e8 s 4.7e8 f 33 c .53e5 t .53e4 ran z 30 32.8 ini d 2300 ran z 30 32.8

21

;-----------------------

;直接底(泥岩)

pro b 9.3e8 s 7.1e8 f 37 c .77e5 t .77e4 ran z 25 30 ini d 2500 ran z 25 30 ;-----------------------

;老底(细砂岩)

pro b 11.8e9 s 10.1e9 f 43 c 1.09e6 t 1.09e5 ran z 0 25 ini d 2500 ran z 0 25

;=============================================

;固定边界

fix x ran x -0.1 0.1 fix x ran x 74.9 75.1 fix y ran y -0.1 0.1 fix y ran y 19.9 20.1 fix z ran z -0.1 0.1

;=============================================

;原岩应力场

ini sxx -25.75e6 grad 0 0 2.5e4 ran z 0 70 ini syy -25.75e6 grad 0 0 2.5e4 ran z 0 70 ini szz -25.75e6 grad 0 0 2.5e4 ran z 0 70 app szz -24e6 ran z 70

;======================================== ;计算设置

set gravity 0 0 -10 set large

;====================================== ;显示垂直应力场分布 pl bcon szz out on

;====================================== ;原岩应力场平衡计算 solve

;====================================== ;调整生成原岩应力场后的节点速度和位移为0 ini xvel 0.0 yvel 0.0 zvel 0.0 ini xd 0.0 yd 0.0 zd 0.0

;====================================== ;开挖巷道

m n ran gro 10 an gr 11 an gr 12 an gr 13 an

;=====定义锚杆，长度2m=================== def inbolts

xtable(1,1)=39.7524 xtable(1,2)=40

xtable(1,3)=36.9437 xtable(1,4)=38.0563 xtable(1,5)=39.0587 xtable(1,6)=35.9413

xtable(1,7)=35

xtable(1,8)=35.2476 xtable(2,1)=41.5544 xtable(2,2)=42

xtable(2,3)=36.4987 xtable(2,4)=38.5013 xtable(2,5)=40.3057 xtable(2,6)=34.6943 xtable(2,7)=33

xtable(2,8)=33.4456

ytable(1,1)=28.5847 ytable(1,2)=27.5 ytable(1,3)=29.9373 ytable(1,4)=29.9373 ytable(1,5)=29.4546 ytable(1,6)=29.4546 ytable(1,7)=27.5 ytable(1,8)=28.5847 ytable(2,1)=29.4525 ytable(2,2)=27.5 ytable(2,3)=31.8872 ytable(2,4)=31.8872 ytable(2,5)=31.0182 ytable(2,6)=31.0182 ytable(2,7)=27.5 ytable(2,8)=29.4525

nn=0 ;锚杆单元序列号 ;隧道锚杆 loop ii (1,19)

loop jj (1,8) nn=nn+1

bx=xtable(1,jj) bz=ytable(1,jj) ex=xtable(2,jj) ez=ytable(2,jj) command

sel cable id nn beg bx ii bz end ex ii ez nseg 4

sel cable pre 60000 prop emod 45e9 xcarea 3.8e-4 gr_per 1 yten 2.5e6 gr_k 17.5e7 gr_c 2e6 endcommand endloop endloop end inbolts

;======================================== ;显示地层移动分布 pl con zdisp out on

22

pl ad sel geo scale 0.005

pl set center 36.7 7.814 28.61 pl set rotation 20 0 20 pl set mag 9.31 pl set dist 217.4 print sel ca le

;======================================== ;巷道开挖平衡计算 solve

print sel ca le sav 4-6.sav

3、FLAC3D计算完毕关机

4、隧道模拟

参数：;问题描述 考虑圣维南原理,取周围岩土的尺寸维隧道尺寸的5~6倍,此处取为15m r;初衬 C50管片外径3.0m，内径2.7m，C50,E＝34.5GPa，V＝0.167,.g;二衬 外径2.7m，内径2.3m,C30,E＝30GPa，V＝0.167（二衬先不考虑）& ;粘质粉土，4.8m，φ＝26°，重度＝19，c＝11KPa，变形模量＝10MPa

;中细纱，1.7m，φ＝28°，重度＝20，c＝0KPa，变形模量＝13Mpa) ;问题描述 考虑圣维南原理,取周围岩土的尺寸维隧道尺寸的5~6倍,此处取为15m

;初衬 C50管片外径3.0m，内径2.7m，C50,E＝34.5GPa，V＝0.167,

;二衬 外径2.7m，内径2.3m,C30,E＝30GPa，V＝0.167

;粘质粉土，4.8m，φ＝26°，重度＝19，c＝11KPa，变形模量＝10MPa+ |& c. g% T P w

;中细纱，1.7m，φ＝28°，重度＝20，c＝0KPa，变形模量＝13Mpa

;-------------------------------前处理 title '沉降计算'

plot set title text 'chenjiang'8 d0 v0 b* o% t plot set rotation 20 0 30 5 e/ \\: b' o. u6 |: Y plot set center 0 0 5 plot set dist 300 plot set mag 0.8

plot add surface green

plot add axes0 @0 G( b c- Q7 D, M, s2 f5 _ plot show

;建立模型---------------------------------4 d gen zone radcyl p0 0 0 0 p1 4.7 0 0 p2 0 50 0 p3 0 0 4.7 & size 10 30 10 4 dim 3 3 3 3 group tunnel gen zone brick p0 4.7 0 0 p1 18 0 0 p2 4.7 50 0 p3 4.7 0 4.7 &size 10 30 5 group brick1! gen zone cshell p0 0 0 0 p1 3 0 0 p2 0 50 0 p3 0 0 3 & dim 2.7 2.7 2.7 2.7 size 1 30 10 6 group segment gen zone cylinder p0 0 0 0 p1 2.7 0 0 p2 0 50 0 p3 0 0 2.7 &*size 6 30 10 group cy, s6 X' K; Z' k/ r( E

gen zone reflect ori 0 0 0 norm 0 0 -1 range z 0 4.73 gen zone brick p0 0 0 4.7 p1 4.7 0 4.7 p2 0 50 4.7 p3 0 0 9.5 &)size 5 30 5 group brick2) X gen zone brick p0 4.7 0 4.7 p1 18 0 4.7 p2 4.7 50 4.7 p3 4.7 0 9.5 &#size 10 30 5 group brick3. gen zone brick p0 0 0 -18 p1 4.7 0 -18 p2 0 50 -18 p3 0 0 -4.7 &$size 5 30 15 group brick4gen zone brick p0 4.7 0 -18 p1 18 0 -18 p2 4.7 50 -18 p3 4.7 0 -4.7 & size 10 30 15 group brick5# gen zon reflect norm 1 0 0 orig 0 0 0

;-赋土参数-------------------------------------------4 model mohr ;E=13Mpa,v=0.3

prop density 2000 bulk 13e6 shear 5eric 28 coh 1e10 tension 1e10 range z -18 4.7

;E=10Mpa,v=0.255 E, s% d) u ?# J( _+ B

prop density 1900 bulk 10e6 shear 4eric 26 coh 1e10 tension 1e10 range z 4.7 9.5$ X ;初始地应力--------------------------------------! ini xdisp 0 ydisp 0 zdisp 0

ini szz -1.85e5 grad 0 0 1.9e4 range z 4.7 9.5 ;19×4.8+20×4.7＝185.2KPa) |' M8 F' c+ E+ d ini sxx -6.67e4 grad 0 0 6.84e3 range z 4.7 9.5 ;侧压力系数取为0.36 （66.67KPa）

ini syy -6.67e4 grad 0 0 6.84e3 range z 4.7 9.5* k) ~8 `0 q% g ini szz -1.85e5 grad 0 0 2.0e4 range z -18 4.7

23

ini sxx -6.67e4 grad 0 0 7.2e3 range z -18 4.7 ini syy -6.67e4 grad 0 0 7.2e3 range z -18 4.7 ;边界条件-------------------------------------

apply nstress -1e4 range z 9.4 9.6* B( h+ C. o! m5 @! D

fix x range x -18.1 -17.9

fix x range x 17.9 18.1/ s1 W( J0 Y& U2 I% c& j fix y range y 49.9 50.1

fix y range y -0.1 0.1# F$ O8 X( Z. V- e1 Y0 n fix z range z -18.1 -17.9 set grav 0 0 -10 small

;solve. u7 |. @$ O9 ^7 r' K0 H' d

;save small.sav0 \\$ R+ S3 l, W$ u8 f\ - |% ]' n# @; `; Q! y ;开挖 求解

model null range group cy/ b r;---------------------------! ini xdisp 0 ydisp 0 zdisp 0 ;---------------------------

model elas range group segment

prop density 2500 bulk 34.5e9 shear 14.78e9 range group segment

;E=13Mpa,v=0.3* n$ {5 D9 j$ q m3 O7 ?+ u$ J: w3 prop density 2000 bulk 13e6 shear 5eric 28 coh 0 range group brick1

prop density 2000 bulk 13e6 shear 5eric 28 coh 0 range group tunnel/ b$ z' m/ V: d6 M& A' t! S

prop density 2000 bulk 13e6 shear 5eric 28 coh 0 range group brick4

prop density 2000 bulk 13e6 shear 5eric 28 coh 0 range group brick5

;E=10Mpa,v=0.251 j; W9 @& N7 o ~+ v3 G

prop density 2000 bulk 13e6 shear 5eric 28 coh 0 range z 4.7 9.5\hist unbal

hist gp xdis 3 0 0& o, B: N! @\ hist gp xdis -3 0 0 hist gp zdis 0 0 9.5

hist gp zdis 0 0 7.51 y1 N) h: N( `' u hist gp zdis 0 0 5.5. t\ hist gp zdis 0 0 4.5 hist gp zdis 0 0 3

hist gp zdis 0 0 -3* Z8 t, O' u5 Z1 @6 z set large- Q1 ?8 p. w7 E) L\ ;solve ;save n.sav

5、显示最大位移点的坐标参考

def get_gp_maxdisp3 R2 R) n) M( Q# f! D$ z

gp0_disp = gp_xdisp(gp_head)*gp_xdisp(gp_head) gp0_disp = gp0_disp + gp_ydisp(gp_head)*gp_ydisp(gp_head) gp0_disp = gp0_disp + gp_zdisp(gp_head)*gp_zdisp(gp_head) gp0_disp = sqrt(gp0_disp) # S/ N+ h( @# C( E& }0 p- O% p_gp=gp_head: ^* E @6 e6 Z$ r ;找最大值* H/ ~4 `. m; ?/ E

loop while p_gp # null8 ]( v9 N* \\9 K) L l; X2 M gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)

gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp) gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp) gp_disp = sqrt(gp_disp) : U* q! w+ w& k5 ?/ g if gp_disp>gp0_disp gp0_disp=gp_disp endif ' Y' r( ` p_gp = gp_next(p_gp)% Z3 i4 b4 t$ H endloop ;找最大值的坐标 p_gp=gp_head ss=0+ ]1 loop while p_gp # null

gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)) I7 D Z0 z gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp) gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp) gp_disp = sqrt(gp_disp)

if gp_disp=gp0_disp. u5 [( r$ H/ i' m ss=ss+1

xtable(1,ss)=gp_xpos(p_gp) ytable(1,ss)=gp_ypos(p_gp) endif

p_gp = gp_next(p_gp) endloop end

get_gp_maxdisp

6、固结问题输出某一节点随时间的孔压变化

方法其实很简单，如果不需要导出数据，直接hist gp pp就可以，如果要导出数据的话就需要使用fish。例如你要记录每隔1hour即3600s后某一节点的孔压，首先你要用fish来控制solve的时间间隔，然后把每个小时后的孔压记录到一个table中，最后plot table就可以了，要把数据导出时，现log on，然后print table(),然后log off即可，具体操作如下,假设求解10小时，所记录节点坐标为0，0，0：

def pplog

sage=3600 ;渗流时间 % p, Y$ C& F4 `- [9 i. A\ pnt=gp_near(0,0,0) n=1

loop while sage <= 36000 ;计算时间以10小时

24

计 , |: command, ` solve age sage\ endcommand

gridpp=gp_pp(pnt) 9 f2 _5 s\x6 L

xtable(1,n) = sage ytable(1,n) = gridpp sage = sage + 3600 n = n + 1

end_loop( t4 b6 |: H0 c4 _+ E1 a end; B6 B) e) P% m5 O% f U pplog

7、初期支护中的“钢拱架”建立

1、钢拱架一般都采用多个beamsel小段来模拟（调整形状），同一榀钢拱架只要保证组成它的beamsel都属同一id就可以了。这一部分好解决。+ 2、设置每个node与网格的关系。我个人认为：实际的钢拱架都与围岩密贴接触（加垫块调整），这和shell模拟的喷混凝土与围岩的接触关系有异曲同工之处，区别仅在于网格不能带动node（钢拱架）旋转。因此node与zone之间的Link就应该只保留x,y,z，而将link中的xrot,yrot,zrot删除。以上是我对钢拱架如何建立的一些浅显想法，有什么不足之处还请各位版主和大牛们提出中肯意见；以起到抛砖引玉的作用，大家共同交流经验，一起进步。 3、如果钢拱架安装之前，已经有shell（喷混凝土）附在网格上。那么在生成beamsel时，两个端点node可以与shell的node设在同一位置（两个node，分属于shell和beam），然后删除node（beam）与zone的Link，最后建立beam node与shell node 的联结（仅是x,y,z,不包括xrot,yrot,zrot），这保证了beam和shell的实际联结情况——密贴但不会随之转动。 你看看下面的程序吧，应该会对你有所启发。 按照这个思路我试着建立了模型，效果还是不错的.我把link的fish处理程序（一个处理Link的小思路）贴出来吧，大家分享一下，可能有的地方不是很合理，欢迎斧正。 restore 111.sav def parmater ;起始位置 z_=0

;榀数（纵向） p=5+ k ;半径

Rbeam=5.85000

Rshell=6.0084(洞室半径)

;首个node的角度位置, _6 l/ M ?( U\

a=0

det(a)=11.0589 ;beam半径稍小于shell（洞室）半径 def operate_link link_id=1000

z=z_\

loop m(1,p)6 x7 I/ d( Q\ loop n(0,13)

x_shellnode=Rshell*cos((a+n*det(a))*degrad)6 T0 y_shellnode=Rshell*sin((a+n*det(a))*degrad)) f8 z_shellnode=z( Q6 W. }9 Q ]0 o1 ^$ u7 n

x_beamnode=Rbeam*cos((a+n*det(a))*degrad) y_beamnode=Rbeam*sin((a+n*det(a))*degrad)8 z_beamnode=z# sn_pointer=nd_near(x_shellnode,y_shellnode,z_shellnode) bn_pointer=nd_near(x_beamnode,y_beamnode,z_beamnode) shellnode_id=nd_id(sn_pointer) beamnode_id=nd_id(bn_pointer) command

sel set link node_tol 0.3 sel

link id link_id beamnode_id target node tgt_num shellnode_id sel link attach xdir=rigid ydir=rigid zdir=rigid range id link_id sel link attach xrdir=free yrdir=free zrdir=free range id link_id6 end_command

link_id=link_id+1 end_loop end_loop end

operate_link

关于user define model

FLAC3D的二次开发环境提供了开放的用户接口，在软件安装文件中包含了软件自带所有本构模型的源代码，且给出了Mohr-Coulomb模型和应变软化模型的编译示例，因此可以方便地进行本构模型的修改与开发。为了方便起见，下面的说明以建立UserModel模型为例。3 k& b. y& (1) 在模型头文件(usermodel.h)中进行新的本构模型派生类的声明，修改模型的ID(为避免与已有模型冲突，一般要求大于100)、名称和版本，修改派生类的私有成员，主要包括模型的基本参数及程序执行过程中主要的中间变量。

25

本文来源：https://www.bwwdw.com/article/0kdd.html