FLAC-3D-学习疑点-详解

1、FLAC3D资料网上FLAC3D学习笔记1. 请问在圆柱体四周如何施加约束条件？可以用fix . ran cylinder end1 end2 radius r1 cylinder end1 end2 radius r2 not，其中r2<r1，其实就是选择range的时候选两个圆柱的差，即得到边界。命令流如下：newgen zone cyl p1 1 0 0 p2 0 4 0 p3 0 0 1 size 12 30 12gen zone refle norm 1 0 0 gen zone refle norm 0 0 1mo ssprop den 2500 bulk 2e8 she 1

2、e8 co 2e6 fric 45 ten 1e6 dil 10prop ftab 1 ctab 2 dtab 3table 1 0 45 0.05 42 0.1 40 1 40table 2 0 2e6 0.05 1e6 0.1 5e5 1 5e5table 3 0 10 0.05 3 0.1 0fix x range end1 1 0 0 end2 1 4 0 rad 1 end1 0.5 0 0 end2 0.5 4 0 rad 0.5回答：加上的约束不对！fix x range end1 1 0 0 end2 1 4 0 rad 1 end1 0.5 0 0 end2 0.5 4 0

3、rad 0.5回答：上面的rad 是不是相差太大了？是不是在1左右来限制了2.关于计算不收敛问题;初始网格gen zone radcyl p0 15 0 -15.75 p1 add 15 0 0 p2 add 0 6 0 p3 add 0 0 15.75 &size 3 4 8 5 dim 2.25 2.25 2.25 2.25 rat 1 1 1 1.1 fillgen zone radtun size 3 4 3 5 dim 0.75 2.25 0.75 2.25 rat 1 1 1 1.1 & p0 15 0 -15.75 p1 add 0 0 -14.25 p

4、2 add 0 6 0 p3 add 15 0 0 fillgen zone reflect dip 90 dd 270 ori 15 0 -15.75; 摩尔库仑模型model mohrmacro niyan 'bulk 3.2825e3 shear 1.5150e3 fric 20 coh 0.07'macro meiceng 'bulk 2.1735e3 shear 8.8917e2 fric 19 coh 0.05'prop niyan range z =-13.5 0prop meiceng range z =-15.5 -13.5prop niyan

5、 range z = -30 -15.5;为开挖定义整体部分range name in z -15.75 0 cyl end1 15 0 -15.75 end2 15 6 -15.75 rad 2.25 range name hangdao x 12.75 17.25 y 0 6 z -16.50 -15.75;进行加载fix x y z range z -30.001 -29.999app sxx=-20 range x -0.001 0.001app sxx=-20 range x 29.999 30.001app syy=-20 range y -0.001 0.001app syy=-

6、20 range y 5.999 6.001app szz=-20 range z -0.001 0.001 solve（1）model null range in model null range hangdaorange name in z -15.75 0 cyl end1 15 0 -15.75 end2 15 6 -15.75 rad 2.25range name hangdaoliangce1 x 12.749 12.751 y 0 6 z -16.50 -15.75range name hangdaoliangce2 x 17.249 17.251 y 0 6 z -16.50

7、-15.75range name hangdao x 12.75 17.25 y 0 6 z -16.501 -16.499apply nstress -0.05 plane dip 90 dd 0 range in ; 给巷道底部加压力app nstress -0.05 plane dip 90 dd 0 range hangdao ; 给巷道直墙部分（左侧）加压力apply nstress -0.05 plane dip 90 dd 0 range hangdaoliangce2 ; 给巷道直墙部分（右侧）加压力 apply nstress -0.05 plane dip 90 dd 0

8、range hangdaoliangce1;对点进行监测hist unbalhist gp xdisp 12.75 0 -15.75:监测左帮hist gp zdisp 15 0 -13.5 ;监测顶板hist gp zdisp 15 0 -16.5 ;监测底板solve问：第一个solve平衡了,是第二个solve不能平衡，这是为什么？如果我把支护力加大，也不能平衡回答：我遇到的不收敛的原因，有两个：1）检查你的模型是不是没有建好，相互之间是不是没有连接好；2）仔细检查你的参数 ，是不是已经发生了破坏。回答：力施加有问题，隧道底部隆起破坏，所以无法计算问：我这个模型就是该模型就是要出现巷道底

9、臌现象，如果按照你的意思，力应该怎样施加才对？ 回答：如果无法收敛，应该用step命令，关键问题是要运行多少步。3.问：命令流如下：g 10 10m edef distortloop i (2,jzones)loop j (2,jzones)x(i,j) =x(i,j)+mag*(urand-0.5)y(i,j)=y(i,j)+mag*(urand-0.5)end_loopend_loopendprop d 1000 sh 1e8 bu 2e8其中：x(i,j) =x(i,j)+mag*(urand-0.5)中的urand是否是随机生成数的命令？回答：urand ：random number

10、drawn from uniform distribution between0.0 and 1.0；grand ：random number drawn from normal distribution: mean = 0.0; standard deviation = 1.04. 问：怎么能把一个PLOT的图像数据导出来以便用其他软件绘图？回答：用set log on 命令，把数据导出来，转到excel里处理一下，然后用surfer或者什么作图软件绘制就行了。5. 问：用命令建立模型后，如何显示点的坐标？回答：使用plo blo gro gpnum on命令6.interface问题问：我

11、加了一个层面，初始平衡就出现了这么大的位移？是怎么原因？;No1gen zone brick p0 0 0 0 p1 10 0 0 p2 0 10 0 p3 0 0 8 &size 8 8 6 rat 1.02 1.02 1.02gen zone brick p0 0 0 9 p1 10 0 9 p2 0 10 9 p3 0 0 11 &size 8 8 2 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,6) range z 0 11group layer1 range z 9 11;interface 1

12、 face range plane norm (0,0,1) origin (0,0,8) dist 0.01ini z add -1.0 range group layer1;No2gen zone brick p0 0 0 11 p1 10 0 11 p2 0 10 11 p3 0 0 13 &size 8 8 2 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,1) range z 11 13group layer2 range z 11 13;interface 2 face range plane norm

13、 (0,0,1) origin (0,0,10) dist 0.01ini z add -1.0 range group layer2;No3gen zone brick p0 0 0 13 p1 10 0 13 p2 0 10 13 p3 0 0 14 &size 8 8 1 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,1) range z 13 14group layer3 range z 13 14;interface 3 face range plane norm (0,0,1) origin (0,0,

14、12) dist 0.01ini z add -1.0 range group layer3;No4gen zone brick p0 0 0 14 p1 10 0 14 p2 0 10 14 p3 0 0 15 &size 8 8 1 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,1) range z 14 15group layer4 range z 14 15interface 4 face range plane norm (0,0,1) origin (0,0,13) dist 0.01ini z add

15、 -1.0 range group layer4;No5gen zone brick p0 0 0 15 p1 10 0 15 p2 0 10 15 p3 0 0 17 &size 8 8 2 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,15) range z 15 17group layer5 range z 15 17interface 5 face range plane norm (0,0,1) origin (0,0,14) dist 0.01ini z add -1.0 range group lay

16、er5;No6gen zone brick p0 0 0 17 p1 10 0 17 p2 0 10 17 p3 0 0 18 &size 8 8 1 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,1) range z 17 18group layer6 range z 17 18interface 6 face range plane norm (0,0,1) origin (0,0,16) dist 0.01ini z add -1.0 range group layer6;No7gen zone brick

17、p0 0 0 18 p1 10 0 18 p2 0 10 18 p3 0 0 20 &size 8 8 2 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,1) range z 18 20group layer7 range z 18 20;interface 7 face range plane norm (0,0,1) origin (0,0,17) dist 0.01ini z add -1.0 range group layer7;No8gen zone brick p0 0 0 20 p1 10 0 20

18、p2 0 10 20 p3 0 0 24 &size 8 8 3 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,17) range z 20 24group layer8 range z 20 24;interface 8 face range plane norm (0,0,1) origin (0,0,19) dist 0.01ini z add -1.0 range group layer8;No9gen zone brick p0 0 0 24 p1 10 0 24 p2 0 10 24 p3 0 0 29

19、 &size 8 8 4 rat 1.02 1.02 1.02gen zone reflect dip 90 dd 90 origin(0,0,0) range z 24 29group layer9 range z 24 29;interface 9 face range plane norm (0,0,1) origin (0,0,23) dist 0.01ini z add -1.0 range group layer9;-边界条件-fix x range x -20.1 -19.9fix x range x 19.9 20.1fix y range y -0.1 0.1fix

20、y range y 9.9 10.1fix z range z -0.1 0.1;-初始条件-apply szz -26e6 range z 27.9 28.1initial szz -1.08e6 grad 0 0 2.6e4 range z 0 28set gravity 0 0 -10initial density 2.6e3initial sxx -45.8e6 grad 0 0 2.6e4 range z 0 28initial syy -45.8e6 grad 0 0 2.6e4 range z 0 28;-model mohr;-材料属性-;-No1-pro bulk 8e9 s

21、hear 7e9 fric 25 coh 7.5e6 ten 2.6e6 range z 0 8 ;-No2-pro bulk 7.8e9 shear 3.5e9 fric 30 coh 6.5e6 ten 2.4e6 range z 8 10;-No3-pro bulk 7e9 shear 3.0e9 fric 22 coh 5.5e6 ten 2.2e6 range z 10 12;-No4-pro bulk 6e9 shear 2.3e9 fric 26 coh 5.5e6 ten 2.0e6 range z 12 13;-No5-pro bulk 6.2e9 shear 2.0e9 f

22、ric 28 coh 5.2e6 ten 1.2e6 range z 13 14;-No6-pro bulk 5.8e9 shear 2.5e9 fric 35 coh 6.0e6 ten 1.3e6 range z 14 16;-No7-pro bulk 9e9 shear 3.2e9 fric 30 coh 10.3e6 ten 1.2e6 range z 16 17;-No8-pro bulk 8.1e9 shear 3.1e9 fric 26 coh 8.0e6 ten 1.5e6 range z 17 19;-No9-pro bulk 7.6e9 shear 2.6e9 fric 2

23、3 coh 8.5e6 ten 1.8e6 range z 19 23;-No10-pro bulk 7.5e9 shear 2.1e9 fric 26 coh 8.0e6 ten 1.9e6 range z 23 28;int 1 prop ks 8e10 kn 8e10 fric 30 dil 6 tension 2e7 bslip=on;int 2 prop ks 8e10 kn 8e10 fric 30 dil 6 tension 2e7 bslip=on;int 3 prop ks 8e10 kn 8e10 fric 35 dil 8 tension 2e7 bslip=onint

24、4 prop ks 9e11 kn 9e11 fric 32 dil 5 tension 2e7 bslip=onint 5 prop ks 9e11 kn 9e11 fric 32 dil 5 tension 2e7 bslip=onint 6 prop ks 8e11 kn 8e11 fric 33 dil 5 tension 2e7 bslip=on;int 7 prop ks 8e10 kn 8e10 fric 33 dil 5 tension 2e7 bslip=on;int 8 prop ks 8e10 kn 8e10 fric 33 dil 5 tension 2e7 bslip

25、=on;int 9 prop ks 8e10 kn 8e10 fric 33 dil 5 tension 2e7 bslip=onhist unbalstep 5000 apply szz -26e6 range z 27.9 28.1initial szz -1.08e6 grad 0 0 2.6e4 range z 0 28set gravity 0 0 -10initial density 2.6e3initial sxx -45.8e6 grad 0 0 2.6e4 range z 0 28initial syy -45.8e6 grad 0 0 2.6e4 range z 0 28i

26、nitial szz -1.08e6 grad 0 0 2.6e4 range z 0 28initial density 2.6e3initial sxx -45.8e6 grad 0 0 2.6e4 range z 0 28initial syy -45.8e6 grad 0 0 2.6e4 range z 0 28这是给的模型加的初始地应力，难道这样加地应力，模型就自动平衡了吗？回答：你用ini平衡，为什么还要计算5000步？问：我的INI是用来施加初始地应力的，难道用INI那样加应力边界条件就不用step 或solve算平衡了吗？ 模型就能自动平衡吗？;-初始条件-apply szz

27、-26e6 range z 27.9 28.1initial szz -1.08e6 grad 0 0 2.6e4 range z 0 28set gravity 0 0 -10initial density 2.6e3initial sxx -45.8e6 grad 0 0 2.6e4 range z 0 28initial syy -45.8e6 grad 0 0 2.6e4 range z 0 28回答：把上面的sxx和syy的剃度2.6e4改为9e3试一下，xx和yy方向的梯度不可能与zz方向的梯度一样的。7.问：prop mat 和zone定义的材料特性有什么区别？回答：前者是对整个

28、模型（如果整个模型只分成一个region的话），后者是对具体的zone的。8.关于ff边界设置的问题自由场地边界（FF边界）在动力分析中经常采用，它的原理是采用粘滞阻尼器来模拟静止边界，将自由场节点的不平衡力加到主体的网格边界上。注意：模型底部的动力边界条件应当在app ff之前，进行app ff以后底部的边界条件就自动转为自由场地边界。现举出一简单例子,比如静力模型的尺寸为10*10*10,有以下几种边界施加方式1）free x y z rang z -.1 .1apply squiet dquiet nquiet rang z -.1 .1apply nstress value hist

29、wave rang z -.1 .1apply ff2）free x y z rang z -.1 .1apply ffapply squiet dquiet nquiet rang z -.1 .1apply nstress value hist wave rang z -.1 .13）free x y z rang z -.1 .1apply squiet dquiet nquiet rang z -.1 .1apply ffapply nstress value hist wave rang z -.1 .1可以验证,以上的施加方式得到的结果是完全一样的.所谓的自由场边界最重要的特征就是

30、复制特征,现对以上各例子做简要分析.1）第二三两句就是在10*10的底面范围内分别设置静止边界和动力条件,当加了ff后,ff的网格的条件即从临近的maingrid的zone进行复制.所以ff的底面的位移边界和动力和maingrid完全一样.2）第二句直接加ff边界,注意这时ff底面是free的,而且动力为0,但是接下来的3.4两句仅仅给定的是z=-.1至.1的范围,导致ff的底面也包括在内了,所以作用方式和命令流1完全相同.3）同2的分析,同样可以得到结论但是对于以上三个命令流中的两句apply squiet dquiet nquiet rang z -.1 .1apply nstress v

31、alue hist wave rang z -.1 .1如果在maingrid内也指定了范围的话,情形就不一样了,比如2命令流改为apply squiet dquiet nquiet rang z -.1 .1 x 0 10 y 0 10apply nstress value hist wave rang z -.1 .1 x 0 10 y 0 10这样就相当于ff的底面为ff且动力为0,而maingrid的底面为quiet且动力为nstress=value*wave.假如先重力平衡，分1）fix z rang bottom情况在加quiet边界和应力波边界的时候需要先free z rang

32、bottom吗，手册上也有app dquiet squiet,后app nvel 0 rang bottom,如果app sxz his wave 是不是就不用free z了2）fix x y z rang bottom情况因为手册上有加quiet边界，反力就会自动加载在quiet边界上，所以是不是根本就不用free x y z呢?像free 与不free的地表响应加速度差别很大呢.手册中有这个例子：free x y z rang z -.1 .1apply squiet dquiet nquiet rang z -.1 .1apply nstress value hist wave rang

33、 z -.1 .1apply ff我的理解是：第一句：释放底面约束，为下面的输入提供条件第二句：在底面施加quiet边界条件第三句：在底面输入外荷载第四句：是在模型的四周施加自由场边界，底面是当然没有自由场边界的。手册讲ff时有个图形，很清楚的显示只在模型的四周和四个边角处有自由边界条件的！另外，很容易理解上面几种情况的计算结果是相同的。因为flac里一般是不讲先后顺序的，只要你在solve之前都设置好了就可以了。不过对于apply ff命令手册里好像提到了一个限制：（K17）The dynamic boundary conditions at the base of the model sh

34、ould be specified before applying the free-field. 即：在施加自由场边界之前应设置好模型底面的动力边界条件。9.复杂模型内部接触面的建立方法最近经常遇到一些朋友在工程项目当中遇到计算模型中接触面方面的问题，大多都很类似，主要是因为在前期计算当中忽略了可能存在的接触面影响，从而在建模中没有考虑到建立接触面。而由于FLAC3D的接触面的建立必须存在两个界面，因此通常的做法是在建模中首先建立两个分开的模型，然后再一个模型表面建立接触面，最后将另外一个模型“移动”过来。而这种方法在已完成的模型当中就不能凑效了，同时因为模型通常都较复杂，建模本身花费的精力

35、就是巨大的，如果修改模型的话，那造成大量不必要的重复劳动。本文推荐一种更加实用的接触面建立方法，通过238版本后的impgrid、pgrid等命令进行模型导入导出，花费很少的时间就可以解决复杂模型的接触面建立问题。下面用一个实例来做一下说明。工况：3*3*3的模型，两个group，模型正中间的一个单元为gropu 2，其他为group 1，如图所示。 现在要在14号单元（即group1）周围加上接触面，模拟实际工程中可能存在的内部结构的接触面问题。下面是工作步骤：（1）将模型文件存为1.sav；（2）将除了需要加接触面的单元以外的所有单元删除，仅保存group 2命令为：del ran gro

36、up 2 not；（3）在group 2 周围建立接触面,命令为：interface 1 face range x 1 y 1 2 z 1 2interface 1 face range x 2 y 1 2 z 1 2interface 1 face range x 1 2 y 1 z 1 2interface 1 face range x 1 2 y 2 z 1 2interface 1 face range x 1 2 y 1 2 z 1interface 1 face range x 1 2 y 1 2 z 2建成后可以看到接触面的效果：plo inter yel ske可以用inter

37、face 1 face命令来默认建立整个模型的外表面接触面，这个将更加方便实用。对于复杂模型，大家可以试一试。（4）这时候需要将上面建好的group 2 和接触面进行保存save 2.sav（5）重新restore先前的模型，并将需要建立接触面的实体删除，再将整个模型用expgrid命令进行导出：rest 1.savdel ran group 2expgrid 1.fac3d下面的工作是这样的：将已建好的group 2和接触面的文件rest，并把其他的模型用impgrid命令进行导入:rest 2.savimpgrid 1.flac3d得到的就是已经建立好接触面的模型 下面进行简单的计算，以验

38、证接触面是否工作。model elaprop bulk 20e8 shear 30e8fix x y z ran z 0ini den 2000set gra 0 0 -10interface 1 prop kn 20e6 ks 20e6 c 10e3 f 15app nstr -200e3 ran x 0 1 y 1 2 z 3solve可以看到接触面的接应力分布结果：10.关于gauss_dev对性质进行高斯正态分布的问题？根据手册上的说明：下面的命令设定一个平均摩擦角为40度，标准方差是±5。则命令如下：prop friction 40 gauss_dev 2问题：请问gaus

39、s_dev 2中的2是如何计算的？如果把±5改为±10，则命令应如何写？40×5211.子程序运行方式问题如下:进行迭代时子程序是如何运行的？下面命令流运行后，红色子程序仅运行一次，而绿色部分则在每次迭带时均调用。能否请高手详细结实一下程序的运行规则？为什么绿色每次调用，而红色只调用一次?newtitle Direct shear testset log onset logfile fa3.4.loggen zone brick size 12 1 10 p0 4 0 6 p1 16 0 6 p2 4 1 6 p3 4 0 11gen zone brick siz

40、e 20 1 10 p1 20 0 0 p2 0 1 0 p3 0 0 5range name bot z 0 5range name top z 6 11interface 1 face range z 5int 1 prop ks 4e4 kn 4e4 fric 30 dil 6 ;tension 1e10 bslip=onini z add -1.0 range top;plo surf lorange interface white axes blackmodel eprop bulk 45e3 sh 30e3fix x y z range z 0fix x range x 0fix

41、x range x 20plot create view_intplot add surfaceplot add interface redplot show;pauseapply nstress -10 range z 10step 0plot contour szz interface white axes black;step 100solve;pausesave dsta.savini xvel 5e-7 range topfix xvel range topdef ini_jdispvalnd = 0.0count = 0.0p_in=i_node_head(i_head)loop

42、while p_in # nullif in_ztarget(p_in) # null thenvalnd = valnd + in_pen(p_in)count = count + 1.0end_ifp_in = in_next(p_in)end_loopnjdisp0 = valnd/countendini_jdispdef sstavvalns = 0.0valss = 0.0valsd = 0.0valnd = 0.0count = 0.0p_in=i_node_head(i_head)loop while p_in # nullif in_ztarget(p_in) # null t

43、henvalns = valns + in_nstr(p_in)*in_area(p_in)valss = valss + in_sstr(p_in,1)*in_area(p_in)valsd = valsd + in_sdisp(p_in,1)valnd = valnd + in_pen(p_in)count = count + 1.0end_ifp_in = in_next(p_in)end_loopsstav = valss/(12.0*1.0)nstav = valns/(12.0*1.0)sjdisp = valsd/countnjdisp = valnd/count - njdis

44、p0endsstavhist ns 1hist sstav nstav sjdisp njdisp ini xdis 0 ydis 0 zdis 0step 2500save dst.savplot his -1 vs -3;pauseplot his -4 vs -3print nstav;pauseret12. reflect问题问：）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 fillp

45、lot surfgen zone reflect dip 0 dd 90 这个是以xy平面作为对称面来对称。那么下面这句话：gen zone reflect dip 90 dd 90按照上面的分析应该按照xz平面对称，可是manual是yz平面，按理说应该是yz面；）关于dd的理解：dip-direction angle, measured in the global xy-planeclockwise from the positive y-axis 这个是从y轴到x顺时针方向所转过的角度；dip dip angle, measured in the negative z-direction

46、 fromthe global xy-plane 这个是xy平面与z轴负方向的夹角。dip表示对称参照面与xy平面的夹角，对称参照面与xy平面的夹角在xy平面的投影是一条射线 ，dd表示y轴正向顺时针到那条射线的夹角。首先应该按照dd的方向大体确定这个面的朝向，dd 指的是从y轴正方向按顺时针（clockwise）方向转向所要确定面的法线方向在xy平面上的投影的夹角，然后再确定dip，dip指的是从xy平面转向所要确定的平面的角gen zone reflect dip 90 dd 90是以YZ面.gen zone reflect dip 90 dd 0 才是以xz为对称面的.dip-direc

47、tion 这个是从y轴到x顺时针方向所转过的角度，如果dd=0,那不就是y轴了吗？所以gen zone reflect dip 90 dd 0 不应是xz面，应该是yz面啊？回答：dip-direction 是倾向，也就是所求面的法向在XY面上的投影的，dd 0表示倾向是0度既是Y轴正方向，所以是XZ面，那dip=0时，即平面与xy的夹角为0,也就是平行于xy面，还要定义dd干吗？dd为任何值不都一样吗？13.第二次国际会议上的一个示例模型def tunnel_datax0 = 0.0y1 = 0.0z0 = 0.0rad = 1.0base = 0.7height = 0.5dx = 2.0

48、dz = 2.0y2 = 5.0y3 = 15.0xr = 4.0zt = 4.0nx = 4nr = 3nu = 4nz1 = 4nxr = 5nzt=5ny1 = 8ny2=4raty = 1.2ratr = 1.5;* DERIVED PARAMETERS *nz2 = 2*nx ; upper-right and top zonesratr2=1/ratr;Compute x-coordinatesx1 = x0+radx2 = x0+dxx3 = x0+base; Compute arc center (x4,z0) to fit arc-points (x1,z0) and (x0

49、+base,z0-height)x4=x0+(rad*rad-base*base-height*height)/(2.0*(rad-base); Compute radius of lower arcrad2=rad-(x4-x0)x7=x0+sqrt(2.0)/2*rad; compute inner points for radtunnel - adjust based on nx,nu zonesx5=(x0*nu+x7*nx)/(nx+nu)x5c=(nu*x0+nx*x7)/(nx+nu)x5b=(nu*x0+2*nx*x7)/(2*nx+nu)x6=(x0*nu+x2*nx)/(n

50、x+nu);Compute z coordinatesz1 = z0+radz2 = z0+dzz4 = z0-dzz5 = z0-heightz6 = (z0+z4)/2z8=z0+sqrt(2.0)/2*radz7=(z0*nu+z8*nx)/(nx+nu)z7c=(nu*z0+nx*z8)/(nx+nu)z7b=(nu*z0+2*nx*z8)/(2*nx+nu);slope constant of lower arcdxz=(x3-x1)/(z5-z0);slope constants of upper arcdxz2=(z8-z0)/(x7-x1)dxz3=(x7-x0)/(z8-z1

51、)enddef make1if y1#y2 thencommand; top and upper-right outsidegen zone radcyl size nx ny1 nz2 nr ratio 1 1 1 ratr group 'g1' &p0 x0 y1 z0 & p1 x2 y1 z0 & p2 x0 y2 z0 & p3 x0 y1 z2 &  p4 x2 y2 z0 &  p5 x0 y2 z2 & p6 x2 y1 z2 & p7 x2 y2 z2 & p8 x1

52、y1 z0 & p9 x0 y1 z1 & p10 x1 y2 z0 & p11 x0 y2 z1        end_commandend_ifif y2#y3 thencommandgen zone radcyl size nx ny2 nz2 nr ratio 1 raty 1 ratr group 'g1' &p0 x0 y2 z0 & p1 x2 y2 z0 & p2 x0 y3 z0 & p3 x0 y2 z2 &

53、60; p4 x2 y3 z0 & p5 x0 y3 z2 & p6 x2 y2 z2 & p7 x2 y3 z2 &  p8 x1 y2 z0 & p9 x0 y2 z1 & p10 x1 y3 z0 &  p11 x0 y3 z1       end_commandend_ifenddef make2if y1#y2 thencommand; lower-right outsidegen zone radcyl size nx ny1

54、nz1 nr ratio 1 1 1 ratr group 'g2' &p0 x4 y1 z0 &  p1 x2 y1 z4 &  p2 x4 y2 z0 &   p3 x2 y1 z0 &  p4 x2 y2 z4 &  p5 x2 y2 z0 &  p6 x2 y1 z6 &  p7 x2 y2 z6 &  p8 x3 y1 z5 &  p9 x1 y1 z0& p10 x3 y2 z5

55、&  p11 x1 y2 z0       end_commandend_ifif y2#y3 thencommand; lower-right outsidegen zone radcyl size nx ny2 nz1 nr ratio 1 raty 1 ratr group 'g2' &p0 x4 y2 z0 & p1 x2 y2 z4 &  p2 x4 y3 z0 & p3 x2 y2 z0 &  p4 x2 y3 z4 &

56、amp; p5 x2 y3 z0 &  p6 x2 y2 z6 & p7 x2 y3 z6 & p8 x3 y2 z5 &  p9 x1 y2 z0 &  p10  x3 y3 z5 &  p11 x1 y3 z0 end_commandend_ifenddef make3if y1#y2 thencommand; bottom outsidegen zone brick size nx ny1 nr ratio 1 1 ratr2 group 'g3' &am

57、p;p0 x0 y1 z4 & p1 x6 y1 z4 & p2 x0 y2 z4 &  p3 x0 y1 z5 & p4 x6 y2 z4 &  p5 x0 y2 z5 &  p6 x5 y1 z5 &  p7 x5 y2 z5       end_commandend_ifif y2#y3 thencommand; bottom outsidegen zone brick size nx ny2 nr ratio 1 raty

58、 ratr2 group 'g3' &p0 x0 y2 z4 & p1 x6 y2 z4 &  p2 x0 y3 z4 &  p3 x0 y2 z5 &  p4 x6 y3 z4 & p5 x0 y3 z5 &  p6 x5 y2 z5 &  p7 x5 y3 z5         end_commandend_ifenddef make4if y1#y2 thencommand; bottom-right outsidegen zone brick size nu ny1 nr ra