6.3. 碎屑流

本算例展示如何使用CuDEM模拟滑坡碎屑流。如 图 6.3 所示，地形通过数字高程模型构建，而滑动体则由滑动前后的数字高程数据之差求得。本例子我们直接提供了相应数据，如果用户需要更换地形等，则需要自行预先处理相应数据。

../_images/debrisflowsetup.png

图 6.3 碎屑流模型构建：地形由12万三角面片组成，滑动体由160万颗粒组成

6.3.1. 导入基础包

导入pycudem和pySudoMath:

1import sys,os,math,random
2import numpy as np
3from dem.core import *
4from dem.cudem import *
5from pySudoMath import *
6# 设置随机种子
7random.seed(1234343)

6.3.2. 设置GUI或命令行运行

为了让脚本既能在GUI中运行，又能直接通过命令行运行，我们做如下设置：

1gui = True
2sim = None
3if gui:
4    app.core.setAnalysisType(AnalysisType="CUDEM")
5    sim=app.core.sim
6else:
7    sim = DEMSim()

6.3.3. 通用设置

1#########general setting##########
2sim.setOutputDir("./")
3sim.setFilePrefix("testmesh_")
4#################Digital elevation mesh
5#put your stl file here for the terrain
6if gui:
7    app.frame.scene.openSTL("./tmp1.stl")

6.3.4. 设置地形模型

  • 设置数字高程

1de = DigitalElevation()
2data = np.loadtxt("dc7_dem_1411slide.txt",  skiprows = 6)  - 1597.4
3data = np.flip(data,0)
4de.loadData2d(data)
5de.setCellSize(1.0)

  • 设置地形材料属性

1demat = Material()
2demat.kn = 1e7
3demat.ks = 1e7
4demat.friction = 0.5
5demat.density = 1000.0
6demat.vdamping = 0.5
7de.material = demat

  • 加入到模拟

1sim.digitalElevation = de

6.3.5. 设置滑动体

  • 定义滑动体颗粒函数

 1def PointSpawn(points, radius = 0.25, spawn = True):
 2    pos = list()
 3    par_num = 0
 4    for p in points:
 5        delta = [-radius, radius]
 6        par_num += 1
 7        if not spawn:
 8            pos.append(Vector3f(p[0],p[1],p[2]))
 9        else:
10            for i in range(2):
11                for j in range(2):
12                    for k in range(2):
13                        pos.append(Vector3f(p[0]+delta[i],p[1]+delta[j],p[2]+delta[k]))
14    return pos

  • 生成滑动体颗粒信息

1psize = 0.125#0.25
2aabb_centers = np.loadtxt("../pos_list.txt")
3aabb_centers += np.array([115.5, 130.5-1,-1597.4]) #shift data
4pos = PointSpawn(aabb_centers, radius = 0.25)
5pos = PointSpawn(pos, radius = 0.125)
6
7rs_p = np.ones(len(pos))*psize
8rs = rs_p #np.hstack((rs_w,rs_p))

6.3.6. 设置模拟域

  • 设置模拟域大小

设置模拟域的AABB大小以及模拟域中分域大小（采用基于cell的算法需要此设置）：

1sim.domain.setAABB(0,230,0,260,0,1714-1597.4)#domain size should be large enough
2sim.domain.setCellSize(psize*4.0)

  • 设置墙作为模拟域边界

1sim.walls.addWall(Vector3f(0,0,0), 1, 2) #bottom
2# sim.walls.addWall(Vector3f(0,0,1714-1597.4), -1, 2) #top
3sim.walls.addWall(Vector3f(0,0,0), 1, 0) #left
4sim.walls.addWall(Vector3f(230,0,0), -1, 0) #right
5sim.walls.addWall(Vector3f(0,0,0), 1, 1) #front
6sim.walls.addWall(Vector3f(0,260,0), -1, 1) #back

  • 设置墙的材料，采用线弹性模型

1wmat = Material()
2wmat.kn = 1e7 #法向刚度
3wmat.ks = 1e7 #切向刚度
4wmat.friction = 0.5 #摩擦系数
5sim.walls.material = wmat #把材料赋值到所有墙

6.3.7. 设置默认材料

颗粒采用默认材料，参考墙的材料：

1sim.material.kn = 1e5
2sim.material.ks = 1e5
3sim.material.friction = 0.5
4sim.material.density = 2000.0
5sim.material.vdamping = 0.5

6.3.8. 求解器设置

  • 设置运动求解器

设置阻尼，重力加速度及 alpha 。其中， alpha 将会影响邻居列表的更新频率。

1sim.integrator.damping = 0.1
2sim.integrator.gravity = Vector3f(0.0, 0.0, -10.0)
3sim.integrator.alpha = 1.2

  • 设置邻居列表属性

颗粒运动超过 delta 倍颗粒粒径后更新邻居列表； 光追接触检测时放大的AABB大小为 ORTextend 。

1sim.nlist.delta = 1.2
2sim.nlist.ORTextend = 2.0*psize*1.2

  • 设置内存动态优化

设置颗粒重新内存从新排列频率 reorderInterval （以时步计）；设置是否开启光追接触计算 useORT ； 设置计算时间步长 dt 。

1sim.reorderInterval = 500 #2000
2sim.useORT = RTflag
3sim.dt = 1e-3

6.3.9. 生成颗粒

  • 生成颗粒位置和大小，本例子采用单一粒径颗粒：

1aabb_centers = np.loadtxt("../pos_list.txt")
2par_num = 0
3aabb_centers += np.array([115.5, 130.5-1,-1597.4]) #shift data
4pos = PointSpawn(aabb_centers, radius = 0.25)
5pos = PointSpawn(pos, radius = 0.125)
6
7rs_p = np.ones(len(pos))*psize
8rs = rs_p #np.hstack((rs_w,rs_p))

  • 根据颗粒位置和半径生成颗粒：

1sim.createParticles(pos,rs)

6.3.10. 运行程序

  • GUI中命令行

1app.run(10000)

  • 终端命令行

1sudosim -i -f packing.py

然后在Python解释器中运行：

1sim.RunPthread(10000)