基于颗粒变形的开源离散元模拟框架YADE的实现及应用

SoftwareX 6（2017）118原始软件出版物YADE中可变形颗粒的离散元模拟Martin Haustein*，Anton Gladkyy，Rüdiger SchwarzeTU Bergakademie Freiberg，Institut für Mechanik und Fluiddynamik，Lampadiusstr.4，09596 Freiberg，德国ar t i cl e i nf o文章历史记录：2016年12月14日收到2017年4月26日收到修订版，2017年保留字：亚德离散元法颗粒材料变形a b st ra ct在本文中，我们描述了开源离散元框架YADE和一个新的变形引擎的实现YADE是一个高度可扩展的软件包，允许使用基于粒子的数值方法模拟颗粒材料领域的当前工业问题在现代DEM框架中，现在可以使用纯粹和简单的离散元方法来描述粉末和颗粒状材料（如金属颗粒）变形是通过扩大球形颗粒的半径来实现的，这取决于它们的重叠，使得材料的体积保持恒定。©2017作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本YADE 2016.06a用于此代码版本的代码/存储库的永久链接http://github.com/ElsevierSoftwareX/SOFTX-D-16-00098GNU通用公共许可证GNU General Public License 2使用git的代码版本控制系统使用C++、Python、Open MPI的编译要求，操作环境依赖cmake，gcc，boost，Qt，freeReader 3，libQGLViewer，Python，numpy，IPython，matplotlib，eigen，gdb，sqlite3，Loki，VTK，CGAL，SuiteSparse，OpenBLAS，Metis如果可用，链接到开发人员文档/手册https://yade-dem.org/doc/如有疑问，请发送支持电子邮件至https://answers.launchpad.net/yade/+addquestion或yade-users@lists.launchpad.net1. 介绍颗粒物质的模拟对于解决科学和工程中的许多问题具有重要意义材料的复杂行为可以描述为离散粒子相互作用的结果这一领域的一个重要数值方法是由Cundall和Strack在岩石工程领域引入的离散单元法（DEM）[1]。该理论的重要应用可以在岩石工程[1，2]，粉末冶金[3]，材料科学[4]，工艺工程[5]以及制药和食品工业[6]等领域中此外，* 通讯作者。电子邮件地址：martin. imfd.tu-freiberg.de（M. Haustein）。http://dx.doi.org/10.1016/j.softx.2017.05.001流体动力学（CFD）和有限元法（FEM）是可能的。最通用的软件包之一是YADE [7]。YADE被世界各地的许多研究人员广泛使用，拥有超过25个贡献者的大型，活跃和不断增长的社区。该软件受益于YADE社区添加的大量功能，例如颗粒-一系列出版物，该软件的使用可以在[16]中找到在这篇文章中，YADE软件将介绍一个例子，通过引入一个新的变形引擎，模拟材料从接触区到自由未压缩空间的移动，在模拟过程中的代码的可扩展性。通过这种方式，可以近似接触区的材料重新分布，而无需计算密集耦合，2352-7110/©2017作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxM. Haustein等人/SoftwareX 6（2017）118119=−333用FEM [12]或Voronoi单元计算DEM [17]。此外，该发动机可以与许多不同的接触定律一起使用。工程和科学中的许多过程涉及球形颗粒的变形，如粉末压实和金属颗粒处理。重要的是要注意，当忽略变形时，颗粒材料的孔隙率和相对密度的计算可能是有缺陷的在YADE中的典型DEM模拟中，球形颗粒填充在一个有限的空间中，如球磨机或多边形的成型压力机。几何图形可以从许多常用的格式导入，如STL或网格，由CAD程序生成通过求解牛顿与分子动力学模拟相反，在DEM中，从接触定律推导出粒子间的力，其中力是这些定律的简单形式是描述材料的粘弹性和塑性特性的弹簧/阻尼器系统的DEM和接触建模的进一步理论在别处描述[18粒子的力和运动主要是通过它们之间的相互作用来计算的。此外，在YADE中，可以应用所谓的发动机来实现重力或旋转。这些发动机 模拟的基本执行单元，并可用于在模拟中集成粒子变形，以获得 更 逼 真 的 压 力 下 粒 子 行 为 。 通 过 在 YADE 或 外 部 程 序 （ 如ParaView [21，22]）中观察所收集数据的轨迹和后处理，可以更深入地了解微观过程。2. 软件描述YADE（Yet another dynamic engine）是一个高度可扩展的开源框架，用于使用离散元方法的数值模型该软件是用C++编写的，具有高灵活性和高性能的Python这使得模拟数百万个粒子相互作用成为 可能。YADE 用于基于Linux 的操作系统 ，可从Ubuntu和Debian等常见发行版的官方存储库获得;源代码可在GitHub上免费获得，并根据GPL 2开源许可证进行许可。YADE的第一个版本于2004年推出从2007年到2010年，代码被重构，并变得更加一致[23]。Python用于仿真设置、控制和后处理。Python绑定是通过Boost- Python库实现的这种解释型编程语言向用户隐藏了它使YADE成为一个适合教学目的的软件，并被许多大学成功地用于向学生介绍DEM基础知识YADE有一个基于IPython的交互式shell此外，可以将Yade作为模块导入Python脚本。YADE使用了大量的元编程技术（在编译时借助宏和模板生成代码），这简化了新接触定律或引擎的编程，但妨碍了代码的调试模板的使用保证了严格类型范式，但增加了编译时间。粒子化是任何基于粒子的仿真软件的一个非常重要的部分。YADE在大多数重要的循环上使用OpenMP（开放式缺点是目前缺少MPI（消息传递接口）功能，这应该在未来开发。模拟可以作为独立应用程序运行，也可以以批处理模式运行，在某些特定的平台上，但减少了软件的代码库，增加了其可读性，减少了代码重复。YADE编译所需的最重要的库是Eigen（用于线性代数的模板头库），QGlviewer（图形界面），VTK（后处理工具），几个Python库，CGAL（几何包，可选）等。内联文档是在编译过程中按需生成的。Sphinx用于这些目的，并以PDF、HTML和EPUB格式生成文档Kozicki等人提出的初始YADE代码是使用C++03标准编写的[24]。自2011-2012年以来主要的编译器是GCC，但Clang也可以与禁用的并行化一起使用，因为Clang中对OpenMP的支持较弱YADE通过为许多Ubuntu和Debian版本构建'yadedaily' Deb-package来使用持续集成一个buildbot，由Univ主办。GrenobleAlpes在每次提交后重新编译代码，并在夜间检测编译失败。这个buildbot还在每次编译后运行单元和功能测试YADE有一个活跃的社区，正在开发代码并帮助用户回答他们的问题。Launchpad基础设施被用作通信平台以及bugTracker [25].3. 变形引擎3.1. 变形算法YADE的功能非常灵活，可以适应许多科学和工程问题。 当前的主题是粉末和金属的压实。描述变形的一种简单方法是假设颗粒的总体积恒定，并且材料从接触区重新分布到颗粒之间的自由空间[26]。这样，半径为r的原始球体变形为半径为r′的更大球体，并减小了一个球冠（见图1）。已知这种方法对于非多孔固体是可接受的[27]。Raji [28]发表了一个类似的方法，但据我们所知，目前还没有在现代开源软件中实现在通常的离散元法中，粒子碰撞后相互重叠，排斥力由粒子间距离δ求得。重叠区域的体积消失，颗粒物质的总体积不被保留（见图1）。 1 a）。在该DEM公式中，重叠区域的体积在颗粒上均匀地重新分布，并且每个接触伙伴的半径增加，使得体积和质量保持恒定（参见图1b）。排斥力由δ′计算，δ′现在也增加在球形颗粒的情况下，重叠体积对应于球冠的体积，其可以通过简单的几何观察来计算变形前球体的体积应等于变形后球体的体积减去球冠的体积。对于帽的高度h，半径在未变形的球体r，和粒子的半径后，变形r′，方程（1）得到。4πr3=4πr′3−πh2（r′−h）（ 1）YADE具有许多外部的、强制的和可选的依赖性。依赖关系使编译过程复杂化简化后，使用几何关系hr′d（见图1）。 1），其中d是到球体中心的距离该软件适用于科学和工程领域的参数化研究120M. Haustein等人/SoftwareX 6（2017）118∑[（）（）]1= − −+ii··IJIJFig. 1. 两个球体接触的示意图，半径为ri，距离接触区为di。在（a）中，通常的DEM公式用重叠区域（深灰色）和重叠距离δ示出。在（b）中，示出了具有变形引擎的DEM的扩展。重叠区域（虚线）中的材料在自由区域（深灰色）中重新分布，导致半径ri′和重叠距离δ′增加。到接触区，方程（1）可以重写为N触点[28]Nr3r ′3r ′d22 r ′d.（2）4i=1球体/壁或球体/刻面接触的距离d是颗粒半径与穿透深度δ之间的差。将如此计算的d添加到等式中的总和。（2）类似于球体/球体接触的d。变形粒子半径r′的三次函数的解可以用数值方法求出，例如用用于计算变形球体的半径的算法可以描述为以下伪代码：对于每个粒子，对于每一个接触-伙伴做：如果球体-球体接触：如果球壁接触，则计算d_Sphere_Sphere，否则计算d_Sphere_Wallendif将d附加到contact_vector末尾，求解Eq. 2对于r该算法同样适用于多粒度颗粒集合。在球体的变形过程中，每个粒子的半径独立地因此，在不均匀变形的情况下总是形成多分散填料。对于宽粒度分布，必须使用小的时间步长，以便算法能够收敛。变形过程是完全可逆的。接触松弛后，变形恢复，直到半径为r的原始球形颗粒恢复。新计算的半径被添加到粒子属性中，并且每个粒子上的力的以下计算是使用变形球体的半径完成对于变形引擎的实现，通过修改State.hpp文件，将用于半径变化的新属性dR添加到每个粒子。在不使用附加属性的情况下更新粒子半径是不可能的，因为必须在 当量（二）、使用dR属性允许灵活的后处理和变形模型的进一步开发。这种实施策略的缺点是，每一个接触法必须 在使用变形引擎之前进行修改。背景下在保持向后相容性的前提下，对粘弹性接触定律（ViscoelasticPM）和粘弹塑性法向接触模型（LudingPM）进行了修正。通过使用预处理器指令，YADE可以在没有变形引擎的情况下编译。YADE旨在让用户可以根据自己的需求包括或排除功能，例如SPH，OpenMP甚至GUI。表1在YADE中进行模拟的参数。物业价值在s 1 10−5中，颗粒密度（kg/m3）2000粒子半径（单位：m 5）10−3kn，单位N/m 1300ks，单位N/m 1300Cnin Ns/m 0.186Csin Ns/m 0.186摩擦系数µ0.5此外，必须考虑的是，接触检测也必须针对变形球体的增加的半径进行调整。这可以在YADE中通过更改IGeomFunctor中的交互检测因子和InsertionSortCollider中的aabbEnlargeFactor来完成。否则，无法检测到球体的接触，仅由于变形而增加的半径而重叠。除了接触定律和粒子状态的文件外，还使用“RealRad”属性扩展了VTKRecorder.cpp。这表示变形球体的半径，可以在ParaView中可视化。在YADE中包含Python允许人们通过使用PyRunner引擎在纯Python中运行复杂的计算，用户可以在高级编程语言中实现因此，修改可以简单快捷地完成，并且可以很容易地修改最终的错误这种方法的缺点是执行速度慢对于执行简单的计算、绘图和测试新的想法来说，它的速度足够快，但是对于复杂的问题，C++代码的实现是不可避免的3.2. 两壁间单个球体的变形作为变形引擎的测试，将直径为10 mm的单个球体放置在两个面对的平行壁之间。两壁之间的距离不断缩小，球体也因此被压缩。模拟的参数如表1所示。为了计算接触伙伴i和j之间的法向力fn，使用线性弹性接触fn=knδij+Cnδstecij，（3）其中kn是法向接触刚度，δij是粒子i和j的法向重叠，δstecij是δij和Cn的时间导数粘性阻尼系数。通过使用变形引擎，该线性模型显示非线性行为，因为重叠是半径的函数颗粒变形越大，抵抗进一步压缩的阻力就越大这意味着一个简单的接触定律可以用于粉末压制弹性体M. Haustein等人/SoftwareX 6（2017）118121=−图二、 直径为10 mm的球体在距离为10 mm、5 mm和3.3 mm的两个壁之间压缩时的变形。显示了半径为ri的原始球体（白色球体）和半径为ri′的变形球体（彩色网格）网格的颜色与变形程度相对应（关于此图例中颜色的参考解释，请读者参考本文的网络版本上壁再次达到零。使用与第3.2节所述相同的模拟参数和接触定律。这些模拟是在使用和不使用变形引擎的情况下完成的。如果允许颗粒变形，则与未变形的球体相比，需要较低的应变，直到达到力极限，如图5所示。此外，如图6所示，颗粒床的孔隙率随着变形而显著降低。这里，应变e被定义为e（L0L）/L0，其中L是从上壁到箱底的距离。L0是从底部到上壁初始位置由于颗粒物质的重新排列，e的值并不终止于零。对于如上所述的立方体盒中400个粒子的测试情况，在没有任何变形的情况下，实现了每秒约1900次迭代（it/s），而变形实现了每秒约1900次迭代（it/s）。图三. 来自[27]的橡胶球变形（无量纲）的实验数据与从YADE获得的模拟数据每个应变e的实验数据用离散符号标记，YADE用连续（彩色）线计算。虚线表示在压缩期间上壁的位置(For对本图图例中所指颜色的解释，读者可参考本文的网络版。）材料.此外，复杂的接触定律，如由Luding [19]提出的粘合剂、弹塑性法向接触模型，可以与变形引擎一起用于描述致密颗粒（如金属粉末或颗粒）的塑性变形将在这一领域进行进一步的此外，在切线方向上，弹簧-阻尼器系统（并联）与摩擦系数为记录变形球体的半径作为施加的应变e的函数。图2示出了处于不同压缩阶段的球体。用YADE软件计算的变形球半径可以与橡胶球在两壁之间变形的实验数据进行比较 如图 3、模拟结果与数据一致[27]。该方法的问题是接触区附近的半径部分高估和自由区的低估原因是球的常曲率不可能变形为椭圆粒子。3.3. 盒子中一组粒子的压缩在另一个测试案例中，将400个直径为10.0 mm的球形颗粒随机分布在80 mm长的立方体盒中。在第一步中，颗粒在重力作用下被压实成松散的随机填料（图4左图）。在第二步中，以0.01 m/s的速度降低上壁，并记录该壁上的力（图11）。4、中间和右边）。上壁的初始位置在达到1000 N的最大力后，以0.01 m/s的速度提升上壁模拟停止时，在C++中作为部分引擎测量1100 it/s。在Python中完全实现并由YADE的PyRunner引擎调用时，每秒只能进行10次迭代。很明显，Python中的纯实现对于真实世界的模拟来说效率不够，但是对于10到30个粒子的测试，计算工作量是足够的。对于某些模拟，仅在每十步更新变形也是足够的，从而可以进一步减少计算时间。4. 影响YADE可以通过使用CAD程序或直接在软件中构建的模型来模拟许多工业过程Python界面允许用户在计算运行时调整模拟参数，如时间步长后处理，如在模拟过程中直接绘制数据，可用于跟踪进度。新的接触定律和引擎的实现最近，新的功能，如光滑粒子流体动力学（SPH）和湿颗粒物质的毛细作用力已经成为可用[29由于采用了体积不变的假设，使得粒子的变形得以实现，因此YADE可以用纯离散元方法来描述橡胶球形状变化的实验数据这似乎是一种很有前途的方法来模拟软球形颗粒，如橡胶，甚至金属颗粒的压缩。这个模型的适用性，可变形泡沫气泡的模拟将很快进行检查据我们所知，用于模拟可变形聚合物和金属泡沫的纯DEM方法尚未发表。大颗粒系综在压力作用下变形的模拟通常颗粒的重叠一定不能大到导致一个真实的模拟大的重叠违反了质量守恒，因为物质在这种状态下会消失。变形引擎允许在YADE中使用物理正确的模型进行DEM模拟，从而获得更可靠的孔隙率，相对密度和颗粒材料的其他属性数据122M. Haustein等人/SoftwareX 6（2017）118图四、 400个球形颗粒在一个立方体盒子里的压实。颜色对应于变形球体的直径（蓝色= 10.0 mm，白色= 10.5 mm，红色=10.5 mm）。11.0 mm）。为了清楚起见，未示出上壁压缩前使用松散的随机填料（左）。之后，上冲头下降（中），直到达到1000 N的最大力（右）。（关于此图例中颜色的参考解释，请读者参考本文的网络版本基于几何考虑的可变形的球形颗粒。在整个模拟过程中，材料的体积保持不变。这一假设允许以适当的精度再现实验数据涉及颗粒变形的工业过程可以在一个完善的软件框架中使用纯DEM方法来近似引用[1] 宾夕法尼亚州坎德尔一种模拟块状岩石系统中渐进的大规模移动的计算机模型。输 入 ：程序Symp. int. Soc. 岩石机械南西. 一九七一年[2] Boon CW，Houlsby GT，Utili S.基于线性规划的岩石分层新方法。ComputGeotech2015;65：12-29.[3] 李文，李文辉，李文辉，等.用离散元法研究粉末冷压成形. 粉末技术2011;208（2）：537-41。图五. 在有和没有变形的情况下，400个颗粒的压实的在z方向上的壁上的力见图6。有和没有变形引擎的400个颗粒的球填料的孔隙率。孔隙率随着颗粒变形而显著降低。所提出的变形引擎仅基于几何考虑，不使用特定于材料的参数。这意味着它可以很容易地用作DEM中可变形球形颗粒的一级近似。它与所使用的接触定律无关，可以灵活地应用于各种材料行为。5. 结论在本文中，我们提出了DEM框架YADE旁边的变形引擎模拟压实过程。这种扩展的YADE代码允许逼真的模拟[4] Nitka M ，Tejchman J.混凝土单轴压缩和 拉伸 行 为 的DEM 建 模。 颗 粒物 质2015;17（1）：145-64。[5] 应用离散元法对滚磨机进行计算机模拟的回顾：第一部分-接触力学。Int JMiner Process2003;71（1 - 4）：73-93.[6] RajiA，Favier J. 农业和食品颗粒材料在压缩载荷下变形的离散元模型。I：理论、模型开发和验证。 食品工程杂志2004;64（3）：359-71.[7] SmilauerV，Emanuele C，Chareyre B，Dorofeenko S，Duriez J，Dyck N，et al Yadedocumentation. 第2版2015年。[8] ChareyreB，Cortis A，Catalano E，Barthélemy E.稠密球形填料中粘性流和诱导力的孔隙尺度模拟。Transp多孔介质2012;92（2）：473-93。[9] Lominé F，Scholtès L，Sibille L，Poullain P.用耦合格子Boltzmann/离散元方法模拟颗粒介质中的流固相互作用：应用于管道侵蚀。Int J Numer AnalMethods Geomech2013;37（6）：577-96.[10] 莫林河推移质输移的最小流体-离散元耦合模型。物理流体2015;27（11）：113302。[11] [10]杨文军，李文军.在离散元法的一般框架中对可变形结构的建模。土工织物和土工膜2016;44（2）：143-56。[12] 弗伦宁湾一种有效的有限元/离散元程序，用于模拟3D颗粒集合体的压缩。计算方法应用机械工程2008;197（49-50）：4266-72.[13] 郭N，赵杰.颗粒介质分层多尺度建模的FEM/DEM耦合方法。Internat J NumerMethods Engrg2014;99（11）：789 - 818.[14] BoonCW，Houlsby GT，Utili S.离散元法中凸多边形与多面体颗粒接触检测的一种新算法。Comput Geotech2012;44：73-82.[15] EliáJ. Simulation of railway ballast using crushable polyhedral particles. 粉末技术2014;264：458-65。[16] https://www.yade-dem.org/doc/publications.html网站。[17] 放大图片创作者：A-S用扩展的截球模型对受限粉末压缩过程进行离散元模拟的实验评价. 粉末技术2015;284：257-64。[18] 放大图片作者：Jing L，Stephansson O.岩石工程离散单元法基础：理论与应用。见：岩土工程发展，第85卷，阿姆斯特丹：牛津，爱思唯尔;2007年。M. Haustein等人/SoftwareX 6（2017）118123[19] 泸定湾离散单元法简介。接触力模型的基础，以及如何进行微观-宏观过渡到连续介质理论。Rev Eur Civil2008;12（7-8）：785-826。[20] 作者：Matuttis H-G，Chen JF.理解离散元方法：颗粒和多体系统的非球形颗粒模拟。Hoboken，NJ：John Wiley& Sons Inc.; 2014年。[21] 张文辉，张文辉，张文辉.ParaView：一个用于大数据可视化的终端用户工具。在：可视化手册. Elsevier;2005年。[22] 阿亚契特大学paraView指南：针对paraView 4.3版进行了更新。全彩色版。Los Alamos：Kitware;2015.[23] Šmilauer V.在yade平台上使用离散单元法的凝聚粒子模型[Ph. D. thesis）。布拉格，格勒诺布尔：捷克技术大学，格勒诺布尔第一大学;2010年。[24] Kozicki J，Donzé FV.一个新的开源软件，用于使用离散建模方法进行数值模拟。计算方法应用机械工程2008;197（49-50）：4429-43.[25] https://launchpad.net/yade网站。[26] ShimaS，Tatara Y，Iio M，Shu C，Lucero CJ. 橡胶球在径向压缩下的大变形：第2部分：多种橡胶材料的实验和理论与实验的比较。JSME Int J Ser A，Mech MaterEng 1993;36（2）：197-205.[27] 林永良，王登明，陆文明，林永生，董国良.软球形颗粒的压缩和变形。化学与工程科学2008;63（1）：195-203。[28] Raji AO。散装农业颗粒变形的离散元建模[博士学位]。thesis）。纽卡斯尔大学;1999年。[29] ScholtèsL，Chareyre B，Nicot F，Darve F.考虑毛细效应的颗粒材料细观力学。Internat J Engrg Sci2009;47（1）：64-75.[30] 格拉茨基不同毛细桥模型在离散元法中的应用比较。颗粒物质2014;16（6）：911 - 20.[31] DuriezJ，Wan R.湿粒状土壤中接触角的力学影响 ActaGeotech2017;12（1）：67-83.