10.1. 一个简单的例子

DEMPM 建模脚本绝大部分与 MPM2D 相同，用户可在熟悉 MPM2D 操作后回到本章。 MPM2D 教程中 双轴示例 将被借用于此，用来展示如何进行多尺度模拟。

 1app.core.setAnalysisType(AnalysisType="MPMDEM")
 2import mpm as task
 3sim=app.core.sim
 4
 5# 此处将会到处DEM RVE求解器
 6import pygodem
 7
 8height = 2.0
 9elemsize = 0.1
10
11sim.grid.setRange(0, 3, 0, 1+height)
12sim.grid.setNumEle(30, 10+int(height/elemsize))
13sim.grid.printInfo()
14
15particle_damping = 0.2
16grid_damping = 0
17picfrac = 1.0
18
19timestep = 1e-3
20confiningStress = 1e5
21
22# sim内部RVE求解器材料设置
23sim.godem.kn = 1e6
24sim.godem.ks = 1e6
25sim.godem.mu = 0.5
26
27# 多尺度模式下仅设置虚材料
28mat = task.Ghostmat()
29mat.rho = 1.
30
31sim.mats=[mat]
32sim.mats[0].setIndex(0)
33
34pic = 2
35left = task.Rectangle(1.0, 1.05, 0.3, 0.3+height, task.MPTTRACSHAPE)
36left.setConsolidation(-confiningStress, -confiningStress)
37left.setTracBc(4, task.NORMALDIR, -confiningStress)
38left.setMatProperty(task.MCMAT, 0, pic)
39left.setDamping(particle_damping, grid_damping, picfrac)
40
41right = task.Rectangle(1.95, 2, 0.3, 0.3+height, task.MPTTRACSHAPE)
42right.setTracBc(2, task.NORMALDIR, -confiningStress)
43right.setConsolidation(-confiningStress, -confiningStress)
44right.setMatProperty(task.MCMAT, 0, pic)
45right.setDamping(particle_damping, grid_damping, picfrac)
46
47bot = task.Rectangle(0.5, 2.5, 0.2, 0.3, task.RIGIDSHAPE)
48bot.setBc(True, 0, True, 0)
49bot.setMatProperty(task.RIGIDMAT, 0, 2)
50
51top = task.Rectangle(0.5, 2.5, 0.3+height, 0.4+height, task.RIGIDSHAPE)
52top.setBc(True, 0, True, -0.02)
53top.setMatProperty(task.RIGIDMAT, 0, 2)
54
55soil = task.Rectangle(1.05, 1.95, 0.3, 0.3+height, task.NONEBC)
56soil.setConsolidation(-confiningStress, -confiningStress)
57soil.setMatProperty(task.MCMAT, 0, pic)
58soil.setDamping(particle_damping, grid_damping, picfrac)
59
60sim.shapes = [left, right, top, bot, soil]
61
62# 读入RVE初始文件
63sim.setRvePath('./myrve.gst')
64
65sim.setTimeDt(30, timestep)
66sim.writer.setDir("./Results", "biax")
67sim.writer.setStartTime(0)
68sim.writer.setEndTime(30)
69sim.writer.setIntervalTime(0.4)
70sim.setMethod(task.GIMP)
71sim.simTypeMasks = task.MANAL
72sim.PreAnalysis()
73sim.setNumThreads(1)

在此脚本中，DEMPM 脚本仅在包的导入，虚材料设置以及 setRvePath 与 MPM2D 不同，其余设置完全相同:

  • import pygodem: 激活 GODEM 的Python上下文环境，此命令执行后，便可进行 sim.godem.xxx 系列操作

  • Ghostmat: 在多尺度环境下，无需 MPM2D 提供材料本构模型，仅有虚材料 mpm.Ghostmat 被调用。

  • setRvePath: 用户需导入 gst 初始化文件进行多尺度模拟，gst文件的生成请参见操作指令详细解读。

10.2. RVE制备

RVE的构成如 图 10.1 所示，RVE由离散元颗粒集合体构成，颗粒在MPM提供的周期边界条件下运动，并反馈给MPM相应的应力以及其他物理场状态量。

../_images/nonsphrve.png

图 10.1 DEM RVE示意图

不同类型RVE特点如下：

  • 圆盘RVE：运行效率高，但无法准确捕获材料的各向异性力学特性。

  • 非球RVE：运行效率较低，可在多尺度模拟中考虑颗粒形状，进而反映颗粒的非共轴性以及各向异性等，以便达到贴合实际工况的目的。

所有与多尺度RVE相关的命令都可通过调用 pygodem 完成，如下脚本展示了如何创建RVE的gst文件：

 1# pySudoMath 封装了一系列数学容器，例如矢量以及矩阵
 2from pySudoMath import *
 3
 4app.core.setAnalysisType(AnalysisType="GODEM")
 5sim=app.core.sim
 6
 7# pySudoMath的Vector2f二维矢量
 8boxsize = Vector2f(10.0,10.0)
 9boxsize *= 0.025
10
11# 颗粒半径范围
12radius_range = Vector2f(0.25e-2, 0.5e-2)
13# 目标生成颗粒数量
14par_num = 400
15
16# 定义颗粒接触的相关材料参数：阻尼，法向刚度，切向刚度
17# 密度，z轴的厚度，摩擦系数，时步
18sim.damping = 0.3
19sim.kn = 1e6
20sim.ks = 1e6
21sim.rho = 2650.0
22sim.z_dim = 0.075
23sim.mu = 0.001
24sim.dt = 1e-5
25
26# 颗粒随即位置和半径生成坐标及半径数组
27pos_r = sim.generatePacking(boxsize, radius_range, par_num)
28
29# 利用坐标和半径数组初始化RVE上下文
30sim.generateRVEs(pos_r, boxsize)
31
32# 设定等向固结边界
33sim.consolSetting(1e5)
34
35# 开始固结
36app.run(10000)
37
38# 固结完毕后保存输出gst文件
39sim.save('sph_100.gst')

依此法便可生成如 图 10.1 (a) 所示的圆盘RVE初始文件。

备注

generatePacking 是 GODEM 为用户提供的在矩形空间内随机生成不重叠的圆盘颗粒的接口，它返回的是一个大小等于目标颗粒数量的三维矢量，该矢量分别包含颗粒坐标x，y以及颗粒半径。用户可自定义三维数组来替换该函数，例如：

def regularPacking(r,w,h):
   pos_r = []
   for j in range(h):
       for i in range(w):
           x = 2.0*r*i + r
           y = 2.0*r*j + r
       pos_r.append(Vector3f(x,y,r))
   return pos_r

代码为均匀正交堆积颗粒生成示例。

在RVE通过 generateRVEs 生成后，用户指定相应围压使得RVE固结进入各向等压状态，在此过程中用户可通过改变摩擦系数 sim.mu 来控制最终试样孔隙比。固结完成后，利用 save 命令保存相应试样，以便于多尺度运行调用。

10.3. RVE求解器设置

在 MPM2D 脚本中，用户可能需要对DEM求解器进行一系列设置，以下为常用命令。

10.3.1. 非粘结接触

  • 接触法向刚度 sim.godem.kn ：

app.core.setAnalysisType(AnalysisType="MPMDEM")
sim=app.core.sim
sim.godem.kn = 1e5

  • 接触切向刚度 sim.godem.ks ：

sim.godem.ks = 1e5

  • 摩擦系数 sim.godem.mu ：

sim.godem.mu = 0.5

10.3.2. 粘结接触

  • 粘结法向强度 sim.godem.bondNormalStrength ：

sim.godem.bondNormalStrength = 1e7

  • 粘结切向强度 sim.godem.bondShearStrength ：

sim.godem.bondShearStrength  = 6e6

  • 胶结强度劣化系数 sim.meltingTs ：

sim.godem.bondShearStrength  = 10

  • 模拟工况设置：

sim.simTypeMasks = mpm.MANAL | mpm.BONDMODEL

表示考虑胶结接触的力学模拟。

10.4. 热力耦合颗粒柱溃散模拟

本节给出多尺度模拟的接触劣化的颗粒柱溃散：

 1app.core.setAnalysisType(AnalysisType="MPMDEM")
 2import mpm as task
 3sim=app.core.sim
 4
 5sim.grid.setRange(0, 0.6, 0, 0.15)
 6sim.grid.setNumEle(120, 30)
 7sim.grid.printInfo()
 8
 9particle_damping = 0.2
10grid_damping = 0
11picfrac = 1.0
12pic = 2
13dt = 5e-5
14
15Mat1 = task.Ghostmat()
16Mat1.rho = 1
17
18Mat = task.Ghostmat()
19Mat.rho = 2650.
20Mat.setThermalProp(0.5, [500, 500])
21Mat.Kn = 1e7
22Mat.Ks = 1e7
23Mat.friction = 0.15
24
25sim.mats=[Mat1, Mat]
26sim.mats[0].setIndex(0)
27sim.mats[1].setIndex(1)
28
29left = task.Rectangle(0., 0.005, 0.01, 0.15, task.RIGIDSHAPE)
30left.setBc(task.FIXMINUS, 0, False, 0)
31left.setMatProperty(task.RIGIDMAT, 0, 1)
32# 多尺度中温度边界诱发接触劣化
33left.setThermalBc(100)
34
35right = task.Rectangle(0.205, 0.210, 0.01, 0.15, task.RIGIDSHAPE)
36right.setBc(task.FIXPLUS, 0, False, 0)
37right.setMatProperty(task.RIGIDMAT, 0, 1)
38
39pts = [0, 0.005, 0.6, 0.005, 0.6, 0.01, 0., 0.01]
40foot = task.DemPolygon(dt, 10, 50, 5e5, 5e5, 0.15)
41foot.setPts(pts)
42foot.usePolygonContact(False)
43foot.addForce2Pt(False)
44foot.setBc(True, 0, True, 0)
45foot.setMatProperty(task.RIGIDMAT, 0, 1)
46
47soil = task.Rectangle(0.005, 0.205, 0.01, 0.11, task.GRAVITYSHAPE)
48soil.setBc(False, 0, True, -9.8)
49soil.setMatProperty(task.MCMAT, 1, pic)
50soil.setDamping(particle_damping, grid_damping, picfrac)
51
52sim.shapes = [left, right, foot, soil]
53
54# 多尺度RVE路径设置，摩擦系数设置以及粘结接触强度设置
55sim.setRvePath('./rve10kpa.gst')
56sim.godem.mu = 0.5
57sim.godem.bondNormalStrength = 1e6
58sim.godem.bondShearStrength = 1e6
59
60sim.setTimeDt(20, dt)
61sim.writer.setDir("./Results", "collapse")
62sim.writer.setStartTime(0)
63sim.writer.setEndTime(20)
64sim.writer.setIntervalTime(0.05)
65
66# 此时仅有力学模拟
67sim.setMethod(task.GIMP)
68sim.simTypeMasks = task.MANAL
69
70sim.PreAnalysis()
71sim.setNumThreads(1)
72
73def gravityHook(steps = 10000):
74    partg = -9.8*sim.iter()/float(steps)-0.098
75    app.print('Current g:', partg)
76    soil.setBc(False, 0, True, partg)
77
78pyhook1 = task.PyHook()
79pyhook1.command = 'gravityHook()'
80pyhook1.reset(iterPeriod = 100, totalRuns = 100)
81pyhook1.dead = False
82sim.hooks = [pyhook1]
83
84app.run(10000)
85
86# 固结完毕后，移除右侧挡板，启动多尺度热力耦合模拟
87right.setBc(False, 0, False, 0)
88# 工况为热力以及粘结接触模拟
89sim.simTypeMasks = task.MANAL | task.BONDMODEL | task.TANAL
90# 设定粘结接触的线性劣化温度
91# 当温度低于10度时强度开始劣化，直至0度完全退化为非粘结接触
92sim.meltingTs = 10.
93soil.setDamping(0.05, 0., 0.05)
94app.Run(10000)

多尺度模拟与MPM模拟相比，多了控制工况为 mpm.BONDMOEDEL ，该工况为RVE粘结接触工况。与之对应，用户需使用 sim.godem.bondNormalStrength 和 sim.godem.bondShearStrength 来指定粘结接触强度，并通过 sim.meltingTs 指定粘结接触劣化开始温度，默认终止温度为0摄氏度，粘结强度在此区间为线性劣化。

10.5. 热力相变诱发边坡失稳模拟

本例子给出热力相变接触劣化的边坡失稳案例，

../_images/slope.png

图 10.2 热力相变诱发边坡失稳示意图

  1app.core.setAnalysisType(AnalysisType="MPMDEM")
  2import mpm as task
  3sim=app.core.sim
  4
  5width = 50.
  6height = 20.
  7elemsize = 0.25
  8numx = int(width/elemsize)
  9numy = int(height/elemsize)
 10
 11sim.grid.setRange(0, width, 0, height)
 12sim.grid.setNumEle(numx, numy)
 13sim.grid.printInfo()
 14
 15pic = 2
 16particle_damping = 1.0
 17grid_damping = 1.0
 18picfrac = 1.0
 19
 20mat1=task.Ghostmat()
 21mat1.rho=1
 22
 23mat2 = task.Ghostmat()
 24mat2.rho = 1000
 25
 26# 设置相变参数-同MPM设置
 27mat2.setInitFraction([0.5697, 0.188735, 0.241565])
 28mat2.setInitDensity([2650., 1000., 1000.])
 29mat2.setThermalProp([900, 2100, 4180], [252288., 252288., 193536., 193536., 48384., 48384.])
 30mat2.freezepoint = 0.
 31mat2.latentratio = 3.33e5
 32mat2.unfrozenlimit = 0.058
 33mat2.unfrozenratio = 0.16
 34
 35# 多尺度接触强度设置
 36sim.godem.kn = 1e6
 37sim.godem.ks = 1e6
 38sim.godem.mu = 0.3
 39sim.godem.bondNormalStrength = 1e7
 40sim.godem.bondShearStrength = 3.25e6
 41
 42sim.mats=[mat1, mat2]
 43sim.mats[0].setIndex(0)
 44sim.mats[1].setIndex(1)
 45
 46# 节点热源，此处使用节点边界NDSHAPE
 47midnodeheat = task.Rectangle(5.4, 5.6, 6, 14.3, task.NDHEATSHAPE)
 48midnodeheat.setThermalBc(50)
 49
 50left = task.Rectangle(elemsize, 2.*elemsize, 2.*elemsize, 16, task.RIGIDSHAPE)
 51left.setBc(1, 0, False, 0)
 52left.setThermalBc(-5)
 53left.setMatProperty(task.RIGIDMAT, 0, 1)
 54left.setDamping(particle_damping, grid_damping, picfrac)
 55
 56right = task.Rectangle(40.+2.*elemsize, 40.+3.*elemsize, 2.*elemsize, 16, task.RIGIDSHAPE)
 57right.setBc(1, 0, False, 0)
 58right.setThermalBc(-5)
 59right.setMatProperty(task.RIGIDMAT,0, 1)
 60right.setDamping(particle_damping, grid_damping, picfrac)
 61
 62bot = task.Rectangle(elemsize, 50, elemsize, 2.*elemsize, task.RIGIDSHAPE)
 63bot.setBc(True, 0, True, 0)
 64bot.setThermalBc(-5)
 65bot.setMatProperty(task.RIGIDMAT, 0, 1)
 66
 67pts2 = [2.*elemsize, 2.*elemsize,
 68        40+2.*elemsize, 2.*elemsize,
 69        40+2.*elemsize, 4.+2.*elemsize,
 70        20+2.*elemsize, 4.+2.*elemsize,
 71        10+2.*elemsize, 14+2.*elemsize,
 72        2.*elemsize, 14+2.*elemsize]
 73soil = task.Polygon(task.GRAVITYSHAPE)
 74soil.setPts(pts2)
 75soil.setMatProperty(task.MCMAT, 1, pic)
 76soil.setDamping(particle_damping, grid_damping, picfrac)
 77
 78sim.shapes = [left, midnodeheat, right, bot, soil]
 79
 80dt = 0.001
 81sim.setTimeDt(8, dt)
 82sim.writer.setDir("./Results", "slope")
 83sim.writer.setStartTime(0)
 84sim.writer.setEndTime(20000)
 85sim.writer.setIntervalTime(2000*dt)
 86
 87sim.setMethod(task.GIMP)
 88# 阶段一使用力学模拟
 89sim.simTypeMasks = task.MANAL
 90
 91sim.PreAnalysis()
 92sim.setNumThreads(1)
 93
 94# 设置边坡初始温度
 95soil.setInitTemperature(-5)
 96
 97# 重力固结分阶段二和阶段三
 98# 阶段二在重力固结30%阶段完成，期间RVE接触为非粘结
 99# 阶段三在重力固结30%-100%，此时RVE激活粘结接触
100addbond = False
101def gravityHookFunc(steps = 10000):
102    global addbond
103    partg = -9.8*sim.iter()/float(steps)
104    app.print('Current g:', partg)
105    soil.setBc(0,0,1,partg)
106    if((not addbond) and abs(partg)/9.8 >= 0.3):
107        addbond = True
108        sim.simTypeMasks = task.MANAL | task.BONDMODEL
109pyhook1 = task.PyHook()
110pyhook1.command = 'gravityHookFunc()'
111pyhook1.reset(iterPeriod = 100, totalRuns = 100)
112pyhook1.dead = False
113sim.hooks = [pyhook1]
114
115# 重力固结完毕后，重力应力有粘结应力的贡献
116app.run(10000)
117# 激活相变和热模拟，相变接触劣化无需设置 meltingTs
118sim.simTypeMasks = task.MANAL | task.PHASEANAL | task.BONDMODEL | task.TANAL
119soil.setDamping(0.1, 0.1, 0.01)
120app.run(100000)