; -----------------------------+; Author: qiaoch@yeah.net ; 2012/2/20; -----------------------------+newdef SetVariables1 array cx0(5),cz0(5),cx1(5),cz1(5) ;store inner and outer x z coordinates of bolts _ySection=0.8 _yspace = 0.8 ; space between rockbolt along y axes _cableSpace = 2.4 ; space between cables _boltRad = 2.5 ;length of rockbolt _cableRad = 6.5 ;length of cable _cabEMod = 200e9 ;Young Model of cable and rockbolt _excaLength = 3*_yspace ;开挖进尺 _beganOnY = 0.0 ; start point on Y axes for excavation _cableSteps = 100 ;steps after cable installed _tunRad = 2.05 ;巷道上半拱半径 _ntunRad = (-1)*_tunRad _sideWHeight = 1.5 _sideWallz = (-1)*_sideWHeight _boltAng = 2*asin(_yspace/2/_tunRad) _cableAng = 2*asin(_cableSpace/2/_tunRad) _rightWall0 = _tunRad - 0.01 ;for liner install _rightWall1 = _tunRad + 0.01 _leftWall0 = (-1)*_tunRad + 0.01 _leftWall1 = (-1)*_tunRad - 0.01endSetVariables1;--------------------------------------------------------------------def SetVariables2 _endony = _beganOnY + _excaLength _beganOnY_01 = _beganOnY - 0.01 _beganOnY_10 = _beganOnY + 0.01 _endony_01 = _endony - 0.01 _endony_10 = _endony + 0.01 _ySection2 = _ySection + 1*_yspace _ySection3 = _ySection + 2*_yspaceendSetVariables2;--------------------------------------------------------------------def RightCoords loop n (1,4) cx0(n) = _tunRad*cos((pi/2)-_boltAng*float(n)) cz0(n) = _tunRad*sin((pi/2)-_boltAng*float(n)) cx1(n) = (_tunRad + _boltRad)*cos((pi/2)-_boltAng*float(n)) cz1(n) = (_tunRad + _boltRad)*sin((pi/2)-_boltAng*float(n)) end_loopend;--------------------------------------------------------------------def LeftCoords loop m (1,4) cx0(m) = _tunRad*cos((pi/2)+_boltAng*float(m)) cz0(m) = _tunRad*sin((pi/2)+_boltAng*float(m)) cx1(m) = (_tunRad + _boltRad)*cos((pi/2)+_boltAng*float(m)) cz1(m) = (_tunRad + _boltRad)*sin((pi/2)+_boltAng*float(m)) end_loopend;---------------------------定位拱顶节点-----------------------------def FindNode1 ;---begin node of cable--- _x0=cx0(n) _y0=_ySection _z0=cz0(n) ;---end node of cable----- _x1=cx1(n) _y1=_y0 _z1=cz1(n)end;---------------------------定位右竖墙节点---------------------------def FindNode2 ;---begin node of cable--- _x0=cx0(n) _y0=_ySection _z0=cz0(n) ;---end node of cable----- if n = 1 then _x1 = _x0 + _boltRad _y1 = _y0 _z1 = _z0 else _x1 = _tunRad + _boltRad*cos(60*degrad) _y1 = _y0 _z1 = _z0 - _boltRad*sin(60*degrad) endifend;---------------------------定位左竖墙节点----------------------------def FindNode3 ;---begin node of cable--- _x0=cx0(n) _y0=_ySection _z0=cz0(n) ;---end node of cable----- if n = 1 then _x1 = _x0 - _boltRad _y1 = _y0 _z1 = _z0 else _x1 = _ntunRad - _boltRad*cos(60*degrad) _y1 = _y0 _z1 = _z0 - _boltRad*sin(60*degrad) endifend;------------------------定位锚索起节点坐标-----------------def FindCableNode ;begin node of cable _x0=cx0(n) _y0=_ySection _z0=cz0(n) ;end node of cable _x1=cx1(n) _y1=_y0 _z1=cz1(n)end;---------------------------安装锚杆----------------------------def InstallBolt command sel cable begin _x0 _y0 _z0 end _x1 _y1 _z1 nseg=5 end_commandend;---------------------------安装锚索----------------------------def InstallCable command sel cable begin _x0 _y0 _z0 end _x1 _y1 _z1 nseg=13 end_commandend;---------------------------rockbolt----------------------------def SetBoltProps command sel cable pro emod _cabEMod xcarea 3.8e-4 gr_per 0.25 & yten 25e4 gr_k 3.8e7 gr_c 4e5 sel cable pretension 8e4 end_commandend;------------------------ cable -----------------------------------def SetCableProps command sel cable pro emod _cabEMod xcarea 2.6e-4 gr_per 0.22 & yten 35e4 gr_k 3.6e7 gr_c 3.5e5 sel cable pretension 8e4 end_commandend;--------------------------------------------------------------------def rWallCables ;--- install the rockbolts on the straight wall of right --- cx0(1) = _tunRad cz0(1) = cz0(5) - 0.7 cx1(1) = _tunRad + _boltRad cz1(1) = cz0(1) cx0(2) = _tunRad cz0(2) = cz0(1) - 0.65 loop n (1,2) FindNode2 InstallBolt SetBoltProps end_loopend;-----------------------------------------------------def lWallCables ;install the cable on the straight wall of left cx0(1) = _ntunRad cz0(1) = cz0(5) - 0.7 cx1(1) = _ntunRad - _boltRad cz1(1) = cz0(1) cx0(2) = _ntunRad cz0(2) = cz0(1) - 0.65 loop n (1,2) FindNode3 InstallBolt SetBoltProps end_loopend;-------------------------------------------------------def InstallBolts ;install right top cables RightCoords loop n (1,4) FindNode1 InstallBolt SetBoltProps end_loop ;install left top cables LeftCoords loop n (1,4) FindNode1 InstallBolt SetBoltProps end_loop ;install the cable on the top of arch _ndd=nd_near(0,_ySection,_tunRad) _nddid=nd_id(_ndd) _x0=nd_pos(_ndd,2,1) _y0=nd_pos(_ndd,2,2) _z0=nd_pos(_ndd,2,3) ;end node of cable _x1=0 _y1=_ySection _z1=_tunRad + _boltRad ;--------------------- InstallBolt SetBoltPropsend;---------------------------------;计算锚索端点坐标def cableCoords cx0(1) = _tunRad*sin(_cableAng) cz0(1) = _tunRad*cos(_cableAng) cx1(1) = (_tunRad + _cableRad)*sin(_cableAng) cz1(1) = (_tunRad + _cableRad)*cos(_cableAng) cx0(2) = 0.0 cz0(2) = _tunRad cx1(2) = 0.0 cz1(2) = _tunRad + _cableRad cx0(3) = (_ntunRad)*sin(_cableAng) cz0(3) = _tunRad*cos(_cableAng) cx1(3) = -1*(_tunRad + _cableRad)*sin(_cableAng) cz1(3) = (_tunRad + _cableRad)*cos(_cableAng)end;--------------------------------------------------;-------------------安装锚索-----------------------def InstallCables cableCoords loop n (1,3) FindNode1 InstallCable SetCableProps end_loopend;---------------------------------------------------def Excavatingcommand mod null range y _beganOnY _endony group upexcav mod null range y _beganOnY _endony group lowexcav ; sel shell id=1 range cyl end1 0 _beganOnY 0 end2 0 _endony 0 rad _tunRad z 0 _tunRad sel shell id=1 range x _rightWall0 _rightWall1 y _beganOnY _endony z _sideWallz 0.1 ; liner on right straight wall sel shell id=1 range x _leftWall0 _leftWall1 y _beganOnY _endony z _sideWallz 0.1 ; liner on left straight wall sel delete shell range cyl end1 0 _endony_01 0 end2 0 _endony_10 0 rad _tunRad ; sel node fix y xr zr range y _beganOnY_01 _beganOnY_10 sel node fix y xr zr range y _endony_01 _endony_10 ;--- Assign the property of Liner --- sel shell id=1 prop iso=(2e9, 0.25) thick=0.3end_commandend;------------------------------------------------------------def InstallBoltOfOneSection if _beganOnY = 0 then command print _beganOnY end_command else _ySection = _beganOnY InstallBolts rWallCables lWallCables endif _ySection = _beganOnY + _yspace InstallBolts rWallCables lWallCables _ySection = _beganOnY + 2*_yspace InstallBolts rWallCables lWallCablesend;------------------------------------------------------------def InstallCableOfOneSection _ySection = _beganOnY + _yspace/2 InstallCablesend;------------------------------------------------------------; Geo model ;------------------------------------------------------------gen zone radcylinder dim _tunRad _tunRad _tunRad _tunRad p0 0 0 0 p1 30 0 0 p2 0 12.0 0 p3 0 0 30 & ratio 1.0 1.0 1.0 1.1 size 6 30 16 20 fill group uptunnel group upexcav range cyl end1 0 0 0 end2 0 12.0 0 rad _tunRad group uptunnel group uptunnel range group upexcav not;------------------------------------------------------------gen zone radtunnel dim _sideWHeight _tunRad _sideWHeight _tunRad p0 0 0 0 p1 0 0 -30 p2 0 12.0 0 & p3 30 0 0 ratio 1.0 1.0 1.0 1.1 size 7 30 6 20 fill group lowtunnel range group uptunnel not group upexcav notgroup lowexcav range x 0 _tunRad z -1.7 0 group lowtunnel;------------------------------------------------------------gen zone reflect norm -1,0,0 ori 0,0,0;------ Assign model property -----------model elpro bulk 10e10 she 10e10 ;------------------------------------------------def DriveInitial _sxxval = _kh * _szzval _sxgrad = _sxxval / (-960) _szgrad = _szzval / (-960) _szztop = (930 * _szzval)/960endset _szzval = -24.5e6 set _kh = 1.0DriveInitial;--- Boundary Condition ---apply szz _szztop range z 29.99 30.01fix z range z -30.01 -29.99 fix x range x -30.01 -29.99fix x range x 29.99 30.01fix y range y -0.01 0.01fix y range y 11.9 12.1 ;//attention!;--- initial condition ---set grav 0 0 -10ini density 2500;-------------------------ini szz _szzval grad 0 0 _szgrad ini sxx _sxxval grad 0 0 _sxgrad ini syy _sxxval grad 0 0 _sxgrad ;-------------------------plot create _view01plot add surfaceplot set mag 5plot show;-------------------------solve; save initialstress.sav;-----------------------------------------------------ini xdisp 0 ydisp 0 zdisp 0;-----------------------------------------------------def derive s_mod1 = y_mod1 / (2.0 * (1.0 + p_ratio1)) b_mod1 = y_mod1 / (3.0 * (1.0 - 2.0 * p_ratio1)) s_mod2 = y_mod2 / (2.0 * (1.0 + p_ratio2)) b_mod2 = y_mod2 / (3.0 * (1.0 - 2.0 * p_ratio2)) s_mod3 = y_mod3 / (2.0 * (1.0 + p_ratio3)) b_mod3 = y_mod3 / (3.0 * (1.0 - 2.0 * p_ratio3))endset y_mod1 = 3.35e9 p_ratio1 = 0.25 set y_mod2 = 4.20e9 p_ratio2 = 0.25set y_mod3 = 4.10e9 p_ratio3 = 0.25derive;----------------------------------------model hoekbrownprop bulk=b_mod1 shear=s_mod1 hbs=0.0054 hbmb=1.203 hbsigci=2.8E7 & hba=0.508 hbs3cv=0 range z -30. -5.prop bulk=b_mod2 shear=s_mod2 hbs=0.0084 hbmb=1.403 hbsigci=3.1E7 & hba=0.508 hbs3cv=0 range z -5. 6.prop bulk=b_mod3 shear=s_mod3 hbs=0.0082 hbmb=1.302 hbsigci=3.0E7 & hba=0.506 hbs3cv=0 range z 6. 30.;pauseprop density 2560;---------------------------------------set _beganOnY = 0SetVariables2ExcavatingInstallBoltOfOneSection;-------------------------------------InstallCableOfOneSectionstep _cableSteps;pause;save excavate01.sav;-------------------------------------set _beganOnY = 2.4SetVariables2ExcavatingInstallBoltOfOneSectionstep _cableSteps;pause;save excavate02.sav;--------------------------------------set _beganOnY = 4.8SetVariables2ExcavatingInstallBoltOfOneSectionInstallCableOfOneSectionstep _cableSteps;pause;save excavate03.sav;-------------------------------------set _beganOnY = 7.2SetVariables2ExcavatingInstallBoltOfOneSectionInstallCableOfOneSectionstep _cableSteps;save excavate04.sav;-------------------------------------set _beganOnY = 9.6SetVariables2ExcavatingInstallBoltOfOneSectionInstallCableOfOneSectionstep _cableSteps;save excavate05.sav;-------------------------------------solvesave excavateFinal.sav
下面是我搜到的一些解释： FLAC3D用户中有较大比例是研究地下结构支护的，由其是在煤矿方面的应用更多。因此，结合自己以前编写的命令流，发一个简单的、供初学者参考的例子。实为供初学者入门参考，为抛砖引玉。
在巷道的模拟当中有这样几个关键问题需要注意：
1）边界条件的设置及地应力的施加
一般巷道模型在z坐标方向都不可能从地表开始建模，因此截取的模型部分不要忘了加一个顶面的应力边界--apply szz；如果模型较大，最好考虑初始应力的渐变。
2）锚杆的安装
锚杆的坐标通过循环语句依据锚杆的长度自动生成，为体现打在不同纵深的锚杆则只需要改变原锚杆坐标的y坐标即可，需要注意的是用于存储锚杆起点和终点坐标的数组需要事先采用array关键词进行声明，而通常fish里面使用的变量一般都不需要显式的声明其变量类型，数组则不同。数组的大小是固定的。
3）围岩本构
本人倾向于采用Hoek-Brown本构模型模拟围岩，可以依据围岩表观特征按照GSI指标对围岩的强度适当折减。高级的应用还可以考虑岩石的峰后效应。
4）应力释放
应力释放确实是一个较难准确考虑的问题。巷道截面一旦开挖就会出现应力的释放，但究竟但支护结构安装前究竟释放了多大的应力很难准确考量。相关讨论可以参考论坛里的帖子，借鉴隧道中的应力释放来考虑。本人倾向于依据施工现场测得或观察到的围岩位移和变形情况对应力释放的大小进行估计。以此确定模拟中支护结构安装前的计算步数。在本示例中没有采用复杂的应力释放方式。
5）、因为喷层使用liner模拟的，所以，显示liner的法向应力的菜单操作：plotitems -> add ->structure element-> liner ->coupling -> stress -> normal
6）、显示锚杆轴力：plot add sel cable force 最终效果，如图：
7）、在建模时要定义好分组，在开挖时是分别开挖upexcav和lowexcav这两个分组，因此，分组的位置和范围对，才能正确的开挖
