FLAC-3D-學習疑點-詳解

1、FLAC3D資料網(wǎng)上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.關(guān)于計算不收斂問題;初始網(wǎng)格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 ; 給巷道直墻部分（左側(cè)）加壓力apply nstress -0.05 plane dip 90 dd 0 range hangdaoliangce2 ; 給巷道直墻部分（右側(cè)）加壓力 apply nstress -0.05 plane dip 90 dd 0

8、range hangdaoliangce1;對點進行監(jiān)測hist unbalhist gp xdisp 12.75 0 -15.75:監(jiān)測左幫hist gp zdisp 15 0 -13.5 ;監(jiān)測頂板hist gp zdisp 15 0 -16.5 ;監(jiān)測底板solve問：第一個solve平衡了,是第二個solve不能平衡，這是為什么？如果我把支護力加大，也不能平衡回答：我遇到的不收斂的原因，有兩個：1）檢查你的模型是不是沒有建好，相互之間是不是沒有連接好；2）仔細檢查你的參數(shù) ，是不是已經(jīng)發(fā)生了破壞。回答：力施加有問題，隧道底部隆起破壞，所以無法計算問：我這個模型就是該模型就是要出現(xiàn)巷道底

9、臌現(xiàn)象，如果按照你的意思，力應該怎樣施加才對？ 回答：如果無法收斂，應該用step命令，關(guān)鍵問題是要運行多少步。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是否是隨機生成數(shù)的命令？回答：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的圖像數(shù)據(jù)導出來以便用其他軟件繪圖？回答：用set log on 命令，把數(shù)據(jù)導出來，轉(zhuǎn)到excel里處理一下，然后用surfer或者什么作圖軟件繪制就行了。5. 問：用命令建立模型后，如何顯示點的坐標？回答：使用plo blo gro gpnum on命令6.interface問題問：我

11、加了一個層面，初始平衡就出現(xiàn)了這么大的位移？是怎么原因？;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定義的材料特性有什么區(qū)別？回答：前者是對整個

28、模型（如果整個模型只分成一個region的話），后者是對具體的zone的。8.關(guān)于ff邊界設置的問題自由場地邊界（FF邊界）在動力分析中經(jīng)常采用，它的原理是采用粘滯阻尼器來模擬靜止邊界，將自由場節(jié)點的不平衡力加到主體的網(wǎng)格邊界上。注意：模型底部的動力邊界條件應當在app ff之前，進行app ff以后底部的邊界條件就自動轉(zhuǎn)為自由場地邊界?，F(xiàn)舉出一簡單例子,比如靜力模型的尺寸為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可以驗證,以上的施加方式得到的結(jié)果是完全一樣的.所謂的自由場邊界最重要的特征就是

30、復制特征,現(xiàn)對以上各例子做簡要分析.1）第二三兩句就是在10*10的底面范圍內(nèi)分別設置靜止邊界和動力條件,當加了ff后,ff的網(wǎng)格的條件即從臨近的maingrid的zone進行復制.所以ff的底面的位移邊界和動力和maingrid完全一樣.2）第二句直接加ff邊界,注意這時ff底面是free的,而且動力為0,但是接下來的3.4兩句僅僅給定的是z=-.1至.1的范圍,導致ff的底面也包括在內(nèi)了,所以作用方式和命令流1完全相同.3）同2的分析,同樣可以得到結(jié)論但是對于以上三個命令流中的兩句apply squiet dquiet nquiet rang z -.1 .1apply nstress v

31、alue hist wave rang z -.1 .1如果在maingrid內(nèi)也指定了范圍的話,情形就不一樣了,比如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時有個圖形，很清楚的顯示只在模型的四周和四個邊角處有自由邊界條件的！另外，很容易理解上面幾種情況的計算結(jié)果是相同的。因為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.復雜模型內(nèi)部接觸面的建立方法最近經(jīng)常遇到一些朋友在工程項目當中遇到計算模型中接觸面方面的問題，大多都很類似，主要是因為在前期計算當中忽略了可能存在的接觸面影響，從而在建模中沒有考慮到建立接觸面。而由于FLAC3D的接觸面的建立必須存在兩個界面，因此通常的做法是在建模中首先建立兩個分開的模型，然后再一個模型表面建立接觸面，最后將另外一個模型“移動”過來。而這種方法在已完成的模型當中就不能湊效了，同時因為模型通常都較復雜，建模本身花費的精力

35、就是巨大的，如果修改模型的話，那造成大量不必要的重復勞動。本文推薦一種更加實用的接觸面建立方法，通過238版本后的impgrid、pgrid等命令進行模型導入導出，花費很少的時間就可以解決復雜模型的接觸面建立問題。下面用一個實例來做一下說明。工況：3*3*3的模型，兩個group，模型正中間的一個單元為gropu 2，其他為group 1，如圖所示。 現(xiàn)在要在14號單元（即group1）周圍加上接觸面，模擬實際工程中可能存在的內(nèi)部結(jié)構(gòu)的接觸面問題。下面是工作步驟：（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得到的就是已經(jīng)建立好接觸面的模型 下面進行簡單的計算，以驗

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可以看到接觸面的接應力分布結(jié)果：10.關(guān)于gauss_dev對性質(zhì)進行高斯正態(tài)分布的問題？根據(jù)手冊上的說明：下面的命令設定一個平均摩擦角為40度，標準方差是±5。則命令如下：prop friction 40 gauss_dev 2問題：請問gaus

39、s_dev 2中的2是如何計算的？如果把±5改為±10，則命令應如何寫？40×5211.子程序運行方式問題如下:進行迭代時子程序是如何運行的？下面命令流運行后，紅色子程序僅運行一次，而綠色部分則在每次迭帶時均調(diào)用。能否請高手詳細結(jié)實一下程序的運行規(guī)則？為什么綠色每次調(diào)用，而紅色只調(diào)用一次?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面；）關(guān)于dd的理解：dip-direction angle, measured in the global xy-planeclockwise from the positive y-axis 這個是從y軸到x順時針方向所轉(zhuǎn)過的角度；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）方向轉(zhuǎn)向所要確定面的法線方向在xy平面上的投影的夾角，然后再確定dip，dip指的是從xy平面轉(zhuǎn)向所要確定的平面的角gen zone reflect dip 90 dd 90是以YZ面.gen zone reflect dip 90 dd 0 才是以xz為對稱面的.dip-direc

47、tion 這個是從y軸到x順時針方向所轉(zhuǎn)過的角度，如果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