面向环境力学的离散元分析软件研发和工程应用

季顺迎　赵金凤　狄少丞　孙珊珊

摘要：离散元计算分析软件对解决环境力学中的离散介质问题有至关重要的作用.针对环境灾害中的非规则颗粒单元，以圆球颗粒为基本单元，分别构造镶嵌组合单元、黏结组合单元、扩展圆盘单元和扩展多面体单元等，并在此基础上开发基于球形颗粒离散元方法的分析软件（Software of Spherical Particlebased Discrete Element Method，SDEM）.该软件可模拟碎冰、岩石和道砟等颗粒材料的力学行为，能直观展现这些力学过程的发生、发展和演化；基于GPU的并行计算实现离散元大规模计算的高效性.对SDEM软件在地质灾害、工程海冰和铁路道床等领域的应用进行介绍.

关键词：离散元方法； 环境力学； 寒区工程； 地质灾害； 有砟铁路道床

中图分类号： O347.7； TP319

文献标志码： B

0前言

美国学者CUNDALL等[1]在20世纪70年代提出的离散元法（Discrete Element Method，DEM） 发展迅速.特别是近年来，随着计算机技术的发展，该方法已广泛地应用于岩土、化工、生物和环境等诸多领域.[23]在DEM的工程应用中，一直伴随着相应软件的不断发展.CUNDALL和STRACK最早研发二维圆盘程序BALL和三维圆球程序TRUBAL.[4]1989年英国Aston大学THORNTON将TRUBAL改进为TRUBALAston版，后定名GRANULE.目前，基于DEM的商业软件主要有美国Itasca公司的PFC2D和PFC3D，英国DEMSolution公司的EDEM.此外，澳大利亚的相关学者也开展出色的工作.

在具有独立知识产权的国产软件研发中，吉林大学基于CAD成功解决机械部件与球体颗粒接触作用的仿真过程.[5]北京大学基于DEM的工程计算软件的前后处理系统，以人机交互方式实现建模、数据传输和动画演示等全部过程，并主要应用于冲击、碰撞和断裂等动力学问题中.[6]最近，大连理工大学面向环境力学问题，自行研发出基于球形颗粒离散元方法的分析软件（Software of Spherical Particlebased Discrete Element Method，SDEM）.该软件通过对环境力学领域中相关工程问题的发生机理、演化过程进行数值模拟和动态显示，为理论分析和工程应用提供数值方法和后处理分析技术.此外，基于GPU并行计算的高性能数值算法，实现离散元方法的大规模计算.目前，SDEM已成功应用到工程海冰、地质灾害和铁路道床等环境力学领域.

1SDEM中离散单元的构造

在自然条件下，铁路道砟、浮冰和块石等散体材料都以非规则形态存在.为此，在SDEM中通过黏结、镶嵌和扩展等不同方式进行非规则颗粒构造，以精确描述颗粒几何形态.

1.1镶嵌式组合颗粒单元

针对岩土材料中块石颗粒的非规则性，可采用具有镶嵌模式的组合颗粒单元进行构造.该组合颗粒由不同粒径的球形颗粒按一定组合方式构成，见图1.对于该组合颗粒的质量、质心和转动惯量采用有限分割法确定，并通过四元数方法进行组合颗粒的动力学分量在整体坐标和局部坐标下的转换.[7]采用镶嵌式组合单元可以根据设定的颗粒形态进行任意组合，应用简便.

1.2黏结式组合颗粒单元

对于可破碎的非规则颗粒单元，可以通过球形颗粒黏结的模式进行构造，见图2.该非规则单元可在一定作用力下按一定的破坏准则发生破碎.单元间的黏结方式可采用接触黏结或平行黏结.[4]在模拟岩石材料的破碎特性时，一般采用平行黏结方式，通过黏结球体单元间的作用力和力矩，确定其法向应力和切向应力.采用该平行黏结模式不仅可以构造非规则颗粒单元，而且可以对连续体材料的破坏过程进行模拟.[89]

1.3扩展多面体单元

扩展多面体单元具有计算简便、精度高的特点，可精确描述自然环境下颗粒材料的几何形态.[1011]在一个基础多面体模型的表面设置若干个等粒径的球体颗粒，所有球体颗粒的球心位于基础多面体的表面.采用该扩展模型可将空间上的一个点扩展为一个球体，将一条线段扩展为一个圆柱及柱端的两个半球体，将任一平面扩展为由圆柱体包络的具有一定厚度的板，将任一多面体扩展为具有光滑棱边及棱角的多面体.采用扩展多面体模型构造的块石颗粒见图3.

类似于扩展多面体单元，还可构造扩展圆盘单元：首先对所需要构造的圆形单元进行函数定义，然后采用无限多个球型颗粒对其填充，以构造出所设定单元形态，见图4.由此，圆盘单元间的相互作用转变为球体单元的接触判断和作用力计算.

为便于离散元方法在环境力学领域的应用，需要研发操作简易、应用灵活的前后处理系统.SDEM可显示散体材料内部力链的分布，直观地展现结构体与颗粒以及颗粒之间的接触状态.

2SDEM在环境力学中的应用

目前，SDEM已应用于工程海冰、铁路道床和地质灾害等不同研究领域.

2.1SDEM在工程海冰中的应用

2.1.1船舶在碎冰区航行的离散元模拟

近年来，随着极地科学考察和油气资源开发以及北极夏季通航可行性研究的开展，冰区船舶结构设计和航运性能分析引起人们更多的重视.针对碎冰区海冰的结构特点，采用三维圆盘单元描述碎冰块，对其与船体结构的相互作用过程进行离散元分析.当海冰单元直径在4.0～6.0 m随机分布时，初始密集度为50%，船舶以4.0 m/s的速度行进，此时海冰与船体的相互作用的整体状况见图5，海冰单元与船体的接触碰撞见图6.[12]在不同冰速、冰厚、冰块尺寸和冰块密集度等海冰条件下，对海冰与船体动力作用过程进行数值分析，可为冰区船舶的安全运行和结构设计提供一定借鉴.

2.1.2海冰与海洋结构相互作用的离散元模拟

在冰区油气开发中，海冰会导致海洋结构的强烈振动，并对结构累积疲劳损伤、上部油气管线和设备以及操作人员健康构成很大威胁.目前，我国渤海冰区的导管架海洋平台主要包括直立、锥体和自升式结构等.采用离散元模型可对海冰与不同类型海洋结构相互作用中的随机破碎过程、冰载荷时程和冰激结构振动响应进行数值计算.

在海冰与直立腿导管架海洋平台结构相互作用的离散元模拟中，将海冰离散为若干个规则排列且具有黏结破碎功能的颗粒单元，并通过海冰单轴压缩试验对单元间的黏结强度进行确定.在此基础上，对直立海洋平台结构的随机作用过程进行数值计算，获得不同桩径下的冰载荷和结构冰振响应.海冰与直立海洋平台结构相互作用过程见图7.由此可以看出颗粒的平均作用力，也反映冰内应力场分布的大致规律[8].

在冰区油气开发中，海洋平台上的锥体结构可以有效降低冰力，避免强烈的冰激振动，是目前渤海油气平台的主要结构形式.采用离散元模型可对海冰与锥体结构的相互作用过程进行数值模拟，对比分析上锥体和下锥体冰力的大小和特征.[9]在海冰与锥体结构的相互作用过程中，海冰主要以弯曲破碎为主，颗粒间发生拉伸和剪切破坏.将锥体平台简化为具有等效的质量刚度阻尼的计算模型，并且不考虑锥体的变形.海冰与锥体结构相互作用过程见图8.计算结果表明，DEM可以有效模拟海冰与锥体结构的相互作用过程，有效确定不同冰厚、冰速和冰强度条件下的冰载荷，并可分析不同锥角和锥径下的冰载荷.

在海冰与海洋结构相互作用的离散元模拟中，结构尺度的增大对计算规模和计算效率提出迫切的需求.为此，采用CUDA编程平台实现基于GPU的大规模离散元并行计算.为模拟大面积海冰分布，

取颗粒单元数目为16万个，其与自升式海洋平台的作用过程见图9.由此可确定海冰与多桩腿结构的冰载荷，并分析海冰与海洋结构相互作用过程中的破碎、重叠和堆积特性.

2.2SDEM在地质灾害中应用

2.2.1山体滑坡过程的离散元模拟

在地质灾害中，边坡的破坏形式有很多种.这里采用SDEM对重庆武隆鸡尾山的垮塌过程进行离散元分析.该垮塌在空间分布上具有明显的崩滑区、铲刮区、滑坡主堆积区和碎屑流堆积区等四部分.[13]针对以上滑坡特性，采用具有破碎特性的黏结颗粒模型模拟岩体材料.

首先对地质体进行有限元网格划分，并提取地质体表面的三角形节点坐标构造地质体模型，然后对关键块体进行颗粒填充，见图10.通过设置颗粒间的不同黏结强度以模拟地质体的力学性质.该山体垮塌过程的离散元模拟见图11.由此可对该山体垮塌的形成机理和演化过程进行数值分析.

2.2.2岩石破坏中声发射现象的离散元模拟

对于岩质边坡，声发射特征是典型的内部破坏状态与宏观可测物理量间的指标.这里采用SDEM软件对岩石试样在单轴压缩破坏过程及其伴随的声发射现象进行数值分析.[14]岩石试样在单轴压缩下的离散元模拟见图12.岩石试样在不同压缩应变下的力链结构和破坏状态见图13.通过离散元模拟得到的应力应变曲线及颗粒间的破碎规律可确定声发射特征.

2.3SDEM在风沙区铁路有砟道床中的应用

为模拟铁路有砟道床的动力特性，分别采用组合颗粒单元和扩展多面体单元模拟道砟碎石颗粒.在此基础上，分别对风沙贯入下的道砟力学特性和往复载荷作用下的动力特性进行离散元分析.

2.3.1含沙铁路道砟变形模量的离散元模拟

针对风沙区有砟铁路道床的结构特性，研究其在细沙贯入下的有效变形模量，有助于理解风沙影

响下的道床动力特性.采用离散单元模型对不同含沙率道砟材料的有效变形模量进行数值分析.

在离散元数值模拟中，细沙和道砟颗粒随机放置于刚性圆筒内，见图14.计算得到的不同含沙率道砟材料变形模量见图15.可知，在低含沙率时，沙石混合体的变形模量基本保持不变；在高含沙率时，变形模量随含沙率的增加呈线性降低.基于离散单元模型的数值模拟，还可在细观尺度上对沙石混合体的力链强度、空间分布及配位数进行分析，揭示沙石混合体有效变形模量随含沙量变化的内在机理.[15]

2.3.2铁路道砟在往复载荷下的离散元模拟

对有砟铁路道床在往复载荷作用下的动力特性研究有助于道床结构的合理设计.采用扩展多面体单元可构造非规则的道砟颗粒，并对道砟材料在往复载荷作用下的动力特性进行离散元分析，研究其沉降量和形变刚度的演变过程.

选用12种不同几何形态的扩展多面体道砟颗粒，对道砟材料在往复载荷下的动力特性进行离散元分析，见图16.道砟颗粒粒径在36～63 mm间随机分布.

在90个加载周期内，轨枕沉降量的变化情况见图17.结果表明，轨枕位移在随往复载荷上下振荡过程中，其有效刚度和累积沉降量在初始阶段波动较大，然后随道砟材料的密实过程不断增加并趋于稳定.该扩展多面体单元能够合理地模拟道砟材料的动力行为，由其细观作用过程揭示相应的宏观演化规律.

3结束语

针对环境力学中的颗粒材料计算力学问题，研制出SDEM软件.该软件可通过构造镶嵌组合单元、黏结组合单元以及扩展的圆盘和多面体单元对不同形态的颗粒单元进行相对精确的描述.此外，还开展基于GPU的大规模离散元数值计算，以更好地解决环境力学问题的计算规模和计算效率.采用SDEM对工程海冰、地质灾害和风沙区有砟道床的动力特性进行离散元数值分析，有效解决相应的工程问题.

在目前SDEM基础上，将针对环境力学中的颗粒材料计算力学问题，进一步开展颗粒渗流、复杂边界条件和颗粒形态等方面的研究.

致谢：SDEM的研发得到美国Clarkson大学Hayley H SHEN教授的指导；美国陆军寒区研究与工程实验室（CREEL） Mark A HOPKINS博士对扩展多面体离散元模型的程序研发给予大力支持，在此深表谢意.

参考文献：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孙其诚， 张凌， 等. 颗粒离散元法研究进展[J]. 力学进展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三维离散元法软件开发研究[J]. 华中科技大学学报： 自然科学版， 2005， 33（S1）： 345348.

[6]郑文刚， 刘凯欣. 离散元法工程计算软件的前后处理系统[J]. 计算机工程与科学， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季顺迎， 狄少丞， 李正， 等. 海冰与直立结构相互作用的离散单元数值模拟[J]. 工程力学， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季顺迎， 李紫麟， 李春花， 等. 碎冰区海冰与船舶结构相互作用的离散元分析[J]. 应用力学学报， 2013， 30（4）： 520526.

[13]冯振， 殷跃平， 李滨， 等. 重庆武隆鸡尾山滑坡视向滑动机制分析[J]. 岩土力学， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]严颖， 狄少丞， 苏勇， 等. 风沙影响下铁路道砟变形模量的离散元数值分析[J]. 计算力学学报， 2012， 29（3）： 439445.

（编辑于杰）

参考文献：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孙其诚， 张凌， 等. 颗粒离散元法研究进展[J]. 力学进展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三维离散元法软件开发研究[J]. 华中科技大学学报： 自然科学版， 2005， 33（S1）： 345348.

[6]郑文刚， 刘凯欣. 离散元法工程计算软件的前后处理系统[J]. 计算机工程与科学， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季顺迎， 狄少丞， 李正， 等. 海冰与直立结构相互作用的离散单元数值模拟[J]. 工程力学， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季顺迎， 李紫麟， 李春花， 等. 碎冰区海冰与船舶结构相互作用的离散元分析[J]. 应用力学学报， 2013， 30（4）： 520526.

[13]冯振， 殷跃平， 李滨， 等. 重庆武隆鸡尾山滑坡视向滑动机制分析[J]. 岩土力学， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]严颖， 狄少丞， 苏勇， 等. 风沙影响下铁路道砟变形模量的离散元数值分析[J]. 计算力学学报， 2012， 29（3）： 439445.

（编辑于杰）

参考文献：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孙其诚， 张凌， 等. 颗粒离散元法研究进展[J]. 力学进展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三维离散元法软件开发研究[J]. 华中科技大学学报： 自然科学版， 2005， 33（S1）： 345348.

[6]郑文刚， 刘凯欣. 离散元法工程计算软件的前后处理系统[J]. 计算机工程与科学， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季顺迎， 狄少丞， 李正， 等. 海冰与直立结构相互作用的离散单元数值模拟[J]. 工程力学， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季顺迎， 李紫麟， 李春花， 等. 碎冰区海冰与船舶结构相互作用的离散元分析[J]. 应用力学学报， 2013， 30（4）： 520526.

[13]冯振， 殷跃平， 李滨， 等. 重庆武隆鸡尾山滑坡视向滑动机制分析[J]. 岩土力学， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]严颖， 狄少丞， 苏勇， 等. 风沙影响下铁路道砟变形模量的离散元数值分析[J]. 计算力学学报， 2012， 29（3）： 439445.

（编辑于杰）