數值模型 PFC 3D 程式碼

def Parameter X0 = 0.0

def ComputeParameter

shDipR = shDip / 180 * pi

2. 設定砂岩塊體參數：

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

;01b_DefineParameter_B.txt

;定義砂岩尺寸、強度

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

def Parameter_B ssB = 1.8 ssH = 3.6 ssW = 2.0 ssOH = 0.4 end

def ComputeParameter_B

ssX1 = (shW / 2) - (ssW / 2)

3. 軟岩分層：

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

;03_R3_LayerShBall.txt

;軟岩分層

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

def LayerShBall layerD = 0.04 layerShift = 0.0 layerSpan = 0.25 initialY = shY1

initialZ = shZ2 + (shY2 - shY1) * tan(shDipR) - layerShift dY = 0.01

dZ = dY * tan(shDipR) searchY1 = initialY searchY2 = initialY + dY searchZ1 = initialZ searchZ2 = initialZ + layerD

i = 1 ;層數計算

m = 1 ;控制外圈開閉

loop while m = 1

layerName = 'ShL' + string(i)

n = 1 ;控制內圈開閉

command

group layerName owner ShBallCenter range x ssX1 ssX2 y searchY1 searchY2 z searchZ1 searchZ2 end_command

loop while n = 1

searchY1 = searchY1 + dY searchY2 = searchY2 + dY searchZ1 = searchZ1 - dZ searchZ2 = searchZ2 - dZ command

group layerName range x ssX1 ssX2 y searchY1 searchY2 z searchZ1 searchZ2 end_command

if searchY2 > boxY2 then n = 0

end_if

if searchZ1 < boxZ1 then n = 0

end_if end_loop

searchY1 = initialY searchY2 = initialY + dY

searchZ1 = initialZ - layerSpan * i searchZ2 = searchZ1 + layerD greenIndex = 1 - i * 0.1 command

plot add ball rgb 0.3 greenIndex 0.9 range group layerName end_command

if searchZ1 < boxZ1 then m = 0

4. 主程式Part 1：

;########################### PART 01 ###########################

new

set gen_error off ;當generate 指令生滿球時，關閉警告

call 01a_R3_DefineParameter_A.txt

;===============================================================建外牆

wall id=1 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, boxZ1) (shX2, boxY1, boxZ1) (shX2, boxY2, boxZ1) (shX1, boxY2, boxZ1) wall id=2 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, boxZ1) (shX1, boxY2, boxZ1) (shX1, boxY2, boxZ2) (shX1, boxY1, boxZ2) wall id=3 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, boxZ1) (shX1, boxY1, boxZ2) (shX2, boxY1, boxZ2) (shX2, boxY1, boxZ1) wall id=4 kn=1e14 ks=1e14 f=0.5 face (shX2, boxY1, boxZ1) (shX2, boxY1, boxZ2) (shX2, boxY2, boxZ2) (shX2, boxY2, boxZ1) wall id=5 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY2, boxZ1) (shX2, boxY2, boxZ1) (shX2, boxY2, boxZ2) (shX1, boxY2, boxZ2)

;===============================================================建中牆

wall id=6 kn=1e14 ks=1e14 f=0.5 face (shX1, shY1, boxZ1) (shX1, shY1, boxZ2) (shX2, shY1, boxZ2) (shX2, shY1, boxZ1)

;===============================================================生球 generate id=1, 19999 rad=0.08, 0.08 x=shX1, shX2 y=shY1, boxY2 z=boxZ1, 10.0 generate id=20000, 39999 rad=0.08, 0.08 x=shX1, shX2 y=shY1, boxY2 z=10.0, boxZ2 def DetermineType

if plH # 0 then command

generate id=40000, 49999 rad=0.08, 0.08 x=shX1, shX2 y=boxY1, shY1 z=boxZ1, boxZ2 end_command

end_if end

DetermineType

property den=20000 fric=0.5 kn=1e8 ks=1e8 ;霣降時採用之球體參數

;===============================================================設定球群組 group ShBall range id=1, 49999

group SsBall range id=50000, 69999

;===============================================================顯示設定

plot set perspective off

;===============================================================霣降

wall id=7 kn=1e14 ks=1e14 f=0.5 face (shX1, shY2, shZ2) (shX2, shY2, shZ2) (shX2, boxY2, shZ1) (shX1, boxY2, shZ1) zv=0.3 wall id=8 kn=1e14 ks=1e14 f=0.5 face (shX1, shY2, shZ2) (shX1, boxY2, shZ1) (shX2, boxY2, shZ1) (shX2, shY2, shZ2) wall id=9 kn=1e14 ks=1e14 f=0.5 face (shX1, shY1, Z0) (shX2, shY1, Z0) (shX2, shY2, shZ2) (shX1, shY2, shZ2) zv=0.3 wall id=10 kn=1e14 ks=1e14 f=0.5 face (shX1, shY1, Z0) (shX1, shY2, shZ2) (shX2, shY2, shZ2) (shX2, shY1, Z0) wall id=11 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, Z0) (shX2, boxY1, Z0) (shX2, shY1, Z0) (shX1, shY1, Z0) zv=0.3 wall id=12 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, Z0) (shX1, shY1, Z0) (shX2, shY1, Z0) (shX2, boxY1, Z0)

wall id=13 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, boxZ2) (shX1, boxY2, boxZ2) (shX2, boxY2, boxZ2) (shX2, boxY1, boxZ2)

;===============================================================抬升7、9、11 號牆 set dt 0.002

set g 0 0 -9.81 cyc 5000

;===============================================================刪掉7、9、11 號牆以上的球，移除6 號牆 del ball range z shZ2 boxZ2

del wall 7 9 11 del wall 6

;===============================================================存檔，軟岩形狀完成 save R3_ShBall_1e8.sav

;===============================================================檔名為球勁度

;########################### END 01 ###########################

5. 主程式Part 2：

;########################### PART 02 ###########################

;===============================================================給軟岩鍵結 def SetPBProperty

PBkn = 1e10

property den=2200 fric=0.5 range group ShBall

property pb_kn=PBkn pb_ks=PBks pb_n=PBns pb_s=PBss pb_r=PBr range group ShBall

;===============================================================跑10000 步平衡

wall id=7 kn=1e14 ks=1e14 f=0.5 face (shX1, shY2, shZ2) (shX2, shY2, shZ2) (shX2, boxY2, shZ1) (shX1, boxY2, shZ1) zv=0.3 wall id=9 kn=1e14 ks=1e14 f=0.5 face (shX1, shY1, Z0) (shX2, shY1, Z0) (shX2, shY2, shZ2) (shX1, shY2, shZ2) zv=0.3 wall id=11 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, Z0) (shX2, boxY1, Z0) (shX2, shY1, Z0) (shX1, shY1, Z0) zv=0.3

;===============================================================抬升7、9、11 號牆 set dt 0.002

set g 0 0 -9.81 cyc 10000

;===============================================================刪掉7、9、11 號牆以上的球 del ball range z shZ2 boxZ2

del wall 7 9 11

wall id=7 kn=1e14 ks=1e14 f=0.5 face (shX1, shY2, shZ2) (shX2, shY2, shZ2) (shX2, boxY2, shZ1) (shX1, boxY2, shZ1) zv=0.3 wall id=9 kn=1e14 ks=1e14 f=0.5 face (shX1, shY1, Z0) (shX2, shY1, Z0) (shX2, shY2, shZ2) (shX1, shY2, shZ2) zv=0.3 wall id=11 kn=1e14 ks=1e14 f=0.5 face (shX1, boxY1, Z0) (shX2, boxY1, Z0) (shX2, shY1, Z0) (shX1, shY1, Z0) zv=0.3

;===============================================================抬升7、9、11 號牆（第二次）

cyc 10000

;===============================================================刪掉7、9、11 號牆以上的球（第二次）

del ball range z shZ2 boxZ2 del wall 7 9 11

;===============================================================刪除上牆，再跑10000 步平衡2，補強鍵結 del wall 8 12

property pb_kn=PBkn pb_ks=PBks pb_n=PBns pb_s=PBss pb_r=PBr range group ShBall cyc 10000

prop xvel 0 yvel 0 zvel 0 range group ShBall prop xdis 0 ydis 0 zdis 0 range group ShBall prop xspin 0 yspin 0 zspin 0 range group ShBall plot show

;===============================================================存檔，軟岩含強度 save 07001200.sav

;===============================================================檔名為日期時間

;########################### END 02 ###########################

6. 主程式Part 3：

;########################### PART 03 ###########################

call 01b_R3_DefineParameter_B.txt

;===============================================================建模板牆

wall id=50 kn=1e12 ks=1e12 f=0.5 face (shX1, ssY1B, ssZ1F) (shX1, ssY2B, ssZ1B) (shX2, ssY2B, ssZ1B) (shX2, ssY1B, ssZ1F) wall id=51 kn=1e12 ks=1e12 f=0.5 face (shX1, ssY2B, ssZ1B) (shX2, ssY2B, ssZ1B) (shX2, ssY2T, ssZ2B) (shX1, ssY2T, ssZ2B) wall id=52 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY1B, ssZ1F) (ssX1, ssY2B, ssZ1B) (ssX1, ssY2T, ssZ2B) (ssX1, ssY1T, ssZ2F) wall id=53 kn=1e12 ks=1e12 f=0.5 face (ssX2, ssY1B, ssZ1F) (ssX2, ssY1T, ssZ2F) (ssX2, ssY2T, ssZ2B) (ssX2, ssY2B, ssZ1B) wall id=54 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY1B, ssZ1F) (ssX1, ssY1T, ssZ2F) (ssX2, ssY1T, ssZ2F) (ssX2, ssY1B, ssZ1F) wall id=55 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY1B, ssZ1F) (ssX2, ssY1B, ssZ1F) (ssX2, ssY2B, ssZ1B) (ssX1, ssY2B, ssZ1B)

;===============================================================建上方延伸牆

wall id 56 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY1T, ssZ2F) (ssX1, ssY2T, ssZ2B) (ssX1, ssY2T, boxZ2) (ssX1, ssY1T, boxZ2) wall id 57 kn=1e12 ks=1e12 f=0.5 face (ssX2, ssY1T, ssZ2F) (ssX2, ssY1T, boxZ2) (ssX2, ssY2T, boxZ2) (ssX2, ssY2T, ssZ2B) wall id 58 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY1T, ssZ2F) (ssX1, ssY1T, boxZ2) (ssX2, ssY1T, boxZ2) (ssX2, ssY1T, ssZ2F) wall id 59 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY2T, ssZ2B) (ssX2, ssY2T, ssZ2B) (ssX2, ssY2T, boxZ2) (ssX1, ssY2T, boxZ2)

;===============================================================群組中央軟岩 group ShBallCenter owner ShBall range x ssX1 ssX2 z Z0 shZ2

;===============================================================分層中央軟岩 call 03_R3_LayerShBall.txt

;===============================================================生砂岩球

generate id=50000, 99999 rad=0.08, 0.08 x=ssX1, ssX2 y=ssY1T, ssY2T z=ssZ2B, boxZ2

;===============================================================砂岩群組、給霣降性質 group SsBall range id=50000, 99999

property den=3700 fric=0.1 kn=1e8 ks=1e8 range group SsBall

;===============================================================設定輸出畫面 plot add ball range group ShBallCenter blue ac off

plot add displacement range group SsBall magenta scale 0.02 plot add displacement range group ShBallCenter blue scale 0.02 plot set rotation (0, 0, 270)

plot set size 10

plot set center (ssX1, Y0, shZ2) set plot jpg quality 2 size 2048 1536

;===============================================================霣降砂岩 set dt 0.01

set g 0 0 -25 cyc 110000

;===============================================================建抬升需要之牆、刪除13 號牆

wall id 60 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY2T, ssZ2B)(ssX1, ssY1T, ssZ2F)(ssX2, ssY1T, ssZ2F)(ssX2, ssY2T, ssZ2B) zv=0.3 wall id 61 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY2T, ssZ2B)(ssX2, ssY2T, ssZ2B)(ssX2, ssY1T, ssZ2F)(ssX1, ssY1T, ssZ2F) wall id 62 kn=1e12 ks=1e12 f=0.5 face (ssX1, ssY2T, boxZ2)(ssX2, ssY2T, boxZ2)(ssX2, ssY1T, boxZ2)(ssX1, ssY1T, boxZ2) del wall 13 del balls range z ssZ2F boxZ2

cyc 10000

;===============================================================存檔，軟岩含強度 save 07001200_ss10.sav

;===============================================================檔名為軟岩檔名_砂岩編號

;########################### END 03 ###########################

7. 主程式Part 4：

;########################### PART 04 ###########################

;===============================================================調整砂岩性質 property den=3700 fric=0.5 kn=1e8 ks=1e8 range group SsBall

;===============================================================補強軟岩鍵結 property pb_kn=PBkn pb_ks=PBks pb_n=PBns pb_s=PBss pb_r=PBr range group ShBall

;===============================================================改牆摩擦力 property pb_kn=1e12 pb_ks=1e12 pb_n=1e12 pb_s=5e11 pb_r=8e-2 range group SsBall cyc 20000

del wall 50 52 53 56 57 58 59 60 61 62 cyc 20000

del wall 54 cyc 10000

clump id 1 range group SsBall del wall 8 10 54 55

;===============================================================跑，前10000 步每1000 步輸出 set plot jpg quality 2 size 2048 1536

set dt 0.0001 set g 0 0 -9.81

;===============================================================設定存檔前綴名 def FileDate

filedate = '225-03292200_4206s_' end

plot hardcopy file fileName plot close

cyc 1000

plot hardcopy file fileName plot close

cyc 5000