Discrete particle simulation of particle–fluid flow model formulations and their applicability

PHYSICAL REVIEW E 72, 041307 ͑2005͒
Characterization of force chains in granular material
J. F. Peters,1,* M. Muthuswamy,2,† J. Wibowo,1 and A. Tordesillas2,‡
DOI: 10.1103/PhysRevE.72.041307
PACS number͑s͒: 83.80.Fg, 45.70.Ϫn, 81.05.Rm, 05.40.Ϫa
I. INTRODUCTION
In a granular material, forces are transmitted from one particle to the next via their contacts. When visualized through stress-induced birefringence within assemblies of photoelastic disks ͓1–3͔, the manner of propagation appears as a complex force network that is highly ramified and undergoes rapid changes in branch morphology throughout the material’s deformation history. Over the past few decades, there have been numerous attempts made toward the characterization of this force network, with techniques such as contact dynamics simulations ͓4,5͔, discrete element simulations ͓6͔, the scalar q model of force transmission ͓7,8͔ and its variants ͓9,10͔, the oriented stress linearity model ͓11͔, and models based on statistical mechanics ͓12͔, just to name a few. Two key characteristics of the force network have been identified ͑see, e.g., ͓13,14͔͒: ͑1͒ The majority of particles ͑typically greater than 60%͒ carry less than the average force ͑i.e., “weak network”͒, and ͑2͒ the number of particles carrying larger than the average force ͑i.e., “strong network”͒ decreases exponentially with increasing contact force magnitude.

∂
∂ε
ε
=
( αε μe
) + ( C ε1 G k - C ε2 ρ g ε)
∂x j
∂x j
k
dt
力视为 F f 。
2 理论介绍
∂
∂
( ρ ε) +
( ρ εu )
∂t g
∂x i g i
dω p
(4)
(5)
根据式(2) 和式(6) ,连续相与离散相之间通过
动量交换实现耦合:
→
- v
F f = K gs( v
离散颗粒相的受力分析和运动轨迹计算。 计算过程充分考虑颗粒体积分数对流体的影响。 通过对比不同直径的
水平管道在相同空气流量和物料流量下的颗粒流动形态和压力损失,发现颗粒在不同直径管道内表现出不同的流
动形态。 压力损失在管道直径为 320 mm 时最小,此时的管道直径正好为保持管底流动形态的最大管道直径。
1 引言
粒顺利输送。 然而在水平管道气力输送过程中,气
流方向与重力方向垂直,物料易在管底沉积。 因此
气力卸船 机 作 为 负 压 气 力 输 送 设 备, 因 其 具
要求水平管道具有比垂直管道更高的气流速度。 而
有输送系 统 封 闭、扬 尘 污 染 小、清 仓 彻 底、 构 造 简
且在同一管道输送系统中,物料往往更容易在水平
体对颗粒 的 作 用 力,包 括 曳 力、上 浮 力、 巴 塞 特 升
力、萨夫曼升力、马格纳斯升力等。 当颗粒密度与
流体密度相差较大时,曳力起主要作用,因此将曳
2. 3 连续相与离散相相互作用
基于 CFD-DEM 耦合方法,连续相采用 RNGk-ε
湍流模型,颗 粒 相 采 用 基 于 牛 顿 运 动 定 律 的 离 散

第 54 卷第 4 期2023 年 4 月中南大学学报(自然科学版)Journal of Central South University (Science and Technology)V ol.54 No.4Apr. 2023离散颗粒抑制热喷流红外辐射的大涡模拟胡峰1，孙文静1, 2，张靖周1，单勇1(1. 南京航空航天大学 能源与动力学院，江苏 南京，210016；2. 中国航天科工飞航技术研究院 北京动力机械研究所，北京，100074)摘要：为了探究气溶胶离散颗粒对飞行器排气喷管热喷流3~5 μm 波段的红外辐射的抑制效果，设计地面状态下气溶胶颗粒投射的仿真环境，采用大涡模拟和颗粒离散相模型对含气溶胶颗粒的飞行器排气喷管尾部气固两相剪切流进行数值模拟研究，系统地分析颗粒的质量流量、粒径和喷射速度对离散颗粒空间分布形态以及热喷流红外辐射抑制的影响规律。

研究结果表明：颗粒的质量流量和粒径对于红外抑制效率的影响较为明显，增加颗粒质量流量对颗粒的空间分布形态影响较小，但能够显著提升红外抑制效率；当颗粒粒径大于1.0 μm 时，颗粒空间分布均匀，红外抑制效率最高；颗粒的喷射速度对于颗粒的空间分布以及红外抑制效率的影响较小。

关键词：红外抑制；高速剪切流；大涡模拟；气溶胶颗粒分布；气固相互作用中图分类号：V231.1 文献标志码：A 开放科学(资源服务)标识码(OSID)文章编号：1672-7207（2023）04-1576-16Large eddy simulation of discrete particles suppressing infraredradiation from thermal jetsHU Feng 1, SUN Wenjing 1, 2, ZHANG Jingzhou 1, SHAN Yong 1(1. College of Energy and Power, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China;2. Beijing Power Machinery Institute, China Aerospace Institute of Science and Technology, Beijing 100074, China)Abstract: To investigate the effect of aerosol discrete particles on the infrared radiation suppression in the 3~5 μm band of thermal jets of aircraft exhaust nozzles, the simulation environment of aerosol particle projection under ground conditions was designed. The large eddy simulation(LES) and particle discrete phase model(DPM) were used to numerically simulate the gas-solid two-phase shear flow at the tail of the aircraft exhaust nozzle containing aerosol particles, and the effect law of particle mass flow, size and jet speed on the spatial distribution of discreteparticles and the suppression of infrared radiation of thermal jets was systematically analyzed. The results show收稿日期： 2022 −05 −14； 修回日期： 2022 −07 −23基金项目(Foundation item)：中国博士后科学基金特别资助(站前)项目(2020TQ0143)；江苏省自然科学基金青年基金资助项目(BK20200448) (Project(2020TQ0143) supported by the Postdoctoral Science Foundation of China; Project(BK20200448) supported by the Youth Fund of Jiangsu Natural Science Foundation)通信作者：孙文静，博士，讲师，从事气固两相湍流、湍流燃烧、流动强化传热研究；E-mail ：**************.cnDOI: 10.11817/j.issn.1672-7207.2023.04.033引用格式： 胡峰, 孙文静, 张靖周, 等. 离散颗粒抑制热喷流红外辐射的大涡模拟[J]. 中南大学学报(自然科学版), 2023, 54(4): 1576−1591.Citation: HU Feng, SUN Wenjing, ZHANG Jingzhou, et al. Large eddy simulation of discrete particles suppressing infrared radiation from thermal jets[J]. Journal of Central South University(Science and Technology), 2023, 54(4): 1576−1591.第 4 期胡峰，等：离散颗粒抑制热喷流红外辐射的大涡模拟that the effect of particle mass flow and size on infrared radiation suppression rate is obvious. With the increase of particle mass flow, its effect on the spatial distribution of particles is small, but the infrared suppression efficiency is significantly improved. When the particle diameter is 1.0 μm, the particle space distribution is uniform and the highest infrared suppression rate is achieved. However, the particle injection speed has less effect on the spatial distribution of particles and infrared radiation suppression efficiency.Key words: infrared suppressing; high-speed shear flow; large eddy simulation; aerosol particle distribution; gas-solid interactions气溶胶红外隐身技术是一种主动型应急红外对抗技术，该技术利用附加的机载引气装置，将细微颗粒喷射在发动机热喷流周围形成气溶胶云，借此对排气喷管热内腔和热喷流的强红外辐射进行遮蔽和散射。

考虑组构演化的散粒体各向异性力学特性数值分析与本构模拟摘要：散粒体材料是一种具有重要应用价值的材料，其力学特性与材料内部结构的组构演化密切相关。

本文从组构演化的角度出发，研究了散粒体材料的各向异性力学特性，采用数值分析和本构模拟的方法对其力学行为进行了研究。

首先介绍了散粒体力学特性的研究现状和存在的问题，然后重点阐述了组构演化与材料力学特性之间的关系。

接着，通过建立数学模型和仿真模拟，探究了散粒体材料的各向异性力学特性变化规律和组构演化的影响。

最后，对本文的研究结果和对散粒体材料力学性质的认识加以分析和总结，提出了未来研究的方向和展望。

关键词：散粒体；组构演化；力学特性；数值分析；本构模拟Abstract: Granular material is a kind of material with important application value. Its mechanical properties are closely related to the structural evolution of the material. Starting from the perspective of structural evolution, this paper studies the anisotropic mechanical properties of granular materials, and uses numerical analysis and constitutive simulation tostudy their mechanical behavior. First, the research status and existing problems of granular material mechanics are introduced, and then the relationship between structural evolution and material mechanicalproperties is elaborated. Then, through the establishment of mathematical models and simulation, the change law of anisotropic mechanical properties of granular materials and the influence of structural evolution are explored. Finally, the research results of this paper and the understanding of the mechanical properties of granular materials are analyzed and summarized, and the future research directions and prospects are put forward.Keywords: granular material; structural evolution; mechanical properties; numerical analysis;constitutive simulationGranular materials are materials made up of small, discrete particles, such as sand, gravel, or powders. They are widely used in many fields, including civil engineering, chemical engineering, and geotechnical engineering, among others. Understanding the mechanical properties of granular materials is a fundamental problem in many applications, such as the design of foundations, embankments, and tunnels.The mechanical properties of granular materials are strongly influenced by their microstructure, which changes continuously during deformation. The structure of a granular material is defined by the packingarrangement of its particles, which depends on the particle size, shape, and the forces acting between them. When a granular material is subjected toexternal forces, such as compression or shear, its particles rearrange and its structure evolves.The structural evolution of granular materials has a significant impact on their mechanical behavior, such as their deformation, strength, and stiffness. For example, under compression, a granular material undergoes densification, accompanied by a decrease in its porosity and an increase in its stiffness. However, this process is not a simple isotropic compaction, but rather involves significant anisotropy in the packing structure, particle orientation, and stress transmission. Similarly, under shear, a granular material experiences strain localization, which isalso influenced by the anisotropic particle arrangements and the development of shear bands.To understand the mechanical properties of granular materials and their structural evolution, numerical analysis and constitutive simulation are widely used. Numerical models, such as the discrete element method (DEM), are used to simulate the behavior of granular materials at the particle scale. Constitutive models, such as the Mohr-Coulomb model and the critical statemodel, are used to describe the macroscopic behaviorof granular materials in terms of their stress-strain relationships and failure criteria.Through numerical analysis and constitutive simulation, researchers have made significant progress in understanding the mechanical properties of granular materials and their structural evolution. For example, the anisotropic mechanical properties of granular materials have been characterized and quantified, and the impact of the packing structure and particle interactions on these properties has been explored. Moreover, the effect of structural evolution on the mechanical behavior of granular materials has been investigated, and new constitutive models that takeinto account this effect have been proposed.In summary, the mechanical properties of granular materials are strongly influenced by their microstructure, which evolves continuously during deformation. Numerical analysis and constitutive simulation are powerful tools for studying thebehavior of granular materials and their structural evolution. The research in this field has important implications for many applications, and future studies should focus on developing more accurate and robustmodels for predicting the mechanical behavior of granular materials in different conditionsContinued:One important application of granular materials is geotechnical engineering, where soil and rock are the primary materials utilized. In order to prevent landslides, control erosion, design foundations and tunnels, and excavate deep underground facilities, an accurate prediction of the mechanical behavior of granular materials is essential. Therefore, researchers in this field need to systematically investigate the influence of various factors on the mechanical behavior of granular materials, such as strain rate, temperature, moisture content, particle size distribution, and confining pressure.Another area where granular materials play an important role is in the manufacturing industry, where powders and granules are used to produce a wide range of products, such as pharmaceuticals, chemicals, ceramics, and food products. In order to optimize the manufacturing process and ensure the quality of the end product, it is crucial to understand the mechanisms underlying the deformation and failure of granular materials during processing. Therefore,researchers in this field need to develop new methods for characterizing the microstructure of granular materials, such as X-ray tomography and micromechanical testing, and integrate these methods into numerical simulations to predict the mechanical behavior of granular materials in different processing conditions.In addition, granular materials are also used in the energy industry, where they play a critical role inoil and gas extraction and storage, carbon capture and storage, and nuclear waste disposal. In these applications, the mechanical behavior of granular materials under extreme conditions is of great significance, such as high pressure, high temperature, high radiation dose, and aggressive chemical environment. Therefore, researchers in this field need to develop new testing techniques and modeling approaches that can capture the effects of these extreme conditions on the structural evolution and mechanical behavior of granular materials.Finally, granular materials are also widely used in natural phenomena, such as landslides, avalanches, and earthquakes. In order to understand and predict the behavior of these phenomena, it is essential to develop new experimental techniques and numericalmodels that can accurately describe the deformation and failure of granular materials under dynamic loading conditions. Therefore, researchers in this field need to conduct experiments with advanced visualization techniques and data acquisition methods and develop new computational methods that can capture the complex interactions between particles and the surrounding environment.In conclusion, the mechanical behavior of granular materials is governed by the microstructure of the material, which continuously evolves during deformation. The research in this field has important implications for various applications, including geotechnical engineering, manufacturing, energy, and natural hazards. Future studies should focus on developing more accurate and robust models for predicting the mechanical behavior of granular materials under a wide range of conditions and investigating the influence of various factors on the microstructure and mechanical behavior of granular materialsIn addition to improving our understanding of the mechanical behavior of granular materials, future research should also focus on developing new techniques for characterizing their microstructure andproperties. Accurate quantification of grain size distribution, shape, and orientation can provide valuable insights into the behavior of granular materials and aid in the development of more realistic models. Techniques such as X-ray tomography, confocal microscopy, and laser scanning can provide three-dimensional information on the microstructure of granular materials, while acoustic and ultrasonic imaging can be used to infer properties such as density, porosity, and elastic moduli.Another important area of research is the study of granular materials under extreme conditions, such as high pressures, temperatures, and strain rates. Such conditions are encountered in many industrial processes, such as metal forging and welding, and in natural disasters, such as earthquakes and landslides. Understanding how granular materials behave under such conditions is essential for the development of safe and efficient manufacturing processes and for the prediction and mitigation of natural hazards.Finally, much can be gained by investigating the role of environmental factors on the mechanical behavior of granular materials. Factors such as moisture, temperature, and chemical composition can alter the mechanical properties of granular materials and theirresponse to external stimuli. For example, the presence of water can cause soils to become more plastic and increase their shear strength, while exposure to certain chemicals can lead to changes in bonding and surface tension of granular particles. Research into the effects of environmental factors on granular materials can provide valuable insights into their behavior in natural and industrial settings.In conclusion, granular materials are ubiquitous in nature and industry, and the study of their mechanical behavior offers many important insights into their properties and response to external stimuli. While our understanding of granular materials has advanced significantly in recent years, much remains to be learned about their microstructure and mechanical properties, particularly under extreme conditions and in the presence of environmental factors. Continued research in this field will have importantimplications for a wide range of applications, including geotechnical engineering, manufacturing, energy, and natural hazardsIn conclusion, granular materials are ubiquitous in nature and technology, and their macroscopic behavior is strongly influenced by their microstructure and response to external stimuli. Research in this fieldhas advanced in recent years, but there is still much to be learned about their mechanical properties and behavior under extreme conditions and in the presence of environmental factors. Further research in thisarea will have important implications for various applications, particularly in geotechnical engineering, manufacturing, energy, and natural hazards。

专利名称：PARTICLE SIMULATION SYSTEM 发明人：MATSUZAWA KAZUYA申请号：JP25489289申请日：19890929公开号：JPH03116385A公开日：19910517专利内容由知识产权出版社提供摘要：PURPOSE:To attain analysis with physically high accuracy while sufficiently utilizing the feature of particle simulation by arranging super-particles in an analysis area where particle simulation is executed, and afterwards, handling the transportation of particles, which are remained in the area without being distributed to the super- particles, according to the transient analysis method of fluid simulation. CONSTITUTION:When the transportation of the plural particles is analyzed, the super- particles are arranged in the analysis area, to which the particle simulation is executed, and afterwards, the transportation of the particles, which are remained in the area without being distributed to the super-particles, is handled according to the transient analysis method of the fluid simulation. When the analysis is executed by using a window method, the transportation for the small number of the particles, which are remained in a window without being distributed to the super-particles, and the particles out of the window is handled by the transient analysis method of the fluid simulation and when the supplying source of the particles such as an ohmic electrode, etc., is included in the analysis area, the transient analysis is executed under a contour condition that the number of the particles is fixed in a contour with the supplying source. Thus, while effectively utilizing the feature of the particle simulation, the simulation can be executed with physically high accuracy.申请人：TOSHIBA CORP 更多信息请下载全文后查看。

Particle Simulations
First Law: Every object continues in its state of rest, or of uniform motion in a straight line, unless compelled to change that state by external forces acted upon it.
Particle Simulations
These particles can reprபைடு நூலகம்sent different entities depending on the simulation. This could be:
atoms molecules dust particles snooker balls asteroids planets galaxies
Particle Simulations
More complex objects will require more complex collision detection and response systems.
For example – hockey pucks bouncing off each other.
Our algorithm to update a particle simulation is as simple as:
for(int i = 0; i < N; i++) { p[i] = p[i] + v[i] * h; }
Particle Simulations
Unfortunately this is a rather boring simulation. Each particle moves in a straight line and will never stop. To make it slightly more interesting we will have to enforce collisions. Collisions occur when a ‘hard’ particle hits into another ‘hard’ object.

COREX熔化气化炉风口回旋区CFD+DEM数值模拟孙俊杰;周恒;罗志国;邹宗树【摘要】COREX熔化气化炉风口回旋区是炉况顺行的基础,在冶炼过程中起着十分重要的作用,为了描述其形状和大小,建立了CFD+ DEM(Computational Fluid Dynamics and Discrete Element Method)耦合模型,对回旋区形成过程及大小进行了颗粒尺度的分析.得到床层高度为0.4m,气体速度11.74 m/s的条件下回旋区颗粒空隙度分布,当吹气时间为0.13s时,气体人口附近有颗粒被吹开,随着时间的推进,气体动能吹开的颗粒增多,0.19～0.21 s时,形成的回旋区开始稳定.对入口处不同气体速度条件下回旋区及其附近颗粒速度进行了计算模拟.模拟结果显示,风口附近颗粒在做回旋运动,并且随着入口气体速度的增大,吹开的颗粒增多,回旋区空腔增大,当入口气体速度为11.74 m/s和16.83 m/s时形成的回旋区较稳定,当入口气体速度大于21.90 m/s时形成的回旋区不太稳定.【期刊名称】《东北大学学报（自然科学版）》【年(卷),期】2013(034)006【总页数】5页(P824-827,844)【关键词】熔化气化炉;风口回旋区;DEM;CFD;回旋区颗粒速度【作者】孙俊杰;周恒;罗志国;邹宗树【作者单位】东北大学材料与冶金学院,辽宁沈阳110819;东北大学材料与冶金学院,辽宁沈阳110819;东北大学材料与冶金学院,辽宁沈阳110819;东北大学材料与冶金学院,辽宁沈阳110819【正文语种】中文【中图分类】TF557风口回旋区的形状和大小决定了高炉及COREX熔化气化炉中煤气的一次分布，反映了焦炭的燃烧状态，是炉况顺行的基础，在冶炼过程中起着十分重要的作用. COREX熔化气化炉风口回旋区的物理化学过程与高炉具有相似之处.关于COREX 熔化气化炉风口回旋区的研究较少，研究COREX时可以借鉴高炉风口回旋区的一些研究成果.Sarkar[1]把焦炭和球团矿简化为连续相，用双欧拉方法对高炉回旋区建模计算，此种方法未能准确描述颗粒的运动情况.Natsui[2]利用DEM，对高炉回旋区区域模拟计算，得到料层高度和风口间隔对死料柱高度的影响，此模型得到较精确的三维高炉回旋区部分参数，因其提前设置了回旋区的边界，但模型没有考虑气体流场对离散颗粒的影响，未实现CFD和DEM的耦合.罗志国等[3-4]建立了熔化气化炉物理模型，通过跟踪示踪颗粒的运动信息，得到观察面板处风口回旋区域的颗粒速度场，利用图像处理手段和分形理论，提供了一种确定熔化气化炉回旋区大小的方法.CFD+DEM耦合模型能够模拟工业中大部分同时存在颗粒和流体的流程，此种方法的最大优势是能够得到详细的单个颗粒信息(轨迹、相间力)，而这些正好是阐明复杂流体行为控制机理的关键信息[5].前人使用CFD+DEM方法对回旋区的研究都是针对高炉，未见有人用这种方法对COREX熔化气化炉的回旋区进行模拟计算.文献[6]建立了CFD+DEM耦合模型对高炉复杂的瞬变现象做了初步的模拟，分析了鼓风气体流量对回旋区处颗粒速度和所受力的影响(其计算相间力时仅考虑了曳力).回旋区内存在鼓风气体和焦炭颗粒，气体和焦炭颗粒的运动参数都会对回旋区形状和大小产生影响.本文以商业软件FLUENT为平台，用CFD+DEM方法对风口回旋区进行研究.利用UDF(user defined function)将DEM程序编译进FLUENT，从而实现了DEM与CFD的耦合，对回旋区内颗粒的运动分布规律模拟研究.1 数学模型描述1.1 颗粒相控制方程COREX熔化气化炉回旋区附近焦炭视为离散相，采用DEM模型计算，颗粒的平动和转动用牛顿第二运动定律描述[7-9]如下.(1)(2)式中：mi，Ii，vi和ωi分别代表颗粒i的质量(kg)，转动惯量(kg·m2)，平动速度(m/s)和转动速度(r/s)；fg代表i颗粒的重力(N)；fc，ij，fd，ij和Ti，j分别为颗粒i和j间的接触力(N)，阻尼力(N)以及转矩(N·m)；t为时间(s).当颗粒i与ki个颗粒同时接触时，其所受颗粒间的作用力和力矩矢量需要叠加在一起，fp-f，i是颗粒和流体的相间力，是颗粒相和流体相互作用力的合力，包括曳力、压力梯度力、虚假质量力、Basset加速度力、Magnus力、Saffman升力、滑移-剪切升力[5，10].1.2 气相控制方程风口处的鼓风气体视为连续相，采用类似于传统双流体模型方法进行描述.气相控制方程为计算单元上局部平均变量的质量和动量守恒形式，其表达式如下. +·(ρfεu)=0，(3)+·(ρfεuu)= -p+·(ετ)+ρfεg-fp-f，i.(4)式中：ρf，u分别代表气体密度(kg/m3)，气体速度(m/s)；τ为气体的应力张量(N/m2)；ε为颗粒孔隙度；p为气体压力(Pa).1.3 颗粒和气相的耦合在整个模型计算中，DEM用来计算每个颗粒的运动，气相流动的模拟在远大于颗粒尺寸的计算网格中进行，利用FLUENT提供的用户自定义函数UDF实现颗粒相与气相的耦合.具体耦合步骤如下.1) 使用DEM将颗粒堆积到模拟所需高度，用FLUENT初始化气体流场；2) 通过DEFINE_ON_DEMOND(name)将计算的颗粒初始结构调入，同时根据式(5)和图1计算网格空隙度：(5)式中：Vi表示颗粒i的体积；V表示单元格体积.图1 气相网格空隙度计算方法Fig.1 Calculation method of grid porosity3) 根据DEM得到的颗粒运动速度、位置和流体速度场计算出每个颗粒所受相间力，通过DEFINE_ADJUST(name，d)更新颗粒所受的力、颗粒速度、位置、转矩等的变化，完成DEM的颗粒计算；4) 通过DEFINE_PROPERTY(name，c，t)将密度ρf，黏度μ变为与空隙度ε有关的ε·ρf和ε·μ来实现对密度和黏度的修正；5) 将每个流体单元内所有颗粒所受相间力求和，得到此单元内的颗粒与流体相间力，将该相间力加到气相动量方程的源项中，进行流体运动计算，得出每个计算单元的速度值，根据网格单元中流体速度、颗粒速度更新流体与颗粒的相间力；6) 将每个颗粒所受相间力加入到DEM中就可以计算下一个时间步长内颗粒运动信息，如此迭代计算就可实现颗粒相和气相的持续耦合计算，计算流程如图2所示.图2 颗粒和气相耦合计算流程图Fig.2 Coupled calculation flow chart forparticle phase and gas phase2 模型建立及参数选择在实际物理实验中颗粒粒径0.002 5 m左右，床层高度0.4 m，一个风口所在区域的颗粒数已达到600 000个，受计算条件的限制，本模型简化为一个0.23m×0.005 m×1.0 m的矩形床，右侧墙壁上有一个0.025 m(5个颗粒厚度)×0.005 m的气体入口，入口中心距底部0.11 m.模型在x方向的网格尺寸为0.023 m，y方向的网格尺寸为0.005 m，z方向的网格尺寸为0.005 m.边界条件设置入口为速度入口，出口为outflow，壁面为无滑移边界条件.颗粒相模型中，颗粒-墙之间的碰撞与颗粒-颗粒之间的碰撞一致，只是墙不会产生位移.模拟参数如表1所示.表1 模拟参数Table 1 Simulation parameters颗粒相流体相颗粒形状球形μ1.8×10-5Pa·s颗粒密度950kg/m3ρf1.205kg/m3弹性模量216kPa11.74m/s 摩擦因数0.3入口速度16.83m/s阻尼因数0.1221.90m/s3 结果与分析3.1 回旋区的形成过程图3为气体喷吹倾角为4°，颗粒直径为0.005 m，床层填充高度为0.4 m，颗粒个数为5 700，气体速度为11.74 m/s条件下不同时刻回旋区颗粒运动及空隙度. 图3 床高0.4 m，气体速度11.74 m/s时回旋区形成过程Fig.3 The formation of raceway with bed height of 0.4 m and gas blowing velocity of 11.74m/s(a)—0.13 s时颗粒分布和空隙度； (b)—0.19 s时颗粒分布和空隙度； (c)—0.21 s时颗粒分布和空隙度.从图上可以看出，0.13 s时，气体入口附近开始有颗粒被吹开，空隙度变大，对比两张图可以看出空隙度变大的范围比颗粒吹开的范围大，这是因为空隙度的计算是以网格为基础，没有颗粒的位置及其附近同属于一个单元网格，导致没有颗粒的位置空隙度不等于1，且其附近位置的空隙度变大.随着时间的推进，气体动能吹开的颗粒增多.0.19～0.21 s时，形成的回旋区开始稳定.3.2 回旋区颗粒速度图4为气体喷吹倾角为4°，颗粒直径为0.005 m，床层填充高度为0.4 m，颗粒个数为5 700，回旋区基本稳定条件下不同气体速度时回旋区及其附近颗粒位置及速度.图4 不同气体速度下回旋区及附近颗粒位置及速度值Fig.4 Particle position and velocity in and around raceway with different gas blowing velocity(a)—11.74 m/s； (b)—16.83 m/s； (c)—21.90 m/s.从图中可以看出，在重力、颗粒间摩擦力和气体动能的共同作用下，风口附近颗粒在做回旋运动，且随着入口气体速度的增大，吹开的颗粒增多，回旋区空腔增大；回旋区内部无颗粒，外围颗粒速度基本为0；当入口气体速度为11.74 m/s和16.83 m/s时形成的回旋区较稳定，当入口气体速度大于21.90 m/s时，床层中的大部分颗粒有向上的速度，形成的回旋区不太稳定，有形成流化床的趋势.4 结论本文通过CFD+DEM耦合方法对风口回旋区的形成进行模拟研究，利用UDF将DEM程序编译进FLUENT，从而实现了DEM与CFD的耦合.对回旋区内颗粒的分布运动规律进行了模拟，得到床层高度为0.4 m，气体速度11.74 m/s的条件下回旋区颗粒的空隙度分布.0.13 s时，气体入口附近开始有颗粒被吹开，随着时间的推进，气体动能吹开的颗粒增多，0.19～0.21 s时，形成的回旋区开始稳定；对入口处不同气体速度条件下回旋区及其附近颗粒速度进行了计算模拟.结果显示，在重力、颗粒间摩擦力和相间力的共同作用下，风口附近颗粒在做回旋运动，且随着入口气体速度的增大，吹开的颗粒增多，回旋区空腔增大；当入口气体速度为11.74 m/s和16.83 m/s时形成的回旋区较稳定，当入口气体速度大于21.90m/s时形成的回旋区不太稳定.参考文献：[1] Sarkar S，Gupta G S，Kitamura S.Prediction of raceway shape and size[J].ISIJ International，2007，47(12)：1738-1744.[2] Natsui S， Ueda S， Oikawa M， et al.Optimization of physical parameters of discrete element method for blast furnace and its application to the analysis on solid motion around raceway[J].ISIJ International，2009，49(9)：1308-1315.[3] 罗志国，孙俊杰，狄瞻霞，等.COREX熔化气化炉风口回旋区内表面积研究[J].东北大学学报：自然科学版，2012，33(6)：840-843.(Luo Zhi-guo，Sun Jun-jie，Di Zhan-xia，et al.Study on surface area of raceway in COREX melter-gasifier[J].Journal of Northeastern University：Natural Science，2012，33(6)：840-843.)[4] Luo Z G， Di Z X， Zou Z S， et al.Application of fractal theory on raceway boundary in COREX melter-gasifier model[J].Ironmaking and Steelmaking，2011，38(6)：417-421.[5] Zhou Z Y，Kuang S B，Chu K W，et al.Discrete particle simulation of particle-fluid flow：model formulations and their applicability[J].Journal of Fluid Mechanics，2010，661：482-510.[6] Zhou Z Y，Zhu H P，Yu A B，et al.Numerical investigation of the transient multiphase flow in an ironmaking blast furnace[J].ISIJInternational，2010，50(4)：515-523.[7] 宋伟刚，陈洪亮，李勤良，等.散状物料转载冲击载荷的DEM仿真[J].东北大学学报：自然科学版，2011，32(11)：1631-1634.(Song Wei-gang，Chen Hong-liang，Li Qin-liang，et al.Simulation of handling impact load of bulk material by DEM method[J].Journal of Northeastern University：Natural Science，2011，32(11)：1631-1634.) [8] Luo Z G，Li H F.Mathematical simulation of burden distribution in COREX melter gasifier by discrete element method[J].Journal of Iron and Steel Research International，2012，19(9)：36-42.[9] Wang X，You C F.Evaluation of drag force on a nonuniform particle distribution with a meshless method[J].Particuology，2011，9(3)：288-297.[10]Zhu H P，Zhou Z Y，Yang R Y，et al.Discrete particle simulation of particulate systems：theoretical developments[J].Chemical Engineering Science，2007，62(13)：3378-3396.。

专利名称：Discrete particle electrolyzer cathode andmethod of making same发明人：Stuart I. Smedley,Martin De TezanosPinto,Stephen R. des Jardins,Donald JamesNovkov,Ronald Gulino申请号：US10424539申请日：20030424公开号：US20040168922A1公开日：20040902专利内容由知识产权出版社提供专利附图：摘要：A system for producing metal particles using a discrete particle electrolyzercathode, a discrete particle electrolyzer cathode, and methods for manufacturing the cathode. The cathode has a plurality of active zones on a surface thereof at least partially immersed in a reaction solution. The active zones are spaced from one another by between about 0.1 mm and about 10 mm, and each has a surface area no less than about 0.02 square mm. The cathode is spaced from an anode also at least partially immersed in the reaction solution. A voltage potential is applied between the anode and cathode. Metal particles form on the active zones of the cathode. The particles may be dislodged from the cathode after they have achieved a desired size. The geometry and composition of the active zones are specified to promote the growth of high quality particles suitable for use in metal/air fuel cells. Cathodes may be formed from bundled wire, machined metal, chemical etching, or chemical vapor deposition techniques.申请人：SMEDLEY STUART I.,PINTO MARTIN DE TEZANOS,JARDINS STEPHEN R. DES,NOVKOV DONALD JAMES,GULINO RONALD更多信息请下载全文后查看。

1 NTNU Dep. Petroleum Technology and Applied Geophysics, N-7491 Trondheim - Norway 2 SINTEF Petroleum Research, N-7465 Trondheim - Norway e-mail: liming.li@iku.sintef.no - runeh@ience and Technology – Rev. IFP, Vol. 57 (2002), No. 5
– Stress-dependency of the acoustic properties of a bonded particle assemblage has been studied.
Résumé — Mécanique du réservoir à l’échelle des grains — Cet article résume quelques développements et applications de la modélisation par particules discrètes, basés sur le code commercial PFC (Particle Flow Code). Nous décrivons ci-après les exemples ayant un rapport avec le comportement géomécanique du réservoir : – Des résultats de simulations numériques et d’expériences sont présentés, illustrant une évolution non-
compaction and stress path for simulated reservoir rock, i.e. materials that are formed under stress. This nonlinearity is often lost in core measurements, as a result of stress release during coring (core damage). – Angular and breakable superparticles have been implemented in the particle model and applied to simulate the formation of a compaction band at high stress level. The model is also used to investigate the formation of a shear band at lower confinement. – Numerical scheme of fluid coupling has been developed in both 2D and 3D models and applied to study stress-dependent permeability.

In view of the importance of soil crushability to stresslevel effects, further studies have been carried out to understand the statistics of individual grain breakages. Results of single-particle crushing tests have been compared with the one-dimensional compression of samples of the same u et al., 2001). They showed
de la vitesse de chargement sur la position des courbes de compression sous un angle pratique : arriver a` mener
des simulations DEM aussi rapidement que possible sans cre´er d’erreurs d’inertie. D’autres chemins de contrainte, incluant celui de compression triaxiale ‘draine´e’, celui a` contrainte moyenne effective constante et celui a` volume constant ont e´galement e´te´ simule´s a` partir de diffe´rentes contraintes initiales le long de la courbe vierge e–log p9. Nous comparons les re´sultats nume´riques de manie`re quantitative avec les donne´es des essais sur le sable siliceux. Ces re´sultats sont e´galement e´tudie´s afin de caracte´riser les sols sujets au broyage lors d’une compression normale, a` la rupture et a` l’e´tat critique. Les simulations de sols sujets a` l’e´crasement, utilisant la me´thode DEM, donnent des informations pre´cieuses sur les origines microme´caniques de la plasticite´ des sols.

不同筛面形式振动筛颗粒运移的DEM模拟吴先进;李圣一;方潘;席仲君;侯勇俊;刘有平【摘要】建立了凸、凹和平直三种筛面形式的振动筛离散元筛分模型,并采用JKR 湿颗粒接触模型模拟了钻井液振动筛上的颗粒运移与透筛特性.研究结果表明,平直筛面的振动筛的筛分能力良好,筛面上分固相颗粒具有运移速度快、分布均匀特点;凹形筛面筛分能力较差,固相颗粒容易在筛面凹部堆积;凸形筛面振动筛的筛分能力也弱于平直筛面的振动筛,固相颗粒也容易在筛面两侧堆积.研究结果对钻井振动筛设计和筛分能力评价具有重要参考意义.【期刊名称】《机械研究与应用》【年(卷),期】2017(030)005【总页数】5页(P45-49)【关键词】振动筛;张紧方式;离散元法;JKR湿颗粒【作者】吴先进;李圣一;方潘;席仲君;侯勇俊;刘有平【作者单位】四川宝石机械专用车有限公司,四川广汉 613800;东方锅炉股份有限公司,四川自贡 643001;西南石油大学机电工程学院,四川成都 610500;川庆钻探工程有限公司川西钻探公司,四川成都 610051;西南石油大学机电工程学院,四川成都 610500;四川宝石机械专用车有限公司,四川广汉 613800【正文语种】中文【中图分类】TE926在钻井液振动筛的筛分过程中，固相颗粒运移是一个复杂的动态过程。

朱维兵[1-3]等用传统的筛分颗粒运移模型，从理论上对筛面上的颗粒运移进行了研究。

姚恒申[4]等针对筛面固相颗粒抛掷规律作了理论研究。

目前研究模型的缺陷只是考虑了筛面上的颗粒运移情况，无法分析不同粒度颗粒的运移和透筛行为。

DEM法是近年来出现的离散元研究方法，非常适合于模拟振动筛分过程颗粒的力学行为。

Delaney, G.W.、Dong, K.J、Fernandez, J.W等[5-9]研究了颗粒在振动平直筛面上运移特性，并建立模型证明了DEM方法对研究颗粒碰撞筛面有良好的模拟效果，以及球型颗粒模拟的可行性。

弯管稀相气力输送CFD-DEM法数值模拟杜俊;胡国明;方自强;范召【摘要】利用 CFD-DEM方法对稀相气固两相流在带有弯管的气力输送管道内的流动特性数值模拟。

通过 CFD求解连续气相流场，使用 DEM求解离散颗粒相运动及受力，为提高仿真速度，模型忽略了空隙度对气相的影响。

仿真结果显示了颗粒在弯管内形成颗粒绳及在垂直管道内颗粒绳分散的过程，并获取了颗粒-颗粒、颗粒-壁面的碰撞信息。

对比弯管上下两部分的碰撞情况，颗粒及壁面在弯管下半部分受到的撞击和磨损更严重。

通过对气力输送各类工作参数的研究发现，其对管道内气固两相流的流动特性和碰撞情况有着不同程度的影响。

气流速度对颗粒绳分散影响甚微，但对弯管内颗粒碰撞强度有明显影响。

随着颗粒质量流量的增加，形成的颗粒绳更紧凑、分散速度更慢，并在弯管中形成了阻碍颗粒-壁面碰撞的防护层。

弯径比的增加也能加强颗粒绳的紧凑度，减缓颗粒绳的分散速度。

%CFD-DEM model was used to simulate the gas-solid flow in dilute pneumatic conveying with bends.In CFD-DEM,the discrete particle phase was obtained by Discrete Element Method (DEM),and the flow of continuum gas phase was determined by Computational Fluid Dynamics (CFD).In order to consume less simulation time,the effect of the particle solid fraction on the gas phase was not taken into account. The results showed the particle rope formation in the bend and its dispersion in the vertical pipe,and obtained the particle-particle and particle-wall collision parison of the collision information in the bend,the collision and abrasion seemed more intensive at the bottom of the bend.It was also found that the geometry and parameters have differentmagnitude effects on the gas-solid flow and collisions in the pipe.The gas velocity was considered to be limited influence on the rope dispersion,but significant effect on the collisions in bend section.With the increase of the solid mass flow,the particle rope seemed stronger and more dispersed at a low rate,and there was a shield formed to impede particle-wall collision at the bend section.The increase of the bend radius ratio also made the particle rope stronger,and the dispersion rate lower.【期刊名称】《国防科技大学学报》【年(卷),期】2014(000)004【总页数】6页(P134-139)【关键词】离散单元法;计算流体力学;气力输送;稀相【作者】杜俊;胡国明;方自强;范召【作者单位】武汉大学动力与机械学院，湖北武汉 430072;武汉大学动力与机械学院，湖北武汉 430072;武汉大学动力与机械学院，湖北武汉 430072;武汉大学动力与机械学院，湖北武汉 430072【正文语种】中文【中图分类】TH48弯管作为气力输送系统中的典型零部件，对输送管道内气固两相流的流动特性有着显著影响，研究弯管内气固两相流的流动特性对设计和生产者优化、改善气力输送系统具有重要意义。

Particle Simulations
Particle Simulations
These particles can represent different entities depending on the simulation. This could be:
atoms molecules dust particles snooker balls asteroids planets galaxies
Particle Simulations
Once a collision has been detected the system must respond to the collision. For our hockey pucks, simply reverse the velocity in the direction of the collision.
Particle Simulations
• main loop – for all particles • ‘move’ particle • if ‘collision’ with boundary – respond to collision – for all particles • for all other particles – if ‘collision’ between particles » respond to collision
Particle Simulations
Second Law: The acceleration a of a body is parallel and directly proportional to the net force F acting on the body, is in the direction of the net force, and is inversely propotional to the mass m of the body.

Applied Physics 应用物理, 2017, 7(4), 97-102Published Online April 2017 in Hans. /journal/app https:///10.12677/app.2017.74014文章引用: 吴达岭, 桂南. 非球形颗粒碰撞过程硬颗粒与软颗粒模型比较研究[J]. 应用物理, 2017, 7(4): 97-102.Comparison of Collision Dynamics of Non-Spherical Particles between the Hard-Particle Model and Soft-Particle ModelDaling Wu 1, Nan Gui 21Zhejiang Academy of Safety Science and Technology, Hangzhou Zhejiang2Institute of Nuclear and New Energy Technology, Tsinghua University, Key Laboratory of Advanced Reactor Engineering and Safety, Ministry of Education, BeijingReceived: Apr. 6th , 2017; accepted: Apr. 20th , 2017; published: Apr. 27th, 2017AbstractNon-spherical particles and complex non-spherical particle systems are commonly encountered in energy chemical engineering, mineral engineering, and safety manufacturing engineering. The collision dynamics of non-spherical particles are fundamental issues in the investigations of com-plex particle systems. Discrete element method is the most prevalent models for particle-particle collision of spherical shapes, which incorporates a spring, a dashpot, and a frictional interaction. In this work, a comparative and analytical study on the mathematical models of collision dynamics for non-spherical particles is carried out. The hard non-spherical particle model and the soft-par- ticle model have been compared. It indicates that agreeable simulation results can be found for the two models. The two models have different advantages or disadvantages in accuracy, compu-tational stability, and feasibility, with particular suitability for different cases.KeywordsEnergy Engineering, Chemical Engineering, Safety Manufacturing, Non-Spherical, Collision Model, Discrete Element Method, Hard Particle, Soft-Particle非球形颗粒碰撞过程硬颗粒与软颗粒模型比较研究吴达岭1，桂 南21浙江省安全生产科学研究院，浙江 杭州吴达岭，桂南2清华大学核能与新能源技术研究院，先进反应堆工程与安全教育部重点实验室，北京收稿日期：2017年4月6日；录用日期：2017年4月20日；发布日期：2017年4月27日摘要在能源化工过程、矿业冶金及安全生产过程中均会涉及复杂的实际非球形颗粒体系，而颗粒的非球形碰撞动力学过程是对复杂颗粒体系进行研究必须要考虑的基本问题。

Advances in Geosciences 地球科学前沿, 2022, 12(12), 1646-1659 Published Online December 2022 in Hans. /journal/ag https:///10.12677/ag.2022.1212160颗粒形状对隧道开挖砂性地层稳定性影响的 力学响应分析高 昆1，高晓耕21中铁隧道集团二处有限公司，河北 廊坊2同济大学土木工程学院地下建筑与工程系，上海收稿日期：2022年11月16日；录用日期：2022年12月22日；发布日期：2022年12月31日摘要颗粒形状是影响土体强度及变形特性的重要因素，现有研究主要聚焦于颗粒形状对土体物理力学性质的影响。

为研究颗粒形状对隧道开挖后砂性地层稳定性的影响，选取典型的颗粒形状表征参数构建离散单元模型模拟隧道开挖，系统分析颗粒形状对隧道开挖后砂性地层应力重分布、地层变形及地层–隧道衬砌相互作用的演化规律。

研究结果表明：隧道开挖后砂性地层会产生一定的土拱效应，土拱效应强弱由地层沉降及地层–隧道相互作用表征，尤其是土体的扁平颗粒比方正颗粒在隧道开挖后更有利于形成稳定的土拱。

关键词隧道工程，土拱效应，颗粒形状，稳定性，力学响应Mechanical Response Analysis Due to the Influence of Particle Shape on the Stability of Sandy Ground after Tunnel ExcavationKun Gao 1, Xiaogeng Gao 21China Railway Tunnel Group No. 2 Co., Ltd., Langfang Hebei2Department of Geotechnical Engineering, Tongji University, ShanghaiReceived: Nov. 16th, 2022; accepted: Dec. 22nd, 2022; published: Dec. 31st, 2022AbstractParticle shape is an important factor affecting the strength and deformation characteristics of soils. The existing researches in this field mainly focus on the influence of particle shape on the高昆，高晓耕physical and mechanical properties of soils. In order to explore the influence of particle shape on the stability of sandy ground after excavation, a series of discrete element models are established simulating tunnel excavation in ground with different representative particle shape descriptors. The influences of particle shape on evolutions of ground stress redistribution, ground settlement and ground-tunnel structure interaction after excavation are comprehensively investigated. The results show that there exists certain soil arching effect in cohesion less ground after tunnel exca-vation. The soil arching effect can be manifested from the ground settlement and ground-tunnel structure interaction. Especially, flat particles are more favorable for the formation of stable arching structures than squared particles.KeywordsTunnel Engineering, Soil Arching, Particle Shape, Stability, Mechanical ResponseThis work is licensed under the Creative Commons Attribution International License (CC BY 4.0)./licenses/by/4.0/1. 引言土体是一种典型的颗粒材料，具有明显的非连续性、不均匀性和各向异性。

本栏目编辑　翟小华破・磨80Solutions，2010：112-114.[4]　Jayasundara C T，Yang R Y，Yu A B. Discrete ParticleSimulation of Particle Flow in the IsaMill Process [J]. American Chemical Society，2006(45)：6349-6359.[5]　Jayasundara C T，Yang R Y，YuA B，et al. Effects of discrotation speed and media loading on particle flow and grinding performance in a horizontal stirred mill [J]. International Journal ofMineral Processing，2010(96)： 27-35.[6]　Jayasundara C T，Yang R Y，Yu A B，et al. Discrete particlesimulation of particle flow in IsaMill-Effect of grinding medium properties [J]. Chemical Engineering Journal，2008(135)： 103-112. □(收稿日期：2012-09-14)(修订日期：2012-10-23)球磨机钢球抛落点切向分速度分析刘基博衡阳有色冶金机械总厂　湖南衡阳　421002摘要：通过对 φ3.6 m×6.0 m 球磨机各球层动球抛落点的切向分速度及其对定球的相对速度和动能进行具体计算，证明了切向分速度产生的动能一部分反馈给了筒体，而其相对速度产生的动能则是对矿物进行研磨的主要动力源。

本文的分析结果可为今后进一步研究提供理论参考。

Eurographics/ACM SIGGRAPH Symposium on Computer Animation(2005)K.Anjyo,P.Faloutsos(Editors)Particle-Based Simulation of Granular MaterialsNathan Bell1,Yizhou Yu1and Peter J.Mucha21University of Illinois at Urbana-Champaign2Georgia Institute of TechnologyAbstractGranular materials,such as sand and grains,are ubiquitous.Simulating the3D dynamic motion of such mate-rials represents a challenging problem in graphics because of their unique physical properties.In this paper we present a simple and effective method for granular material simulation.By incorporating techniques from physical models,our approach describes granular phenomena more faithfully than previous methods.Granular material is represented by a large collection of non-spherical particles which may be in persistent contact.The particles represent discrete elements of the simulated material.One major advantage of using discrete elements is that the topology of particle interaction can evolve freely.As a result,highly dynamic phenomena,such as splashing and avalanches,can be conveniently generated by this meshless approach without sacrificing physical accuracy.We generalize this discrete model to rigid bodies by distributing particles over their surfaces.In this way,two-way coupling between granular materials and rigid bodies is achieved.Categories and Subject Descriptors(according to ACM CCS):I.3.5[Computer Graphics]:Computational Geometry and Object Modeling—Physically based modeling I.3.7[Computer Graphics]:Three-Dimensional Graphics and Realism—Animation I.6.8[Simulation and Modeling]:Types of Simulation—Animation1.IntroductionGranular materials,such as sand,soil,powders and grains, present an indispensable part of the real world.They are sim-ply very large ensembles of macroscopic particles[JNB96]. Given their ubiquity,there is an inevitable need for granu-lar material simulation for graphics-related applications such as building virtual environments and animating natural phe-nomena.In fact,there have been a few inspirational tech-niques on this topic in graphics[LM93,CLH96,SOH99, ON03].However,these techniques usually do not have suffi-cient spatial or physical accuracy to describe important gran-ular phenomena.More mature techniques for3D environ-ments still need to be developed.In this paper,we aim to faithfully perform3D dynamic simulations of granular ma-terials.In such simulations,we not only need to simulate the interactions,such as collisions and friction,among the parti-cles,but also the interactions between granular particles and other larger-scale objects.In addition to their role in graphics,granular materials have profound industrial significance.Indeed,the grinding of particles and ores alone accounts for an estimated1.3% of U.S.electrical power consumption[JGD94].Hence there is considerable interest in developing models to better un-derstand granular behavior.Different fromfluids and other materials,granular materials exhibit unique physical proper-ties which demand novel simulation techniques.For exam-ple,a sand pile at rest with a slope lower than the angle of repose behaves like a solid,i.e.,the material remains at rest even though gravitational forces create stresses on its surface and in its interior.If the pile is tilted a few degrees above the angle of repose,grains start toflow.However,the main flow exists in a boundary layer at the pile’s surface with only limited motion in the interior.As another example,when a granular material is held in a cylindrical container,the pres-sure at the bottom of the container does not necessarily—as influids—increase indefinitely as the height of the material is increased[dG99].Instead,for a sufficiently tall column, the pressure saturates to a maximum value independent of the height.Because of the contact forces among grains and static friction with the sides of the container,the container walls support the extra weight.It is this feature that allows the sand in an hourglass toflow at a nearly constant rate. Clearly,an accurate model for granular materials must be able to reliably reproduce such phenomena.c The Eurographics Association2005.Figure1:Friction forces produce a constant rate offlow through the hourglass.In this paper,we present a simple and effective simula-tion method built upon both theoretical and experimental re-sults in physics.We adopt the practice of modeling gran-ular material with discrete elements represented by parti-cles[CS79,HL98],while selecting interparticle interactions based in part on computational cost.This approach permits faithful reproduction of a wide range of both static and dy-namic granular phenomena.Furthermore,the discrete nature of our method obviates the need for an explicit surface repre-sentation,thereby allowing the simulation of sparse events. Interactions among particles are governed by a molecular-dynamics based contact model.Contact forces are obtained explicitly from the relative velocity and overlap between par-ticles.The formulations of these forces are derived from elasticity theory and experimental results.One major advan-tage of this contact model is that large numbers of bodies in persistent contact are handled efficiently.We extend this contact model to rigid bodies by covering their surfaces with particles.A signed distance representa-tion is used to produce a uniform distribution of particles and allow particles to be placed at an offset from the original mesh.The force and torque accumulated from all these sur-face particles are used to integrate the object through rigid body motion.This produces true two-way coupling between granular particles and rigid bodies.As a by-product,interac-tions between different rigid bodies can be simulated using the same particle-based pelling demonstra-tions for these types of interactions have been produced.1.1.Related Work in GraphicsUsing particles to simulate natural phenomena has a long history in computer graphics.Early work,such as[Ree83],developed techniques to animate and render irregular objects with particle systems.However,it was the introduction of particle systems with pairwise interactions that made it pos-sible to simulatefluids and deformable objects.[MP89]sim-ulated viscousfluids with particle interactions based on the Lennard-Jones potential force.[TPF89]paired particles to better simulate deformable objects.[Ton91]improved parti-cle motion by adding additional particle interactions based on heat transfer among particles.There have also been suc-cessful applications of more accurate particle-basedfluid models in graphics.[DC99,MCG03]adapted smoothed par-ticle hydrodynamics(SPH)[Mon94]to compute the motion of deformable substances and liquids.[HMT01]modified the formulation of SPH to simulate hair-hair interactions by modeling hair as afluid-like continuum.Recently,another particle method,Moving Particle Semi-Implicit(MPS),was also developed[YKO96]along with its application to in-compressiblefluid simulation in graphics[PTB∗03]. Although being effective forfluids and deformable ob-jects,these techniques are not appropriate for granular ma-terials for the following reasons.First,the formulations for particle collision were not specifically designed according to the real material properties of granular particles.Some of them were designed to reflectfluid properties,such as zero divergence(though the divergence in granular materials is typically small).Second,static friction among particles is missing from these methods.Only cohesive forces were used to keep particles together.However,static friction and cohe-sion give rise to different effects.Fluid viscosity in some of the methods only resembles dynamic friction. Nevertheless,there have been a few inspirational tech-niques focusing on granular materials in graphics[LM93, CLH96,SOH99,ON03].These techniques usually do not have sufficient spatial or physical accuracy.The techniques in[LM93,CLH96]were derived from physics,but they model soil as a heightfield which precludes the simulation of many interesting3D effects,such as splashing.While the methods in[SOH99,ON03]can produce interesting granular behavior for footprints or sandboxes,their empirical nature prevents them from faithfully simulating high-speed scenar-ios where significant energy is exchanged among particles and objects.While the proposed method is in many ways similar to [LHM95],there are a number of crucial differences.Firstly, our model for particle interaction is more rigorously devel-oped than the damped linear spring force used in their pa-per.Secondly,the non-spherical particles in their method are assembled using stiff springs whereas we consider each grain as a rigid object.Thirdly,we extend the particle-based method to handle two-way coupling with rigid bodies.While their method can produce compelling animations using a small number of particles,we remark that there is no natural extension of their multi-scale aspect to three dimensions,c The Eurographics Association2005.since particles cannot be assumed to be in contact with a ground plane.1.2.Physical Models of Granular MaterialsThere are a number of computational models used to inves-tigate granular materials.Event-driven(ED),or hard-sphere algorithms,are based on the calculation of changes from dis-tinct collisions between particles(typically spheres or poly-hedra).The efficiency and utility of ED methods rely on the assumptions of binary collisions,short contact durations, and relatively long intervals between consecutive events. Hence such methods are well suited to simulate rapidflows, whose behavior resembles that of a gas.However in many cases of interest some or all of these assumptions are vio-lated.In systems with strong dissipation the time between collisions may vanish,resulting in so-called inelastic col-lapse[HL98].Likewise,ED methods are inappropriate for systems with persistent contacts and frequent collisions such as heaps,hoppers,and shaken materials.Since modeling individual particles explicitly is often in-feasible(e.g.fine powders),there is strong interest in con-tinuum models for granular media.Avalanches,landslides, and other large-scale phenomena can in principle be sim-ulated using continuum techniques[AT01].Unfortunately these techniques are relatively new and have not yet repro-duced the rich small-scale details available in the results of their discrete counterparts.Molecular Dynamics(MD)or soft-sphere methods are an-other common way to simulate granular materials.Like ED, MD methods are discrete,with particles also represented by spheres or polyhedra.In ED methods,contacts are resolved through an instantaneous change in the momentum of col-liding particles,thereby enforcing non-interpenetration.In contrast,MD methods permit small overlap between parti-cles to simulate the elastic deformation of material elements. This overlap is then used along with other physical quantities to determine an appropriate contact force.This distinction allows MD methods efficiently model systems with persis-tent contacts and frequent collisions.Moreover,MD meth-ods have demonstrated a wide variety of both static and dy-namic granular phenomena such as heap formation,patterns in vibrating layers,force chains,and particle size segrega-tion[Lee94,SvHSvS04,DR99].Therefore,the simulation technique in this paper is based on MD methods.2.Contact ForcesMD based models are characterized by their handling of con-tact forces.Since our objects will be composed of collections of round particles,it suffices to describe the contact force be-tween a single pair of spheres,computing the total force on a given sphere as the sum over its pairwise interactions.We first define the relevant physical quantities used in the con-tact models and then examine a few common forceschemes.Figure2:Quantities used in collision modeling.We divide the contact force−→F=−→F n+−→F t into normal and tangential components.Let x1and x2be the centers of particles of radius R1and R2respectively.The forces described here are applied to the particle centered at x1.An equal and opposite force is ap-plied to the particle centered at x2.Define the overlapξof the particles and the normalized line of centers−→N as follows:ξ=max(0,R1+R2− x2− x1 ),(1)−→N=x2−x1x2−x1.(2)Let v1and v2denote the velocity of each particle at the contact point and define the relative velocity−→V,relative ve-locity in the normal direction˙ξ,and tangential velocity−→V t:−→V= v1− v2,(3)˙ξ=−→V·−→N,(4)−→Vt=−→V−˙ξ−→N.(5) 2.1.Normal ForcesIn this section we describe a number of frequently used nor-mal forces and their properties.We restrict our attention to force rules of the form−→Fn=f n−→N,(6)f n+k dξα˙ξ+k rξβ=0,(7) with viscous damping coefficient k d and elastic restoration coefficient k r.Particle stiffness is controlled by k r and dis-sipation during collisions by k d.In the simplest case we chooseα=0,β=1and obtain a damped linear spring force,f n+k d˙ξ+k rξ=0.(8) Here k r and k d can be chosen to match experimental data, or manually tuned to produce the desired response.Given a coefficient of normal restitutionεn and contact duration t c the corresponding coefficients can be determined:k d=2m eff−lnεt c,(9)c The Eurographics Association2005.k r=m efft c2(ln2ε+π2),(10)1 m eff =1m1+1m2.(11)Here m eff represents the reduced mass,given particle masses of m1and m2respectively.While the linear model is adequate for many purposes,more accurate options ex-ist.Drawing from elasticity theory,we can chooseβ= 3/2which corresponds to the Hertz contact force between spheres[LL86],f n+k dξα˙ξ+k rξ3/2=0.(12) Hertz theory relates the stiffness coefficient k r to the mate-rial properties of the spheres in contact.Let E1and E2be the Young’s moduli andµ1andµ2the Poisson ratios of the two spheres.Now define the effective radius R eff,Young’s modulus E eff,and corresponding stiffness coefficient k r:1 R eff =1R1+1R2,(13)1 E eff =1−µ21E1+1−µ22E2,(14)k r=43E eff R eff.(15)So far,we have leftαunspecified in the non-linear model.In the linear model(α=0,β=1),the coefficient of nor-mal restitutionεn was constant for a given choice of pa-rameters.However both experimental and theoretical re-sults demonstrate a dependence ofεn on the impact veloc-ity[RPBS99,SP98].These results suggest a relationship like (1−εn)∝v1/5[SDW96];qualitativelyεn decreases with increasing impact velocity.In the non-linear model,differ-ent values ofαgive rise to different dependencies.Figure3demonstrates this relationship for common parameter pairs.Constantεn are obtained in both the linear model and thenon-linear model withα=1/ingα=0,β=3/2hasthe unnatural effect of producing more elastic collisions withincreasing velocity.Of the four,the model that best agreeswith experimental results isα=1/2,β=3/2.This choiceis comprehensively studied in[SP98]and[RPBS99].2.2.Shear ForcesIf desired,one can introduce a force to account for shearfriction between particles.The most basic approach simplydefines a force opposing the tangential velocity.−→Ft=−k t −→Vt(16)Here k t is a viscous damping term not unlike k d.This rule can be problematic since −→F t is not limited by the normal force.A contrasting option accounts for this Coulomb lawFigure3:The coefficient of normal restitution is indepen-dent of the velocity in both the linear modelα=0,β=1and non-linear model withα=1/4.The non-linear model with α=1/2becomes more dissipative with increasing velocity, in agreement with experimental results.However choosing α=0,β=3/2produces the opposite behavior.In all tests k r=1.0E6and k d={0,50,100,250,500}.for friction coefficientµand normal force f n,−→Ft=−µf n−→Vt−→Vt.(17)The numerical difficulties encountered by the discontinuity at −→V t =0can be alleviated by combining(16)with(17):−→Ft=−min(µf n,k t−→Vt)−→V t−→V t .(18) Here the magnitude of k t affects oblique impacts and should be large enough that the behavior of force(18)resembles that of(17)[SDW96].The shear friction forces discussed above can only slow movement in the tangent direction,not stop or reverse it. These rules are inadequate for applications that require truly static friction,such as heap formation,angles of repose,and the quasi-static stages of avalanches.In such situations,there is a threshold force below which the grains do not move at all,opposed by static friction originating in the even smaller scale surface interactions of particle contact.However,im-plementation of a simple threshold rule is history dependent and computationally intensive.A common solution is to in-troduce a virtual spring when particlesfirst touch[LH93]. While the particles remain in contact,the total displacement −→D from the initial contact point is used to compute the shearc The Eurographics Association2005.force:−→F t =−min (µf n ,k t −→D )−→D −→D.(19)Another approach purely relies on particle geometry,us-ing non-spherical particles to create the necessary rough-ness for static friction.Non-spherical particles are often con-structed of polygons [PB95]or spheres connected through springs [PB93],though the latter can become costly for sufficiently stiff springs.The spatial configuration of these spheres or polygons represents the microstructure of non-spherical particles.No explicit static friction force need be used for these compound non-spherical objects.When two non-spherical particles are in contact with each other,the normal forces between their rough structures produce a re-sistance to relative tangential motion,thus creating actual effects of static friction.In order to draw upon the existing theory for spherical particles,we model non-spherical parti-cles with collections of spheres,avoiding overly stiff inter-actions by replacing internal springs with rigid motion con-straints.With static friction obtained via normal forces be-tween compounds of spheres,we can then limit our attention to dynamic shear forces as expressed in (18).3.Granular Material SimulationWe directly model a granular material as a large collection of granular particles which may be in persistent contact when the material is static.However,new external forces may cause some or all of the particles to have relative motion,leading to momentum and energy exchange through colli-sions.The aforementioned contact forces are used for com-puting interactions between both static particle configura-tions and dynamic assemblies.3.1.Granular ParticlesOur model of granular particles is inspired by the non-spherical particles in [PB93].However,instead of using springs connecting round particles,we more simply model a grain as multiple round particles that are constrained to move together as a singlerigid body.Such rigid collectionsFigure 4:Non-spherical particles are used to accurately model static granular phenomena.Here tetrahedron and cube shapes are constructed from spheres.of spheres have been considered in discrete element simu-lations before [WB93],but were unable with then-presentcomputers to simulate the broad scope of granular phenom-enon which can now be simulated.The positions of the round particles are given as offsets from the body’s center of mass.Let o i denote the offset of a given particle,R the current orientation matrix,and x the center of mass of the grain.The current position of the particle is simply x +R o i .Dissipation in the contact force model depends on the relative velocity at the point of con-tact p .We can determine the velocity at the contact point from the velocity v and angular velocity ωof the grain usingv +ω×( p − x ).A contact force −→F applied at a contact pointp contributes a torque ( p − x )×−→F to the body.The net force on a grain is simply the sum of the individual contact forces applied to its particles.The motivations for such a non-spherical grain model are as follows.First,non-spherical particles produce a consid-erable angle of repose and exhibit stick-slip behavior,in ac-cordance with experimental results [PB93].Spherical par-ticles can produce neither the same angle of repose nor stick-slip behavior without some explicit static friction force rule.While one can partially remedy this problem by apply-ing equation (19)to round particles [LH93],non-spherical granular particles offer an arguably more natural solution.Second,we find that the bookkeeping associated with the virtual-springs in (19)is computationally expensive.Third,the idea of approximating objects with round particles ex-tends naturally to general rigid bodies.3.2.Contact DetectionIn this section we briefly discuss an efficient method for identifying particle contacts.Since all our objects are com-pounds of spheres,contact detection reduces to finding over-lapping spheres.Clearly the naive O (n 2)method of considering all pairs of particles is impractical for even small simulations.Assuming all particles are of similar size,we can apply spatial hashing to more efficiently resolve contacts.A voxel size two times the maximum particle radius guarantees that only particles belonging to adjacent voxels can be in contact.Hence the problem of finding all particle contacts reduces to exploring the 27voxels surrounding a given particle.To exploit tempo-ral coherence,the hash structure is kept between time steps and updated whenever a particle moves between voxels.While this scheme provides O (n )lookup of all potential contacts,performing 27individual lookups per particle can dominate the overall computational cost.Further exploiting temporal coherence,we can reduce the cost to one lookup by inserting the particle into all 27neighboring voxels.This approach biases the cost of maintaining the hash structure to updates,rather than queries.Justifying this preference,we note that particles typically move short distances during acThe Eurographics Association 2005.Figure5:A steel ball collides with a sand pile producing a large splash.single adaptively-integrated time step.Moreover,it is clear that particles in persistent contact in a slowly-moving collec-tion are restricted by numerical stability to moving a fraction of a radius in a single time step.Hence it is very likely for a particle to remain within a voxel through many time steps. When using this method we observe a significant reduction in compute time.In simulations with large stagnant regions, the cost of contact detection is negligible.3.3.Interaction with Rigid BodiesIn many situations,granular materials need to interact with rigid bodies.For example,they are excavated by rigid tools and transported in rigid containers.In a landslide,they de-stroy large structures while at other times,we may simply throw small objects at them to see a splash.Realistic inter-action between granular materials and rigid bodies is thus an integral part of granular material simulation.In our method, we generate particles to completely cover the surface of a rigid object and simulate the interaction between granular particles and these surface particles using the contact forces described in Section2.With this method we achieve effi-cient two-way coupling between granular particles and rigid bodies.We create our granular particles from spheres,prescrib-ing a few positions to create the desired shape.In the case of an arbitrary shape requiring hundreds or thousands of parti-cles,such manual placement is not an option.We require a general method for representing shapes with particles.Given their ubiquity in graphics,we introduce a way to utilize bod-ies defined by triangle meshes.3.3.1.Surface Particle Generation for Triangle Meshes While one could place particles using the triangles di-rectly[Tur92],we prefer tofirst construct a signed distance representation and then apply ideas from Witkin and Heck-bert’s[WH94]method for sampling implicits with particles. This approach has the advantage that one can easily place particles completely inside,outside,or at an arbitrary offset from the original mesh.The signed distancefield is com-puted on a voxel grid of resolution comparable to the radius of the particles to be placed.Construction of the signed dis-tancefield follows the algorithm described in[BA02].Since only a narrow band of signed distances near the surface is required,this procedure is fast.The signed distance values are then stored in a hash table.An optional offset is then taken and an iso-surface ex-tracted using Marching Tetrahedra[TPG99].This offset tri-angle mesh is only used for initializing particle positions. The initial positions are chosen randomly within each trian-gle and the number of particles is controlled by the trian-gle’s area A,the particle radius R,and a user defined den-sity value D.Wefirst place DA/(πr2) particles inside the triangle and then one additional particle with probability DA/(πr2)− DA/(πr2) .Once initialized,the particles are then permitted tofloatFigure6:Rigid bodies defined by triangle meshes are cov-ered with particles.The signed distance representation per-mits placement of particles at an arbitrary offset from the mesh surface.Here a camel model is represented by4147 spheres.c The Eurographics Association2005.Figure 7:Structures swept away by anavalanche.Figure 8:A bulldozer clears the remains of the avalanche.along the signed distance field according to a simplified ver-sion of the interaction in [WH94].Each of these "floaters"has a position −→P and a velocity −→V =−→V F +−→V R .Floaters are constrained to the zero-isosurface using a feedback term,−→V F =−φ(−→P )−→N ,where φ(−→P )and −→N are the signed dis-tance and normal at −→P respectively.These values are ob-tained from the signed distance field through trilinear inter-polation.Since a uniform distribution of particles is desired,a simple repulsion model is also used.−→V R =∑P i ∈S \P K ( −→P −−→P i )−→P −−→P i −→P −−→P i .(20)Here K is a kernel function with compact support and S isthe set of floaters where K is non-zero.The floaters are inte-grated through time until the convergence criteria are satis-fied.The final positions of the floaters are then fixed on the rigid body during subsequent simulations with granular par-ticles.We find that this approach produces natural coverings for a wide range of models.Figure 6shows the result of this procedure applied to a camel model.3.4.Time IntegrationIntegration through time is accomplished through rigid body evolution [Bar97].Mass properties for granular particles are obtained from their constituent spheres.In the case of rigid bodies,the defining triangle mesh is used to derive the rel-evant quantities [Mir96].Since the time step required for numerical stability can vary greatly during simulation,an adaptive ODE integrator is a necessity.We choose to solve the equations of motion using an adaptive Runge-Kutta-Fehlburg method (RKF45)[Lam91].This choice would per-haps seem unnatural,given the discontinuities in the lower-order derivatives of motion when sphere pairs transition into and out of contacts.Nevertheless,test simulations on differ-ent flows found that RKF45was more efficient than lower-order methods at typically-specified error tolerances,per-haps because these discontinuities are at isolated points in both space and time.Before computing contact forces,the hash structure described in 3.2is updated to reflect the cur-rent particle positions.After forces have been computed therelevant derivatives are returned to the integrator.4.Results and DiscussionWe have applied our method to several large-scale simula-tions.Figure 5shows a steel ball colliding with a sand pile,represented by 45,494tetrahedron-shaped particles.Blocked by a glass wall,the wave of granular material is halted,and produces a heap near the base of the wall.Here the im-portance of accurate modeling of both static and dynamic phenomena is made clear as the granular material transi-tions rapidly from static to dynamic and then back to static regimes.The two-way coupling between granular materials and rigid bodies is highlighted in Figures 7and 8.In the first sequence an avalanche destroys structures composed of rigid bodies.At first the bodies resist the oncoming mass,produc-ing small splashes against the forward structures and divert-ing the flow.However the structures are soon overwhelmed,and become part of the flow themselves.After the avalanche has laid waste to this virtual environment,a bulldozer passes through to push aside the remnants.Rigid bodies and gran-ular material are swept away by the bulldozer’s blade.In all sequences the ground surface is covered with particles us-ing the same technique as the rigid bodies.In the splash and hourglass,an additional layer of smaller particles are added to increase the roughness of the surface.Since our treatment of rigid bodies also handles their mu-tual interactions,we produced a sequence with 1000rings to mimic the simulation in [GBF03].Two sequences were created to demonstrate varying quantities of friction and ing low friction and dissipation causes the rings to scatter,collecting near the edges of the containing box.When greater friction and dissipation are used,the rings clump closer to the rods at the center of the box.To this point we have distinguished between rigid bodies and gran-ular particles,however we note that in such large-scale sim-ulations we may regard individual rigid bodies as complex grains of granular material.Indeed,in many such cases the distinction is only a matter of scale and quantity.For ex-cThe Eurographics Association 2005.。