PFC落石模拟

Landslide/Rockfall simulation

---

山体滑坡/落石模拟

https://www.youtube.com/watch?v=WSa3909qYmI

        模拟的目的在于通过导入团块的对象文件产生团块的二进制输出。

具体措施：

(i) 使用导入的几何体形成团块模板
(ii) 使用Taghavi(2011)定义的BubblePack算法来生成卵石分布
(iii) 使用Surfcalculate关键字从表面描述中计算惯性属性，假设密度均匀。
(iv) 理解为比率和距离提供的数值所产生的差异。

clump template + create keyword

①bubblepack关键字
通过Taghavi（2011）的Bubble Pack算法，指定创建团块模板卵石。

        距离fdistance：根据Taghavi（2011）的定义，与光滑度的角度测量相对应的距离fdistance为0 < fdistance < 180。fdistance越大，卵石分布越平滑。

        比率fratio ：块状模板中保留的最小和最大卵石的比率，0 < fratio < 1。

②surfcalculate关键字

表面计算，假设密度均匀，从面的描述中计算出惯性属性。

clump template + export s <keyword ...>
将块状模板s导出到文件。默认的文件名是{"clump.p2clp"（2D）；"clump.p3clp"（3D）}。如果没有指定nothrow，那么如果一个文件存在，就会被覆盖，否则就会出现错误。这是一个二进制文件，可以用import关键字重新加载。

        skip-errors：表示如果存在指定的名为s的前一个文件，那么它将被覆盖。

model domain condition periodic

应用周期性边界条件。当球或团块的中心点落在模型域外时，它们会被传回模型的另一侧。

model domain condition destroy

销毁，删除落在模型域之外的球、团块和壁面，这些球、团块和壁面的单个外延（即紧紧包括物体的轴对齐的边界框）。

clump distribute 

①porosity关键字
最终的孔隙率被设置为fporos。默认情况下，fporos={二维为0.16；三维为0.359}，对应于{二维为盘子；三维为球体}的紧密单分散包装的孔隙率。

②bin ibin
指定bin ibin的分布属性。可以指定任何数量的分布的半径范围、体积分数和分布类型。所有分布的体积分数之和必须为1。

        volume-fraction fvfrac：这个分布中的团块的体积分数。所有分布的体积分数之和必须是1。

clump distribute porosity 0.15 diameter number-bins 6 box -30 30 -30 30 -30 30 ...
        bin 1 size 16    15  volume-fraction 0.03  group 'clump1' template 'clump1' ...
        bin 2 size 15    14  volume-fraction 0.12  group 'clump2' template 'clump2' ...
        bin 3 size 14    13  volume-fraction 0.30  group 'clump3' template 'clump3' ...
        bin 4 size 13    12  volume-fraction 0.4   group 'clump4' template 'clump4' ...
        bin 5 size 12    11  volume-fraction 0.12  group 'clump5' template 'clump5' ...
        bin 6 size 11    10  volume-fraction 0.03  group 'clump6' template 'clump6'

; Start with a new file and provide title to the project
; Set a random number so that after multiple runs of the same model, produced results are similar.
model new
model random 10001
model title 'Debris_Flow' 
 
;Aim 1: Produce binary outputs of clumps by importing the object files of clumps

;Action 1: Import object files from the directory
    geometry import 'geo1.stl'
    geometry import 'geo2.stl'
    geometry import 'geo3.stl'
    geometry import 'geo4.stl'
    geometry import 'geo5.stl'
    geometry import 'geo6.stl'


clump template create name 'trial1' geometry 'geo1' surfcalculate bubblepack ratio 0.3 distance 100
clump template create name 'trial2' geometry 'geo1' surfcalculate bubblepack ratio 0.3 distance 50
clump template create name 'trial3' geometry 'geo1' surfcalculate bubblepack ratio 0.3 distance 150

clump template create name 'trial4' geometry 'geo1' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'trial5' geometry 'geo1' surfcalculate bubblepack ratio 0.2 distance 50
clump template create name 'trial6' geometry 'geo1' surfcalculate bubblepack ratio 0.2 distance 150

model save 'state1'
;Action 3: Find a suitable set and create clump templates
clump template create name 'clump1' geometry 'geo1' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'clump2' geometry 'geo2' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'clump3' geometry 'geo3' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'clump4' geometry 'geo4' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'clump5' geometry 'geo5' surfcalculate bubblepack ratio 0.2 distance 100
clump template create name 'clump6' geometry 'geo6' surfcalculate bubblepack ratio 0.2 distance 100


;Action 4: Export each clump template to a binary file and issue nothrow/skip-errors keyword to overwrite pre-esxiting files.
    clump template export 'clump1' nothrow 
    clump template export 'clump2' nothrow 
    clump template export 'clump3' nothrow 
    clump template export 'clump4' nothrow 
    clump template export 'clump5' nothrow 
    clump template export 'clump6' nothrow 
model save 'state2'
;Aim 1 Achieved

model new
model random 10001
;Aim 2: Use brick logic to from a brick from clumps that makes the reproduction and computation fast.

;Action 1:Import the clump templates from the binary files
    clump template import 'clump1' 
    clump template import 'clump2' 
    clump template import 'clump3' 
    clump template import 'clump4' 
    clump template import 'clump5' 
    clump template import 'clump6' 

;put them in a periodic cubic box and allow them to stablize
model domain extent -40 40 condition periodic

;Decide a contact model for future contact detection
contact cmat default model 'linear' ...
                         method deformability   emod 2e9 kratio 2.0 ...
                         property fric '0.1' dp_nratio 0.2

;The error was reported because the volume was typed as avolume
;Provide density and initial damping to Produce a mass of the clumps and provide the fractional content cumulative to 1
clump distribute porosity 0.15 diameter number-bins 6 box -30 30 -30 30 -30 30 ...
        bin 1 size 16    15  volume-fraction 0.03  group 'clump1' template 'clump1' ...
        bin 2 size 15    14  volume-fraction 0.12  group 'clump2' template 'clump2' ...
        bin 3 size 14    13  volume-fraction 0.30  group 'clump3' template 'clump3' ...
        bin 4 size 13    12  volume-fraction 0.4   group 'clump4' template 'clump4' ...
        bin 5 size 12    11  volume-fraction 0.12  group 'clump5' template 'clump5' ...
        bin 6 size 11    10  volume-fraction 0.03  group 'clump6' template 'clump6'
clump attribute density 3000.0 damp 0.9


;To form a grain size distribution curve create a measure region inside the box
;tabulate the particles according to their size in their respective bins with minimum and maximum extent
; Dump the table to form a table chart
measure create id 1 radius 29 ...
                        bins 100 12 20
measure dump id 1 table 'numerical' 

; Create a model history to analuze the disturbance in the model specimen.
model history id 1 mechanical ratio-average

;use timescaling to achieve equillibrium quickly after it set timestep to auto mode
model mechanical timestep scale
model cycle 250 calm 100
model solve ratio-average 1e-3
model cycle 250 calm 50
model solve ratio-average 1e-4
model cycle 250 calm 25
model solve ratio-average 1e-5
model mechanical timestep auto

; Create and export a brick from the stablized clump heap. It can be utilized at later stage in replication 
brick make id 1
brick export id 1 filename 'brick1' nothrow
model save 'state3'


model new
model random 10001
model domain extent -1100 1100 -1100 1100 -300 900 condition destroy destroy destroy 

contact cmat default model 'linear' ...
                    method deformability emod 2e10 kratio 1.0   ...
                    property fric '0.5' dp_nratio 0.2

;Import the slope/mountain as a wall
    geometry import 'slope.stl'
    wall import from-geometry 'slope' skip-errors id 100
    model gravity 9.81

;Import the brick and assemble it, the origin is the bottom-most point from where the brick is replicated in axial directions.
brick import id 1 filename 'brick1' 
brick assemble id 1 origin (-800,700,575) size (2,2,2)
brick assemble id 1 origin (25,625,250) size (2,2,2)


;Group the clumps which are not required
clump group 'del1' slot '2' range plane origin (-671.673,780.009,620.539) normal (0.220853,-0.560137,0.798418) below
clump group 'del2' slot '2' range plane origin (-745.333,735.558,620.036) normal (0.419921,-0.366438,0.830295) below
clump group 'del3' slot '2' range plane origin (-671.673,780.009,680) normal (0,0,1) above

;Remove some of the clumps
clump delete range group 'del1' slot '2'
clump delete range group 'del2' slot '2'
clump delete range group 'del3' slot '2'
clump group 'debris' slot '2' range position-x -1100 1100

;Form boxe to contain and stablize the initial debris
wall generate box (-850,-600) (650,900) (550,750) one-wall 
wall generate box (5,215) (600,850) (150,500)  one-wall

;To avoid model interference and better results negate the body based local damping
clump attribute damp 0.0
model clean all 

;Settle the clumps in the box
model history id 1 mechanical ratio-average
model mechanical timestep scale
model cycle 500 calm 100
model solve ratio-average 1e-3
model cycle 750 calm 50
model solve ratio-average 1e-4
model cycle 1000 calm 25
model solve ratio-average 1e-5
model mechanical timestep auto
model save 'state4'
;Change the frictional properties of the slope
wall delete walls range set id 101 102
wall property 'fric' 0.2
    
wall property 'fric' 0.4 range position-y -1000 280
wall property 'fric' 0.4 range position-x -230 2000
    
wall property 'fric' 0.6 range position-y -1000 -150
wall property 'fric' 0.6 range position-x  260 2000

wall property 'fric' 0.4 range id 2595
wall property 'fric' 0.4 range id 2678
wall property 'fric' 0.4 range id 2772
wall property 'fric' 0.4 range id 1552
wall property 'fric' 0.4 range id 1494
wall property 'fric' 0.4 range id 1442
wall property 'fric' 0.4 range id 1392
wall property 'fric' 0.4 range id 1355
wall property 'fric' 0.4 range id 1314
wall property 'fric' 0.4 range id 1289
wall property 'fric' 0.4 range id 1270
wall property 'fric' 0.4 range id 1252
wall property 'fric' 0.4 range id 4886
wall property 'fric' 0.4 range id 4888

;Add the particle traces
clump trace id 1 (-610.013,738.36,585.659)
clump trace id 2 (-698.949,845.282,679.247)
clump trace id 3 (-602.018,871.206,669.083)

clump trace id 4 (156.198,629.542,381.03)
clump trace id 5 (160.236,743.022,409.412)
clump trace id 6 (43.7976,687.653,398.846)

model save 'initiation'

; Use result logic to save miniature saved state which can be used later to generate image of plots.
model result interval mechanical 0.2 warn off
clump result active on add-attribute group add-attribute velocity add-attribute energy-kinetic
wall result active on add-attribute group
trace result active on
model cycle 1
model save 'initiation' 

model cycle 1
model solve
model save 'deposition'

program return

fish define Make_Movie(dur,inc,name)
  i = 1
  curv = inc
  dur = dur
  namefile = name
  loop while (curv <= dur)
    i = i + 1
    local fname = string.build('%1_%2.png',namefile,i)
        command
          model result import @i skip-fish
          plot 'gro2' export bitmap size 1920 1080 filename @fname
        endcommand
    curv = curv + inc
  endloop
end
@Make_Movie(284.6,0.2,'velocity')