Thus, the three vectors *v*1, *v*2, *v*3 are the three unit eigenvectors corresponding to the three largest eigenvalues of matrix *S*. This step is a principal component analysis (PCA) applied to the barycenters of the clusters. As a result, we choose the *d* three-dimensional axes of the 3*D* star coordinate system as *ai* = (*v*1*i*, *v*2*i*, *v*3*i*), *i* = 1, . . . , *d*.

Obviously, we can also project into 2D coordinates in the same way. However, when comparing and evaluating projections to 2D and 3D visual space (Poco et al., 2011), a quantitative analysis confirms that 3D projections outperform 2D projections in terms of precision. Moreover, a user study indicates that certain tasks can be more reliably and confidently answered with 3D projections. Nonetheless, as 3D projections are displayed on 2D screens, interaction is more difficult.

After having computed the projected clusters, we can display them using star coordinates by rendering a point primitive for each projected data point. A less cluttered and more beautiful display is to render the boundary of the clusters. Considering the cluster that is described by the set of points {*pi* = (*xi*, *yi*, *zi*) : *i* = 1, . . . , *m*} after being projected into the 3D space. In order to compute the boundary of this group of points,we need to have a continuous representation of the group. Therefore, we consider the function

$$\,\_{1}f\_{h}(p) = \sum\_{i=1}^{m} K\left(\frac{p - p\_{i}}{h}\right), \, p \in \mathbb{R}^{3},\tag{20}$$

where *K* is a kernel function and *h* is the bandwidth. Then, we can reconstruct the field over a regular grid and render the boundary set of the points by using standard isosurface extraction methods to extract the boundary surface of the set *S*(*h*, *c*) = {*p* ∈ **R**3 : *fh*(*p*) ≥ *c*}, where *c* is an isovalue. We choose parameter *h* and *c* to guarantee that *S*(*h*, *c*) is connected and has a volume of minimum extension. The kernel function should be sufficiently smooth and have a small compact support. For example, we can choose *K*(*p*)=(1 − ||*p*||2)2 for ||*p*|| ≤ 1 and *K*(*p*) = 0 otherwise and the bandwidth *h* to be equal to the longest length of the minimum spanning tree of these *m* points. In Figure 14(b) we show the visualization of the clusters by rendering such boundary surfaces, where it can be shown that for the chosen kernel isovalue *c* = 9 16 is appropriate. In order to visualize all clusters of the cluster tree, we render the surfaces in a semi-transparent fashion. The resulting visualization shows sequences of nested surfaces, where the inner surfaces represent higher density levels. Figure 14(b) shows the nested density cluster visualization with respect to the cluster tree in Figure 14(a).

#### **7.2 Coordinated views**

18 Will-be-set-by-IN-TECH

for any *d*-dimensional points *x* and *y*, the distance of the images of two *d*-dimensional points is lower than or equal to the distance of the points with respect to the *L*1-norm. Therefore, two points in the multi-dimensional space are projected to 3*D* star coordinates preserving the similarity properties of clusters (at least with respect to the *L*1-norm). In other words, the mapping of *d*-dimensional data to the 3D visual space does not break clusters. The second property that our projection from multi-dimensional feature space into three-dimensional star coordinate systems should fulfill is that separated clusters should not be projected into the same region. The projection into star coordinates may cause severe cluttering of clusters when not carefully choosing the axes (*a*1,..., *ad*). To alleviate the problem of overlapping clusters we introduce a method which chooses a "good" coordinate system. Assume that a hierarchy of high density clusters have *q* mode clusters, which do not contain any higher level densities. Let *mi* be the barycenter of the points within the *i*th cluster, *i* = 1, . . . , *q*. We want to choose a projection that maintains best the distances between clusters. Let {*v*1, *v*2, *v*3} be an orthonormal basis of the candidate three-dimensional space of projections. The desired choice of a 3*D* star coordinate layout is to maximize the distance of the *q* projected barycenters *VTmi*

*T*, i.e. to maximize the objective function

Thus, the three vectors *v*1, *v*2, *v*3 are the three unit eigenvectors corresponding to the three largest eigenvalues of matrix *S*. This step is a principal component analysis (PCA) applied to the barycenters of the clusters. As a result, we choose the *d* three-dimensional axes of the 3*D*

Obviously, we can also project into 2D coordinates in the same way. However, when comparing and evaluating projections to 2D and 3D visual space (Poco et al., 2011), a quantitative analysis confirms that 3D projections outperform 2D projections in terms of precision. Moreover, a user study indicates that certain tasks can be more reliably and confidently answered with 3D projections. Nonetheless, as 3D projections are displayed on

After having computed the projected clusters, we can display them using star coordinates by rendering a point primitive for each projected data point. A less cluttered and more beautiful display is to render the boundary of the clusters. Considering the cluster that is described by the set of points {*pi* = (*xi*, *yi*, *zi*) : *i* = 1, . . . , *m*} after being projected into the 3D space. In order to compute the boundary of this group of points,we need to have a continuous

> *p* − *pi h*

where *K* is a kernel function and *h* is the bandwidth. Then, we can reconstruct the field over a regular grid and render the boundary set of the points by using standard isosurface extraction

*S* = ∑ *i*<*j*

∑ *i*<*j*

star coordinate system as *ai* = (*v*1*i*, *v*2*i*, *v*3*i*), *i* = 1, . . . , *d*.

representation of the group. Therefore, we consider the function

*fh*(*p*) =

*m* ∑ *i*=1 *K*

2D screens, interaction is more difficult.

(*mi* − *mj*)(*mi* − *mj*)*T* . (19)

, *p* ∈ **R**3 , (20)

Since it can be shown that

with *V* = [*v*1, *v*2, *v*3]

with

Generating all clusters and displaying them in star coordinates allows for further analysis of the detected clusters. The simplest interaction method is to select individual clusters by just clicking at the boundary surface. When a cluster is selected, intra-cluster variability is visualized using parallel coordinates, see Figure 14(b) and (c). In both pictures the relation between the selected cluster with the dimension can be observed.

Moreover, we visualize the coordinated view in physical space, which exhibits the spatial location of the selected feature. The rendering in physical space can be preformed by just plotting all particles that belong to the selected feature or by extracting a boundary surface of that feature, i.e., a surface that separates all particles that belong to the feature from all particles that do not belong to the feature. Figure 15 shows an attribute-space rendering of the detected clusters in 3D optimized star coordinates (a), a color-coded object-space rendering of the clustered particles (b), and a separation surface of clusters in object space (c).

Fig. 15. (a) Seven-dimensional attribute space visualization of SPH data set using optimized 3D star coordinates. (b) Object space visualization of cluster distribution. (c) Object space visualization of a separating surface.

For the visualization of enclosing surfaces in attribute as well as in object space, we looked into an alternative approach of enclosing surfaces for point clusters using 3D discrete Voronoi diagrams (Rosenthal & Linsen, 2009). Our system provides three different types of enclosing surfaces. By generating a discrete distance field to the point cluster and extracting an isosurface from the field, an enclosing surface with any distance to the point cluster can be

Fig. 17. Different visualizations of two point clusters (colored red and blue) from the 2008 IEEE Visualization Design Contest data. The clusters were found using density-based clustering of multidimensional feature space and were projected to a 3D visual space using a linear projection. Additionally to the cluster points (a), three types of enclosing surfaces are shown. (b) Isosurface extraction from distance field computed using a 3D discrete Voronoi diagram of resolution 256 × 256 × 256. (c) Hull of the cluster computed from the isosurface

SmoothViz: Visualization of Smoothed Particles Hydrodynamics Data 23

Fig. 18. Coordinated views allow for selecting clusters in cluster tree and investigating properties in attribute space (using parallel coordinates) as well as locations in physical

**8. Interactive visual system for exploration of multiple scalar and flow fields**

Our research results are combined in the SmoothViz software system that is offered to the SPH community via our website (http://vcgl.jacobs-university.de/software). Not all presented features are included yet. Currently, the system consists of three modules responsible

space.

of the distance field. (d) Isosurface extraction from distance field to hull.

generated. As a second type of enclosing surfaces, a hull of the point cluster is extracted. The generation of the hull uses a projection of the discrete Voronoi diagram of the point cluster to an isosurface to generate a polygonal surface. Generated hulls of non-convex clusters are also non-convex. The third type of enclosing surfaces can be created by computing a distance field to the hull and extracting an isosurface from the distance field. This method exhibits reduced bumpiness and can extract surfaces arbitrarily close to the point cluster without losing connectedness. Figure 16 shows the idea of the different approaches starting from an isosurface from the distance field to the point cluster (a), connecting the neighbors that contribute to the surface in (a) to form a non-convex hull (b), and computing surfaces that are equidistant to the computed non-convex hull (b). Figure 17 shows a comparison of the different enclosing surfaces when applied to a cluster of points when projected into optimized star coordinates.

Fig. 16. (a) Extracting an isosurface from the distance field to the point cluster. Voronoi regions on the isosurface induce neighborhoods. (b) Neighbors are connected to form a hull. The image also shows an isosurface extracted from the distance field to the hull.

We extended our work on interactivity by explicitly encoding the cluster hierarchy in a tree that is visually encoded in a radial layout. Coordinated views between cluster tree visualization and parallel coordinates as well as object-space visualizations allow for an interactive analysis of multi-field SPH data (Linsen et al., 2009). The cluster tree allows for the selection of detected clusters, the parallel coordinate plots show the properties of the selected clusters, and object-space visualizations in form of extracted surfaces or particle distributions exhibit the location of the respective clusters in physical space. Figure 18 shows such a visual analysis set-up when applied to the IEEE Visualization Contest data (Rosenthal et al., 2008). We also proposed a method to integrate the parallel coordinates into the cluster tree visualization. The MultiClusterTree approach (Long & Linsen, 2011) uses circular parallel coordinates for the embedding into the radial hierarchical cluster tree layout, which allows for the analysis of the overall cluster distribution. This visual representation supports the comprehension of the relations between clusters and the original attributes. The combination of the 2D radial layout and the circular parallel coordinates is used to overcome the overplotting problem of parallel coordinates when looking into data sets with many records. Figure 19 shows how integrated circular coordinates can provide a good overview of the cluster distribution.

20 Will-be-set-by-IN-TECH

generated. As a second type of enclosing surfaces, a hull of the point cluster is extracted. The generation of the hull uses a projection of the discrete Voronoi diagram of the point cluster to an isosurface to generate a polygonal surface. Generated hulls of non-convex clusters are also non-convex. The third type of enclosing surfaces can be created by computing a distance field to the hull and extracting an isosurface from the distance field. This method exhibits reduced bumpiness and can extract surfaces arbitrarily close to the point cluster without losing connectedness. Figure 16 shows the idea of the different approaches starting from an isosurface from the distance field to the point cluster (a), connecting the neighbors that contribute to the surface in (a) to form a non-convex hull (b), and computing surfaces that are equidistant to the computed non-convex hull (b). Figure 17 shows a comparison of the different enclosing surfaces when applied to a cluster of points when projected into optimized

a) b)

We extended our work on interactivity by explicitly encoding the cluster hierarchy in a tree that is visually encoded in a radial layout. Coordinated views between cluster tree visualization and parallel coordinates as well as object-space visualizations allow for an interactive analysis of multi-field SPH data (Linsen et al., 2009). The cluster tree allows for the selection of detected clusters, the parallel coordinate plots show the properties of the selected clusters, and object-space visualizations in form of extracted surfaces or particle distributions exhibit the location of the respective clusters in physical space. Figure 18 shows such a visual analysis set-up when applied to the IEEE Visualization Contest data (Rosenthal et al., 2008). We also proposed a method to integrate the parallel coordinates into the cluster tree visualization. The MultiClusterTree approach (Long & Linsen, 2011) uses circular parallel coordinates for the embedding into the radial hierarchical cluster tree layout, which allows for the analysis of the overall cluster distribution. This visual representation supports the comprehension of the relations between clusters and the original attributes. The combination of the 2D radial layout and the circular parallel coordinates is used to overcome the overplotting problem of parallel coordinates when looking into data sets with many records. Figure 19 shows how integrated circular coordinates can provide a good overview of the

Fig. 16. (a) Extracting an isosurface from the distance field to the point cluster. Voronoi regions on the isosurface induce neighborhoods. (b) Neighbors are connected to form a hull.

The image also shows an isosurface extracted from the distance field to the hull.

star coordinates.

cluster distribution.

Fig. 17. Different visualizations of two point clusters (colored red and blue) from the 2008 IEEE Visualization Design Contest data. The clusters were found using density-based clustering of multidimensional feature space and were projected to a 3D visual space using a linear projection. Additionally to the cluster points (a), three types of enclosing surfaces are shown. (b) Isosurface extraction from distance field computed using a 3D discrete Voronoi diagram of resolution 256 × 256 × 256. (c) Hull of the cluster computed from the isosurface of the distance field. (d) Isosurface extraction from distance field to hull.

Fig. 18. Coordinated views allow for selecting clusters in cluster tree and investigating properties in attribute space (using parallel coordinates) as well as locations in physical space.

#### **8. Interactive visual system for exploration of multiple scalar and flow fields**

Our research results are combined in the SmoothViz software system that is offered to the SPH community via our website (http://vcgl.jacobs-university.de/software). Not all presented features are included yet. Currently, the system consists of three modules responsible

(a) (b)

(c) (d)

field (shown as splatted isosurface).

**9. Conclusion**

Fig. 20. Screenshots of SmoothViz software system for SPH data exploration: (a)

Three-dimensional particle distribution modeling a White Dwarf passing close to a Black Hole. (b) Maximal intensity projection plot of the density field of a White Dwarf with user defined transfer function; (c) Several density isosurfaces of two White Dwarfs in point-based representation. (d) Interplay of a velocity field (shown with streamlines) and a temperature

SmoothViz: Visualization of Smoothed Particles Hydrodynamics Data 25

We have presented approaches for visualization of SPH data. All methods operate directly on the particles that are distributed in a highly adaptive and irregular manner and that do not have any connectivity. Operating on the particles avoids the introduction of errors that occur when resampling to a grid. Our visualizations focus on surface extractions from such data. We first presented an isosurface extraction from any scalar field of the SPH data. It exploits a fast navigation through a *k*d-tree via an indexing structure and allows for fast isosurface extraction of high quality. Because of approximations made during simulation, it is desirable to add a smoothing term to the isosurface extraction method. This is achieved by the use of level-set methods. Again, the method operates on the particles only. We have presented several ways on how to accelerate the computations including a narrow-band approach, a local variational approach, and a signed distance function computation to any isosurface representation. Extracted isosurfaces are given in form of point clouds. We presented how they can be rendered using an image-space point cloud rendering approach that avoids any pre-computation and thus can immediately applied to any extracted surface. Shadows and transparency are supported at interactive rates. We further extended the work to the extraction of boundary surfaces of features in multi-field data. The attribute space of the multi-field data is being explored using clustering and cluster visualization

Fig. 19. Integrated circular parallel coordinates in clusters tree visualization for data set with hierarchical clusters.

for time-varying data manipulation, scalar field exploration, and flow field visualization. An intuitive graphical user interface (GUI) allows for easy processing and interaction. Additional functionalities and visualizations that are common in the SPH community have been included.

First, the user can load SPH data containing time-varying particle positions and time-varying multiple scalar and vector field values sampled at the particles. A 3D view of the particle distribution at a chosen time step allows the user to adjust the viewing parameters using arbitrary rotation and translation of camera. Loading of successive or preceding time steps from the time-varying series of data sets is as easy as play or rewind in a standard media player. Extracted pathlines can show evolution in time of an individual particle or sets of particles. Figure 20(a) shows the GUI and a particle distribution plot for a chosen time step.

There are two options to represent the structure of a selected scalar field: Maximal intensity projection plots can render any of the scalar fields using one of the build-in color maps and allowing for manually modifying the transfer function. Figure 20(b) shows the GUI for the transfer function modification and the respective maximum intensity plot of a chosen scalar field. Alternatively, isosurfaces can be extracted for interactively selected isovalues and shown using a point splatting technique or a dense point cloud rendering. Figure 20(c) shows a number of nested isosurfaces using point cloud renderings.

Finally, a specified number of streamlines can be computed with respect to the vector field chosen by the user. Combined views are possible to explore multiple fields simultaneously, e.g. multiple isosurfaces together with stream- or pathlines. Figure 20(d) shows an isosurface rendering using point splatting combined with a rendering of selected streamlines. For more details on the system, we refer to the user manual that comes with the SmoothViz software package.

Fig. 20. Screenshots of SmoothViz software system for SPH data exploration: (a) Three-dimensional particle distribution modeling a White Dwarf passing close to a Black Hole. (b) Maximal intensity projection plot of the density field of a White Dwarf with user defined transfer function; (c) Several density isosurfaces of two White Dwarfs in point-based representation. (d) Interplay of a velocity field (shown with streamlines) and a temperature field (shown as splatted isosurface).

#### **9. Conclusion**

22 Will-be-set-by-IN-TECH

Fig. 19. Integrated circular parallel coordinates in clusters tree visualization for data set with

for time-varying data manipulation, scalar field exploration, and flow field visualization. An intuitive graphical user interface (GUI) allows for easy processing and interaction. Additional functionalities and visualizations that are common in the SPH community have

First, the user can load SPH data containing time-varying particle positions and time-varying multiple scalar and vector field values sampled at the particles. A 3D view of the particle distribution at a chosen time step allows the user to adjust the viewing parameters using arbitrary rotation and translation of camera. Loading of successive or preceding time steps from the time-varying series of data sets is as easy as play or rewind in a standard media player. Extracted pathlines can show evolution in time of an individual particle or sets of particles. Figure 20(a) shows the GUI and a particle distribution plot for a chosen time step. There are two options to represent the structure of a selected scalar field: Maximal intensity projection plots can render any of the scalar fields using one of the build-in color maps and allowing for manually modifying the transfer function. Figure 20(b) shows the GUI for the transfer function modification and the respective maximum intensity plot of a chosen scalar field. Alternatively, isosurfaces can be extracted for interactively selected isovalues and shown using a point splatting technique or a dense point cloud rendering. Figure 20(c) shows a

Finally, a specified number of streamlines can be computed with respect to the vector field chosen by the user. Combined views are possible to explore multiple fields simultaneously, e.g. multiple isosurfaces together with stream- or pathlines. Figure 20(d) shows an isosurface rendering using point splatting combined with a rendering of selected streamlines. For more details on the system, we refer to the user manual that comes with the SmoothViz software

number of nested isosurfaces using point cloud renderings.

hierarchical clusters.

been included.

package.

We have presented approaches for visualization of SPH data. All methods operate directly on the particles that are distributed in a highly adaptive and irregular manner and that do not have any connectivity. Operating on the particles avoids the introduction of errors that occur when resampling to a grid. Our visualizations focus on surface extractions from such data. We first presented an isosurface extraction from any scalar field of the SPH data. It exploits a fast navigation through a *k*d-tree via an indexing structure and allows for fast isosurface extraction of high quality. Because of approximations made during simulation, it is desirable to add a smoothing term to the isosurface extraction method. This is achieved by the use of level-set methods. Again, the method operates on the particles only. We have presented several ways on how to accelerate the computations including a narrow-band approach, a local variational approach, and a signed distance function computation to any isosurface representation. Extracted isosurfaces are given in form of point clouds. We presented how they can be rendered using an image-space point cloud rendering approach that avoids any pre-computation and thus can immediately applied to any extracted surface. Shadows and transparency are supported at interactive rates. We further extended the work to the extraction of boundary surfaces of features in multi-field data. The attribute space of the multi-field data is being explored using clustering and cluster visualization

Gingold, R. A. & Monaghan, J. J. (1977). Smoothed particle hydrodynamics - Theory and

SmoothViz: Visualization of Smoothed Particles Hydrodynamics Data 27

Grossman, J. P. & Dally, W. J. (1998). Point sample rendering, *Proceedings of 9th Eurographics*

Hopf, M. & Ertl, T. (2003). Hierarchical splatting of scattered data, *Visualization Conference,*

Hopf, M., Luttenberger, M. & Ertl, T. (2004). Hierarchical splatting of scattered 4d data, *IEEE*

Linsen, L., Long, T. V. & Rosenthal, P. (2009). Linking multi-dimensional feature space cluster

Linsen, L., Long, T. V., Rosenthal, P. & Rosswog, S. (2008). Surface extraction from multi-field

Linsen, L., Müller, K. & Rosenthal, P. (2007). Splat-based ray tracing of point clouds, *Journal of*

Long, T. V. (2010). *Visualizing High Density Clusters in Multidimensional Data*, PhD thesis, Jacobs

Long, T. V. & Linsen, L. (2011). Visualizing high density clusters in multidimensional data using optimized star coordinates, *Journal of Computational Statistics (to appear)* . Lucy, L. B. (1977). A numerical approach to the testing of the fission hypothesis, *Astronomical*

Maur, P. & Kolingerová, I. (2001). Post-optimization of delaunay tetrahedrization, *SCCG '01:*

Molchanov, V., Rosenthal, P. & Linsen, L. (2010). Non-iterative second-order approximation of

Molchanov, V., Rosenthal, P. & Linsen, L. (2011). Variational level-set detection of local

Navratil, P. A., Johnson, J. L. & Bromm, V. (2007). Visualization of cosmological particle-based

Pav, S. E. & Walkington, N. J. (2004). Robust three dimensional delaunay refinement, *13th*

Phong, B. T. (1975). Illumination for computer generated pictures, *Commun. ACM* 18: 311–317. Poco, J., Etemadpour, R., Paulovich, F. V., Long, T. V., Rosenthal, P., de Oliveira, M. C. F.,

Price, D. (2007). SPLASH: An interactive visualisation tool for smoothed particle

Rosenberg, I. D. & Birdwell, K. (2008). Real-time particle isosurface extraction, *Proceedings of*

with 3d projections, *Computer Graphics Forum* 30(3): 1111–1120.

*Proceedings of the 17th Spring conference on Computer graphics*, IEEE Computer Society,

signed distance function for any isosurface representation, *Computer Graphics Forum*

isosurfaces from unstructured point-based volume data, *Schloss Dagstuhl Scientific*

*International Meshing Roundtable*, Sandia National Laboratories, SAND 2004-3765C,

Linsen, L. & Minghim, R. (2011). A framework for exploring multidimensional data

hydrodynamics simulations, *Publications of the Astronomical Society of Australia*

*the 2008 symposium on Interactive 3D graphics and games*, I3D '08, ACM, New York, NY,

visualization to surface extraction from multi-field volume data, *IEEE Computer*

particle volume data using multi-dimensional cluster visualization, *IEEE Transactions*

181: 375–389.

*IEEE* .

*Workshop on Rendering*, pp. 181–192.

*Graphics and Applications* 29(3): 85–89.

*WSCG* 15(1–3): 51–58.

*Journal* 82: 1013–1024.

29(3): 1211–1220.

pp. 145–156.

24: 159–173.

USA, pp. 35–43.

Washington, DC, USA, p. 31.

datasets, 13(6): 1712–1718.

*Visualization Workshop 2009 Follow-up, to appear*.

University.

*Computer Graphics and Applications* 24: 64–72.

*on Visualization and Computer Graphics* 14(6): 1483–1490.

application to non-spherical stars, *Monthly Notices of the Royal Astronomical Society*

methods. Coordinated or integrated views to parallel or circular coordinates, respectively, allow for further visual analysis of the properties of the extracted clusters. Coordinated views to object space allow for the investigation of the spatial distribution of detected features. Enclosing surfaces show the cluster boundaries. The presented functionality has partially been incorporated into the SmoothViz software package including further features such as geometric flow visualization. It allows for interactive exploration and integrated analysis of multiple fields of SPH data.

#### **10. Acknowledgments**

This work was supported by the German Research Foundation (DFG) under grant number LI 1530/6-1.

#### **11. References**

24 Will-be-set-by-IN-TECH

methods. Coordinated or integrated views to parallel or circular coordinates, respectively, allow for further visual analysis of the properties of the extracted clusters. Coordinated views to object space allow for the investigation of the spatial distribution of detected features. Enclosing surfaces show the cluster boundaries. The presented functionality has partially been incorporated into the SmoothViz software package including further features such as geometric flow visualization. It allows for interactive exploration and integrated analysis of

This work was supported by the German Research Foundation (DFG) under grant number LI

Biddiscombe, J., Graham, D. & Maruzewski, P. (2008). Visualization and analysis of SPH data,

Borouchaki, H., Hecht, F., Saltel, E. & George, P. (1995). Reasonably efficient delaunay

CGAL (2011). Computational geometry algorithms library (CGAL), http://www.cgal.org/. Cha, D., Son, S. & Ihm, I. (2009). Gpu-assisted high quality particle rendering, *Computer*

Co, C. S. & Joy, K. I. (2005). Isosurface Generation for Large-Scale Scattered Data Visualization,

Courant, R., Friedrichs, K. & Lewy, H. (1928). Über die partiellen differenzengleichungen der

Dobrev, P., Rosenthal, P. & Linsen, L. (2010a). An image-space approach to interactive

Dobrev, P., Rosenthal, P. & Linsen, L. (2010b). Interactive image-space point cloud rendering

Du, Q. & Wang, D. (2006). Recent progress in robust and quality delaunay mesh generation,

Fraedrich, R., Auer, S. & Westermann, R. (2010). Efficient high-quality volume rendering of sph data, *IEEE Transactions on Visualization and Computer Graphics* 16: 1533–1540. Fraedrich, R., Schneider, J. & Westermann, R. (2009). Exploring the "millennium run" - scalable

George, P. L., Hecht, F. & Saltel, E. (1991). Automatic mesh generator with specified boundary,

mathematischen physik, *Mathematische Annalen* 100(1): 32 – 74. Delaunay, B. N. (1934). Sur la sphere vide, *Bull. Acad. Sci. USSR* 7: 793–800.

based mesh generator in 3 dimensions, *4th International Meshing Roundtable*, Sandia

*in* G. Greiner, J. Hornegger, H. Niemann & M. Stamminger (eds), *Proceedings of Vision, Modeling, and Visualization 2005*, Akademische Verlagsgesellschaft Aka

point cloud rendering including shadows and transparency, *Computer Graphics and*

with transparency and shadows, *in* V. Skala (ed.), *Communication Papers Proceedings of WSCG, The 18th International Conference on Computer Graphics, Visualization and Computer Vision*, UNION Agency – Science Press, Plzen, Czech Republic, pp. 101–108.

rendering of large-scale cosmological datasets, *IEEE Transactions on Visualization and*

multiple fields of SPH data.

*ERCOFTAC Bulletin* 76(9-12).

National Laboratories, pp. 3–14.

*Graphics Forum* 28(4): 1247–1255.

*J. Comput. Appl. Math.* 195(1): 8–23.

*Computer Graphics* 15(6): 1251–1258.

*Comput. Methods Appl. Mech. Eng.* 92(3): 269–288.

GmbH, pp. 233–240.

*Geometry* 12(3): 2–25.

**10. Acknowledgments**

1530/6-1.

**11. References**

**0**

**2**

1*USA* 2*UK*

Donghong Gao1 and Jin Sun2

*Edinburgh, Edinburgh EH9 3JL, Scotland*

**Using DEM in Particulate Flow Simulations**

2*Institute for Infrastructure and Environment, School of Engineering, The University of*

There is no doubt that the dynamics of solid particles is the center of interest in the mineral processing industry, including crushing, grinding, classification, mineral separation, and leaching, just to name a few. But most processes involve fluid as a carrier and/or media, making the study of fluid dynamics and coupling between fluid dynamics and particle motion essential part for such a particulate flow. Various modeling and coupling approaches capable of considering particle behaviors, fluid dynamics and coupling effects, have been actively

By name, particulate flows usually include one or more continuous fluid phase, and one or more type of particles, or say generally, discrete phases. The discrete particles/bubbles/ droplets are often dispersed in a continuous phase, so a discrete phase is also called a disperse phase in continuum multiphase modeling. There may be strong interactions between discrete phases and continuous phase, and strong interactions among discrete phases for dense particulate flows. The coupling physics pose a huge challenge to researchers, since coupling physics between fluid dynamics and particle motion requires coupling numerical modeling approaches. There is no one-fits-all solution for all applications especially after

As computer technology for hardware and software advances so rapidly, it also push scientific and engineering simulations to high standards of requirement with respects to accuracy, fidelity, efficiency. There is increasing research activity of using Discrete Element Method

Discrete Element Method (DEM) (Cundall & Strack, 1979; Landry et al., 2003; Walton, 1992) is a Lagrangian model and is well accepted nowadays to model solid particle behavior. In principle, the DEM is based on the concept that individual particles, each of which is usually assumed to be semi-rigid, are considered to be separate and are connected at boundaries by appropriate contact laws. DEM naturally captures characteristics of each particle, therefore further dynamics like breakage and wear can be modeled locally at the small scale. Using DEM to track dynamics of particles, although the computing cost is high, eliminates the need of modeling fluid dynamics of particle phase, therefore improves fidelity of simulations. Interactions among discrete phases can be addressed more accurately inside DEM, thank to

considering limitations, accuracy, computational costs of various numerical models.

that microscopic physics has been clearly understood and described at most times.

**1. Introduction**

pursued, researched and developed in recent years.

(DEM) in particulate flow simulations.

1*Optimization Services, Metso Minerals, Colorado Springs, CO 80903*

### **Using DEM in Particulate Flow Simulations**

#### Donghong Gao1 and Jin Sun2

1*Optimization Services, Metso Minerals, Colorado Springs, CO 80903* 2*Institute for Infrastructure and Environment, School of Engineering, The University of Edinburgh, Edinburgh EH9 3JL, Scotland* 1*USA* 2*UK*

#### **1. Introduction**

26 Will-be-set-by-IN-TECH

28 Hydrodynamics – Optimizing Methods and Tools

Rosenthal, P. (2009). *Direct Surface Extraction from Unstructured Point-based Volume Data*, Shaker Verlag, Aachen, Germany (Ph.D. Thesis Jacobs University, Bremen, Germany). Rosenthal, P. & Linsen, L. (2006). Direct isosurface extraction from scattered volume data,

Rosenthal, P. & Linsen, L. (2008a). Image-space point cloud rendering, *Proceedings of Computer*

Rosenthal, P. & Linsen, L. (2008b). Smooth surface extraction from unstructured point-based

Rosenthal, P. & Linsen, L. (2009). Enclosing surfaces for point clusters using 3d discrete

Rosenthal, P., Long, T. V. & Linsen, L. (2008). "Shadow Clustering": Surface extraction

Rosenthal, P., Molchanov, V. & Linsen, L. (2010). A narrow band level set method for surface

Rosenthal, P., Rosswog, S. & Linsen, L. (2007). Direct surface extraction from smoothed

Rosswog, S. (2009). Astrophysical smooth particle hydrodynamics, *New Astronomy Reviews*

Rosswog, S., Ramirez-Ruiz, E. & Hix, W. R. (2009). Tidal Disruption and Ignition of White Dwarfs by Moderately Massive Black Holes, *Astrophysical Journal* 695: 404–419. Sapidis, N. S. & Perucchio, R. (1991). Domain delaunay tetrahedrization of arbitrarily shaped

Schindler, B., Fuchs, R., Biddiscombe, J. & Peikert, R. (2009). Predictor-corrector schemes

Walker, R., Kenny, P. & Miao, J. (2005). Visualization of Smoothed Particle Hydrodynamics

voronoi diagrams, *Computer Graphics Forum* 28(3): 999–1006.

*Visualization - EuroVis 2006*, pp. 99–106,367.

*Workshop*, Lehmanns Media - LOB, pp. 50–61.

*Visualization and Computer Graphics* 15: 1243–1250.

(Electronic version http://diglib.eg.org). URL: *http://www.cs.kent.ac.uk/pubs/2005/2230*

New York, NY, USA, pp. 465–480.

*Graphics International*, pp. 136–143.

Visualization Design Contest.

14(6): 1531–1546.

53(4-6): 78 – 104.

*in* B. S. Santos, T. Ertl & K. I. Joy (eds), *Eurographics / IEEE VGTC Symposium on*

volume data using PDEs, *IEEE Transactions on Visualization and Computer Graphics*

from non-equidistantly sampled multi-field 3D scalar data using multi-dimensional cluster visualization, VisWeek 08 Conference Compendium, Winner of IEEE

extraction from unstructured point-based volume data, *in* V. Skala (ed.), *Proceedings of WSCG, The 18th International Conference on Computer Graphics, Visualization and Computer Vision*, UNION Agency – Science Press, Plzen, Czech Republic, pp. 73–80.

particle hydrodynamics simulation data, *Proceedings of the 4th High-End Visualization*

curved polyhedra defined in a solid modeling system, *SMA '91: Proceedings of the first ACM symposium on Solid modeling foundations and CAD/CAM applications*, ACM Press,

for visualization ofsmoothed particle hydrodynamics data, *IEEE Transactions on*

for Astrophysics, *in* L. Lever & M. McDerby (eds), *Theory and Practice of Computer Graphics 2005*, Eurographics Association, University of Kent, UK, pp. 133–138. There is no doubt that the dynamics of solid particles is the center of interest in the mineral processing industry, including crushing, grinding, classification, mineral separation, and leaching, just to name a few. But most processes involve fluid as a carrier and/or media, making the study of fluid dynamics and coupling between fluid dynamics and particle motion essential part for such a particulate flow. Various modeling and coupling approaches capable of considering particle behaviors, fluid dynamics and coupling effects, have been actively pursued, researched and developed in recent years.

By name, particulate flows usually include one or more continuous fluid phase, and one or more type of particles, or say generally, discrete phases. The discrete particles/bubbles/ droplets are often dispersed in a continuous phase, so a discrete phase is also called a disperse phase in continuum multiphase modeling. There may be strong interactions between discrete phases and continuous phase, and strong interactions among discrete phases for dense particulate flows. The coupling physics pose a huge challenge to researchers, since coupling physics between fluid dynamics and particle motion requires coupling numerical modeling approaches. There is no one-fits-all solution for all applications especially after considering limitations, accuracy, computational costs of various numerical models.

As computer technology for hardware and software advances so rapidly, it also push scientific and engineering simulations to high standards of requirement with respects to accuracy, fidelity, efficiency. There is increasing research activity of using Discrete Element Method (DEM) in particulate flow simulations.

Discrete Element Method (DEM) (Cundall & Strack, 1979; Landry et al., 2003; Walton, 1992) is a Lagrangian model and is well accepted nowadays to model solid particle behavior. In principle, the DEM is based on the concept that individual particles, each of which is usually assumed to be semi-rigid, are considered to be separate and are connected at boundaries by appropriate contact laws. DEM naturally captures characteristics of each particle, therefore further dynamics like breakage and wear can be modeled locally at the small scale. Using DEM to track dynamics of particles, although the computing cost is high, eliminates the need of modeling fluid dynamics of particle phase, therefore improves fidelity of simulations. Interactions among discrete phases can be addressed more accurately inside DEM, thank to that microscopic physics has been clearly understood and described at most times.

neighbors based on the theory of integral interpolants using kernel functions which can be differentiated without use of the grid. SPH, as a Lagrangian particle–based method, has its particular characteristics. It has some special advantages over conventional grid–based CFD. The most significant one is the meshfree feature. SPH does not require a pre-defined mesh to provide connections between particles when solving the governing equations. The SPH

Using DEM in Particulate Flow Simulations 31

From a numerical implementation point of view, DEM–SPH, a Lagrangian–Lagrangian model, is the best incorporation for particulate flows because it can totally eliminate the need for a volume mesh. The meshfree feature is very attractive to the mineral processing industry,

However, SPH is not without problems. While Eulerian–Eulerian approaches, modeling naturally discrete particles as continua, have difficulty to give correct constitutive equations, similarly, SPH, modeling naturally continuous media as particles, compromises accuracy in some aspects. It resolves the dissipative term poorly in comparison with grid–based methods. SPH has a limited ability to deal with steep density gradient or other large property changes. Boundary conditions do not fit naturally in the particle approach, so they are difficult to implement in SPH. It is hard to capture fluid dynamics where complex boundary conditions

In summary, there is no single one-fits-all solution. Every model has strengths on some aspects and weaknesses on others, especially considering accuracy and cost factors. People have to be able to, based on preliminary understanding of the physical characteristics of a system of interest, pick up the right models, combine them together, and develop/use appropriate models for the specific system to capture major phenomena to discover and investigate the controlling mechanisms behind them. In this work we present three numerical coupling approaches to capture the physics of interest: one-way coupling with CFD, DEM-CFD

A one-way coupling is basically to run fluid dynamic solver separately from DEM, then import fluid flow solution to DEM, where the fluid effect on particles is considered. The one-way approach is practically important for industry applications, because at complex situations full coupling modelings are hard to converge, if not impossible. It has the advantages of using commercial package. This advantage may become very attractive in industry applications where the flow condition is very complex and density of particles is not very high. A one-way coupling can be extended to so-called 1.5-way coupling if multiphase fluid solver is used instead of single phase solver. We applied one-way coupling to a slurry pump, where solid particles and fluid are well mixed so that it is appropriate to treat slurry as a kind of single phase mixture, and the FLUENT is used to solve the flow field. But the CFD solver cannot give direct answers to our concerns: wear effect of particles on pump structure and particle breaking probability, therefore CFD results are imported into our DEM code to

For the strong coupling physics, the full coupling of DEM-CFD or DEM-SPH is necessary. In the DEM-CFD coupling, we employ a lot of widely accepted Eulerian–Eulerian multi-fluid models that have been intensively studied in the continuum multiphase fluid dynamics. Convergence of the coupling models is usually a huge challenge. The numerical methods are discussed in each section of model descriptions. The segregation of different sizes of particles in a fluidization bed is controlled by both particles motion and fluid dynamics. Due to the

simple geometry of a bed, DEM-CFD is the best candidate for this application.

particles themselves are adaptive to geometry and free surface confinement.

are of critical importance.

coupling, DEM-SPH coupling.

simulate the detailed behavior of individual particles.

where geometry is complex and free surface flow is typical in many applications.

For fluid dynamics modeling, the terminology Computational Fluid Dynamics (CFD) is well-known dedicated to it. A conventional CFD is usually based on continuum mechanics principles and control volume methods. After decades of intensive research, conventional CFD is a well-developed technology with a series of mature well-defined numerical and physical models for single phase, turbulence and multiphase flows. It is basically a grid-based Eulerian model and is usually computationally efficient, especially for single phase flows.

In practice, using CFD is not as easy as expected, since most mineral processing applications involve complex geometry and free surface flows. Generating appropriate volume mesh for complex geometry is a challenge even with help of commercial programs. Free surface flow also adds cost and difficulty to simulations. In the commonly used VOF method, the entire possible physical domain, even the space that is to be occupied by fluid occasionally, has to be meshed, and interface capturing and reconstruction scheme have to be implemented (Gao et al., 2003).

Among CFD approaches, the one named Eulerian–Eulerian model, or say, multi-fluid model, has been extensively studied, implemented in MFIX (Syamlal, 1998; Syamlal et al., 1993), CFX, and FLUENT, and applied to simulations of fluidized beds (Gera et al., 2004; Sun & Battaglia, 2006a). The formulation of this model is essentially based on the continuum fluid dynamics. It considers both the fluid phase and solid particles to be interpenetrating continua whose dynamics are governed by the Navier-Stokes equations (Goldhirsch, 2003; Huilin et al., 2003; Savage, 1998). The particle mixture can be divided into several disperse phases with different properties. Closure of the model requires formulation of constitutive equations for each phases and inter-phase momentum transfer models, where often the most difficulties are encountered and approximations are made (Jenkins & Savage, 1983; Srivastava & Sundaresan, 2003). From the difficulty of building continuum models for granular flows (Gao et al., 2006; Savage, 1998), people realized that many of the physics–based governing equations work well at small scale, but non-linear physics makes derivations of those equations at a larger scale based on simplifications and assumptions no longer valid. For particulate flows with a wide property distribution, modeling errors are easily accumulated and computation costs are largely amplified.

The major drawback of the Eulerian–Eulerian approach is that it cannot capture essential characteristics of individual solid particles regarding size and shape, and thus cannot effectively identify influence of these characteristics on process performances. Contact of individual particles with structure is often the major source of wear and erosion. Sizeand shape-change processes, such as breakage and chemical reaction of individual particles, usually are core features of the mineral processing industry.

However, the extensive Eulerian–Eulerian research laid solid foundation for coupling DEM with CFD. Lots of ideas and equations can be adopted in DEM-CFD full coupling. DEM-CFD assumes the high theoretical fidelity, since each phase is kept to have its natural properties. We treat fluid as continuous continuum, and particles as discrete entities. The basic concepts of interpenetrating phases for multiphase flows still hold, although we only need to model and compute fluid phase. Instead of modeling several disperse particle phases in kinetic theory (Savage, 1998), we derive particle motions directly from DEM, therefore improves modeling accuracy.

Smooth Particle Hydrodynamics (SPH) has been used to simulate fluid dynamics for years (Monaghan, 1988; 1994; Morris et al., 1997). In SPH, a fluid field is represented by particles, each of which is associated with a mass, density, velocity, viscosity, pressure and position. Particles are moved by averaging (smoothing) their interaction with spatial 2 Will-be-set-by-IN-TECH

For fluid dynamics modeling, the terminology Computational Fluid Dynamics (CFD) is well-known dedicated to it. A conventional CFD is usually based on continuum mechanics principles and control volume methods. After decades of intensive research, conventional CFD is a well-developed technology with a series of mature well-defined numerical and physical models for single phase, turbulence and multiphase flows. It is basically a grid-based Eulerian model and is usually computationally efficient, especially for single phase flows. In practice, using CFD is not as easy as expected, since most mineral processing applications involve complex geometry and free surface flows. Generating appropriate volume mesh for complex geometry is a challenge even with help of commercial programs. Free surface flow also adds cost and difficulty to simulations. In the commonly used VOF method, the entire possible physical domain, even the space that is to be occupied by fluid occasionally, has to be meshed, and interface capturing and reconstruction scheme have to be implemented (Gao

Among CFD approaches, the one named Eulerian–Eulerian model, or say, multi-fluid model, has been extensively studied, implemented in MFIX (Syamlal, 1998; Syamlal et al., 1993), CFX, and FLUENT, and applied to simulations of fluidized beds (Gera et al., 2004; Sun & Battaglia, 2006a). The formulation of this model is essentially based on the continuum fluid dynamics. It considers both the fluid phase and solid particles to be interpenetrating continua whose dynamics are governed by the Navier-Stokes equations (Goldhirsch, 2003; Huilin et al., 2003; Savage, 1998). The particle mixture can be divided into several disperse phases with different properties. Closure of the model requires formulation of constitutive equations for each phases and inter-phase momentum transfer models, where often the most difficulties are encountered and approximations are made (Jenkins & Savage, 1983; Srivastava & Sundaresan, 2003). From the difficulty of building continuum models for granular flows (Gao et al., 2006; Savage, 1998), people realized that many of the physics–based governing equations work well at small scale, but non-linear physics makes derivations of those equations at a larger scale based on simplifications and assumptions no longer valid. For particulate flows with a wide property distribution, modeling errors are easily accumulated and computation costs

The major drawback of the Eulerian–Eulerian approach is that it cannot capture essential characteristics of individual solid particles regarding size and shape, and thus cannot effectively identify influence of these characteristics on process performances. Contact of individual particles with structure is often the major source of wear and erosion. Sizeand shape-change processes, such as breakage and chemical reaction of individual particles,

However, the extensive Eulerian–Eulerian research laid solid foundation for coupling DEM with CFD. Lots of ideas and equations can be adopted in DEM-CFD full coupling. DEM-CFD assumes the high theoretical fidelity, since each phase is kept to have its natural properties. We treat fluid as continuous continuum, and particles as discrete entities. The basic concepts of interpenetrating phases for multiphase flows still hold, although we only need to model and compute fluid phase. Instead of modeling several disperse particle phases in kinetic theory (Savage, 1998), we derive particle motions directly from DEM, therefore improves

Smooth Particle Hydrodynamics (SPH) has been used to simulate fluid dynamics for years (Monaghan, 1988; 1994; Morris et al., 1997). In SPH, a fluid field is represented by particles, each of which is associated with a mass, density, velocity, viscosity, pressure and position. Particles are moved by averaging (smoothing) their interaction with spatial

usually are core features of the mineral processing industry.

et al., 2003).

are largely amplified.

modeling accuracy.

neighbors based on the theory of integral interpolants using kernel functions which can be differentiated without use of the grid. SPH, as a Lagrangian particle–based method, has its particular characteristics. It has some special advantages over conventional grid–based CFD. The most significant one is the meshfree feature. SPH does not require a pre-defined mesh to provide connections between particles when solving the governing equations. The SPH particles themselves are adaptive to geometry and free surface confinement.

From a numerical implementation point of view, DEM–SPH, a Lagrangian–Lagrangian model, is the best incorporation for particulate flows because it can totally eliminate the need for a volume mesh. The meshfree feature is very attractive to the mineral processing industry, where geometry is complex and free surface flow is typical in many applications.

However, SPH is not without problems. While Eulerian–Eulerian approaches, modeling naturally discrete particles as continua, have difficulty to give correct constitutive equations, similarly, SPH, modeling naturally continuous media as particles, compromises accuracy in some aspects. It resolves the dissipative term poorly in comparison with grid–based methods. SPH has a limited ability to deal with steep density gradient or other large property changes. Boundary conditions do not fit naturally in the particle approach, so they are difficult to implement in SPH. It is hard to capture fluid dynamics where complex boundary conditions are of critical importance.

In summary, there is no single one-fits-all solution. Every model has strengths on some aspects and weaknesses on others, especially considering accuracy and cost factors. People have to be able to, based on preliminary understanding of the physical characteristics of a system of interest, pick up the right models, combine them together, and develop/use appropriate models for the specific system to capture major phenomena to discover and investigate the controlling mechanisms behind them. In this work we present three numerical coupling approaches to capture the physics of interest: one-way coupling with CFD, DEM-CFD coupling, DEM-SPH coupling.

A one-way coupling is basically to run fluid dynamic solver separately from DEM, then import fluid flow solution to DEM, where the fluid effect on particles is considered. The one-way approach is practically important for industry applications, because at complex situations full coupling modelings are hard to converge, if not impossible. It has the advantages of using commercial package. This advantage may become very attractive in industry applications where the flow condition is very complex and density of particles is not very high. A one-way coupling can be extended to so-called 1.5-way coupling if multiphase fluid solver is used instead of single phase solver. We applied one-way coupling to a slurry pump, where solid particles and fluid are well mixed so that it is appropriate to treat slurry as a kind of single phase mixture, and the FLUENT is used to solve the flow field. But the CFD solver cannot give direct answers to our concerns: wear effect of particles on pump structure and particle breaking probability, therefore CFD results are imported into our DEM code to simulate the detailed behavior of individual particles.

For the strong coupling physics, the full coupling of DEM-CFD or DEM-SPH is necessary. In the DEM-CFD coupling, we employ a lot of widely accepted Eulerian–Eulerian multi-fluid models that have been intensively studied in the continuum multiphase fluid dynamics. Convergence of the coupling models is usually a huge challenge. The numerical methods are discussed in each section of model descriptions. The segregation of different sizes of particles in a fluidization bed is controlled by both particles motion and fluid dynamics. Due to the simple geometry of a bed, DEM-CFD is the best candidate for this application.

assuming soft particles:

criterion |**F**t| ≤ *μ*|**F**n|.

here,

a triangle.

boundary triangle as:

**F**n,*ij* = *k*nδn,*ij* − *λ*n *m*eff **u**n,*ij* , (4) **F**t,*ij* = *k*tδt,*ij* − *λ*t *m*eff **u**t,*ij* , (5)

δn,ij = [(*Ri* + *Rj*) − *rij*]**n***ij*, (6)

**r***ij* = **r***i* − **r***j*, *rij* = |**r***ij*| , (7)

*e*n = exp(−*λ*n*t*c/2), (10)

n/4)−1/2 . (11)

*C*wr,*iE*sh,*i*, (12)

**n***ij* = **r***ij*/*rij* , (8) **n***ij* × **t***ij* = 1 . (9)

where *m*eff = (*mi mj*)/(*mi* + *mj*), subscript *n* denotes the normal direction and *t* the tangential direction, *k* is the spring coefficient, *λ* the damping coefficient, and δ the elastic displacement for a contact, respectively. The damping effect in the tangential direction can be neglected. The tangential displacement δt between particles is obtained by integrating surface relative velocities over time during deformation of the contact. Actually just this history dependent feature makes the computation, especially parallel computation, more expensive. The magnitude of δt is truncated as necessary to satisfy a local Coulomb yield

Using DEM in Particulate Flow Simulations 33

In the above DEM model, the normal compression between two particles is easily written as:

The value of the spring constant should be large enough to avoid particle interpenetration, but not too large to require an unreasonably small time step Δ*t*, since an accurate simulation typically requires Δ*t* ∼ *t*c/50, here *t*c is the characteristic contact time during a collision process between particles. The amount of energy lost in collisions is characterized by the inelasticity through the coefficient of restitution *e*. For the linear spring-dashpot model, the

following relations can be taken as guidance to find the damping coefficient *λ*n

*t*c = *π*(*k*n/*m*eff − *λ*2

Δ*h A* = ∑ *i*

After contact force is calculated, the equations of motion, which are ordinary differential equations, can be numerically integrated to get the particle trajectories. The boundary surfaces are represented by triangles. Any meshing tool generating surface triangles can be used here. In comparison with particle–particle contact, the only difference is the geometry resolution for the particle–triangle contact. The overlap *δ*n is equal to a particle radius minus the distance to

The advantage of DEM is that it can capture behaviors of individual particles, and collect detailed contact information, such as velocity, contact force, shear and impact energy spectrum, so that it can model wear more accurately. Most severe wear happens in particulate flows when solid particles collide on a surface. We model the particle wear effect on a

The normal direction and tangential direction are defined as:

In the DEM-SPH coupling, a multiphase DEM-SPH model is proposed and described in detail. The coupling between the solid particles and the fluid phase are considered via volume fraction, pressure and drag force. In mills, particles are so dense that the behaviors of solid particles dominate the overall physics, and separation of solid particles from fluid is apparent. But fluid dynamics cannot be neglected because fluid damping/drag may change particle–particle contacts and particle–boundary contacts. Furthermore, the rotating nature, presence of a free surface flow, and liner geometry make CFD simulations of mills prohibitively expensive. Numerical scheme is simple at this stage, time stepping scheme could be either explicit Euler or lower order Runge-Kutta for more accuracy. The DEM-SPH is applied to a preliminary study of a mill.

#### **2. DEM model with flow effect**

The DEM simulation is based on a 3D soft particle model (Cundall & Strack, 1979; Silbert et al., 2001; Walton & Braun, 1986) where small deformations and multiple contacts on a particle are allowed, and friction and rotation are also taken into account. Each particle has six degrees of freedom of motion. For simplification of the description of the DEM model, spherical particles are assumed, although in our real application and implementation, we use spherical, tetrahedral and irregular convex shape particles together. The complex shapes definitely add difficulty to the computational geometric and solving for particle rotation, but the basic theory behind it is the same as for all spherical particles. Although our program is fully parallelized for a distributive memory machine, this work is not intended to cover the parallelization scheme and readers who are interested are referred to Plimpton (1995), which is the framework we follow.

The movement of a particle with mass *m*, moment of inertia *I* can be described by Newton's law and the kinematic relation:

$$m\_i \frac{\mathbf{d} \,\mathbf{u}\_i}{\mathbf{d} \,t} = \sum\_j \mathbf{F}\_{\mathbf{c},ij} + \mathbf{F}\_{\mathbf{s}\mathbf{f},i} + m\_i \,\mathbf{g}\_{\mathbf{v}} \tag{1}$$

$$I\_i \frac{\mathbf{d} \,\omega\_i}{\mathbf{d} \,t} = -\sum\_j \mathbf{R}\_i \, \mathbf{n}\_{ij} \times \mathbf{F}\_{\mathbf{c},jj} \,\,\,\,\,\tag{2}$$

$$\frac{\mathbf{d}\,\mathbf{r}\_{l}}{\mathbf{d}\,t} = \mathbf{u}\_{l}.\tag{3}$$

where subscripts *i* and *j* are for identifying particles, **u** is the velocity of the mass center, ω the angular velocity, **r** the position, *R* the radius of the particle, **F**c,*ij* the contact force of particle *j* acting on *i*, **F**sf,*i* the fluid–solid interaction which is assumed to act at the mass center of particles, and **n***ij* the normal direction of the contact pointing to particle *i* from *j*.

The governing equations are simply ordinary differential equations in time. Each particle is evolved by integrating the governing equations and applying the initial condition. The major task of modeling and simulation thus becomes actually formulating and calculating the force terms.

The implementation of contact forces is essentially a reduced version of that employed by Walton & Braun (1986), developed earlier by Cundall & Strack (1979). Contact force **F**c,*ij* = {**F**n,*ij*, **F**t,*ij*} is first calculated from the deformation through the spring-dashpot model, assuming soft particles:

$$\mathbf{F}\_{\mathbf{n},ij} = k\_{\mathbf{n}} \delta\_{\mathbf{n},ij} - \lambda\_{\mathbf{n}} m\_{\mathbf{eff}} \mathbf{u}\_{\mathbf{n},ij} \,. \tag{4}$$

$$\mathbf{F}\_{\mathbf{t},\mathbf{i}\mathbf{j}} = k\_{\mathbf{t}} \delta\_{\mathbf{t},\mathbf{i}\mathbf{j}} - \lambda\_{\mathbf{t}} \, m\_{\mathbf{eff}} \, \mathbf{u}\_{\mathbf{t},\mathbf{i}\mathbf{j}} \, \mathbf{} \tag{5}$$

where *m*eff = (*mi mj*)/(*mi* + *mj*), subscript *n* denotes the normal direction and *t* the tangential direction, *k* is the spring coefficient, *λ* the damping coefficient, and δ the elastic displacement for a contact, respectively. The damping effect in the tangential direction can be neglected. The tangential displacement δt between particles is obtained by integrating surface relative velocities over time during deformation of the contact. Actually just this history dependent feature makes the computation, especially parallel computation, more expensive. The magnitude of δt is truncated as necessary to satisfy a local Coulomb yield criterion |**F**t| ≤ *μ*|**F**n|.

In the above DEM model, the normal compression between two particles is easily written as:

$$\delta\_{\mathbf{n},\mathbf{i}\mathbf{j}} = [(R\_{\mathbf{i}} + R\_{\mathbf{j}}) - r\_{\mathbf{i}\mathbf{j}}] \mathbf{n}\_{\mathbf{i}\mathbf{j}\prime} \tag{6}$$

here,

4 Will-be-set-by-IN-TECH

In the DEM-SPH coupling, a multiphase DEM-SPH model is proposed and described in detail. The coupling between the solid particles and the fluid phase are considered via volume fraction, pressure and drag force. In mills, particles are so dense that the behaviors of solid particles dominate the overall physics, and separation of solid particles from fluid is apparent. But fluid dynamics cannot be neglected because fluid damping/drag may change particle–particle contacts and particle–boundary contacts. Furthermore, the rotating nature, presence of a free surface flow, and liner geometry make CFD simulations of mills prohibitively expensive. Numerical scheme is simple at this stage, time stepping scheme could be either explicit Euler or lower order Runge-Kutta for more accuracy. The DEM-SPH

The DEM simulation is based on a 3D soft particle model (Cundall & Strack, 1979; Silbert et al., 2001; Walton & Braun, 1986) where small deformations and multiple contacts on a particle are allowed, and friction and rotation are also taken into account. Each particle has six degrees of freedom of motion. For simplification of the description of the DEM model, spherical particles are assumed, although in our real application and implementation, we use spherical, tetrahedral and irregular convex shape particles together. The complex shapes definitely add difficulty to the computational geometric and solving for particle rotation, but the basic theory behind it is the same as for all spherical particles. Although our program is fully parallelized for a distributive memory machine, this work is not intended to cover the parallelization scheme and readers who are interested are referred to Plimpton (1995), which

The movement of a particle with mass *m*, moment of inertia *I* can be described by Newton's

**F**c,*ij* + **F**sf,*i* + *mi* **g** , (1)

*Ri* **n***ij* × **F**c,*ij* , (2)

d *t* = **u***i* , (3)

*mi* d **u***i* d *t* = ∑ *j*

*Ii* d ω*i*

d **r***i*

d *t* = − ∑

particles, and **n***ij* the normal direction of the contact pointing to particle *i* from *j*.

*j*

where subscripts *i* and *j* are for identifying particles, **u** is the velocity of the mass center, ω the angular velocity, **r** the position, *R* the radius of the particle, **F**c,*ij* the contact force of particle *j* acting on *i*, **F**sf,*i* the fluid–solid interaction which is assumed to act at the mass center of

The governing equations are simply ordinary differential equations in time. Each particle is evolved by integrating the governing equations and applying the initial condition. The major task of modeling and simulation thus becomes actually formulating and calculating the force

The implementation of contact forces is essentially a reduced version of that employed by Walton & Braun (1986), developed earlier by Cundall & Strack (1979). Contact force **F**c,*ij* = {**F**n,*ij*, **F**t,*ij*} is first calculated from the deformation through the spring-dashpot model,

is applied to a preliminary study of a mill.

**2. DEM model with flow effect**

is the framework we follow.

terms.

law and the kinematic relation:

$$\mathbf{r}\_{\mathbf{i}\mathbf{j}} = \mathbf{r}\_{\mathbf{i}} - \mathbf{r}\_{\mathbf{j}\prime} \qquad r\_{\mathbf{i}\mathbf{j}} = |\mathbf{r}\_{\mathbf{i}\mathbf{j}}|\,,\tag{7}$$

The normal direction and tangential direction are defined as:

$$\mathbf{n}\_{\mathbf{i}\mathbf{j}} = \mathbf{r}\_{\mathbf{i}\mathbf{j}} / r\_{\mathbf{i}\mathbf{j}} \,\tag{8}$$

$$\mathbf{n}\_{\vec{l}\vec{l}} \times \mathfrak{k}\_{\vec{l}\vec{j}} = 1.\tag{9}$$

The value of the spring constant should be large enough to avoid particle interpenetration, but not too large to require an unreasonably small time step Δ*t*, since an accurate simulation typically requires Δ*t* ∼ *t*c/50, here *t*c is the characteristic contact time during a collision process between particles. The amount of energy lost in collisions is characterized by the inelasticity through the coefficient of restitution *e*. For the linear spring-dashpot model, the following relations can be taken as guidance to find the damping coefficient *λ*n

$$\mathbf{e}\_{\mathbf{n}} = \exp(-\lambda\_{\mathbf{n}} \mathbf{t}\_{\mathbf{c}}/2) \,\,\,\,\,\,\tag{10}$$

$$t\_{\rm c} = \pi (k\_{\rm n}/m\_{\rm eff} - \lambda\_{\rm n}^2/4)^{-1/2} \,. \tag{11}$$

After contact force is calculated, the equations of motion, which are ordinary differential equations, can be numerically integrated to get the particle trajectories. The boundary surfaces are represented by triangles. Any meshing tool generating surface triangles can be used here. In comparison with particle–particle contact, the only difference is the geometry resolution for the particle–triangle contact. The overlap *δ*n is equal to a particle radius minus the distance to a triangle.

The advantage of DEM is that it can capture behaviors of individual particles, and collect detailed contact information, such as velocity, contact force, shear and impact energy spectrum, so that it can model wear more accurately. Most severe wear happens in particulate flows when solid particles collide on a surface. We model the particle wear effect on a boundary triangle as:

$$
\Delta \hbar A = \sum\_{i} C\_{\text{wr},i} E\_{\text{sh},i\prime} \tag{12}
$$

where *C*d is drag coefficient, and the particle Reynolds number Re is defined as

*θ*2.65

terminal velocity of a single particle. The following form for *V*r is used

*A* = *θ*4.14

*θ*s = 1 *V*c

*K*(*ξ*) =

paper to the particle volume fraction calculation. Please refer to that for details.

bandwidth. In this study, the following kernel function is used:

∑*Np*

∑*Np*

where the subscript *c* denotes a cell and *i* denotes a particle. *Np* is the total number of particles in the system. *K* is a kernel function, which should be bell-shape function with local support, such as Gaussian function *K*(*ξ*) = exp[−(*ξ*/*w*)2], where *ξ* = (|**x***i* − **x**c|)/*w* and *w* is the

This function is very close to Gaussian function (with scaled bandwidth *w* = 0.45), and it is more efficient to compute than Gaussian function, because it only requires to loop over the particles in the neighbor cells. Recently, Xiao & Sun (2011) dedicated a big portion of their

Similarly, the solid phase velocity is obtained via the following coarse-graining procedure:

*i*=1 *K*(|**x***i* − **x**c|)*Vi*

*i*=1 *K*(|**x***i* − **x**c|)

1 − *ξ*24 if |*ξ*| < 1

*B* = 0.8 *θ*1.28

*A* − 0.06 Re +

*V*r = 0.5

characteristic length) particles.

through the averaging processes:

Re = ¯*d*p|**u**s − **u**f|*ρ*

Using DEM in Particulate Flow Simulations 35

where ¯*d*p is the average diameter of particles, and *V*r is the ratio of the falling velocity to the

From the above governing equations and drag force calculation, we can see volume fractions appear allover everywhere. This is the core modeling concept for interpenetrating disperse flow, whereas the exact particle shape has not been followed, rather the volume/mass conservation must be maintained. This multiphase modeling approach is a perfect compromise for large scale particulate flows with a huge number of small (relative to flow

The volume fraction and velocity of the solid phase are needed for each cell. It sounds simple at the first glance. It may be simple for regular mesh, not for irregular mesh and irregular particles. As we mentioned, particle shape is not of importance at modeling particle fluid interaction, nor should be mesh geometry. We compromise the geometry details for the efficiency of numerical calculation. Many researchers use kernel function to calculate the volume fraction to avoid handling complex mesh and particle geometry. The volume fraction of each cell is obtained from the particle spatial distribution and the volume of each particle

(0.06 Re)2 + 0.12 Re(2B − A) + A2

f if *θ*f ≤ 0.85

f (21)

f if *θ*f > 0.85. (22)

0 if |*ξ*| ≥ 1 (24)

*μ* , (19)

, (20)

(23)

where subscript *i* denotes the particle groups, Δ*h* is the thickness to be worn off the surface, *A* is the triangle area, *E*sh,*i* is the cumulative shear energy over the time period of interest, and *C*wr,*i* is the wear coefficient that is a function of the triangle and particle material properties and the particle size.

Similarly, we model the particle breakage as

$$B\mathbf{r}\_{\rm i} = \mathbb{C}\_{\rm br,i} E\_{\rm imp,i} / m\_{\rm i} \tag{13}$$

where *B*r*i* is the breakage percentage (probability) of the group *i* of particles when passing through a process, *E*imp,*i* is the cumulative impact energy over the retention time of the process, and *C*br,*i* is the breakage coefficient that is a function of the material properties and the particle size. The constants *C*wr,*i* and *C*br,*i* are to be determined from the calibration with experiment or operational data in this work. See the reference for details and recent developments about the wear model (Hollow & Herbst, 2006; Qiu et al., 2001) and particle breakage model (Herbst & Potapov, 2004; Potapov et al., 2007).

#### **3. CFD coupling with DEM**

The basic conceptual theory for CFD coupling with DEM comes from the Eulerian–Eulerian multi-fluid model, where the fluid phase and the solid particle mixture are described as interpenetrating continua. The particle mixture can be divided into a discrete number of phases, each of which can have different physical properties. Generally, *n* sets of governing equations have to be solved for a multiphase flow with *n* phases, and an exponentially increasing number of constitutive equations are required for closure of the model.

The approach of CFD coupling with DEM is proposed to overcome the closure difficulty. The Eulerian control-volume multiphase CFD governing equation is used to describe the fluid dynamics, while the DEM is used to model the solid mixture dynamic behaviors.

The governing equations for incompressible flow, continuity and momentum equations, for the fluid phase, are:

$$\frac{\partial}{\partial t}(\rho \theta\_\mathbf{f}) + \nabla \cdot \rho \theta\_\mathbf{f} \mathbf{u}\_\mathbf{f} = \mathbf{0},\tag{14}$$

$$\frac{\partial}{\partial t}(\rho \theta\_\mathbf{f})\mathbf{u}\_\mathbf{f} + \nabla \cdot (\rho \theta\_\mathbf{f} \mathbf{u}\_\mathbf{f} \mathbf{u}\_\mathbf{f}) = -\rho \nabla P + \mathbf{F}\_\mathbf{d} + \nabla \cdot \mathbf{T} + \rho \theta\_\mathbf{f} \mathbf{g}\_\mathbf{f} \tag{15}$$

where *ρ* is the fluid density, **u**f the fluid velocity, *θ*f the fluid volume fraction, and *P* the pressure, and *F*d the drag force, and **T** the viscous stress tensor.

The drag force **F**d between particles and fluid phase is generally defined as

$$\mathbf{F\_d} = \beta (\mathbf{u\_s} - \mathbf{u\_f}) \,\tag{16}$$

where *β* is the drag force coefficient, **u**s the velocity vector of solid particles, and (**u**s − **u**f) is the slip velocity between the two phases. The volume fraction *θ*s and velocity fields **u**s of the solid phase are obtained through averaging particle data from DEM simulation, thus *θ*f can be obtained from *θ*f = 1 − *θ*s.

Without losing generality, in the DEM-CFD coupling, the drag correlation by Syamlal et al. (1993) is adopted:

$$\beta = \frac{3}{4} \frac{\mathbf{C\_d}}{V\_r^2} \frac{\rho |\mathbf{u\_s} - \mathbf{u\_f}|}{d\_p} \theta\_\mathbf{s} \theta\_\mathbf{f\_f} \tag{17}$$

$$\mathcal{C}\_{\rm d} = \left( 0.63 + 4.8\sqrt{V\_{\rm r}/\text{Re}} \right)^2,\tag{18}$$

6 Will-be-set-by-IN-TECH

where subscript *i* denotes the particle groups, Δ*h* is the thickness to be worn off the surface, *A* is the triangle area, *E*sh,*i* is the cumulative shear energy over the time period of interest, and *C*wr,*i* is the wear coefficient that is a function of the triangle and particle material properties

where *B*r*i* is the breakage percentage (probability) of the group *i* of particles when passing through a process, *E*imp,*i* is the cumulative impact energy over the retention time of the process, and *C*br,*i* is the breakage coefficient that is a function of the material properties and the particle size. The constants *C*wr,*i* and *C*br,*i* are to be determined from the calibration with experiment or operational data in this work. See the reference for details and recent developments about the wear model (Hollow & Herbst, 2006; Qiu et al., 2001) and particle

The basic conceptual theory for CFD coupling with DEM comes from the Eulerian–Eulerian multi-fluid model, where the fluid phase and the solid particle mixture are described as interpenetrating continua. The particle mixture can be divided into a discrete number of phases, each of which can have different physical properties. Generally, *n* sets of governing equations have to be solved for a multiphase flow with *n* phases, and an exponentially

The approach of CFD coupling with DEM is proposed to overcome the closure difficulty. The Eulerian control-volume multiphase CFD governing equation is used to describe the fluid

The governing equations for incompressible flow, continuity and momentum equations, for

where *ρ* is the fluid density, **u**f the fluid velocity, *θ*f the fluid volume fraction, and *P* the

where *β* is the drag force coefficient, **u**s the velocity vector of solid particles, and (**u**s − **u**f) is the slip velocity between the two phases. The volume fraction *θ*s and velocity fields **u**s of the solid phase are obtained through averaging particle data from DEM simulation, thus *θ*f can

Without losing generality, in the DEM-CFD coupling, the drag correlation by Syamlal et al.

*ρ*|**u**s − **u**f| *dp*

0.63 + 4.8*V*r/Re

2

(*ρθ*f) + ∇ · *ρθ*f**u**f = 0 , (14)

(*ρθ*f)**u**f + ∇ · (*ρθ*f**u**f**u**f) = −*ρ*∇ *P* + **F**d + ∇ · **T** + *ρθ*f**g** , (15)

**F**d = *β*(**u**s − **u**f), (16)

*θ*s*θ*f , (17)

, (18)

increasing number of constitutive equations are required for closure of the model.

dynamics, while the DEM is used to model the solid mixture dynamic behaviors.

*B*r*i* = *C*br,*iE*imp,*i*/*mi*, (13)

and the particle size.

**3. CFD coupling with DEM**

the fluid phase, are:

*∂ ∂t*

*∂ ∂t*

be obtained from *θ*f = 1 − *θ*s.

(1993) is adopted:

Similarly, we model the particle breakage as

breakage model (Herbst & Potapov, 2004; Potapov et al., 2007).

pressure, and *F*d the drag force, and **T** the viscous stress tensor.

The drag force **F**d between particles and fluid phase is generally defined as

*β* = 3 4 *C*d *V*2 *r*

*C*d =  where *C*d is drag coefficient, and the particle Reynolds number Re is defined as

$$\text{Re} = \frac{\vec{d}\_{\text{p}} |\mathbf{u}\_{\text{s}} - \mathbf{u}\_{\text{f}}| \rho}{\mu} \,\tag{19}$$

where ¯*d*p is the average diameter of particles, and *V*r is the ratio of the falling velocity to the terminal velocity of a single particle. The following form for *V*r is used

$$V\_r = 0.5 \left( A - 0.06 \,\text{Re} + \sqrt{(0.06 \,\text{Re})^2 + 0.12 \,\text{Re}(2\text{B} - \text{A}) + \text{A}^2} \right), \tag{20}$$

$$A = \theta\_{\rm f}^{4.14} \tag{21}$$

$$B = \begin{cases} 0.8 \,\theta\_{\text{f}}^{1.28} & \text{if} \quad \theta\_{\text{f}} \le 0.85\\ \theta\_{\text{f}}^{2.65} & \text{if} \quad \theta\_{\text{f}} > 0.85. \end{cases} \tag{22}$$

From the above governing equations and drag force calculation, we can see volume fractions appear allover everywhere. This is the core modeling concept for interpenetrating disperse flow, whereas the exact particle shape has not been followed, rather the volume/mass conservation must be maintained. This multiphase modeling approach is a perfect compromise for large scale particulate flows with a huge number of small (relative to flow characteristic length) particles.

The volume fraction and velocity of the solid phase are needed for each cell. It sounds simple at the first glance. It may be simple for regular mesh, not for irregular mesh and irregular particles. As we mentioned, particle shape is not of importance at modeling particle fluid interaction, nor should be mesh geometry. We compromise the geometry details for the efficiency of numerical calculation. Many researchers use kernel function to calculate the volume fraction to avoid handling complex mesh and particle geometry. The volume fraction of each cell is obtained from the particle spatial distribution and the volume of each particle through the averaging processes:

$$\theta\_{\mathbf{s}} = \frac{1}{V\_{\mathbf{c}}} \frac{\sum\_{i=1}^{N\_p} K(|\mathbf{x}\_i - \mathbf{x}\_{\mathbf{c}}|) V\_i}{\sum\_{i=1}^{N\_p} K(|\mathbf{x}\_i - \mathbf{x}\_{\mathbf{c}}|)} \tag{23}$$

where the subscript *c* denotes a cell and *i* denotes a particle. *Np* is the total number of particles in the system. *K* is a kernel function, which should be bell-shape function with local support, such as Gaussian function *K*(*ξ*) = exp[−(*ξ*/*w*)2], where *ξ* = (|**x***i* − **x**c|)/*w* and *w* is the bandwidth. In this study, the following kernel function is used:

$$K(\xi) = \begin{cases} \left[1 - \xi^2\right]^4 & \text{if } |\xi| < 1\\ 0 & \text{if } |\xi| \ge 1 \end{cases} \tag{24}$$

This function is very close to Gaussian function (with scaled bandwidth *w* = 0.45), and it is more efficient to compute than Gaussian function, because it only requires to loop over the particles in the neighbor cells. Recently, Xiao & Sun (2011) dedicated a big portion of their paper to the particle volume fraction calculation. Please refer to that for details.

Similarly, the solid phase velocity is obtained via the following coarse-graining procedure:

**4. SPH coupling with DEM**

*θ* and actual fluid density *ρ* by

values of *γ* and *c*0. The speed of sound is:

SPH equations are obtained by interpolating fluid dynamics governing equations over disordered mass points in the influence range of an interpolation kernel function (Monaghan, 1988). The kernels are analytical functions which can be differentiated without using a mesh. Although control volume CFD can be tuned to get accurate solutions of physical problems, it requires tremendous work, including generating the mesh, in order to couple with DEM to account for the multiphase flows. On the other hand, SPH is a method which gives reasonable

Using DEM in Particulate Flow Simulations 37

Monaghan & Kocharyan (1995) originally built SPH multiphase models for interpenetrating multi-fluids. In the later work, Monaghan (1997) improved the multiphase SPH solver by using an implicit drag technique. In this work, we will modify the interpenetrating multiphase SPH model to couple SPH for fluid dynamics with DEM for solid particles. SPH

d *t* = −*ρ*ˆ ∇ · **u** , (28)

*ρ*ˆ + **g** , (29)

(31)

*ρ*ˆ = *ρ θ* . (30)

0 (*ρ* − *ρ*0), (32)

*ρ* (33)

(**u**s − **u**) + ∇ · **T**

accuracy and couples well with the particle method DEM without requiring a mesh.

is only used to model fluid phase, and particles are represented and evolved by DEM The governing equations, continuity and momentum equations, for the fluid phase are:

> + *β ρ*ˆ

where *β* is the drag force coefficient, **u**s the velocity vector of solid particles and **T** the viscous stress tensor. The fluid density *ρ*ˆ in the multiphase model is related to the fluid volume fraction

The equation of state has to be defined in order to fully describe the dynamics of the fluid. The actual equation of state of incompressible flow is very stiff, requiring extremely small time steps. In SPH, the fluid pressure is an explicit function of local fluid density. Therefore, it is necessary to use a quasi-compressible equation of state and an artificial speed of sound as a reference value. Monaghan (2000) used the equation of state similar to that defined for water:

> *ρ ρ*0

where *ρ*0 is the reference density, *c*0 the speed of sound at *ρ*0, and *γ* a constant with physical meaning of the ratio of specific heat for ideal gases. Monaghan (2000) took *γ* = 7 for incompressible flows. The choice *γ* = 7 results in large changes in pressure from small changes or perturbations in density. We confirmed that Equation 31 works well for single phase flow SPH. However, in this work considering multiphase free surface flow, the density changes could also be exaggerated by changes or errors in calculation of volume fraction (see below). We also follow suggestions of Morris et al. (1997) regarding choosing appropriate

> *c* = *γP*

*γ* − 1

*P* = *ρ*0*c*2 0 *γ*

≈ *c*2

d *ρ*ˆ

d **u**

d *t* = −∇ *P ρ*

$$\mathbf{u}\_{\mathbf{s}} = \frac{\sum\_{i=1}^{N\_p} K(|\mathbf{x}\_i - \mathbf{x}\_{\mathbf{c}}|) V\_i \mathbf{u}\_i}{\sum\_{i=1}^{N\_p} K(|\mathbf{x}\_i - \mathbf{x}\_{\mathbf{c}}|) V\_i} \tag{25}$$

On the DEM side, two major solid–fluid interaction terms acts on a particle. In Equation 1 force *F*sf,*i* equals the sum of buoyancy force and drag force acting on particle *i*:

$$\mathbf{F\_{sf,i}} = -V\_{\mathbf{i}} \,\nabla P - F\_{\mathbf{d,i}\,\mathbf{i}}\,\tag{26}$$

The drag of Equation 16 is the aggregate drag force acting on all the particles in a cell. For the motion equation 1 of a solid particle, the drag force on each particle is needed. There are many ways of redistributing the aggregate force back to individual particles. This can be obtained from the cell whose center is nearest to the particle, or from the neighboring cells whose contribution portion is determined by distances between particles and cell centers. Each particle shares the drag force proportional to its surface area, or proportional to its volume, or as functions of other properties. The uncertainty also occurs with buoyancy force. There is no sure answer to the question of where pressure values are used in calculation of buoyancy force on a particle. Here in the application of DEM-CFD full coupling, we take the simple scheme by distributing the nearest cell force according to surface areas of particles. For the one-way coupling applications, particles are small and well mixed with fluid phase, so that one-way coupling is justified, and Wen-Yu drag relation (Li & Kuipers, 2002; Rong &

$$\mathbf{F}\_{\rm sf,i} = -V\_{\rm i} \nabla P - \frac{3}{4} \frac{\mathbf{C}\_d |\mathbf{u}\_i - \mathbf{u}\_{\rm f}| \rho \theta\_{\rm f}^{-1.65}}{D\_{\rm i}} \left(\mathbf{u}\_i - \mathbf{u}\_{\rm f}\right) \tag{27}$$

where *Vi* is the volume, and *Di* the hydraulic diameter of particle *i*. For the drag coefficient *C*d, please refer to DEM–SPH coupling section.

Horio, 1999) is used. Thus, the drag force calculation can be as simple as:

For DEM-CFD coupling, the fluid equations are solved using a solver provided by a library in the OpenFOAM (Open Field Operation and Manipulation) toolbox (Rusche, 2002). This solver is used in the current study with modifications to accommodate the fact that only the fluid phase is solved and the disperse phase is tracked in the Lagrangian framework. Finite volume method is used to discretize the equations on an unstructured mesh. For the time integration, Euler implicit scheme is used, which has only first order accuracy but is unconditionally stable. The convection and diffusion terms are discretized with a blender of central differencing (second-order accurate) and upwind differencing (first order accurate). The advantage of blended differencing is that high accuracy is achieved while the boundness of the solution is ensured. A sophisticate stepping control and interpolation over time was brought up by (Xiao & Sun, 2011) to enhance accuracy and convergence of the DEM–CFD coupling solver.

The velocity-pressure coupling is handled with the modified PISO algorithm (Rhie & Chow, 1983). In this algorithm, momentum equation is first solved to get a predicted velocity field, and then the pressure equation (obtained by combining momentum equation into continuity equation) is solved for a corrected velocity field. This process is repeated until the velocity field satisfies continuity equation. The PISO algorithm prevents from the decoupling of pressure-velocity and the oscillation in the solution, eliminating the necessity of a stagger grid. Therefore, a collocated grid is used in the models, where all the variables are stored in the cell centers, thus it is a significant simplification over a stagger grid.

#### **4. SPH coupling with DEM**

8 Will-be-set-by-IN-TECH

On the DEM side, two major solid–fluid interaction terms acts on a particle. In Equation 1

The drag of Equation 16 is the aggregate drag force acting on all the particles in a cell. For the motion equation 1 of a solid particle, the drag force on each particle is needed. There are many ways of redistributing the aggregate force back to individual particles. This can be obtained from the cell whose center is nearest to the particle, or from the neighboring cells whose contribution portion is determined by distances between particles and cell centers. Each particle shares the drag force proportional to its surface area, or proportional to its volume, or as functions of other properties. The uncertainty also occurs with buoyancy force. There is no sure answer to the question of where pressure values are used in calculation of buoyancy force on a particle. Here in the application of DEM-CFD full coupling, we take the simple scheme by distributing the nearest cell force according to surface areas of particles. For the one-way coupling applications, particles are small and well mixed with fluid phase, so that one-way coupling is justified, and Wen-Yu drag relation (Li & Kuipers, 2002; Rong &

*i*=1 *K*(|**x***i* − **x**c|)*Vi***u***i*

(25)

*i*=1 *K*(|**x***i* − **x**c|)*Vi*

*Cd*|**u***i* − **u**f|*ρθ*−1.65

where *Vi* is the volume, and *Di* the hydraulic diameter of particle *i*. For the drag coefficient

For DEM-CFD coupling, the fluid equations are solved using a solver provided by a library in the OpenFOAM (Open Field Operation and Manipulation) toolbox (Rusche, 2002). This solver is used in the current study with modifications to accommodate the fact that only the fluid phase is solved and the disperse phase is tracked in the Lagrangian framework. Finite volume method is used to discretize the equations on an unstructured mesh. For the time integration, Euler implicit scheme is used, which has only first order accuracy but is unconditionally stable. The convection and diffusion terms are discretized with a blender of central differencing (second-order accurate) and upwind differencing (first order accurate). The advantage of blended differencing is that high accuracy is achieved while the boundness of the solution is ensured. A sophisticate stepping control and interpolation over time was brought up by (Xiao & Sun, 2011) to enhance accuracy and convergence of the DEM–CFD

The velocity-pressure coupling is handled with the modified PISO algorithm (Rhie & Chow, 1983). In this algorithm, momentum equation is first solved to get a predicted velocity field, and then the pressure equation (obtained by combining momentum equation into continuity equation) is solved for a corrected velocity field. This process is repeated until the velocity field satisfies continuity equation. The PISO algorithm prevents from the decoupling of pressure-velocity and the oscillation in the solution, eliminating the necessity of a stagger grid. Therefore, a collocated grid is used in the models, where all the variables are stored in

the cell centers, thus it is a significant simplification over a stagger grid.

*Di*

f *Vi*

(**u***i* − **u**f), (27)

**F**sf,*i* = −*Vi* ∇*P* − *F*d,i , (26)

**u**s = ∑*Np*

force *F*sf,*i* equals the sum of buoyancy force and drag force acting on particle *i*:

Horio, 1999) is used. Thus, the drag force calculation can be as simple as:

4

**F**sf,*i* = −*Vi* ∇*P* − 3

*C*d, please refer to DEM–SPH coupling section.

coupling solver.

∑*Np*

SPH equations are obtained by interpolating fluid dynamics governing equations over disordered mass points in the influence range of an interpolation kernel function (Monaghan, 1988). The kernels are analytical functions which can be differentiated without using a mesh. Although control volume CFD can be tuned to get accurate solutions of physical problems, it requires tremendous work, including generating the mesh, in order to couple with DEM to account for the multiphase flows. On the other hand, SPH is a method which gives reasonable accuracy and couples well with the particle method DEM without requiring a mesh.

Monaghan & Kocharyan (1995) originally built SPH multiphase models for interpenetrating multi-fluids. In the later work, Monaghan (1997) improved the multiphase SPH solver by using an implicit drag technique. In this work, we will modify the interpenetrating multiphase SPH model to couple SPH for fluid dynamics with DEM for solid particles. SPH is only used to model fluid phase, and particles are represented and evolved by DEM The governing equations, continuity and momentum equations, for the fluid phase are:

$$\frac{\mathbf{d}\,\rho}{\mathbf{d}\,t} = -\rho \,\nabla \cdot \mathbf{u}\_{\prime} \tag{28}$$

$$\frac{\mathbf{d}\,\mathbf{u}}{\mathbf{d}\,t} = -\frac{\nabla\,P}{\rho} + \frac{\beta}{\beta}(\mathbf{u}\_s - \mathbf{u}) + \frac{\nabla\cdot\mathbf{T}}{\rho} + \mathbf{g}\_{\prime} \tag{29}$$

where *β* is the drag force coefficient, **u**s the velocity vector of solid particles and **T** the viscous stress tensor. The fluid density *ρ*ˆ in the multiphase model is related to the fluid volume fraction *θ* and actual fluid density *ρ* by

$$
\boldsymbol{\mathfrak{\rho}} = \boldsymbol{\mathfrak{\rho}} \,\,\boldsymbol{\mathfrak{\theta}}\,.\,\tag{30}
$$

The equation of state has to be defined in order to fully describe the dynamics of the fluid. The actual equation of state of incompressible flow is very stiff, requiring extremely small time steps. In SPH, the fluid pressure is an explicit function of local fluid density. Therefore, it is necessary to use a quasi-compressible equation of state and an artificial speed of sound as a reference value. Monaghan (2000) used the equation of state similar to that defined for water:

$$P = \frac{\rho\_0 c\_0^2}{\gamma} \left[ \left( \frac{\rho}{\rho\_0} \right)^{\gamma} - 1 \right] \tag{31}$$

$$\approx c\_0^2 \left(\rho - \rho\_0\right),\tag{32}$$

where *ρ*0 is the reference density, *c*0 the speed of sound at *ρ*0, and *γ* a constant with physical meaning of the ratio of specific heat for ideal gases. Monaghan (2000) took *γ* = 7 for incompressible flows. The choice *γ* = 7 results in large changes in pressure from small changes or perturbations in density. We confirmed that Equation 31 works well for single phase flow SPH. However, in this work considering multiphase free surface flow, the density changes could also be exaggerated by changes or errors in calculation of volume fraction (see below). We also follow suggestions of Morris et al. (1997) regarding choosing appropriate values of *γ* and *c*0. The speed of sound is:

$$
\sigma = \sqrt{\frac{\gamma P}{\rho}} \tag{33}
$$

The kernel function is the commonly used cubic spline function

*πh*3

⎧ ⎪⎪⎪⎨

1 − 3 2 *q*2 + 3 4

The Gidaspow drag correlation, which combines the Wen-Yu relation and the Ergun equation, is commonly used in CFD multiphase modeling (Gera et al., 1998; 2004; Li & Kuipers, 2002;

Using DEM in Particulate Flow Simulations 39

1.75|**u***a* − **u***j*| *θaDj*

1 4 *q*3 if 0 ≤ *q* < 1

*ρaθaθj θj* ≤ 0.2

*ρaθaθj θj* > 0.2

*aj* (42)

�

(38)

(39)

(40)

(41)

(2 − *q*)3 if 1 ≤ *q* < 2 0 otherwise

⎪⎪⎪⎩

|**u***a* − **u***j*|*C*d*θ*−2.65 *a Dj*

+

0.44 Re ≥ 1000

Re (1 + 0.15*Re*0.687) 10−4 < *Re* < 1000

where *Dj* is solid particle effective hydraulic diameter. The *ρa*, *θa* and *θj* terms are intentionally grouped together as factors, because these terms will be canceled when *βaj* is substituted back

From the comparison of governing equations of DEM–CFD to those of DEM–SPH, one apparent difference is that DEM–CFD calculates aggregate coupling force first, and distributes back to particles, while the DEM–SPH calculates coupling force for individual particles and collects the total force on SPH particle. Due to the high cost the drag force calculation, the cost of DEM–SPH may be much higher than DEM–CFD. Further modification to the DEM–SPH

No matter which way of coupling, we have to obtain the collective volume fractions *θj* and *θa*. Using kernel function for volume fraction calculation is the most natural way in the SPH framework. Please refer to Monaghan (1997); Monaghan & Kocharyan (1995) for more theory support. The most simple form of fluid fraction in the SPH type of kernel functions can be

*j*

where *W*∗ is a kernel function that can be different from *W*, and a different smoothing length should be used. The smoothing length should be at least twice the maximum particle size. After that, *θj* can be obtained by smoothing over SPH particles and using the same smoothing

One can use the same way as in the above DEM–CFD, Equation 23, for volume fraction calculation, however, using a imaginary sphere as *V*c, because there is no real mesh in SPH. The imaginary sphere radius should be at least twice the maximum particle size. In this work,

*VjW*∗

*θa* = 1 − ∑

150(1 − *θa*)*μa θ*2 *aρaD*2 *j*

*W*(*r*, *h*) = 1

where *q* = *r*/*h* and *r* is the distance between particles.

Rong & Horio, 1999) and it is used here:

⎧ ⎪⎪⎪⎨

3 4

�

⎪⎪⎪⎩

� 24

Re = *ρaθa*|**u***a* − **u***j*|*Dj μa*

by adopting DEM–CFD approach is the future work of authors.

*βaj* =

*C*d =

to Equations 36 and 37.

defined as follows:

length as in SPH equations.

SPH is based on the theory of integral interpolants. If the kernel functions are some types of delta functions, then a field variable *A*(**r**) can be approximated by the weighted averaging over a limited range of neighboring particles as:

$$A(\mathbf{r}) = \int A(\mathbf{r'}) \mathcal{W}(\mathbf{r} - \mathbf{r'}, h) \mathbf{d} \, \mathbf{r'} \simeq \sum\_{a} A\_{a} \frac{m\_{a}}{\rho\_{a}} \mathcal{W}(\mathbf{r} - \mathbf{r}\_{a}, h) \, \, \, \tag{34}$$

where *ma* denotes mass of SPH particle *a* at the position **r***a*, and similar notations for *ρa*, *Aa*. *W*(**r** − **r***a*, *h*) is the kernel which is a function of smoothing length *h* and distance between positions **r** and **r***a*.

For clarification of the description, the subscripts *a* and *b* are used for the SPH fluid particles, and *i* and *j* for the DEM solid particles. Integrating the governing equations and simplifying the integrals as above, we can get the overall SPH governing equations:

$$\frac{\mathbf{d}\,\hat{\boldsymbol{\theta}}}{\mathbf{d}\,\mathbf{t}} = \sum\_{b} m\_{b} \mathbf{u}\_{ab} \cdot \nabla\_{a} W\_{ab} \tag{35}$$

$$\begin{split} \frac{\mathbf{d}\,\mathbf{u}\_{a}}{\mathbf{d}\,t} &= -\sum\_{b} m\_{b} \left( \frac{P\_{a}\theta\_{a}}{\hat{\rho}\_{a}^{2}} + \frac{P\_{b}\theta\_{b}}{\hat{\rho}\_{b}^{2}} \right) \nabla\_{a} W\_{ab} \\ &- \sum\_{j} \frac{P\_{a}V\_{j}}{\hat{\rho}\_{a}} \nabla\_{a} W\_{aj} \\ &+ \frac{1}{3} \sum\_{j} \frac{\beta\_{aj}V\_{j}}{\hat{\rho}\_{a}\theta\_{j}} \left( \frac{\mathbf{u}\_{ja} \cdot \mathbf{r}\_{ja}}{r\_{ja}} \right) W\_{ja} \frac{\mathbf{r}\_{ja}}{r\_{ja}} \\ &+ \mathbf{r}\_{ab} + \mathbf{g}\_{a} \end{split} \tag{36}$$

where *Vj* is the volume of DEM particle *j*, *θj* the solid phase volume fraction at the position of DEM particle *j*, *θa* the fluid phase volume fraction at the position of SPH particle *a*, and τ*ab* the viscous term. The solid–fluid inter-phase interaction force is represented by the second and third group of terms on the right hand side of Equation 36. The second term is the pressure gradient on a solid particle, and third term is the drag force. The drag force between two particles acts along their center line, working like dashpot damping in the DEM models. Thus the inter-phase interaction in Equation 1 can be written as:

$$\mathbf{F}\_{\rm sf,j} = -\sum\_{a} m\_{a} \frac{P\_{a}}{\hat{\rho}\_{a}} \nabla\_{\hat{\jmath}} \mathcal{W}\_{\rm ja} - \frac{1}{3} \sum\_{a} m\_{a} \frac{\mathcal{\beta}\_{a\dot{\jmath}} V\_{\dot{\jmath}}}{\hat{\rho}\_{a} \theta\_{\dot{\jmath}}} \left(\frac{\mathbf{u}\_{ja} \cdot \mathbf{r}\_{ja}}{r\_{ja}}\right) \mathcal{W}\_{\rm ja} \frac{\mathbf{r}\_{ja}}{r\_{ja}}.\tag{37}$$

In the above equations, we have used the notation

$$\mathbf{u}\_{aj} = \mathbf{u}\_a - \mathbf{u}\_{j'i'}$$

for a vector, and

$$\begin{aligned} \mathcal{W}\_{ab} &= \mathcal{W}(r\_{ab}) \\ \nabla\_a \mathcal{W}\_{ab} &= \frac{\mathbf{r}\_{ab}}{r\_{ab}} \frac{\mathbf{d} \, \mathcal{W}}{\mathbf{d} \, r}(r\_{ab}) = \frac{\mathbf{r}\_{ab}}{r\_{ab}} \left(\frac{\mathbf{d} \, \mathcal{W}}{\mathbf{d} \, r}\right)\_{ab} \end{aligned}$$

for a kernel and a kernel derivative, where *rab* is the distance between particle *a* and *b*. Similar notation is used for other terms.

10 Will-be-set-by-IN-TECH

SPH is based on the theory of integral interpolants. If the kernel functions are some types of delta functions, then a field variable *A*(**r**) can be approximated by the weighted averaging

where *ma* denotes mass of SPH particle *a* at the position **r***a*, and similar notations for *ρa*, *Aa*. *W*(**r** − **r***a*, *h*) is the kernel which is a function of smoothing length *h* and distance between

For clarification of the description, the subscripts *a* and *b* are used for the SPH fluid particles, and *i* and *j* for the DEM solid particles. Integrating the governing equations and simplifying

� � ∑*a*

*Aa ma ρa*

*mb***u***ab* · ∇*aWab* (35)

∇*aWab*

**u***ja* · **r***ja rja*

 *Wja* **r***ja rja*

. (37)

 *Wja* **r***ja rja* +τ*ab* + **g***a* , (36)

*W*(**r** − **r***a*, *h*), (34)

over a limited range of neighboring particles as:

 *A*(**r** �

)*W*(**r** − **r** � , *h*)d **r**

the integrals as above, we can get the overall SPH governing equations:

d *ρ*ˆ d *t* = ∑ *b*

d **u***a*

the inter-phase interaction in Equation 1 can be written as:

*ma Pa Vj ρ*ˆ*a*

In the above equations, we have used the notation

**F**sf,*j* = − ∑*a*

for a vector, and

notation is used for other terms.

d *t* = − ∑

*b mb Paθa ρ*ˆ2 *a* + *Pbθb ρ*ˆ2 *b*

> *Pa Vj ρ*ˆ*a*

> > *βajVj ρ*ˆ*aθj*

where *Vj* is the volume of DEM particle *j*, *θj* the solid phase volume fraction at the position of DEM particle *j*, *θa* the fluid phase volume fraction at the position of SPH particle *a*, and τ*ab* the viscous term. The solid–fluid inter-phase interaction force is represented by the second and third group of terms on the right hand side of Equation 36. The second term is the pressure gradient on a solid particle, and third term is the drag force. The drag force between two particles acts along their center line, working like dashpot damping in the DEM models. Thus

> 3 ∑*a ma βajVj ρ*ˆ*aθj*

**u***aj* = **u***a* − **u***j* ,

d *r* (*rab*) = **r***ab*

for a kernel and a kernel derivative, where *rab* is the distance between particle *a* and *b*. Similar

*rab*

d *W* d *r*

 *ab* ,

∇*jWja* − 1

*Wab* = *W*(*rab*),

*rab*

d *W*

∇*aWab* = **r***ab*

∇*aWaj*

**u***ja* · **r***ja rja*

− ∑ *j*

+ 1 3 ∑ *j*

*A*(**r**) =

positions **r** and **r***a*.

The kernel function is the commonly used cubic spline function

$$\mathcal{W}(r,h) = \frac{1}{\pi h^3} \begin{cases} 1 - \frac{3}{2}q^2 + \frac{3}{4}q^3 & \text{if} \quad 0 \le q < 1\\ \frac{1}{4}(2-q)^3 & \text{if} \quad 1 \le q < 2\\ 0 & \text{otherwise} \end{cases} \tag{38}$$

where *q* = *r*/*h* and *r* is the distance between particles.

The Gidaspow drag correlation, which combines the Wen-Yu relation and the Ergun equation, is commonly used in CFD multiphase modeling (Gera et al., 1998; 2004; Li & Kuipers, 2002; Rong & Horio, 1999) and it is used here:

$$\beta\_{a\dot{j}} = \begin{cases} \frac{3}{4} \frac{|\mathbf{u}\_{a} - \mathbf{u}\_{\dot{j}}| \mathbb{C}\_{\rm d} \theta\_{a}^{-2.65}}{D\_{\dot{j}}} \rho\_{a} \theta\_{a} \theta\_{\dot{j}} & \theta\_{\dot{j}} \le 0.2\\ \left[ \frac{150(1 - \theta\_{a}) \mu\_{a}}{1 - \theta\_{\dot{j}}} + \frac{1.75|\mathbf{u}\_{a} - \mathbf{u}\_{\dot{j}}| }{} \right] \rho\_{a} \theta\_{\dot{a}} \theta\_{\dot{i}} & \theta\_{\dot{i}} > 0.2 \end{cases} \tag{39}$$

$$\left\{ \begin{array}{ll} \left\lfloor \frac{150^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\cdots7}\right\}}\right\}}}}}}}}+\cdots+\left\{\frac{150^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\cdots7}\right\}\right\}}\right\}}}}}\right\}}\right\}}\\\theta\right\}}\\\theta\right)^{\left\{\frac{150^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\left\{120^{\cdots7}\right\}\right\}}\right\}}}}}}}\right\}}\right\}}\right\}}\right\}}\right\}}\right}\\}\right}\\$$

$$\mathbf{C\_d} = \begin{cases} \frac{\Delta \mathbf{r}}{\text{Re}} (1 + 0.15 \text{Re}^{0.687}) & 10^{-4} < \text{Re} < 1000\\ 0.44 & \text{Re} \ge 1000 \end{cases} \tag{40}$$

$$\text{Re} = \frac{\rho\_a \theta\_a |\mathbf{u}\_a - \mathbf{u}\_j| D\_j}{\mu\_a} \tag{41}$$

where *Dj* is solid particle effective hydraulic diameter. The *ρa*, *θa* and *θj* terms are intentionally grouped together as factors, because these terms will be canceled when *βaj* is substituted back to Equations 36 and 37.

From the comparison of governing equations of DEM–CFD to those of DEM–SPH, one apparent difference is that DEM–CFD calculates aggregate coupling force first, and distributes back to particles, while the DEM–SPH calculates coupling force for individual particles and collects the total force on SPH particle. Due to the high cost the drag force calculation, the cost of DEM–SPH may be much higher than DEM–CFD. Further modification to the DEM–SPH by adopting DEM–CFD approach is the future work of authors.

No matter which way of coupling, we have to obtain the collective volume fractions *θj* and *θa*. Using kernel function for volume fraction calculation is the most natural way in the SPH framework. Please refer to Monaghan (1997); Monaghan & Kocharyan (1995) for more theory support. The most simple form of fluid fraction in the SPH type of kernel functions can be defined as follows:

$$\theta\_a = 1 - \sum\_j V\_j W\_{aj}^\* \tag{42}$$

where *W*∗ is a kernel function that can be different from *W*, and a different smoothing length should be used. The smoothing length should be at least twice the maximum particle size. After that, *θj* can be obtained by smoothing over SPH particles and using the same smoothing length as in SPH equations.

One can use the same way as in the above DEM–CFD, Equation 23, for volume fraction calculation, however, using a imaginary sphere as *V*c, because there is no real mesh in SPH. The imaginary sphere radius should be at least twice the maximum particle size. In this work,

values at *t*

order accuracy.

**5. Applications**

parameters for the same size of DEM particle.

**5.1 One-way coupling for a pump**

*n*+1 can be obtained using the trapezoidal rule. This method can achieve second

Boundary conditions for SPH particles are similar to that for solid DEM particles: spring-dashpot contact models are applied to the particles assuming spherical particles with diameter equal to the initial separation distance. The contact coefficients are taken from DEM

Using DEM in Particulate Flow Simulations 41

For slurry in pump operation, solid particles are well mixed with fluid, the particle volume fraction is usually not very high, velocity is high but streamlines do not interweave together, so that the particle–particle collision probability is low, thus the particle–boundary collision is the main interest here since particle–boundary interaction is the major source of wear and the reason for particle breakage under these conditions. DEM simulations were also performed to estimate the breakage percentage and liberation percentage of the particles during the operation of the pump at different conditions. The commercial software FLUENT is used

The CFD simulation of a pump is not easy and straightforward like pipe flow. First of all, for high fidelity we use the manufacturer's geometry files in IGES format. These geometry files need a lot of cleaning, merging, and smoothing work in Gambit before appropriate CFD meshes can be generated. Moreover, due to the rotating part, the multiple reference frame technique has been used. Due to the unsymmetrical geometry, time dependent solutions have

Fluid field solutions including pressure and velocity fields are imported into our DEM program, where DEM tracer particles are created and evolved as described above. During the simulations, particles are first set at the inlet segment of the pump, then particles flow along with fluid to enter the pump. Particles leaving the pump from the outlet are re-inserted at the inlet as in semi-periodic boundary condition, so that particle dynamics gradually reaches steady state. Once steady state is reached, shear and impact energy data are collected and

We performed the wear study on the Metso slurry pump HM300. Ultrasonic measurements of lining thickness as a function of location number is marked on Fig. 1. One hundred thousand (100,000) ore particles from size 0.5 mm to 3 mm are generated for DEM. Shear energy spectra are used in wear analysis as described by Equation 12. The prediction error was found to be less than 4%. The comparison of the wear prediction with measurements is shown in Fig. 2. We performed the ore breakage and valuable particle liberation study on Metso slurry pump MM300. We employed 5 different sizes (1.4 mm, 5.2 mm, 8.17 mm, 10.6 mm and 15 mm) of DEM particles to represent valuable particles, respectively. There are a total of 400,000 particles for ore, and 1000 particles for each of the 5 groups of valuable particles. Based on the impact energy spectra, the particle breakage percentage can be calculated according to Equation 13. The liberation percentage is proportional to be the breakage rate of the ores. The valuable particle breakage percentages and liberation percentages are estimated as shown in Table 1. The valuable particle liberation percentages are under 5% and the valuable particle breakage percentages are very small (< 0.1%) at the pump conditions under consideration. The trend of liberation percentages changing with flow rate and with solid concentration meet the qualitative expectation. As the volume flowrate increases, although the increased velocity enhances the collision of particles, but the retention time decreases, therefore, the liberation

to obtain the fluid dynamics solutions assuming single phase mixture flow.

to be pursued, thus the sliding mesh method has been used.

recorded for at least three revolutions of the pump.

we use a combination of both approaches. We first obtain the intermediate volume fraction in a way like Equation 23, then smooth it over SPH particles in the normal way as

$$
\theta\_b = \sum\_a \theta\_a \frac{m\_a}{\rho\_a} \mathcal{W}\_{ab} \, . \tag{43}
$$

We have to go back to formulate the viscous term in order to complete the set of modeling equations. This work employs a SPH viscous diffusion model used by Morris et al. (1997):

$$\begin{split} \boldsymbol{\tau}\_{ab} &= \sum\_{b} \frac{m\_{b}(\mu\_{a} + \mu\_{b}) \mathbf{r}\_{ab} \cdot \nabla\_{a} W\_{ab}}{\rho\_{a} \rho\_{b} \, r\_{ab}^{2}} \mathbf{u}\_{ab} \\ &= \sum\_{b} \frac{m\_{b}(\mu\_{a} + \mu\_{b}) \mathbf{u}\_{ab}}{\rho\_{a} \rho\_{b} \, r\_{ab}} \left( \frac{\mathbf{d} \, W}{\mathbf{d} \, r} \right)\_{ab} \end{split} \tag{44}$$

where *μ* is the dynamic viscosity. This expression uses only the first order kernel derivative. It conserves translation momentum accurately, while angular momentum is only approximately conserved. For applications with low fluid velocities like in this work, this formulation is appropriate.

In this work, several practical approaches suggested in a series of SPH publications by Monaghan (1989; 1994; 2000) are employed in order to build a stable, robust solver. First, the so-called XSPH velocity correction is employed:

$$\frac{\mathbf{d}\,\mathbf{r}\_a}{\mathbf{d}\,t} = \mathbf{u}\_{\text{xsph},a} = \mathbf{u}\_a - 0.5 \sum\_b \frac{m\_b}{0.5(\rho\_a + \rho\_b)} W\_{ab} \mathbf{u}\_{ab} \,. \tag{45}$$

The XSPH variant is used to move the particles and also in the continuity equation for consistency. The adjustment is important for free surface flows and useful for high speed flows. It basically keeps the particles more ordered and moves particles in a velocity similar to the average velocity in their neighborhood to prevent fluid penetration. Second, artificial pressure is introduced into the momentum equation to avoid SPH particle clustering. See Monaghan (2000) for details.

The overall time advancing scheme for a DEM and SPH coupling system is: 1) calculate the solid and fluid coupling terms based on the old field values, 2) calculate the fluid–fluid interactions and integrating SPH particles. 3) calculate the solid–solid contact forces and move solid particles.

Two different time stepping schemes for SPH particle integration have been applied and tested in this work. One scheme is analogous to the explicit method of control volume CFD for incompressible flow, where the source term of pressure and force terms are calculated based on the old velocity field at *t n*, then the pressure equation is solved, and finally the new pressure field is used to update the velocity to *t n*+1. Similarly in the SPH solver, at the first sweep, the density changes and force terms are calculated based on fields at *t n*; at the second sweep, the pressure and density are updated; at the third sweep, the velocity fields are updated and the particles are moved. The method is enhanced by incorporation of the leap-frog approach: velocities are updated at intervals midway, *t n*+1/2, between time steps *t n* and *t n*+1.

The other scheme is the simple predictor-corrector (Monaghan & Kocharyan, 1995), or, the second order Runge-Kutta method. First, values (velocity, density, position) at *t n*+1/2 are predicted from *t n* and *t n*−1/2. Then, force and other changes under the predicted conditions are calculated. After that, field values at *t n*+1/2 are corrected using the new changes. Finally, values at *t n*+1 can be obtained using the trapezoidal rule. This method can achieve second order accuracy.

Boundary conditions for SPH particles are similar to that for solid DEM particles: spring-dashpot contact models are applied to the particles assuming spherical particles with diameter equal to the initial separation distance. The contact coefficients are taken from DEM parameters for the same size of DEM particle.

#### **5. Applications**

12 Will-be-set-by-IN-TECH

we use a combination of both approaches. We first obtain the intermediate volume fraction in

*θa ma ρ*ˆ*a*

We have to go back to formulate the viscous term in order to complete the set of modeling equations. This work employs a SPH viscous diffusion model used by Morris et al. (1997):

> *mb*(*μa* + *μb*)**u***ab ρaρb rab*

where *μ* is the dynamic viscosity. This expression uses only the first order kernel derivative. It conserves translation momentum accurately, while angular momentum is only approximately conserved. For applications with low fluid velocities like in this work, this formulation is

In this work, several practical approaches suggested in a series of SPH publications by Monaghan (1989; 1994; 2000) are employed in order to build a stable, robust solver. First,

*b*

The XSPH variant is used to move the particles and also in the continuity equation for consistency. The adjustment is important for free surface flows and useful for high speed flows. It basically keeps the particles more ordered and moves particles in a velocity similar to the average velocity in their neighborhood to prevent fluid penetration. Second, artificial pressure is introduced into the momentum equation to avoid SPH particle clustering. See

The overall time advancing scheme for a DEM and SPH coupling system is: 1) calculate the solid and fluid coupling terms based on the old field values, 2) calculate the fluid–fluid interactions and integrating SPH particles. 3) calculate the solid–solid contact forces and move

Two different time stepping schemes for SPH particle integration have been applied and tested in this work. One scheme is analogous to the explicit method of control volume CFD for incompressible flow, where the source term of pressure and force terms are calculated based

pressure and density are updated; at the third sweep, the velocity fields are updated and the particles are moved. The method is enhanced by incorporation of the leap-frog approach:

The other scheme is the simple predictor-corrector (Monaghan & Kocharyan, 1995), or, the

second order Runge-Kutta method. First, values (velocity, density, position) at *t*

*mb*(*μa* + *μb*)**r***ab* · ∇*aWab ρaρb r*2 *ab*

*Wab* . (43)

*Wab***u***ab* . (45)

*n*; at the second sweep, the

*n*+1.

*n*+1/2 are

*n* and *t*

(44)

**u***ab*

 *ab*

d *W* d *r*

*mb* 0.5(*ρa* + *ρb*)

*n*, then the pressure equation is solved, and finally the new pressure

*n*+1/2, between time steps *t*

*n*−1/2. Then, force and other changes under the predicted conditions

*n*+1. Similarly in the SPH solver, at the first sweep, the

*n*+1/2 are corrected using the new changes. Finally,

a way like Equation 23, then smooth it over SPH particles in the normal way as

τ*ab* = ∑ *b*

the so-called XSPH velocity correction is employed:

d **r***a*

Monaghan (2000) for details.

on the old velocity field at *t*

field is used to update the velocity to *t*

velocities are updated at intervals midway, *t*

*n* and *t*

are calculated. After that, field values at *t*

solid particles.

predicted from *t*

appropriate.

= ∑ *b*

d *t* = **u**xsph,*a* = **u***a* − 0.5∑

density changes and force terms are calculated based on fields at *t*

*θb* = ∑*a*

#### **5.1 One-way coupling for a pump**

For slurry in pump operation, solid particles are well mixed with fluid, the particle volume fraction is usually not very high, velocity is high but streamlines do not interweave together, so that the particle–particle collision probability is low, thus the particle–boundary collision is the main interest here since particle–boundary interaction is the major source of wear and the reason for particle breakage under these conditions. DEM simulations were also performed to estimate the breakage percentage and liberation percentage of the particles during the operation of the pump at different conditions. The commercial software FLUENT is used to obtain the fluid dynamics solutions assuming single phase mixture flow.

The CFD simulation of a pump is not easy and straightforward like pipe flow. First of all, for high fidelity we use the manufacturer's geometry files in IGES format. These geometry files need a lot of cleaning, merging, and smoothing work in Gambit before appropriate CFD meshes can be generated. Moreover, due to the rotating part, the multiple reference frame technique has been used. Due to the unsymmetrical geometry, time dependent solutions have to be pursued, thus the sliding mesh method has been used.

Fluid field solutions including pressure and velocity fields are imported into our DEM program, where DEM tracer particles are created and evolved as described above. During the simulations, particles are first set at the inlet segment of the pump, then particles flow along with fluid to enter the pump. Particles leaving the pump from the outlet are re-inserted at the inlet as in semi-periodic boundary condition, so that particle dynamics gradually reaches steady state. Once steady state is reached, shear and impact energy data are collected and recorded for at least three revolutions of the pump.

We performed the wear study on the Metso slurry pump HM300. Ultrasonic measurements of lining thickness as a function of location number is marked on Fig. 1. One hundred thousand (100,000) ore particles from size 0.5 mm to 3 mm are generated for DEM. Shear energy spectra are used in wear analysis as described by Equation 12. The prediction error was found to be less than 4%. The comparison of the wear prediction with measurements is shown in Fig. 2.

We performed the ore breakage and valuable particle liberation study on Metso slurry pump MM300. We employed 5 different sizes (1.4 mm, 5.2 mm, 8.17 mm, 10.6 mm and 15 mm) of DEM particles to represent valuable particles, respectively. There are a total of 400,000 particles for ore, and 1000 particles for each of the 5 groups of valuable particles. Based on the impact energy spectra, the particle breakage percentage can be calculated according to Equation 13. The liberation percentage is proportional to be the breakage rate of the ores. The valuable particle breakage percentages and liberation percentages are estimated as shown in Table 1. The valuable particle liberation percentages are under 5% and the valuable particle breakage percentages are very small (< 0.1%) at the pump conditions under consideration. The trend of liberation percentages changing with flow rate and with solid concentration meet the qualitative expectation. As the volume flowrate increases, although the increased velocity enhances the collision of particles, but the retention time decreases, therefore, the liberation

Solid mass frac- Flowrate valuable-part valuable-part tion in slurry m3/s liberate-% break-% 64% 0.278 3.79% 0.06% 72% 0.278 3.68% 0.06% 80% 0.278 3.62% 0.07% 72% 0.250 4.64% 0.06% 72% 0.278 3.68% 0.06% 72% 0.306 2.97% 0.06%

Using DEM in Particulate Flow Simulations 43

Table 1. Valuable particle liberation and breaking analysis in slurry pump.

< *h*p >= ∑*N*p

To reduce computational cost while still keeping the essential physics of the process, a smaller bed as shown in Fig. 3, thus a smaller number of particles, is used, while the aspect ratio of

Two types of particles with the same density but different sizes were perfectly mixed initially. The simulation was conducted for 10 seconds. The results in Fig. 4 show the close to zero

The initial bed configuration as well as those at *t* = 6 s are presented in Fig. 5. The left panel shows that initially the particles are perfectly mixed and randomly packed. At *t* = 6 s, the two types of particles are still well mixed without visible segregation. The results are better