9.1. 一个简单的例子

MPM2D 底层使用C++语言编写，使用Python语言封装，以便于用户更为快捷的学习和使用。本章假定用户已完成 SudoSimGUI 的安装。若用户对 SudoSimGUI 的安装仍有疑问，请参阅 SudoSimGUI 安装章节或访问我们的官方网址寻求帮助。

more info SudoSimLab官网: https://www.sudosimlab.com/

MPM2D 的安装及使用均可基于 SudoSimGUI ，本章将围绕 MPM2D 的程序架构及Python脚本建模展开，不再赘述 SudoSimGUI 的操作教程，如有需要请参见 此处。

列表 9.1 砂土自由落体
 1app.core.setAnalysisType(AnalysisType = "MPM")
 2import mpm as task
 3sim=app.core.sim
 4
 5grid = sim.grid
 6grid.setRange(0, 2.6, 0, 2.6)
 7grid.setNumEle(26, 26)
 8grid.printInfo()
 9
10mat = task.Mc()
11mat2.rho = 2650
12mat2.E = 1e7
13mat2.nu = 0.2
14mat2.c = 0.01
15mat2.phi = 20.
16mat2.psi = 0.1
17mat2.tension = 0.01
18mat2.infinistate = False
19sim.mats = [mat]
20sim.mats[0].setIndex(0)
21
22soil = task.Rectangle(0.85, 1.75, 0.3, 2.3, task.GRAVITYSHAPE)
23soil.setBc(False, 0, True, -9.8)
24soil.setMatProperty(task.MCMAT, 0, 2)
25sim.shapes = [soil]
26
27sim.setNumThreads(1)
28sim.setTimeDt(15., 1e-3)
29sim.setMethod(task.GIMP)
30sim.simTypeMasks = task.MANAL
31
32sim.PreAnalysis()
33app.run(1)

上文代码演示了如何使用 MPM2D 模拟砂土的自由落体运动：

  • 创建一个名为 sandgravity.py 文件，并将上述内容拷贝至该文件内

  • 运行 SudoSimGUI，并按照教程将 sandgravity.py 导入到 SudoSimGUI

  • 在 SudoSimGUI 右下角的命令行内，输入如下命令

列表 9.2 SudoSimGUI - sandgravity.py
>>> app.run(100)

9.2. 模型文件核心成员

文件 sandgravity.py 给出了 MPM2D 运行所必需的架构，

MPM2D
  |-- grid
  |-- mats
  |-- shapes
  |-- sim
       |-- NumThreads
       |-- TimeDt
       |-- Method
       |-- TypeMasks

如上所示的 MPM2D 核心成员均为Python对象，它们由底层C++对象抽象而来。

注意

用户无需关注纷繁复杂的C++对象构造、内存分配以及生命周期管理，只需了解对象的Python构造函数及其在 MPM2D 求解中所承担职责即可完成MPM模拟

MPM2D 核心成员：

  • grid

    存储背景网格的几何信息，包括网格单元和节点坐标、数量以及尺寸。

  • mats

    存储所有参与运算的材料模型及边界条件抽象材料，例如弹性本构模型、摩尔-库伦本构模型以及位移或温度边界抽象材料。

  • shapes

    用于不同材料或不同边界条件物质点的分类管理。MPM2D 所有物质点均属于特定的 shape 对象，该对象负责物质点的初始化生成、物质点边界条件的动态控制以及物质点间相互作用。

  • sim

    MPM2D 核心对象，为所有上述成员的管理者。sim 对象负责物质点全局模拟设置，包括模拟时步、形函数设置、模拟类型（力，热，热力或相变）以及并行设置。

9.3. 模型文件核心指令

本节将对 sandgravity.py 文件展开解析。

9.3.1. GUI 和 CLI 初始化

MPM2D 同时支持 GUI 模式与 CLI 模式。

>>> app.core.setAnalysisType(AnalysisType = "MPM")
>>> import mpm as task

如上为二者初始化代码：

  • GUI 模式下，该命令将建立 SudoSimGUI 与 MPM2D 通信，完成 MPM2D 在 SudoSimGUI 中的前后处理以及内存分配等。通常该步运行完毕后，用户可在 SudoSimGUI 中进行 MPM 文件导入 或者 直接建模。

  • CLI 模式下，用户通过直接导入 mpm 包来初始化 MPM2D。

>>> import mpm as task

9.3.2. 主控成员: sim

无论是 GUI 模式，还是 CLI 模式，用户均需导入 MPM2D 的主控制成员 sim 来完成建模及程序运行。

  • GUI 模式的 sim 初始化如下：

>>> sim=app.core.sim

  • CLI 模式的 sim 初始化如下：

>>> sim=task.Analysis()

9.3.3. 程序元信息

MPM2D 程序除主控成员 sim 外仍需要一系列元信息来辅助建模，元信息包括：

  • 从C++封装的成员构造函数信息：例如 grid 、mats 以及 shapes 等。

  • 从C++封装的成员类方法信息：例如 grid.setRange 等。用户可直接在Python端完成对C++底层程序的调用。

  • 材料/边界/模拟类型的枚举类信息：例如 task.MCMAT 表示摩尔-库伦弹塑性本构的枚举，task.GRAVITYSHAPE 表示重力边界条件的枚举，task.MANAL 表示力学模拟工况的枚举。

  • 还有其他方便用户前后处理，数据输出的枚举信息，详细请参见 此处。

用户可随时使用 help 命令在 SudoSimGUI 或 cmd 中查询主控成员 sim 以及元信息。例如，输入如下命令查询 grid 信息：

>>> help(sim.grid)

用户可在 SudoSimGUI 返回的帮助信息中获取 grid 的构造方法以及成员函数和成员变量。

Help on Grid in module mpm object:

class Grid(pybind11_builtins.pybind11_object)
 |  Method resolution order:
 |      Grid
 |      pybind11_builtins.pybind11_object
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(...)
 |      __init__(self: mpm.Grid) -> None
 |
 |  getNumEle(...)
 |      getNumEle(self: mpm.Grid) -> tuple
 |
 |      get the number of elements along x and y directions
 |
 |  printInfo(...)
 |      printInfo(self: mpm.Grid) -> None

9.3.4. 背景网格: grid

用户可通过Python赋值命令获取 sim.grid 。

grid = sim.grid

MPM2D 使用二维结构网格，所有网格具有相同的x轴或y轴尺寸，网格在x、y轴的覆盖范围由 grid.setRange(XSTART, XEND, YSTART, YEND) 来设置。

grid.setRange(0, 2.6, 0, 2.6)

网格的尺寸可通过指定x、y轴的单元数量间接设定，MPM2D 将根据 grid.setNumEle(XNUM, YNUM) 输入的网格数量自动计算求得网格尺寸。

grid.setNumEle(26, 26)

网格的坐标和尺寸设置完毕后，用户可通过 grid.printInfo 方法来打印 grid 内部信息，以此来检查模型设置是否正确。

grid.printInfo()

备注

事实上，绝大部分 MPM2D 的Python对象都被封装了 printInfo 方法，用户可尝试调用不同对象的 printInfo 来获取设置信息或者获取程序运行时信息。

9.3.5. 材料本构模型: mats

MPM2D 提供了一系列材料本构模型，包括线弹性本构、摩尔库伦本构、Drucker-Prager本构，修正剑桥本构等，本例子中使用的是摩尔库伦本构。

mat = task.Mc()
mat2.rho = 2650
mat2.E = 1e7
mat2.nu = 0.2
mat2.c = 0.01
mat2.phi = 20.
mat2.psi = 0.1
mat2.tension = 0.01
mat2.infinistate = False

注意

所有材料本构均以默认构造函数生成，用户可通过 help 函数查询不同本构模型的参数并进行相应设置。 MPM2D 额外提供了材料的 isfinistate 属性用于决定是否材料的应变更新采用 小应变 理论还是 有限应变 理论。前者采用欧拉应变表征材料应变，多用于小旋转变形问题，求解效率高；后者采用格林应变来表征应变，多适用于大旋转变形问题，求解效率较低。

若 infinistate = True 表明使用 有限应变，反之则为 小应变 。

sim.mats = [mat]

当所有材料初始化完毕后，需被添加到 sim.mats 方可被主程序调用，用户直接使用Python列表命令来向 MPM2D 添加材料。

sim.mats[0].setIndex(0)

用户需根据材料的添加顺序为每个材料设置 唯一 的索引，该索引将被 sim.shapes 对象用于绑定物质点与对应材料模型。

9.3.6. 形状: shapes

MPM2D 通过 sim.shapes 来管理所有物质点； shapes 负责定义物质点初始区域、物质点材料本构模型以及物质点和节点的边界条件施加。 MPM2D 包含基本的几何形状，例如矩形、圆盘、环形以及任意多边形，用户可通过几何形状来对具有不同属性或者边界的物质点进行分组。

本节例子中使用了矩形形状来生成物质点，同时该形状被施加 重力 边界条件。

soil = task.Rectangle(0.85, 1.75, 0.3, 2.3, task.GRAVITYSHAPE)

边界条件可分别在x和y轴单独施加，函数形参为 shape.setBc(isXSET, VALUEX, isYSET, VALUEY)。

soil.setBc(False, 0, True, -9.8)

矩形形状内物质点的材料为摩尔库伦本构，参数使用索引为0的材料对象，同时形状生成物质点的线密度为2。

soil.setMatProperty(task.MCMAT, 0, 2)

MPM2D 在执行计算前，将会统一为每个 shape 生成对应的物质点，物质点的密度取决于用户输入，若无输入则默认为1。该例子中，用户指定物质点线密度为2，这意味着程序将会对每个网格单元以 2 x 2 的正交格式进行剖分，并逐一判断所有子单元的中心是否处于 soil 的矩形形状内，并将所有被矩形包含的物质点添加到 MPM2D 内存中。

物质点的体积均等于子单元的面积 (0.0025 m 2)，物质点的质量由索引为0的材料对象的密度参数确定 (6.625 kg)，物质点的初始速度、初始应力等状态量由边界条件管理对象通过 shape 施加，详情参见 边界条件 。

sim.shapes = [soil]

注意

所有 shape 对象定义和设置完毕后，用户直接通过Python列表操作将所有 shape 对象添加到 sim.shapes ，无需同 sim.mats 一样额外设置索引。

9.3.7. 模拟设置

当 MPM2D 求解所需的网格、物质点以及边界条件设定完毕后，在计算开始前，用户需进行模拟全局参数及工况设置。

  • MPM2D 支持 OpenMP 并行计算，用户通过 setNumThreads 指定计算线程 (至少为1)。

sim.setNumThreads(1)

  • 设定模拟的总时间以及时步。

sim.setTimeDt(15., 1e-3)

  • 设定插值形函数，目前 MPM2D 支持 task.LINEAR 、 task.GIMP 以及 task.CPDI 三种插值格式，分别对应 线性 、GIMP 以及 LINEARCPDI 形函数。

sim.setMethod(task.GIMP)

  • 设定模拟工况，目前 MPM2D 支持 task.MANAL 、 task.TANAL 以及 task.PHASEANAL 三种模拟工况，分别对应 力学 、热学 以及 相变 模拟。

sim.simTypeMasks = task.MANAL

注意

本节中的四项设置均不可忽略！

9.3.8. 计算循环

在 MPM2D 计算前，用户需执行:

sim.PreAnalysis()

完成所有的对象内存分配、物质点生成以及边界和初始条件初始化。

  • GUI 模式运行：

>>> app.run()
>>> app.run(100)

SudoSimGUI 的程序计算通过 run() 或者 run(100) 触发，前者将运行程序至所指定的总时间，而后者运行给定的迭代步 (100)。

  • CLI 模式运行：

>>> sim.run()
>>> sim.run(100)

与 GUI 模式同理，仅有 run 的执行者变为 sim 而非 SudoSimGUI 的成员 app 。

警告

不执行 sim.PreAnalysis() 可能导致未定义行为！

9.4. 材料本构类型

MPM2D 目前支持如下材料本构模型：

  • 线弹性模型 mpm.Elastic

  • Mohr Coulomb 模型 mpm.Mc

  • Drucker Prager 模型 mpm.DruckerPrager

  • Modified Cambridge Clay 模型 mpm.CamClay

  • Neo Hookean 模型 mpm.NeoHookean

  • Newtonian 流体模型 mpm.NewtonLiquid

  • 边界条件虚模型 mpm.GhostMat

所有的材料本构模型的定义均为 object = mpm.Name() ，例如线弹性模型定义为:

myelatic = mpm.Elastic()

材料的参数分为:

  • 共同参数: 密度 rho，接触法向及切向刚度 kn 和 ks , 接触摩擦系数 friction ，有限应变或小应变布尔值 infinistate ，平面应变或平面应力布尔值 planestrain 。

  • 特异参数: 每个本构模型自身特有的状态参量，用户可通过 help(Name) 来获取，例如通过：

help(mpm.Mc)

获取 Mohr Coulomb 模型 特有参数，用户将得到以下结果：

class Mc(Material)
|  Data descriptors defined here:
|  E
|      Young's modulous
|  c
|      cohesion
|  nu
|      Poisson's ratio
|  phi
|      friction angle
|  psi
|      dilation angle
|  tension
|      tension cut-off

综上，材料的声明及其参数定义有如下格式：

myelastic = mpm.Elastic()
myelastic.rho = 1000
myelastic.E = 1e4
myelastic.nu = 0.25
myelastic.infinistate = False
myelastic.planestrain = True

myelastic.printInfo()

提示

用户可通过调用 printInfo 方法来查询材料被正确设置。

材料的 infinistate 默认值为 False ，即默认为小应变模式； planestrain 默认值为 True ，即为平面应变模式。

用户通过添加材料到 sim.mats 完成材料对象在主控成员 sim 中的注册，若有多个材料参与注册则有:

mymat1 = mpm.Elastic()
mymat1.rho = 1000
mymat1.E = 1e4
mymat1.nu = 0.25

mymat2 = mpm.Mc()
mymat2.E = 1e5
mymat2.nu = 0.3
mymat2.c = 1000
mymat2.phi = 15
mymat2.psi = 0.01
mymat2.tension = 0.01

sim.mats = [mymat1, mymat2]
sim.mats[0].setIndex(0)
sim.mats[1].setIndex(1)

用户需把每个材料添加到 sim.mats 并为其指定顺序且唯一索引。

注意

边界条件虚模型 mpm.Ghostmat 是一种特殊的材料本构模型，该本构模型不参与任何应力-应变计算，仅作为 刚性物质点 或者多尺度模拟 DEM RVE 的占位符。刚性物质点是一种仅用于提供速度边界或者温度边界的特殊物质点，本身不携带任何其他物理量； DEM RVE是一种利用离散元方法代替传统本构模型的求解技术，本身不携带 MPM 中密度以外的任何其他物理量。由于程序实现的原因，二者皆需通过 mpm.Ghostmat 完成必要的求解初始化。

mpm.Ghostmat 的声明及注册与其他材料相似：

mygho = mpm.Ghostmat()
mygho.rho = 1.

sim.mats = [mygho]

每个材料都被映射为一个 枚举变量 以便于用户更为快捷的绑定材料到特定形状，枚举全部为 大写字母，它们是：

  • mpm.Elastic: mpm.ELASTICMAT

  • mpm.Mc: mpm.MCMAT

  • mpm.DruckerPrager: mpm.DPMAT

  • mpm.CamClay : mpm.CAMCLAYMAT

  • mpm.NeoHookean: mpm.NEOMAT

  • mpm.NewtonLiquid: mpm.NEWTONMAT

  • mpm.Ghostmat : mpm.GHOSTMAT

9.5. 形状类型

MPM2D 的 形状 模型最核心的任务为辅助物质点或物质点集合的生成。在模拟较为复杂的问题中，通常求解域不再为简单的矩形或者圆盘，可能是多种形状的组合形状；出于此考虑，MPM2D 提供了通用简单形状以及可灵活扩展的形状，具体如下：

  • 矩形 mpm.Rectangle

  • 圆盘 mpm.Sphere

  • 圆环 mpm.Ring

  • 多面体 mpm.Polygon

  • DEM驱动圆盘 mpm.DemSphere

  • DEM驱动多面体 mpm.DemPolygon

每个形状的构造函数同时定义了形状的几何约束以及边界条件类型。在 sandgravity.py 例子中：

soil = task.Rectangle(0.85, 1.75, 0.3, 2.3, task.GRAVITYSHAPE)
soil.setBc(False, 0, True, -9.8)
soil.setMatProperty(task.MCMAT, 0, 2)

矩形的x范围为0.85到1.75，y范围为0.3到2.3，边界类型为重力边界。

函数不仅需被指定范围及边界，还需被指定材料的类型及材料对象索引，材料类型参见 材料枚举 ，而材料索引则为材料在 sim.mats 中的顺序号。

备注

形状的几何约束将被程序内部用作检测特定的物质点是否处于形状内，检测的标准为形状的线密度。若形状线密度为2，则程序将会以每个背景网格单元包含4个物质点的准则来循环扫略所有物质点，并将满足条件的物质点记录用以物质点的生成。

9.5.1. 矩形 mpm.Rectangle

矩形形状构造函数如下：

mpm.Rectangle(xmin, xmax, ymin, ymax, bctype)
参数:

  • xmin (float) – x轴起始坐标

  • xmax (float) – x轴终止坐标

  • ymin (float) – y轴起始坐标

  • ymax (float) – y轴终止坐标

  • bctype (int) – 边界条件类型

9.5.2. 圆盘 mpm.Sphere

圆盘形状构造函数如下：

mpm.Sphere(centerx, centery, radius, bctype)
参数:

  • centerx (float) – 圆盘中心x坐标

  • centery (float) – 圆盘中心y坐标

  • radius (float) – 圆盘半径

  • bctype (int) – 边界条件类型

9.5.3. 圆环 mpm.Ring

圆环形状构造函数如下：

mpm.Sphere(centerx, centery, innerradius, externalradius, bctype)
参数:

  • centerx (float) – 圆环中心x坐标

  • centery (float) – 圆环中心y坐标

  • innerradius (float) – 圆环内半径

  • externalradius (float) – 圆环外半径

  • bctype (int) – 边界条件类型

9.5.4. 多面体 mpm.Polygon

多面体状构造函数如下：

mpm.Polygon(bctype)
参数:

bctype (int) – 边界条件类型

注意

多面体的构造函数相较规则形状较为特殊，仅包含边界类型。MPM2D 的多面体由一系列外边界轮廓点来定义，方法为 setPts 。

若用户已知 N 个多面体外边界点坐标，且用户保证尾部点不与首部点重合，则有如下 mpm.Polygon 声明及定义方法：

mypoly = mpm.Polygon(mpm.NONEBC)
inputpts = [x0, y0,
            x1, y1,
            ...
            xn-1, yn-1]
mypoly.setPts(inputpts)

9.5.5. DEM圆盘 mpm.DEMSphere

所有具有DEM前缀的形状均为边界条件或接触形状，此类形状与 材料定义 处的 mpm.Ghostmat 类似，并不参与物质点的生成，而是作为刚体实体直接参与到计算当中。若用户模拟圆盘刚体或者可自由移动的圆形DEM颗粒与物质点颗粒的相互作用，则可使用 mpm.DEMSphere 。

mpm.DemSphere(timestep, searchrange, searchinterval, kn, ks, friction)
参数:

  • timestep (float) – 求解时步

  • searchrange (float) – 搜索物质点范围

  • searchinterval (float) – 搜索物质点频率

  • kn (float) – 接触法向刚度

  • ks (float) – 接触切向向刚度

  • friction (float) – 接触摩擦系数

用户通过如下命令声明 mpm.DemSphere ：

mydemsph = mpm.DemSphere(1e-3, 10, 100, 1e5, 7e4, 0.5)
px = 0.5
py = 0.7
radius = 1
mydemsph.setPts([px, py, radius])

用户定义了圆心在 (0.5, 0.7) 半径为1的 mpm.DemSphere，该圆盘每隔100迭代步扫略10倍半径内的区域，并将接触的物质点添加到内存中进行接触求解，接触求解的法向刚度为1e5，切向刚度为7e4，接触摩擦为0.5。

9.5.6. DEM多面体 mpm.DemPolygon

若用户模拟更为复杂的刚体形状，用户可考虑将该形状离散为多面体，并使用 mpm.DemPolygon 来模拟 。mpm.DemPolygon 的构造函数与 mpm.DemSphere 完全相同，setPts 与 Polygon 相同。

mydempoly = mpm.DemPolygon(1e-3, 10, 100, 1e5, 7e4, 0.5)
inputpts = [x0, y0,
            x1, y1,
            ...
            xn-1, yn-1]
mydempoly.setPts(inputpts)

注意

DEM前缀形状既可协同位移边界条件来表征刚体，亦可被赋予质量以及转动惯量来表征结构物实体，DEM形状和MPM物质点的相互作用由 SudoSim 的非球接触算法驱动。

9.6. 物质点类型

9.6.1. 常规物质点/非刚体物质点

默认状态下，shapes 会生成常规物质点或者非刚体物质点，此类物质点具有所有的物理量：质量、动量、应力和应变等。在 MPM2D 中用户无法直接创建并添加物质点到主控成员 sim 中，而是通过 shapes 的托管方式来完成。但用户可通过如下方法读取非常规物质点信息：

outnrpts = sim.nrpts

该命令将会导出所有的非刚体物质点信息到 outnrpts 变量，用户可将其视作包含物质点指针的一维数组。用户可进一步从 sim.nrpts 来提取单个物质点的只读信息：

ptid = 0
targetpt = sim.nrpts[ptid]

pos = targetpt.getPos()
strain = targetpt.getStrain()
stress = targetpt.getStress()
temperature = targetpt.getTemperature()
linkednodes = targetpt.getNodeList()

上文代码展示了如何获取指定索引物质点的位置、应力、应变、温度以及映射节点列表。

9.6.2. 边界物质点/刚体物质点

当用户需要施加移动边界条件，或者用户需要施加变边界条件时，可考虑使用边界物质点/刚体物质点。用户无法直接添加刚体物质点，而是通过 shapes 来代理添加：

ghomat = mpm.Ghostmat()
sim.mats = [ghomat]

parbc = mpm.Rectangle(1, 5, 1, 2, mpm.RIGIDSHAPE)
parbc.setBc(True, 0, True, 0)
#parbc.setBc(False, 0, True, 5.)
bot.setMatProperty(mpm.RIGIDMAT, 0, 1)
sim.shapes = [parbc]

上文在矩形区域内添加了刚体物质点，矩形区域被施加 mpm.RIGIDSHAPE 枚举表示该区域不再生成非刚体物质点而是刚体物质点。刚体物质点的边界条件在x和y轴分别独立并由 setBc 指定，True 表示施加速度边界，紧随其后为速度边界数值。刚体形状的材料本构为 mpm.Ghostmat ，对应材料索引为 Ghostmat 在 sim.mats 内索引。

用户可通过如下方法读取刚体物质点信息：

outrpts = sim.rigidpts

用户可进一步从 sim.nrpts 来提取单个刚体物质点的只读信息：

ptid = 0
targetpt = sim.rigidpts[ptid]

pos = targetpt.getPos()
vel = targetpt.getVel()

注意

  • 刚体物质点索引和非刚体物质点 互相独立。

  • 刚体物质点的边界值可以为非零值，若为非零值，刚体物质点在运行状态会以被施加的边界速度值进行移动，以此来表征 移动边界 。

9.7. 边界条件类型

MPM2D 提供多种边界条件以适用于不同工况的大变形模拟：

  • 速度边界：通过刚体物质点或者节点来施加

  • 体力边界/重力边界：直接施加于非刚体物质点

  • 面力边界：直接施加于非刚体物质点

  • DEM边界：通过定义DEM前缀形状来施加

  • 物质点初值速度：直接施加于非刚体物质点

边界的枚举变量如下：

  • 刚体物质点速度边界： mpm.RIGIDSHAPE

  • 节点速度边界： mpm.NDSHAPE

  • 体力边界/重力边界： mpm.GRAVITYSHAPE

  • 面力边界： mpm.MPTTRACSHAPE

  • DEM边界： mpm.DemPolygon / mpm.DemSphere

  • 物质点初值速度： mpm.MPTVELSHAPE

  • 无边界条件： mpm.NONEBC

9.7.1. 速度边界

用户通过指定 mpm.RIGIDSHAPE 或者 mpm.NDSHAPE 来指定边界形状

parbc = mpm.Rectangle(1, 5, 1, 2, mpm.RIGIDSHAPE)
parbc.setBc(False, 0, True, -0.5)
ndbc = mpm.Rectangle(0.9, 5.1, 0.9, 2.1, mpm.NDSHAPE)
parbc.setBc(True, 0, True, 0.0)

注意

  • 刚体物质点边界将被施加到刚体物质点线性形函数映射的所有节点，而当使用 mpm.NDSHAPE 时，矩形区域会扫掠所有节点而非物质点，所以节点的形状坐标应略大于刚体物质点范围，以便程序正确捕获被施加速度边界的目标。

  • 尽管节点速度边界也可取值为非零值，但考虑到移动的非刚体物质点形函数变化，通常建议用户的节点边界只用作固定边界。

9.7.2. 体力边界/重力边界

用户通过指定 mpm.GRAVITYSHAPE 来指定边界形状

bodybc = mpm.Rectangle(1, 5, 1, 2, mpm.GRAVITYSHAPE)
bodybc.setBc(False, 0, True, -9.8)
bodybc.setBc(True, -9.8*np.sin(np.pi/6.), True, -9.8*np.cos(np.pi/6.))

用户可通过改变体力/重力边界值实现旋转坐标模拟。

9.7.3. 面力边界

面力边界仅对目标非刚体物质点施加单侧压或拉应力而非体力。面力边界需要被提供施加拉/压应力方向、被施加的面以及力的数值：

rightbc = task.Rectangle(1.7, 1.8, 0.3, 2.3, task.MPTTRACSHAPE)
right.setTracBc(2, task.NORMALDIR, q) # 面索引、力方向、力的数值

物质点的面索引以底部为1开始，沿着逆时针方向依次为 1,2,3,4。物质点的力方向枚举为 mpm.XDIR、mpm.YDIR、mpm.NORMALDIR 以及 mpm.TANGENTDIR。方向分别表示力被施加于x方向、y方向、面的法向以及切向。当物质点的力施加方向为法向时，负值力表示压应力，正值表示拉应力；当方向为切向时，负值表示使得物质点发生顺时针转动的力，正值表示使物质点发生逆时针转动的力。

9.8. 多物理场模拟设置

MPM2D 默认情况下运行力学模拟，本章旨在为用户提供可进一步自定义的热学或者热力学模拟。

9.8.1. 热学模拟

MPM2D 材料的热学共同参数的设置，包括比热容以及热传导系数设置：

myelastic = mpm.Elastic()
myelastic.setThermalProp(1, [0.1])
myelastic.setThermalProp(1, [0.1, 0.2])
myelastic.setThermalProp(1, [0.1, 0.2, 0.3, 0.4])
myelastic.thermalExpCoeff = 1.

setThermalProp 接受 float 以及 float[] 参数分别用作比热容及热传导系数张量，热传导张量列表的大小可分别为1，2，4。thermalExpCoeff 表征材料热响应对力响应的影响因子。

温度边界 ：

mytempbc = mpm.Rectangle(0, 1, 0, 1, mpm.RIGIDSHAPE)
mytempbc.setBc(True, 0, True, 0)
mytempbc.setThermalBc(100)

setThermalBc 用于设置温度边界，可与力学边界 setBc 共同作用。

主控成员 sim 设置：

sim.simTypeMasks = mpm.MANAL | mpm.TANAL

若要进行热力学模拟，只需在 sim.simTypeMasks 中以 BIT:OR 同时添加力及热的工况标记；若只需单独进行热或者力学模拟，则只保留对应的标记。

9.8.2. 相变模拟

MPM2D 相变模拟基于热学模拟，但材料参数的赋值有所不同。基于饱和混合介质假定，材料最多可能出现固相、液相以及冰相，所以三者的密度、初始体积分数以及热参数需被指定：

myelastic = mpm.Elastic()
# solid, ice, water
myelastic.setInitFraction([2650, 1000, 1000])
myelastic.setInitDensity([0.5, 0.25, 0.25])
# heat specific capacity: solid ice water
# heat conductivity: solid1 solid2 ice1 ice2 water1 water2
myelastic.setThermalProp([800, 2000, 4000], [9, 9, 2, 2, 0.6, 0.6])

液相与冰相的转换阈值及相变潜热系数设置参见如下命令：

myelastic.freezepoint = 0.
myelastic.latentratio = 100.

MPM2D 考虑了细颗粒较粗颗粒介质在相变中的不同力学响应的影响（未冻水），用户应指定材料是否为粗或细颗粒：

iscoarse = True
myelastic.coarsegrain = iscoarse
if(iscoarse):
    myelastic.Tres = 0.2
else:
    myelastic.unfrozenlimit = 0.05
    myelastic.unfrozenratio = 0.5

粗粒土的相变假定瞬间完成，但由于数值实现，需指定略大于 freezepoint 的 Tres 以保证数值收敛。细粒土的相变具有冷/热锋滞后性，滞后曲线以 unfrozenlimit 和 unfrozeonratio 表征。用户可通过 myelastic.usecontentfrozen 来控制细粒土未冻水是否由含水量还是饱和度来控制。最后设定 myelastic.thermalExpCoeff 完成相变材料设定。

myelastic.thermalExpCoeff = 1.

主控成员 sim 设置：

  • 若要进行相变模拟，需同时指定 PHASEANAL 和 TANAL：

sim.simTypeMasks = mpm.PHASEANAL | mpm.TANAL

  • 若要进行热力相变模拟，则需指定：

sim.simTypeMasks = mpm.MANAL | mpm.PHASEANAL | mpm.TANAL

9.9. 运行时及后处理

GUI 模式下，用户可直接使用 SudoSimGUI 进行实时结果提取以及渲染，请参见 相应教程 。

CLI 模式下，MPM2D 提供两种方式用于用户的计算结果保存：

  • 二进制压缩文件：二进制压缩文件保存了所有物质点的所有信息，包括当前位置、位移、速度、应力、应变、温度等详细状态量信息；该文件较大，为降低文件保存对程序计算耗时的影响，文件以压缩格式保存为 result.xxx 文件。文件后缀为运行迭代步数，例如 result.0 和 result.1000，用户可通过官网获取二进制解析程序将二进制文件解析为 VTK 文件。

  • ASCII文本文件：ASCII文本文件为用户关注的、占用空间较小的运行时文件，例如用户可能关系某个边界物质点在运行过程中的接触反力，或者关注某一组物质点运行过程中的平均速度或者平均位移，用户可通过ASCII模式快速保存并可用任意文本编辑器随时访问。

  1. 用户通过 sim.writer 进行二进制文件保存：

sim.writer.setDir("Results", "mycase")
sim.writer.setStartTime(0)
sim.writer.setEndTime(25)
sim.writer.setIntervalTime(0.5)

如上命令设置了 sim.writer 保存路径为 ./Results/mycase/ 文件夹，保存时间范围为0到25秒，保存频率为0.5秒。若要解析二进制文件，假设用户已下载到名为 exdata 的解析程序：

./exdata ./result.*

exdata 将自动解析当前文件夹所有二进制文件，并将其对应 VTK 文件输出至当前路径。

  1. 用户通过 mpm.PyHook 来进行ASCII文件动态保存：

outdata = []
def savemyascii():
   global outdata
   shapeid = 0
   forcedir = 1
   force = sim.getShapeReactForce(shapeid, forcedir)
   outdata.append([sim.iter(), force])
   with open('./mydata.dat', 'w') as fw:
      for d in outdata:
         fw.write("%.5f\t%.5f\n"%(d[0], d[1]))

myhook = mpm.PyHook()
myhook.command = 'savemyascii()'
# Run to 50000
#myhook.reset(iterPeriod = 500, totalRuns = 100)
# Infinite execution
myhook.reset(iterPeriod = 500)
myhook.dead = False
sim.hooks = [myhook]

mpm.PyHook 接受一个字符串的Python函数，会以固定周期，例如示例中500迭代步来执行目标函数；若用户指定总运行次数100，mpm.PyHook 会在50000步后停止运行；否则，mpm.PyHook 将会持续运行直到程序退出。用户也可以通过 mpm.PyHook.dead=True 来显式中止运行。用户可以通过 mpm.PyHook 和 sim.writer 组合的形式达到更为灵活的结果保存效果。

9.10. 双轴剪切模拟

../_images/biax.png

图 9.1 双轴剪切示意图

双轴剪切如 图 9.1 所示，刚体物质点表征的刚性加载板将以用户设定速度进行剪切，试样左右两侧被施加100 kPa围压。在模拟双轴剪切前，首先要对物质点进行分区以便确定形状及边界条件设置。物质点分布如 图 9.1 (b) 所示，顶部及底部区域将被设置矩形范围的刚体物质点，非刚体区域的左右两侧一个网格单元宽度的矩形将被设置面力边界，其余区域仅生成物质点，无边界条件施加。

建模脚本如下：

  1# 初始化 SudoSimGUI 及 mpm 主控成员 sim
  2app.core.setAnalysisType(AnalysisType="MPMDEM")
  3import mpm as task
  4sim=app.core.sim
  5
  6# 试样高度及单元大小
  7height = 2.0
  8elemsize = 0.1
  9
 10# 设置网格尺寸及单元数量
 11sim.grid.setRange(0, 3, 0, 1+height)
 12sim.grid.setNumEle(30, 10+int(height/elemsize))
 13sim.grid.printInfo()
 14
 15# 阻尼设置： 颗粒阻尼、网格阻尼以及PIC/FLIP系数
 16particle_damping = 0.2
 17grid_damping = 0
 18picfrac = 1.0
 19
 20# 时步及围压
 21timestep = 1e-3
 22confiningStress = 1e5
 23
 24# 材料设定：材料1为刚体虚材料，材料2为非刚体材料
 25mat1 = task.Ghostmat()
 26mat1.rho = 1.
 27
 28mat2 = task.Mc()
 29mat2.rho = 2650.
 30mat2.E = 1e5
 31mat2.nu = 0.3
 32mat2.c = 1000.
 33mat2.phi = 20.
 34mat2.psi = 10.
 35mat2.tension = 100
 36
 37sim.mats=[mat1, mat2]
 38sim.mats[0].setIndex(0)
 39sim.mats[1].setIndex(1)
 40
 41# 形状设定
 42# 线密度为2，每个单元包含4个物质点
 43pic = 2
 44# 该形状具有面力边界
 45left = task.Rectangle(1.0, 1.05, 0.3, 0.3+height, task.MPTTRACSHAPE)
 46# setConsolidation用于为材料提供应力初值：(sxx, syy)
 47left.setConsolidation(-confiningStress, -confiningStress)
 48# 4表示物质点左侧表面
 49# NORMALDIR表示围压会始终施加于面法向，即使物质点发生旋转
 50# 面力边界仍会追踪旋转矩阵，保证面力正确施加
 51# 负值面力表示压应力
 52left.setTracBc(4, task.NORMALDIR, -confiningStress)
 53# 左侧区域使用索引为1的材料：Mohr Coulomb本构
 54left.setMatProperty(task.MCMAT, 1, pic)
 55# 阻尼被施加于形状
 56left.setDamping(particle_damping, grid_damping, picfrac)
 57
 58right = task.Rectangle(1.95, 2, 0.3, 0.3+height, task.MPTTRACSHAPE)
 59# 右侧边界与左侧边界均为面力边界，仅有面力索引由4变为2，表示右侧边界
 60right.setTracBc(2, task.NORMALDIR, -confiningStress)
 61right.setConsolidation(-confiningStress, -confiningStress)
 62right.setMatProperty(task.MCMAT, 1, pic)
 63right.setDamping(particle_damping, grid_damping, picfrac)
 64
 65# 底部刚性板由刚体物质点表示，被约束所有自由度，材料为虚材料
 66bot = task.Rectangle(0.5, 2.5, 0.2, 0.3, task.RIGIDSHAPE)
 67bot.setBc(True, 0, True, 0)
 68bot.setMatProperty(task.RIGIDMAT, 0, 2)
 69
 70# 顶部刚性板由刚体物质点表示，具有向下的移动速度边界，材料为虚材料
 71top = task.Rectangle(0.5, 2.5, 0.3+height, 0.4+height, task.RIGIDSHAPE)
 72top.setBc(True, 0, True, -0.02)
 73top.setMatProperty(task.RIGIDMAT, 0, 2)
 74
 75# 试样中心无任何边界条件施加
 76soil = task.Rectangle(1.05, 1.95, 0.3, 0.3+height, task.NONEBC)
 77soil.setConsolidation(-confiningStress, -confiningStress)
 78soil.setMatProperty(task.MCMAT, 1, pic)
 79soil.setDamping(particle_damping, grid_damping, picfrac)
 80
 81# 注册所有形状到 sim
 82sim.shapes = [left, right, top, bot, soil]
 83
 84# 设置模拟总时间30秒，时步0.001
 85sim.setTimeDt(30, timestep)
 86# 设置压缩文件存储路径、时间及频率
 87sim.writer.setDir("./Results", "biax")
 88sim.writer.setStartTime(0)
 89sim.writer.setEndTime(30)
 90sim.writer.setIntervalTime(0.4)
 91# 设置形函数GIMP，工况为力学模拟
 92sim.setMethod(task.GIMP)
 93sim.simTypeMasks = task.MANAL
 94# 单线程模拟，以及程序初始化
 95sim.PreAnalysis()
 96sim.setNumThreads(1)
 97# 动态保存函数，用户获取双轴剪切过程中应力-应变关系
 98ss_st=[]
 99def savess():
100    global ss_st
101    fyy = sim.getShapeReactForce(3, 1)
102    styy = sim.iter()*sim.dt
103    ss_st.append([styy, fyy])
104    with open('./sscurve.dat', 'w') as fw:
105        for d in ss_st:
106            fw.write(f'\t{d[0]}\t{d[1]}\n')
107# 钩子函数用于用户自定义动态输出
108pyhook1 = task.PyHook()
109pyhook1.command = 'savess()'
110pyhook1.reset(iterPeriod = 500, totalRuns = 60)
111pyhook1.dead = False
112sim.hooks = [pyhook1]

9.11. 颗粒柱溃散模拟

../_images/collapse.png

图 9.2 颗粒柱溃散示意图

本小节将会给出颗粒柱溃散例子，在该例子中将会引入镜像边界以及摩擦边界，脚本如下：

  1# 初始化 SudoSimGUI 及 mpm 主控成员 sim
  2app.core.setAnalysisType(AnalysisType="MPMDEM")
  3import mpm as task
  4sim=app.core.sim
  5
  6sim.grid.setRange(0, 0.6, 0, 0.15)
  7sim.grid.setNumEle(120, 30)
  8sim.grid.printInfo()
  9
 10# 阻尼及时步设定
 11particle_damping = 0.2
 12grid_damping = 0
 13picfrac = 1.0
 14pic = 2
 15dt = 5e-5
 16
 17# 虚材料
 18Mat1 = task.Ghostmat()
 19Mat1.rho = 1
 20
 21# Mohr Coulomb 材料
 22Mat = task.Mc()
 23Mat.E = 8.4e5
 24Mat.nu = 0.3
 25Mat.phi = 18.9
 26Mat.psi = 0.1
 27Mat.tension = 0.1
 28Mat.rho = 2650.
 29# 接触参数激活用于与DEM底板的摩擦接触求解
 30Mat.Kn = 5e5
 31Mat.Ks = 5e5
 32Mat.friction = 0.15
 33
 34sim.mats=[Mat1, Mat]
 35sim.mats[0].setIndex(0)
 36sim.mats[1].setIndex(1)
 37
 38# 左侧镜像边界: FIXMINUS 表示左侧挡板只约束向左的位移，而不约束向右的位移
 39left = task.Rectangle(0., 0.005, 0.01, 0.15, task.RIGIDSHAPE)
 40left.setBc(task.FIXMINUS, 0, False, 0)
 41left.setMatProperty(task.RIGIDMAT, 0, 1)
 42
 43# 右侧镜像边界: FIXPLUS 表示右侧挡板只约束向右的位移，而不约束向左的位移
 44right = task.Rectangle(0.205, 0.210, 0.01, 0.15, task.RIGIDSHAPE)
 45right.setBc(task.FIXPLUS, 0, False, 0)
 46right.setMatProperty(task.RIGIDMAT, 0, 1)
 47
 48# DEM驱动的形状用于模拟摩擦底板
 49# 多面体的节点坐标由Python列表顺序给出
 50pts = [0, 0.005, 0.6, 0.005, 0.6, 0.01, 0., 0.01]
 51foot = task.DemPolygon(dt, 10, 50, 5e5, 5e5, 0.15)
 52foot.setPts(pts)
 53# 若usePolygonContact为True，则物质点根据当前变形矩阵的多面体形状计算与底板接触
 54# 若为False，物质点始终以圆盘形状参与接触计算
 55foot.usePolygonContact(False)
 56# 若addForce2Pt为True则以体力方式对物质点施加接触力
 57# 若为False，则以面力方式，通过节点积分的方式施加接触力
 58foot.addForce2Pt(False)
 59foot.setBc(True, 0, True, 0)
 60foot.setMatProperty(task.RIGIDMAT, 0, 1)
 61
 62# 非刚体材料被施加重力边界
 63soil = task.Rectangle(0.005, 0.205, 0.01, 0.11, task.GRAVITYSHAPE)
 64soil.setBc(False, 0, True, -9.8)
 65soil.setMatProperty(task.MCMAT, 1, pic)
 66soil.setDamping(particle_damping, grid_damping, picfrac)
 67
 68sim.shapes = [left, right, foot, soil]
 69
 70# 运行及结果保存时间设置
 71sim.setTimeDt(20, dt)
 72sim.writer.setDir("./Results", "collapse")
 73sim.writer.setStartTime(0)
 74sim.writer.setEndTime(20)
 75sim.writer.setIntervalTime(0.05)
 76
 77# 插值形函数及工况设置
 78sim.setMethod(task.GIMP)
 79sim.simTypeMasks = task.MANAL
 80
 81# 程序初始化及并行设置
 82sim.PreAnalysis()
 83sim.setNumThreads(1)
 84
 85# 分段加载重力，保证准静态的重力固结
 86def gravityHook(steps = 10000):
 87    partg = -9.8*sim.iter()/float(steps)-0.098
 88    app.print('Current g:', partg)
 89    soil.setBc(False, 0, True, partg)
 90
 91# 利用钩子函数完成重力固结
 92pyhook1 = task.PyHook()
 93pyhook1.command = 'gravityHook()'
 94pyhook1.reset(iterPeriod = 100, totalRuns = 100)
 95pyhook1.dead = False
 96sim.hooks = [pyhook1]
 97
 98# 重力固结
 99app.run(10000)
100# 颗粒柱开始溃散
101# 移除右侧边界
102right.setBc(False, 0, False, 0)
103# 减小颗粒阻尼
104soil.setDamping(0.05, 0., 0.05)
105app.Run(10000)

9.12. 二维圆盘热模拟

../_images/thermal.png

图 9.3 二维圆盘热传导模拟

 1# 初始化 SudoSimGUI 及 mpm 主控成员 sim
 2app.core.setAnalysisType(AnalysisType="MPMDEM")
 3import mpm
 4sim=app.core.sim
 5
 6length = 1.0
 7eleNum = 20
 8elemsize = length / eleNum
 9pic = 2
10
11# 网格设置
12xrange = [0, length+2.0*elemsize]
13yrange = xrange
14sim.grid.setRange(xrange[0],xrange[1],yrange[0],yrange[1])
15sim.grid.setNumEle(eleNum+2, eleNum+2)
16sim.grid.printInfo()
17
18# 设置热力学材料参数
19myg = mpm.Ghostmat()
20myg.rho = 1
21
22mye=mpm.Elastic()
23mye.rho = 1
24mye.setThermalProp(1.0, [1.0,1.0])
25sim.mats=[myg, mye]
26sim.mats[0].setIndex(0)
27sim.mats[1].setIndex(1)
28
29# 左侧及底部设置刚体物质点，提供温度为100摄氏度的边界条件
30left = mpm.Rectangle(0*elemsize, 1.*elemsize, 1*elemsize, (eleNum+1)*elemsize, mpm.RIGIDSHAPE)
31left.setThermalBc(100.0)
32left.setMatProperty(mpm.RIGIDMAT, 0, pic)
33
34bot = mpm.Rectangle(0*elemsize, xrange[1] -0*elemsize, 0*elemsize, 1.*elemsize, mpm.RIGIDSHAPE)
35bot.setThermalBc(100.0)
36bot.setMatProperty(mpm.RIGIDMAT, 0, pic)
37
38soil = mpm.Rectangle(1*elemsize, xrange[1] - 1.*elemsize, 1.*elemsize, (eleNum+1)*elemsize, mpm.NONEBC)
39soil.setMatProperty(mpm.ELASTICMAT, 1, pic)
40
41sim.shapes = [left, bot, soil]
42
43# 设置运行时长、时步以及保存时间
44sim.setTimeDt(10, 1e-4)
45sim.writer.setDir("./Results", "thermal2d")
46sim.writer.setStartTime(0)
47sim.writer.setEndTime(10)
48sim.writer.setIntervalTime(0.1)
49
50# 设置插值形函数以及热学模拟
51sim.setMethod(mpm.GIMP)
52sim.simTypeMasks = mpm.TANAL
53
54# 内存初始化以及模拟
55sim.PreAnalysis()
56sim.setNumThreads(1)
57
58app.run()

9.13. 一维无限长冰柱融化模拟

../_images/thaw.png

图 9.4 一维冰柱融化模拟

  1# 初始化 SudoSimGUI 及 mpm 主控成员 sim
  2app.core.setAnalysisType(AnalysisType="MPMDEM")
  3import mpm
  4sim=app.core.sim
  5
  6# 背景网格设置
  7elemsize = 0.01
  8width = 5
  9height = 0.03
 10sim.grid.setRange(0, width, 0, height)
 11sim.grid.setNumEle(500, 3)
 12sim.grid.printInfo()
 13
 14# 固、冰及液相参数
 15soilc = 0.1
 16icec = 0.1
 17waterc = 0.1
 18
 19soilk = 0.1
 20icek = 0.1
 21waterk = 0.1
 22
 23soilrho = 1
 24icerho = 1
 25waterrho = 1
 26
 27hcs = [soilc, icec, waterc]
 28conducts = [soilk, soilk, icek, icek, waterk, waterk]
 29
 30# 初始体积分数、密度
 31initfracs = [0, 1., 0.]
 32initrhos = [soilrho, icerho, waterrho]
 33
 34mg = mpm.Ghostmat()
 35mg.rho = 1.
 36
 37mat=mpm.Elastic()
 38mat.E = 1.0
 39mat.nu = 0.25
 40
 41# 以下参数用以确定临界时步
 42mat.rho = min(initrhos)
 43mat.setThermalProp(min(hcs), [max(conducts), 0., 0., max(conducts)])
 44mat.thermalExpCoeff = thermalExpAlpha
 45mat.setInitFraction(initfracs)
 46mat.setInitDensity(initrhos)
 47mat.setThermalProp(hcs, conducts)
 48# 相变温度及潜热系数
 49mat.freezepoint = 0.
 50mat.latentratio = 1
 51# 粗粒土相变设置
 52mat.Tres = 0.19
 53mat.coarsegrain = True
 54
 55sim.mats=[mg, mat]
 56sim.mats[0].setIndex(0)
 57sim.mats[1].setIndex(1)
 58
 59# 形状
 60# 恒定温度边界
 61left = mpm.Rectangle(0, elemsize, elemsize, elemsize*2, mpm.RIGIDSHAPE)
 62left.setThermalBc(1)
 63left.setMatProperty(mpm.RIGIDMAT, 0, 1)
 64
 65# 无边界
 66soil = mpm.Rectangle(elemsize, width-elemsize, elemsize, elemsize*2, mpm.NONEBC)
 67soil.setMatProperty(mpm.ELASTICMAT, 1, 1)
 68soil.setDamping(0., 0., 0.)
 69
 70sim.shapes = [left, soil]
 71
 72# 程序自动计算临界时步并设置模拟时长
 73dt = sim.getThermalCriticalTimeStep()
 74sim.setTimeDt(0.5, dt)
 75
 76# 模拟保存设置
 77sim.writer.setDir("./Results", "coarse")
 78sim.writer.setStartTime(0)
 79sim.writer.setEndTime(10)
 80sim.writer.setIntervalTime(200*dt)
 81
 82# 工况设置
 83sim.setMethod(mpm.GIMP)
 84sim.simTypeMasks = mpm.PHASEANAL | mpm.TANAL
 85
 86# 初始化
 87sim.PreAnalysis()
 88sim.setNumThreads(1)
 89
 90# 函数用以计算热锋前沿位置
 91def getthawfront():
 92   maxpx = 0
 93   for pt in sim.nrpts:
 94       # 检测热锋前沿位置
 95       if(pt.getPhaseState() == 2):
 96           if(pt.getPos()[0]>maxpx):
 97               maxpx = pt.getPos()[0]
 98   return maxpx
 99
100data = []
101def recordThawFront():
102    global data
103    thisruntime = sim.dt*sim.iter()
104    data.append([thisruntime, getthawfront()])
105    with open('./coarse_thawfront.dat', 'w') as fw:
106        for d in data:
107            fw.write('%.4f\t%.4f\n'%(d[0], d[1]))
108
109totrun = 200000
110pyhook1 = mpm.PyHook()
111pyhook1.command = 'record()'
112pyhook1.reset(iterPeriod = int(totrun/30.), totalRuns = 30)
113pyhook1.dead = False
114sim.hooks = [pyhook1]
115
116# 设定初始温度值
117soil.setInitTemperature(-0.05)
118# 程序计算
119app.run(totrun)