How Should the Discrete Element Method Be Applied in Industrial Systems ? : A Review †

In this paper, we describe an industrial application of the discrete element method (DEM). The DEM has been applied to various powder systems thus far and therefore appears to be an established approach. However, it cannot be applied to many industrial systems because of several critical problems such as modeling of large-scale simulations, complexly shaped wall boundaries and free surface fluid flow. To solve these problems, novel models were developed by our group. A coarse-grain DEM model was developed for large-scale simulations, a signed distance function-based wall boundary model was developed for complexly shaped walls and a DEM-moving particle semi-implicit method was developed for solid-liquid flow involving a free surface. The adequacy of these models was demonstrated through verification and validation tests. Our approach shows promise in industrial


Introduction
Granular flows usually differ from fluid flows.They are regarded as complex because constitutive equations are influenced by many parameters such as coefficient of restitution, coefficient of friction, particle size distribution, humidity, etc. Modeling of granular flow is extremely difficult.It may not be possible to model even simple phenomena in industries.Many requirements for numerical technologies should be applied in industrial design.The question arises as to how granular flow should be simulated in industrial systems.Granular flow could be simulated when individual particle behavior could be modeled.The discrete element method (DEM) (Cundall and Strack, 1979) was developed based on this idea and it became a successful model in the simulation of granular flow.The DEM is a Lagrangian approach and solid particle behavior is computed based on the Newton's second law of motion.The contact force acting on a solid particle is calculated using the Voigt model, namely, a spring-dashpot system.Recent improvements in computer hardware have allowed for an increase in number of solid particles and for com-plex phenomena such as solid-fluid coupling problems to be simulated.When solid-fluid coupling problems are computed, the DEM is coupled with computational fluid dynamics (CFD), where the local volume average technique (Anderson and Jackson, 1967) is used in the governing equations and semi-empirical momentum exchange coefficients (Wen and Yu, 1966;Ergun, 1952) are often used.This approach is referred to as the DEM-CFD method (Tsuji et al., 1993).Tsuji's contribution was significant to simulate solid-fluid coupling problems (Tsuji, 2007) and the DEM-CFD method became a standard global approach.Tsuji also investigated the application of a linear spring model, devised an efficient calculation setup (Tsuji et al., 1993) and proposed how physical properties should be set in a non-linear spring system (Tsuji et al., 1992).His research was exceptional and his proposal has become the standard approach.
Unfortunately, some critical problems exist for industrial application of the DEM.Although it has been applied to various engineering fields, the number of calculated particles is extremely small compared with the actual number in industrial processes.A maximum of one million calculated particles can be used on a single personal computer, whereas over one billion particles are required.Hence, the number of calculated particles becomes extremely small in industrial applications.In addition, modeling of solid-liquid flow involving a free surface is also a significant problems.Limited simulations of a wet ball mill and slurry transportation have been performed thus far because accurate simulation of free surface fluid flow is difficult, and this problem is common in a single flow.Recently, mesh-free Lagrangian CFD (Koshizuka et al., 1996;Monaghan, 1988) such as smoothed-particle hydrodynamics or the moving particle semi-implicit, volumeof-fluid (Yokoi, 2007;Xiao et al., 2011;Ohta et al., 2011) and front tracking methods (Vu et al., 2015) have been used extensively in simulations of free surface fluid flow.Indeed, solid-liquid flow could not be simulated in industrial systems.Perhaps coupling the DEM with these models for free surface fluid flow is a promising approach to simulate solid-liquid flow involving a free surface.Currently, models for solid-liquid flow that include a free surface are required.Finally, modeling of an arbitrary shaped wall is also important in industrial systems.In previous studies, a combination of simple surface equations (Mio et al., 2002;Liu et al., 2011;Sakai et al., 2005), the location of DEM particles and meshes (Shigeto and Sakai, 2011;Mateo-Ortiz et al., 2014) was used to create a wall boundary in the DEM.The combination of simple surface equations was difficult to create an arbitrarily shaped domain.In addition, the wall boundary model composed of the DEM particles exhibited difficulties in terms of accuracy for application to industrial systems.Although a mesh-based wall was often used to model a geometrically complex wall, collision detection may be difficult in programming the algorithm.
Hence, the DEM exhibits several problems when it was applied to industrial systems.To solve these problems, the authors have developed original models, namely, the coarse-grain DEM (Sakai and Koshizuka, 2009;Sakai et al., 2010;Sakai et al., 2012b;Sakai et al., 2014) for largescale powder systems, the DEM-MPS method (Sakai et al., 2012a;Yamada and Sakai, 2013;Sun et al., 2014) for solid-liquid flow involving a free surface, and the signed distance function (SDF)-based wall boundary model (Shigeto and Sakai, 2013;Sakai et al., 2015) for granular flow in a geometrically complex device.These models are described in this paper.

Large-scale discrete element simulation 2.1 Existing problems
Although the DEM seems to be well established, it is problematic in that the number of calculated particles is restricted when the simulation needs to be completed in a reasonable time on a single personal computer.This is especially critical because a number of DEM simulations are performed without using a supercomputer and the calculation time increases significantly on a single personal computer when an excessive number of particles are used.This problem cannot be solved even with the latest multicore processors (Liu et al., 2014;Yue et al., 2014).Many industries require application of the DEM to large-scale powder systems on a single personal computer.

Coarse-grain DEM model
In order to solve this problem, the authors have developed a physics model, which is termed the coarse-grain DEM model (Sakai and Koshizuka, 2009;Sakai et al., 2010;Sakai et al., 2012b;Sakai et al., 2014).The coarsegrain model was developed to model contact, drag, gravitational and cohesive forces as the scaling law model.As illustrated in Fig. 1, l 3 original particles exist in the coarse-grain particle whose size is l times larger than the original particle.Here, the coarse-grain ratio l was 2.0.In the coarse-grain model, the translational motion of the coarse-grain particle is assumed to be equivalent to be the average of that of the original particles.Hence, the velocity and displacement of the coarse-grain particle are assumed to be averages of the original particles, namely, , respectively.Besides, the original particles in the coarse-grain particle are assumed to rotate around each mass center with equal angular velocity, where the angular velocity is averaged as illustrated in Fig. 1(b).The contact force acting on the coarse-grain particle was estimated assuming that the total energy of the coarse-grain particle is consistent with that of the original particles.When a binary collision of the coarsegrain particles occurs, binary collisions between each of the original particles (i.e., l 3 binary collisions due to the original particles) are assumed to occur simultaneously.The contact force acting on the coarse-grain particle is simulated using springs, dash-pots and a friction slider, as for the existing DEM.From the assumption of the translational motion of the coarse-grain model, the drag force and gravitational force acting on a coarse-grain particle become l 3 times larger than that of the original particle.Consequently, the following relationship is obtained between the coarse-grain and original particles: where m, v, l, F f , V, p, F C and F g are the solid mass, solid velocity, coarse-grain ratio, drag force, particle volume, pressure, contact force and gravitational force, respectively.Subscripts CGM and O refer to the coarse-grain model and original particle, respectively.An important point of the coarse-grain model is that original physical properties such as density, coefficient of friction and coefficient of restitution can be applied directly.Now, modeling of the contact force is described in the coarse-grain model.The contact force can be divided into a normal and tangential component.Both components are functions of relative particle position and velocity.The normal component of the contact force is given by: where k, δ and η are the stiffness, the overlap and the damping coefficient, respectively.As described above, the translational motion of the coarse-grain particle is assumed to agree with that of the original particles.Hence, the CGM subscript can be replaced by O in Eq. ( 2).In the coarse-grain model, the normal component of the displacement is modeled by overlap as is the case with the existing DEM.Hence, the potential energy, i.e., elastic energy, is agreed between the coarse-grain particle and group of original particles by the modeling.
In the same manner as for the normal component, the tangential component of the contact force is given by: where µ is the coefficient of friction.When the solid particle does not slip on the surface of solid particle or wall, upper equation is applied.The translational motion of the coarse-grain particle is assumed to agree with that of the original particles.
From the assumption of the translational motion of the coarse-grain model, namely, agreement of kinetic energy between the coarse-grain and original particles, the drag force acting on a coarse-grain particle becomes l 3 times larger than that of the original particle.The drag force F f acting on the coarse-grain particle is given by: where ε, u and β are the void fraction, the fluid velocity and an interphase momentum transfer coefficient, respectively.
The following relationship is obtained for the rotational motion, namely, angular velocity and torque, in the original and coarse-grain particles: where ω, I and T are the angular velocity, the inertial momentum of the particle and the torque, respectively.Eq. ( 5) also ensures that the rotational energy between the coarsegrain particle and group of original particles agrees.The governing equations of the gas phase, namely, the continuity and the Navier-Stokes equations are the same as those of the original DEM-CFD method even for the coarse-grain model.
The solid particle motion, namely, the velocity, position, angular velocity and angle, could be updated using time integration schemes (Fraige and Langston, 2004;Danby et al., 2013;Kruggel-Emden et al., 2008).A semi-implicit finite volume method that uses a staggered grid was used to discretize the incompressible fluid flow.A hybrid scheme (Patankar, 1980), which is a combina-tion of a first-order upwind scheme and a second-order central scheme was used in the convection term.The pressure-velocity coupling was based on the fractional step method for incompressible fluid flows.

Numerical example of coarse-grain model
Adequacy of the coarse-grain model with high coarsegrain ratio is shown through the validation tests (Sakai et al., 2014).The validation results were described in this paper.
The calculation domain was 50 mm wide, 20 mm long and 200 mm high.The domain was discretized using 2.5 mm structural grids in the x, y and z directions.The grid size was chosen to be large relative to the solid particle diameter, where the fluid flow characteristics could be simulated accurately.The number of grids was 20, 80 and 8, in the x, y and z directions, respectively.The number of calculated particles was 230,400, which corresponded to the actual mass of powder used in the experiment.Initially, the spherical particles were packed randomly.Air was injected from the base, where the superficial velocity was 0.070 m/s.The density, stiffness, coefficient of restitution and coefficient of friction were 2,500 kg/m 3 , 10 N/m, 0.9 and 0.3, respectively.The gas was air, where the density and the viscosity were 1.0 kg/m 3 and 1.8 × 10 -5 Pa s, respectively.
Fig. 2 illustrates typical snapshots obtained from simulation and experimental results of the fluidized state at quasi-steady state.In the simulation and experiment, splashes and bubbles were observed in the quasi-steady state.From the appearance, macroscopic flow properties agreed qualitatively between the simulation and experiment.
The pressure drop and bed height were measured in the experiment and simulation for quantitative comparison.Fig. 3 compares the pressure drop between the simulation and experiment, which agreed quantitatively.The bed height was measured by visual observation to compare the calculated and experimental results fairly.The simulated bed height had an error of 7.8 % against the experimental height, which is an allowable margin of error from an engineering perspective.
The coarse-grain model can therefore simulate the bulk bed state by setting a higher coarse-grain ratio.This implies that the particle-particle and particle-fluid interactions were modeled precisely in the coarse-grain model.We can easily predict that the coarse-grain model may simulate other phenomena related to bulk states, such as bubbling, channeling and slugging.The coarse-grain model has a sufficient potential to be used with much high coarse grain ratio.Hereafter, the possibility will be investigated by various validation tests.

Existing problems
Solid-liquid flows involving a free surface occur in chemical engineering applications such as wet ball milling and slurry transport.Whereas numerical simulation of the process should be used for the design and optimization of operational conditions, such simulations have not yet been implemented extensively.This situation has arisen because of difficulties in modeling and computationally expensive calculations.
Currently, there is no broadly used numerical approach for modeling such solid-liquid flows, although some simulations have been reported previously.These approaches exhibit some problems, namely, some difficulties in mod- eling large deformations of the fluid free surface, difficulties in modeling hydrodynamic forces and difficulties in modeling contact force.Therefore, the existing methodology cannot be used for many solid-liquid flow systems in chemical engineering applications.The development of a new method is desired to solve the above difficulties or problems.Many industrial solid-liquid flows involve large deformations and/or a high volume fraction of solid phase, where the contact force acting on a solid particle becomes important.Lagrangian method such as the SPH or MPS method makes it possible to simulate free surface fluid flow precisely.These methods have been applied to various engineering fields thus far.Lagrangian-Lagrangian approaches are applicable to these important industrial flows.Hence, the MPS method and the DEM are combined to simulate solid-liquid flow involving free surfaces (Sakai et al., 2012a;Yamada and Sakai, 2013;Sun et al., 2014).The DEM-MPS method whose algorithm was fully explicit is described in this paper.

DEM-MPS method
For a locally averaged description of continuous fluid motion, the mass and momentum conservation equations can be written as: where ρ, p, ν, f and g are fluid density, pressure, kinematic viscosity, drag force and gravitational acceleration, respectively.
The governing equation for the solid particles is given by Newton's second law of motion: In a Lagrangian particle method such as the MPS method, the continuum phase is modeled as an assembly of fluid particles.In the MPS method, fluid particles are not connected and therefore so-called weight functions are used to define discretization schemes on the particle positions.
Weakly compressible and incompressible methods were developed in the MPS method.In the weakly compressible approach, fluid pressure is calculated by the modified Tait equation of state, which agrees well with the real behavior of water under typical circumstances.
where c is the speed of sound in a fluid at reference density and γ is a constant often set to 7.0 for water.For numerical purposes within a weakly compressible method, the speed of sound c can be tuned artificially to relax the restriction on time increments because of stability conditions.The c should be 10 times higher than the maximum main stream velocity to constrain density fluctuations up to 1 % and guarantee a nearly divergence-free velocity field.
The pressure gradient and viscous terms are discretized by the original MPS method.

Numerical example of DEM-MPS method
In order to illustrate the adequacy of the DEM-MPS method, a validation test was performed in our previous study (Sun et al., 2014).A cylindrical tank was used with inner diameter of 120 mm and depth of 100 mm.Water and glass beads were placed in the tank.The average diameter of the glass beads was ~2.7 mm, and their density was 2500 kg/m 3 .In the laboratory experiment, the tank was revolved at 63 rpm and 1.25 kg of glass beads were inserted.Water was injected afterwards and adjusted to fill approximately half of the tank.
2D and 3D simulations were performed in this study.The diameter of the DEM particles was 2.7 mm and their density was 2,500 kg/m 3 .The stiffness was set to be 10,000 N/m, coefficient of restitution was 0.9 and coefficient of friction was 0.3.For the liquid phase, the physical properties are equal to those of water, namely, the density was 1000 kg/m 3 and the kinematic viscosity was 0.89 mm 2 /s.The fluid particle size, which was defined as the initial interval between two fluid particles, was chosen to be 2.7 mm.
Fig. 5(a) and (b) shows a snapshot of the 2D and 3D calculation results at quasi steady state.In this figure, liquid phase was visualized by using meshes.These observations have been confirmed by our experiment.The bed width and height were measured in the validation tests as shown in Fig. 5(c), where 2D and simulation results are given by blue and red line respectively.The bed width and height obtained from the 2D and 3D simulation results agreed quantitatively with the experimental data, though bed shape obtained from the 3D simulation was more precise than the 2D simulation.Hence, adequacy of the DEM-MPS method was demonstrated through the validation tests.

Arbitrary shape boundary wall model 4.1 Existing problems
In the previous DEM simulations, the calculation domain, i.e., the wall boundary, was modeled using surface equations (Sakai et al., 2005) or meshes (Shigeto and Sakai, 2011;Mateo-Ortiz et al., 2014).The use of surface equations is an easy approach to create simple regimes.Modeling of geometrically complex walls is often required in industrial systems.These domain shapes are substantially more difficult to express with combinations of the equations.Hence, arbitrary domain shapes were modeled using meshes.The meshes are more adaptable for arbitrary shapes than the equations, but there is a difficulty in introducing the algorithm of collision detection between particles and walls.To solve these problems, a simple boundary model was proposed, where the wall boundary is expressed by a scalar field, which is referred to as the SDF.The SDF model was based on the level set method (Osher and Fedkiw, 2002), which was originally proposed in the CFD.The uniqueness of the model is that collisions can be detected without complicated procedures.We developed the SDF model, which includes conservation of energy in non-dissipative systems.The DEM/ SDF can also be applied to cohesive particle and rotating vessel systems (Shigeto and Sakai, 2013).

DEM/SDF model
As mentioned above, we developed the SDF model to retain energy conservation in a non-dissipative system.This idea is important to simulate cohesive and contact forces precisely in the DEM simulations.The contact force model based on the DEM/SDF is addressed now.The contact force is modeled based on elastic potential energy when the solid particles interact with walls.The energy is related to displacement by: where ϕ and r are the SDF and minimal distance between the wall and particle, respectively.The gradient of the energy P gives the elastic force between the particle and the SDF based wall boundary: Hence, the normal component of the contact force is given by: As a matter of course, our SDF model is applicable to cohesive particle systems that are affected by van der Waals force.In the same manner as for the elastic force, the van der Waals force is modeled based on the potential energy, namely: where h SDF is the intersurface distance.
The gradient of Eq. ( 13) yields the van der Waals force between the particle and the SDF based wall boundary: The conservation of energy is guaranteed theoretically in our SDF model.Hence, this scheme may stabilize the calculation.In addition, the algorithm is simpler the existing wall boundary composed of meshes.We expected that the SDF will become a standard wall boundary model the DEM simulations.Hereafter the DEM/SDF is going to be applied to investigate lots of phenomena in a geometrically complex boundary.

Screw conveying
In order to illustrate the adequacy of the DEM/SDF method, a verification test was performed in a screw conveyor (Shigeto and Sakai, 2013).In this test, simulations were carried out with the mesh-and new SDF-based models.Initially, 62,500 particles was located randomly at the bottom of the calculation domain.The particles have a uniform diameter of 200 μm and density of 2,500 kg/m 3 .Fig. 6 designates the scalar value distribution of the SDF model of the inner screw blade.The sign becomes posi-tive inside the calculation domain and vice versa.The screw was rotated at 30 rpm.In the verification tests, cohesive particles were used, where the Hamaker constant was 5.0 × 10 -20 J.
Fig. 7 shows typical snapshots at 0.5, 1.0 and 1.5 s in the verification tests.Similar particle distributions, agreement of velocity distribution of the solid particles and agreement of transported mass were obtained in the SDF and mesh systems.Besides, the load torque on the screw agreed well between the SDF and mesh models.Consequently, these results illustrate that the accuracy of the SDF model is equivalent to the existing mesh model in the reproduction of granular flows in the screw conveying system.

Twin-screw kneader
Application of the DEM/SDF method to a twin-screw kneader is described according to past study (Sakai et al., 2015).In order to show the adequacy of the DEM/SDF method, a validation test was performed in a twin-screw kneader where the focus was on dry powder mixing.Simulation and experimental results were compared for the apparent solid particle location and particle bed height.For the validation system, the height, width and length were 106 mm, 186 mm and 46 mm, respectively.Four paddles were used in the depth direction.The paddles were rotated clockwise, and the paddle rotational velocity was 60 rpm in the validation tests.Glass beads (~1.0 mm) were used in the tests.The total amount of powder was 0.524 kg.Calculation conditions were set to compare the simulation and experimental results fairly.Mono-sized solid particles with 1.0 mm diameter were used in the calculations.A simple linear contact model was used in this study.The spring constant, coefficient of friction and coefficient of restitution were set at 1,000 N/m, 0.3 and 0.9, respectively.The number of solid particles in the calculation was 400,000 and their mass corresponded to that in the experiment.The time step was 1.0 × 10 -5 s and the number of iterations was 1.6 × 10 6 .The wall boundary was modeled by the SDF.Fig. 8 shows the SDF for four cross-sectional views of the paddle.Once again, in the figure, positive and negative regions exist inside and outside the paddle, respectively.Distance and sign were given for the entire computational domain, i.e., the two paddles and the vessel.Fig. 9 shows the validation test results.The simulation and experiment agreed qualitatively in terms of the spatial distribution of particles, namely, more particles were located on the left side.The quantitative bed height was compared between the simulations and the experiments.The bed heights were estimated to be 101.81 and 101.24 mm, respectively, in the simulation and experiment.Thus, differences between the simulation and experiment were quite small.The validation test result indicates that the DEM/SDF simulates the particle location precisely.

Conclusions
New modeling of the flow simulations were described for industrial applications.At first, problems of the existing DEM was designated from a view point of the industrial application such as modeling of large-scale simulations, complexly shaped wall boundaries and free surface fluid flow.In order to solve these problems, novel models were developed by our group.A coarse-grain DEM model was developed for large-scale simulations, a SDF-based wall boundary model was developed for complexly shaped walls and a DEM-MPS method was developed for solid-liquid flow involving a free surface.The adequacy of these models was demonstrated through verification and validation tests.Combination of these models becomes a key technology in applying the DEM to actual powder process.