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

Abstract

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 applications.

1. 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 complex 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.

At present, the DEM appears to be an established approach. It has been applied in various engineering fields and many DEM applications exist. With respect to single flow such as granular flow, the DEM has been applied to solid mixing (Radeke et al., 2010; Chandratilleke et al., 2012; Halidan et al., 2014; Eitzlmayr et al., 2014; Huang and Kuo, 2014), die filling (Guo et al., 2009; Wu, 2008; Mateo-Ortiz et al., 2014; Yue et al., 2014), powder transportation (Tsuji et al., 1992; Mio et al., 2009; Kuang et al., 2013), and granulation (Washino et al., 2013), among others. Since the DEM-CFD method was developed, the DEM has been applied to new fields along with fluidized beds and cyclones. The DEM-CFD method was applied to simulations of coarse particles that belong to the Geldart B and D group (Kawaguchi et al., 1998; Wang et al., 2013; Hernández-Jiménez et al., 2014; Olaofe et al., 2014; Patil et al., 2014) and fine particles such as Geldart A (Pandit et al., 2006; Kobayashi et al., 2013; Hong et al., 2015). The DEM-CFD method has been validated in fluidized (Alobaid et al., 2013; Wu et al., 2013) or spouted (Jajcevic et al., 2013; Salikov et al., 2014) beds. The DEM-CFD method can be used in simulations of geometrically complex fluidized beds such as circulating fluidized beds (Xu et al., 2012; Lu et al., 2014). Hence, DEM technologies are regarded as a reliable 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, volume-of-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 large-scale 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.

2. 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 multi-core 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.

2.2 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 coarse-grain model was developed to model contact, drag, gravitational and cohesive forces as the scaling law model. As illustrated in Fig. 1, l³ 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, v̄_O = v_CGM and δ̄_O = δ_CGM, 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 coarse-grain particles occurs, binary collisions between each of the original particles (i.e., l³ 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³ times larger than that of the original particle. Consequently, the following relationship is obtained between the coarse-grain and original particles:   

$$\begin{array}{ll}{m}_{\text{CGM}}{\mathit{v}}_{\text{CGM}}\hfill & ={\mathit{F}}_{{\text{f}}_{\text{CGM}}}-V_{\text{CGM}}\nabla p+\sum {\mathit{F}}_{{\text{C}}_{\text{CGM}}}+{\mathit{F}}_{{\text{g}}_{\text{CGM}}}\hfill \\ \hfill & =l^3{\overline{\mathit{F}}}_{{\text{f}}_{\text{O}}}-l^3{V}_{\text{O}}\nabla p+l^3\sum {\overline{\mathit{F}}}_{{\text{C}}_{\text{O}}}+l^3{\overline{\mathit{F}}}_{{\text{g}}_{\text{O}}}\hfill \end{array}$$(1)

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.

Fig. 1

Image of coarse-grain model (coarse-grain ratio l = 2.0). Reprinted with permission from Ref. (Sakai et al., Chem. Eng. J., 2014). Copyright: (2014) Elsevier B.V.

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:   

$$\begin{array}{ll}{\mathit{F}}_{{\text{C}}_{{\text{n}}_{\text{CGM}}}}\hfill & =l^3\left(-k{\mathit{\delta }}_{{\text{n}}_{\text{CGM}}}-\eta {\mathit{v}}_{{\text{n}}_{\text{CGM}}}\right)\hfill \\ \hfill & =l^3\left(-k{\overline{\mathit{\delta }}}_{{\text{n}}_{\text{O}}}-\eta {\overline{\mathit{v}}}_{{\text{n}}_{\text{O}}}\right)\hfill \end{array}$$(2)

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:   

$$\begin{array}{ll}{\mathit{F}}_{{\text{C}}_{{\text{t}}_{\text{CGM}}}}\hfill & =\{\begin{array}{l}{l}^3\left(-k{\mathit{\delta }}_{{\text{t}}_{\text{CGM}}}-\eta {\mathit{v}}_{{\text{t}}_{\text{CGM}}}\right)\\ -l^3\mu \left|{\mathit{F}}_{{\text{C}}_{{\text{n}}_{\text{CGM}}}}\right|\left|{\mathit{v}}_{{\text{t}}_{\text{CGM}}}\right|/\left|{\mathit{v}}_{{\text{t}}_{\text{CGM}}}\right|\end{array}\hfill \\ \hfill & =\{\begin{array}{l}{l}^3\left(-k{\overline{\mathit{\delta }}}_{{\text{t}}_{\text{O}}}-\eta {\overline{\mathit{v}}}_{{\text{t}}_{\text{O}}}\right)\\ -l^3\mu \left|{\overline{\mathit{F}}}_{{\text{C}}_{{\text{n}}_{\text{O}}}}\right|\left|{\overline{\mathit{v}}}_{{\text{t}}_{\text{O}}}\right|/\left|{\overline{\mathit{v}}}_{{\text{t}}_{\text{O}}}\right|\end{array}\hfill \end{array}$$(3)

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³ times larger than that of the original particle. The drag force F_f acting on the coarse-grain particle is given by:   

$$\begin{array}{ll}{\mathit{F}}_{{\text{f}}_{\text{CGM}}}\hfill & =\frac{\beta }{1-\epsilon }(\mathit{u}-{\mathit{v}}_{\text{CGM}})V_{\text{CGM}}\hfill \\ \hfill & =l^3\frac{\beta }{1-\epsilon }(\mathit{u}-{\overline{\mathit{v}}}_{\text{O}})V_{\text{O}}\hfill \end{array}$$(4)

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:   

$$\begin{array}{ll}{\dot{\mathit{\omega }}}_{\text{CGM}}\hfill & =\frac{{\mathit{T}}_{\text{CGM}}}{I_{\text{CGM}}}=\frac{r_{\text{CGM}}\times {\mathit{F}}_{{\text{C}}_{{\text{CGM}}_{\text{t}}}}}{I_{\text{CGM}}}\hfill \\ \hfill & =\frac{l{r}_{\text{O}}\times l^3{\overline{\mathit{F}}}_{{\text{C}}_{{\text{O}}_{\text{t}}}}}{l^5{I}_{\text{O}}}=\frac{l^4{\overline{\mathit{T}}}_{\text{O}}}{l^5{I}_{\text{O}}}=\frac{1}{l}{\overline{\dot{\mathit{\omega }}}}_{\text{O}}\hfill \end{array}$$(5)

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 coarse-grain 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 combination 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.

2.3 Numerical example of coarse-grain model

Adequacy of the coarse-grain model with high coarse-grain 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³, 10 N/m, 0.9 and 0.3, respectively. The gas was air, where the density and the viscosity were 1.0 kg/m³ 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.

Fig. 2

Particle spatial distribution at quasi-steady state in the validation test.

Reprinted with permission from Ref. (Sakai et al., Chem. Eng. J., 2014). Copyright: (2014) Elsevier B.V.

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.

Fig. 3

Comparison of pressure drop between calculated and experimental results in a fluidized bed.

Reprinted with permission from Ref. (Sakai et al., Chem. Eng. J., 2014). Copyright: (2014) Elsevier B.V.

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.

3. Modeling of solid-liquid flow involving a free surface

3.1 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 modeling 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.

3.2 DEM-MPS method

For a locally averaged description of continuous fluid motion, the mass and momentum conservation equations can be written as:   

$$\frac{\text{D}(\epsilon \rho )}{\text{D}t}+\epsilon \rho \nabla \cdot \mathit{u}=0$$(6)

  

$$\frac{\text{D}\mathit{u}}{\text{Dt}}=-\frac{1}{\rho }\nabla p+v{\nabla }^2\mathit{u}-\frac{1}{\epsilon \rho }\mathit{f}+\mathit{g}$$(7)

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:   

$$m\frac{\text{d}\mathit{v}}{\text{dt}}={\mathit{F}}_{\text{C}}+{\mathit{F}}_{\text{f}}-V\nabla p+m\mathit{g}$$(8)

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.   

$$p=\{\begin{array}{cc}\frac{{\rho }_0{c}^2}{\gamma }\left({\left(\frac{N}{n_0}\right)}^{\gamma }-1\right)& n\ge \epsilon n_0\\ 0& n\ge \epsilon n_0\end{array}$$(9)

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.

Fig. 4

Initial configuration of fluid particles. Reprinted with permission from Ref. (Sun et al., Chem. Eng. J., 2014). Copyright: (2014) Elsevier B.V.

The pressure gradient and viscous terms are discretized by the original MPS method.

3.3 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³. 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³. 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³ and the kinematic viscosity was 0.89 mm²/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.

Fig. 5

Quasi-steady solution of circulating flow by the DEM-MPS method.

Reprinted with permission from Ref. (Sun et al., Chem. Eng. J., 2014). Copyright: (2014) Elsevier B.V.

4. 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).

4.2 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:   

$$P_{\text{elastic}}^{\text{SDF}}=\frac{1}{2}k{\left|{\mathit{\delta }}_{\text{n}}^{\text{SDF}}\right|}^2=\frac{1}{2}k{\left(r-\varphi \right)}^2$$(10)

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:   

$$\begin{array}{ll}{\mathit{F}}_{\text{elastic}}^{\text{SDF}}\hfill & =\nabla P_{\text{elastic}}^{\text{SDF}}\hfill \\ \hfill & =-k(r-\varphi )\cdot \nabla \varphi \hfill \\ \hfill & =-k{\mathit{\delta }}^{\text{SDF}}\cdot \left|\nabla \varphi \right|\hfill \end{array}$$(11)

Hence, the normal component of the contact force is given by:   

$${\mathit{F}}_{{\text{C}}_{\text{n}}}^{\text{SDF}}=-k{\mathit{\delta }}_{\text{n}}^{\text{SDF}}\left|\nabla \varphi \right|-\eta {\mathit{v}}_{\text{n}}$$(12)

which conserves the energy in non-dissipative systems because of the correction term |∇ϕ|.

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:   

$$P_{\text{vdw}}^{\text{SDF}}=-\frac{H_{\text{A}}\cdot r}{6h^{\text{SDF}}}$$(13)

where h^SDF is the intersurface distance.   

$$h^{\text{SDF}}=\varphi -r$$(14)

The gradient of Eq. (13) yields the van der Waals force between the particle and the SDF based wall boundary:   

$${\mathit{F}}_{\text{vdw}}^{\text{SDF}}=\nabla P_{\text{vdw}}^{\text{SDF}}=-\frac{H_{\text{A}}\cdot r}{6{\left(h^{\text{SDF}}\right)}^2}\cdot \nabla \varphi .$$(15)

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 in the DEM simulations. Hereafter the DEM/SDF is going to be applied to investigate lots of phenomena in a geometrically complex boundary.

4.3 Numerical examples of DEM/SDF model

4.3.1 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³. Fig. 6 designates the scalar value distribution of the SDF model of the inner screw blade. The sign becomes positive 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⁻²⁰ J.

Fig. 6

Cross-sectional view of SDF in a screw conveyer.

Reprinted with permission from Ref. (Shigeto and Sakai, Chem. Eng. J., 2013). Copyright: (2013) Elsevier B.V.

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.

Fig. 7

Typical snapshots of granular flow in a screw conveying system.

Reprinted with permission from Ref. (Shigeto and Sakai, Chem. Eng. J., 2013). Copyright: (2013) Elsevier B.V.

4.3.2 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⁶. 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.

Fig. 8

Cross-sectional view of SDF in a twin-screw kneader in a twin-screw kneader.

Reprinted with permission from Ref. (Sakai et al., Chem. Eng. J., 2015). Copyright: (2015) Elsevier B.V.

Fig. 9

Validation test result in a twin-screw kneader.

Reprinted with permission from Ref. (Sakai et al., Chem. Eng. J., 2015). Copyright: (2015) Elsevier B.V.

5. Conclusions

New modeling of the granular 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.

Acknowledgments

My heartfelt gratitude and appreciation goes to Prof Seiichi Koshizuka whose enthusiastic tutelage was of inestimable value for my life. I’m sincerely grateful to Prof. Jusuke Hidaka and to Prof. Yutaka Tsuji for their encouragement and instructions. I would like to thank Dr. Yoshinori Yamada, Dr. Yusuke Shigeto and Dr. Xiaosong Sun for their insightful discussion on developing the numerical models.

Nomenclature

ω

angular velocity (rad/s)

CGM

coarse-grain model

l

coarse-grain ratio (−)

μ

coefficient of friction (−)

F_C

contact force (N)

η

damping coefficient (−)

δ

displacement (m)

F_f

drag force (N)

f

drag force (N/m³)

ρ

fluid density (kg/m³)

u

fluid velocity (m/s)

g

gravitational acceleration (m/s²)

F_g

gravity (N)

H_A

Hamaker constant (J)

I

inertia momentum of particle (kg m²)

β

inter-phase momentum transfer coefficient (kg/m²s)

h

intersurface distance (m)

ν

kinematic viscosity (m²/s)

O

original

m

particle mass (kg)

PW

particle-to-wall

v

particle velocity (m/s)

V

particle volume (m³)

p

pressure (N/m²)

SDF

signed distance function (−)

ϕ

signed distance function (−)

c

speed of sound (m/s)

k

stiffness (N/m)

T

torque (N m)

ε

void fraction (−)

Author’s short biography

Mikio Sakai

Dr Mikio Sakai is an Associate Professor at The University of Tokyo. He earned his Ph.D. degree from The University of Tokyo in 2006 and became an Assistant Professor in 2007. He then became an Associate Professor in 2008. His interests include modeling of multi-phase flow such as gas-solid flows, solid-liquid flows and gas-liquid-solid flows, the arbitrary shaped wall boundary model for granular flows and parallel computing techniques. He is now an associate editor of Chemical Engineering Science.

References

•  Alobaid  F.,  Ströhle  J.,  Epple  B., Extended CFD/DEM model for the simulation of circulating fluidized bed, Adv. Powder Technol., 24 (2013) 403–415.
•  Anderson  T.B.,  Jackson  R., Fluid Mechanical Description of Fluidized Beds. Equations of Motion, Fluid Mech. Descr. Fluid. Beds. Equations Motion, 6 (1967) 527–539.
•  Chandratilleke  G.R.,  Yu  a.B.,  Bridgwater  J., A DEM study of the mixing of particles induced by a flat blade, Chem. Eng. Sci., 79 (2012) 54–74.
•  Cundall  P.A.,  Strack  O.D.L., A discrete numerical model for granular assemblies, Geotechnique, 29 (1979) 47–65.
•  Danby  M.,  Shrimpton  J.,  Palmer  M., On the optimal numerical time integration for DEM using Hertzian force models, Comput. Chem. Eng., 58 (2013) 211–222.
•  Eitzlmayr  A.,  Koscher  G.,  Reynolds  G.,  Huang  Z.,  Booth  J.,  Shering  P.,  Khinast  J., Mechanistic Modeling of Modular Co-Rotating Twin-Screw Extruders, Int. J. Pharm., 474 (2014) 157–176.
•  Ergun  S., Fluid flow through packed columns, Chem. Eng. Progr., 48 (1952) 89–94.
•  Fraige  F.Y.,  Langston  P.a, Integration schemes and damping algorithms in distinct element models, Adv. Powder Technol., 15 (2004) 227–245.
•  Guo  Y.,  Kafui  K.D.,  Wu  C.Y.,  Thornton  C.,  Seville  J.P.K., A coupled DEM/CFD analysis of the effect of air on powder flow during die filling, AIChE J., 55 (2009) 49–62.
•  Halidan  M.,  Chandratilleke  G.R.,  Chan  S.L.I.,  Yu  A.B.,  Bridgwater  J., Prediction of the mixing behaviour of binary mixtures of particles in a bladed mixer, Chem. Eng. Sci., 120 (2014) 37–48.
•  Hernández-Jiménez  F.,  Li  T.,  Cano-Pleite  E.,  Rogers  W.,  Acosta-Iborra  A., Characterization of the particle-wall frictional forces in pseudo-2D fluidized beds using DEM, Chem. Eng. Sci., 116 (2014) 136–143.
•  Hong  K.,  Chen  S.,  Wang  W.,  Li  J., Fine-grid two-fluid modeling of fluidization of Geldart A particles, Powder Technol., 2015.
•  Huang  A.N.,  Kuo  H.P., Developments in the tools for the investigation of mixing in particulate systems—A review, Adv. Powder Technol., 25 (2014) 163–173.
•  Jajcevic  D.,  Siegmann  E.,  Radeke  C.,  Khinast  J.G., Large-scale CFD-DEM simulations of fluidized granular systems, Chem. Eng. Sci., 98 (2013) 298–310.
•  Kawaguchi  T.,  Tanaka  T.,  Tsuji  Y., Numerical simulation of two-dimensional fluidized beds using the discrete element method (comparison between the two- and three-dimensional models), Powder Technol., 96 (1998) 129–138.
•  Kobayashi  T.,  Tanaka  T.,  Shimada  N.,  Kawaguchi  T., DEM-CFD analysis of fluidization behavior of Geldart Group A particles using a dynamic adhesion force model, Powder Technol., 248 (2013) 143–152.
•  Koshizuka  S.,  Nobe  A.,  Oka  Y., Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Eng., 123 (1996) 421–434.
•  Kruggel-Emden  H.,  Sturm  M.,  Wirtz  S.,  Scherer  V., Selection of an appropriate time integration scheme for the discrete element method (DEM), Comput. Chem. Eng., 32 (2008) 2263–2279.
•  Kuang  S.B.,  Li  K.,  Zou  R.P.,  Pan  R.H.,  Yu  A.B., Application of periodic boundary conditions to CFD-DEM simulation of gas-solid flow in pneumatic conveying, Chem. Eng. Sci., 93 (2013) 214–228.
•  Liu  H.,  Tafti  D.K.,  Li  T., Hybrid parallelism in MFIX CFD-DEM using OpenMP, Powder Technol., 259 (2014) 22–29.
•  Liu  P.Y.,  Yang  R.Y.,  Yu  A.B., Dynamics of wet particles in rotating drums: Effect of liquid surface tension, Phys. Fluids, 23 (2011).
•  Lu  L.,  Xu  J.,  Ge  W.,  Yue  Y.,  Liu  X.,  Li  J., EMMS-based discrete particle method (EMMS-DPM) for simulation of gas-solid flows, Chem. Eng. Sci., 120 (2014) 67–87.
•  Mateo-Ortiz  D.,  Muzzio  F.J.,  Méndez  R., Particle size segregation promoted by powder flow in confined space: The die filling process case, Powder Technol., 262 (2014) 215–222.
•  Mio  H.,  Higuchi  R.,  Ishimaru  W.,  Shimosaka  A.,  Shirakawa  Y.,  Hidaka  J., Effect of paddle rotational speed on particle mixing behavior in electrophotographic system by using parallel discrete element method, Adv. Powder Technol., 20 (2009) 406–415.
•  Mio  H.,  Kano  J.,  Saito  F.,  Kaneko  K., Effects of rotational direction and rotation-to-revolution speed ratio in planetary ball milling, Mater. Sci. Eng. A, 332 (2002) 75–80.
•  Monaghan  J.J., An introduction to SPH, Comput. Phys. Comm., 48 (1988) 89–96.
•  Ohta  M.,  Kikuchi  D.,  Yoshida  Y.,  Sussman  M., Robust numerical analysis of the dynamic bubble formation process in a viscous liquid, Int. J. Multiph. Flow, 37 (2011) 1059–1071.
•  Olaofe  O.O.,  Patil  A.V.,  Deen  N.G.,  van der Hoef  M.A.,  Kuipers  J.A.M., Simulation of particle mixing and segregation in bidisperse gas fluidized beds, Chem. Eng. Sci., 108 (2014) 258–269.
•  Osher  S.,  Fedkiw  R., Level Set Methods and Dynamic Implicit Surfaces, Springer, 2002.
•  Pandit  J.K.,  Wang  X.S.,  Rhodes  M.J., On Geldart Group A behaviour in fluidized beds with and without cohesive interparticle forces: A DEM study, Powder Technol., 164 (2006) 130–138.
•  Patankar  S.V., Numerical heat transfer and fluid flow, Routledge, 1980.
•  Patil  A.V.,  Peters  E.aJ.F.,  Kolkman  T.,  Kuipers  J.aM., Modeling bubble heat transfer in gas-solid fluidized beds using DEM, Chem. Eng. Sci., 105 (2014) 121–131.
•  Radeke  C.A.,  Glasser  B.J.,  Khinast  J.G., Large-scale powder mixer simulations using massively parallel GPU architectures, Chem. Eng. Sci., 65 (2010) 6435–6442.
•  Sakai  M.,  Abe  M.,  Shigeto  Y.,  Mizutani  S.,  Takahashi  H.,  Viré  A.,  Percival  J.R.,  Xiang  J.,  Pain  C.C., Verification and validation of a coarse grain model of the DEM in a bubbling fluidized bed, Chem. Eng. J., 244 (2014) 33–43.
•  Sakai  M.,  Koshizuka  S., Large-scale discrete element modeling in pneumatic conveying, Chem. Eng. Sci., 64 (2009) 533–539.
•  Sakai  M.,  Shibata  K.,  Koshizuka  S., Effect of Nuclear Fuel Particle Movement on Nuclear Criticality in a Rotating Cylindrical Vessel, J. Nucl. Sci. Technol., 42 (2005) 267–274.
•  Sakai  M.,  Shigeto  Y.,  Basinskas  G.,  Hosokawa  A.,  Fuji  M., Discrete element simulation for the evaluation of solid mixing in an industrial blender, Chem. Eng. J., 279 (2015) 821–839.
•  Sakai  M.,  Shigeto  Y.,  Sun  X.,  Aoki  T.,  Saito  T.,  Xiong  J.,  Koshizuka  S., Lagrangian-Lagrangian modeling for a solid-liquid flow in a cylindrical tank, Chem. Eng. J., 200–202 (2012a) 663–672.
•  Sakai  M.,  Takahashi  H.,  Pain  C.C.,  Latham  J.-P.,  Xiang  J., Study on a large-scale discrete element model for fine particles in a fluidized bed, Adv. Powder Technol., 23 (2012b) 673–681.
•  Sakai  M.,  Yamada  Y.,  Shigeto  Y.,  Shibata  K.,  Kawasaki  V.M.,  Koshizuka  S., Large-scale discrete element modeling in a fluidized bed, Int. J. Numer. Meth. Fluids, 64 (2010) 1319–1335.
•  Salikov  V.,  Antonyuk  S.,  Heinrich  S.,  Sutkar  V.S.,  Deen  N.G.,  Kuipers  J.aM., Characterization and CFD-DEM modelling of a prismatic spouted bed, Powder Technol., 270 (2014) 622–636.
•  Shigeto  Y.,  Sakai  M., Arbitrary-shaped wall boundary modeling based on signed distance functions for granular flow simulations, Chem. Eng. J., 231 (2013) 464–476.
•  Shigeto  Y.,  Sakai  M., Parallel computing of discrete element method on multi-core processors, Particuology, 9 (2011) 398–405.
•  Sun  X.,  Sakai  M.,  Sakai  M.-T.,  Yamada  Y., A Lagrangian-Lagrangian coupled method for three-dimensional solid-liquid flows involving free surfaces in a rotating cylindrical tank, Chem. Eng. J., 246 (2014) 122–141.
•  Tsuji  Y., Multi-scale modeling of dense phase gas-particle flow, Chem. Eng. Sci., 62 (2007) 3410–3418.
•  Tsuji  Y.,  Kawaguchi  T.,  Tanaka  T., Discrete particle simulation of two dimensional fluidized bed, Powder Technol., 77 (1993) 79–87.
•  Tsuji  Y.,  Tanaka  T.,  Ishida  T., Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol., 71 (1992) 239–250.
•  Vu  T.V.,  Tryggason  G.,  Homma  S.,  Wells  J.,  Takakura  H., Front Tracking Computation of Trijunction Solidification with Volume Change, Procedia IUTAM, 15 (2015) 14–17.
•  Wang  J.,  van der Hoef  M.a.,  Kuipers  J.aM., Particle granular temperature of Geldart A, A/B and B particles in dense gas-fluidized beds, Chem. Eng. Sci., 97 (2013) 264–271.
•  Washino  K.,  Tan  H.S.,  Hounslow  M.J.,  Salman  a.D., A new capillary force model implemented in micro-scale CFD-DEM coupling for wet granulation, Chem. Eng. Sci., 93 (2013) 197–205.
•  Wen  C.Y.,  Yu  Y.H., Mechanics of fluidization, Chem. Eng. Progr. Symp. Ser., 62 (1966) 100–111.
•  Wu  C.Y., DEM simulations of die filling during pharmaceutical tabletting, Particuology, 6 (2008) 412–418.
•  Wu  G.,  Ouyang  J.,  Yang  B.,  Li  Q., Use of compromise-based local porosities for coarse grid DEM simulation of bubbling fluidized bed with large particles, Adv. Powder Technol., 24 (2013) 68–78.
•  Xiao  F.,  Ii  S.,  Chen  C., Revisit to the THINC scheme: A simple algebraic VOF algorithm, J. Comput. Phys., 230 (2011) 7086–7092.
•  Xu  M.,  Chen  F.,  Liu  X.,  Ge  W.,  Li  J., Discrete particle simulation of gas-solid two-phase flows with multi-scale CPU-GPU hybrid computation, Chem. Eng. J., 207–208 (2012) 746–757.
•  Yamada  Y.,  Sakai  M., Lagrangian-Lagrangian simulations of solid-liquid flows in a bead mill, Powder Technol, 239 (2013) 105–114.
•  Yokoi  K., Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm, J. Comput. Phys., 226 (2007) 1985–2002.
•  Yue  X.,  Hao  Z.,  Chunhai  K.,  Congshu  L.,  Shi  S.,  Yuanqiang  T.,  Chunsheng  F., A GPU-based Discrete Element Modelling Code and its Application in Die Filling, Comput. Fluids, 110 (2014) 235–244.