Original Research Papers
Discrete Element Method (DEM) for Industrial Applications: Comments on Calibration and Validation for the Modelling of Cylindrical Pellets
2015 Volume 32 Pages 236-252
Details
2015 Volume 32 Pages 236-252
As a consequence of increasing computer power and more readily useable commercial and open source codes, Discrete Element Method (DEM) is becoming widely used across a range of applications to simulate increasingly complex processes. This is exacerbating the challenge of setting up simulations for industrial applications. The literature on input parameter selection is divided. A number of papers report methods for their direct measurement. Others, by contrast, propose a “calibration” approach where the particle properties are derived as adjustable parameters by quantitative comparison of experimental and simulation results. This paper reports on the calibration and validation of the input parameters for the specific case of cylindrical tablets represented by conjoined spheres. The initial steps are to not only assign and optimise the DEM input parameters but also optimise the shape representation; what degree of linearity of the edges and angularity of the corners are required to accurately reflect the cylindrical shape. The model was used to simulate two configurations of a rotating drum: an “attrition tester” with a single longitudinal baffle and an un-baffled drum. The results were compared qualitatively and quantitatively with experimental data. While the qualitative comparison was good in most cases, detailed quantitative comparison fared less well, with some significant errors. This study highlights some of the key issues for a wider-spread of industrial applications for DEM.
In recent years, advances in computing speed and power as well as improvements in programming have opened the way to model complex granular flow using the Discrete Element Method (DEM), originally presented in (Cundall and Strack, 1979). The particulate system is modelled as an assembly of singular discrete and interacting particles. It has been used to investigate a number of complex particulate systems, as a result of considerable scientific advances in the development of models for capturing the information and detailed interactions at the particle level. The DEM technique provides much needed insight into the mechanisms governing particle flow. The discrete nature of the outputs for such numerical simulations can enhance fundamental understanding of the granular motion and can then help to improve the design and operation of systems involving particulate material (Cleary, 2000). A key application of DEM is the possibility to model equipment of complex geometry and the complex kinematics for these geometries. The developers of commercial DEM packages, such as EDEM (DEM-Solutions, 2014) or open source codes such as LIGGGHTS (CFDEM Project, 2014) have placed emphasis on integration of DEM with computer aided design CAD packages.
The value of DEM hitherto is demonstrated by the broad variety of applications reported in the literature and the wide range of industries where it has been applied, such as chemicals, pharmaceuticals, ceramics, metal, food and agricultural. Many DEM simulations have been published in the literature modelling diverse granular processes such as comminution (Djordjevic et al., 2003; Mori et al., 2004; Gudin et al., 2006; Powell et al., 2008; Tavares and Carvalho, 2009), granulation (Gantt and Gatzke, 2005; Hassanpour et al., 2009; Nakamura et al., 2009), flow in a hopper (Cleary and Sawley, 2003; Anand et al., 2008; Ketterhagen et al., 2009; Ketterhagen and Hancock, 2010), die filling for tableting (Wu, 2008; Guo et al. 2010), fracture of agglomerates (Kafui et al., 1994; Thornton et al., 1996; Ning et al., 1997; Liu et al., 2010), packing of particles (Matuttis et al., 2000; Zhang et al., 2001; Dutt et al., 2005; Lochmann et al, 2006; Fu et al, 2006; Aste et al., 2007), bulk compression of particles (Foo et al., 2004; Hassanpour and Ghadiri, 2004; Samimi et al., 2005; Markauskas and Kacianauskas, 2006; Mehrotra et al., 2009), flow in screw extruders and conveyers (Moysey and Thompson, 2005; Owen and Cleary, 2009), vibratory screening, filling of dragline bucket, conveyor belt design, earth-mover bulldozer plate design (Cleary, 2000, 2010; Zhang et al., 2008). Another significant application for DEM simulations has been to study the mixing of granular materials, which is still poorly understood and is almost certainly carried out in non-optimal fashion in most applications (Roberts, 1998). Granular mixing in blenders has also been extensively studied to compare mixer designs or look at the effect of operating variables on mixing, segregation and flow patterns (Muguruma et al., 1997; Arratia et al., 2006; Chaudhuri et al., 2006; Yang et al., 2008, Xu et al., 2010; Marigo et al., 2013; Alizadeh et al., 2014).
1.1 Barriers to the use of DEM for industrial applicationsDue to computing considerations there are some limitations in the current state of the art for DEM software for industrial applications. In particular there are restrictions on the number, size and shape of particles that can realistically be modelled. These limitations will be overcome in due course by code development and parallel computing using high performance clusters. The main constraint and difficulty for the application of DEM in an industrial context will be related to the determination (or calibration), in a reasonable time scale, of the input parameters necessary for the simulations. These limitations lead to uncertain accuracy of the DEM modelling and often DEM simulations results can be considered only qualitative, or at best semi-quantitative instead of being quantitative; indicative rather than predictive.
The representation of real particle shapes is still one of the key issues along with the size and number of particles that can be modelled. The majority of DEM simulations consider only a small number of particles with diameters in the order of mm and a maximum number of particles in the order few thousands to achieve a reasonable simulation time (Lim, 2002; Lemieux et al., 2008). The simple spherical shape and the low number of particles that DEM can readily handle at the moment is generally not comparable to real industrial systems which typically consist of millions of small, irregular particles often with a particle size distribution. Therefore, modelling capabilities are mainly limited by high computational requirements that results from the explicit time integration scheme based on sequential calculations over the desired modelled time period, with a low time-step. For example, the simulation of one million non-spherical particles in a complex 3D geometry is not reasonably performed by a single processor and this has been addressed by parallelisation of computers using high performance cluster (HPC) or graphic processing unit (GPU) (Radeke et al., 2010; Stoltz et al., 2013). Consider however powder blending in a real industrial mixing operations some of the segregation phenomena that occur over a long period of time: in the order of second or minutes. However, typical times that can currently be simulated in DEM are usually much shorter: in the order of fractions of a second.
Another major challenge for DEM simulations is to simulate non-spherical particles in an efficient manner. The effect of particle shape can be important and it can be a dominant factor in many cases. In dynamic systems, for example mixing operations or hopper discharge, differences in particle shapes can have a significant effect since particle shape contributes to one of the mechanisms that can promote segregation (Cleary and Sawley, 2002; Li et al., 2004; Roskilly et al., 2010; Metzger et al., 2011). The unconfined yield strength, a factor that influences the bulk flow properties of a powder, is also affected by particle shape. The particle shape controls the number of contacts per adjacent particle and the direction of those contacts relative to a line passing through the centre of two adjacent particles (Johanson, 2009). Particle shape, in particular circularity, affects the packing characteristics in terms of void fraction and bed height (Miyajima et al., 2001).
Different methods of modelling complex shapes have been presented in literature, for example complex particles can be represented in DEM codes as clusters of inter-penetrating or non-penetrating spheres, polyhedrons, ellipses, ellipsoids or polygons. Clusters of inter-penetrating or non-penetrating spherical particles with different or constant diameters are rigidly connected to form an approximation of any desired shape by using the multi-sphere (MS) method (Favier et al., 1999; Kruggel-Emden et al., 2008; Kodam et al., 2010; Bharadwaj et al., 2010). The development of non-spherical shapes to represent the real particle can be complicated since more sophisticated particle contact models should be developed rather than the usual simple sphere contact model. The detection of the contact and subsequent calculation of the forces and torques for edge-edge, edge-corner, corner-corner contacts can be complicated and computationally expensive in case of polygons and polyhedrons (Favier et al., 1999). However, by using the multi-sphere method it is possible to ensure computational efficiency for contact detection and force calculation and many DEM packages implement such methods (Kodam et al., 2010). A drawback of this method is that to precisely represent a complex shape requires clustering of many spherical particles which increases computational memory requirements proportionately.
Probably the most important issue in the application of DEM for real industrial applications is, however, the specification of the material or particle properties required as input parameters by DEM; whether they are evaluated by measurement or by “calibration”. This is particularly true in the case of real, non-ideal particles, especially when they are of complex shape and small (sub-millimetre). The simulation input parameters are in fact often not measured and the values are sometimes assumed without proper justification since there is no robust method for their determination (Asaf et al., 2007). Whilst DEM is being increasingly used, experimental quantitative validations of the simulation results are rarely reported (Chung and Ooi, 2008). It is often not clear how the input parameters are selected or when they are calibrated or measured and the final simulation it is often not validated. Thinking about predictivity of the resulting model, in order to interpret the simulation results confidently correct input parameters should be selected in addition to validation of the simulation results against experimental data (Dong and Moys, 2003).
Two main approaches for the definition of the DEM input parameters can be found in the literature. The first approach is widely employed and involves the determination of DEM simulation parameters by comparing numerical results with some bulk experiments and “calibrating” (more realistically described as fitting) the input parameters to achieve a combination thereof such that the simulations and experiments give comparable macroscopic results. For example, an automated optimisation algorithm to select the ideal set of parameters that best match a DEM simulations with experiments in the shortest time frame possible has been developed by DEM-Solutions Ltd for their EDEM software. The idea is to conduct a laboratory scale experiment that recreates the typical flow regime and bulk behaviour seen in an idealised application and run a series of simulations varying the combination of input parameters that allow a global fitting between DEM simulations and experiments (Johnstone, 2010). This was a sand-pile test where the angle of repose for the experimental pile was matched to the one generated by DEM simulations as functions of various combinations of input parameters. Since the only output that is checked is the angle and the input parameters can have multiple combinations many sets of these would give the same macroscopic output (Curry et al., 2010). A similar DEM parametric approach has been carried out for the operation of a Freeman powder rheometer to understand how the measured force and torque values are affected by key particle properties such as particle size and shape (Asaf et al., 2007). Numerical simulations based on the Hertz Mindlin contact model were compared with experimental results whereby it was shown that particle shape and material parameters such as friction coefficients (sliding and rolling) were important factors determining the force and the torque on the rheometer blade. It was observed interestingly that the DEM fitted input parameters for glass beads were different to those reported in the literature. It was also noted that good agreement with experimental results could be obtained for various combinations of particle shape, sliding and rolling friction coefficients.
These two examples suggest that since the parameters are identified using macroscopic comparison, the parameters for the contact laws seem to be poorly linked to measurable physical parameters of the powder material and different combinations of input parameters can lead to similar end results (Pizette et al., 2010). They emphasise that parameter estimation via a fitting route where the number of parameters exceeds the number of experimental variables is a fundamentally ill-posed problem and thus inevitably has no unique solution. These observations infer that given the general approach to determination of these input parameters they do not necessarily have a physical meaning but rather are simply treated as adjustable parameters, irrespective of the fact that the model is composed to give them a distinct physical meaning.
The second approach involves the development of systematic methodology for the direct experimental determination of the required input parameters. Some of the necessary input parameters are easily measurable while others are not, and this strongly depends, for example on the size and shape of the particles. The influence of particle size, shape, solid density and porosity is of course important and their determination is relatively straightforward. For example for the particle shapes, particle descriptors can be generated from digital image segmentation technique and implemented as particle shape parameters to generate corresponding irregular shaped DEM particles (Williams et al., 2014). Conversely, measurement of single-particle properties such as the coefficient of restitution, friction, Young’s modulus of elasticity and Poisson’s ratio are possible with available tests but are not always practical (Grima and Wypych, 2011). Considering the other particle contact properties, attempts are presented in the literature, but mainly applied to particles with diameters in the order of mm. The determination of DEM input parameters has been reported for pharmaceutical tablet processing operations. Friction was measured using a rheometer used as “a Pin-on-disk” where the measurement of the friction coefficient between particles and aganst a rotating disk was measured for relative large tablets (Suzzi et al., 2012). Rolling friction was measured for large pharmaceutical tablets and several common materials, such as glass beads and steel ball bearings using a slope with a given height and by measuring the distance travelled by the particle (Ketterhagen et al., 2010). The coefficient of restitution has been measured extensively in the literature using high-speed video, this is calculated as the ratio between the final velocity after the impact and the initial velocity before the impact for spherical particles. For non-spherical particles however the impact will result in a more complex trajectory involving particle rotation and the description of such impact is thus significantly more complicated. The three-dimensional motion of irregular shaped particles to determine the resulting coefficient of restitution, using polyethylene pellets, showed the wide spread for the measured coefficient of restitution (Hastie, 2013). Elastic and elastic-plastic macroscopic mechanical properties for mm size zeolite granules were used to calibrate the mechanical properties of the modelled granule in DEM using the force-displacement behaviour measured by the diametrical compression test (Müller and Tomas, 2014). Overall, it is clear that obtaining the input parameters for the real materials is a challenge, especially if the particle size is small and often even without considering other important attributes of each particle (e.g. surface asperities, particle shape, stiffness, hardness and size). It is not clear that this is a realistic objective for heterogeneous materials with irregularities in particle shape and size, non-uniform surface asperities, inconsistency in moisture content and non-linearity of friction (Suzzi et al., 2012).
A further complication is that particles are all single entities within a population. Commonly, only average material values are used in DEM simulations as input parameters. Particles can have a size distribution and some (most?) of the input parameters in truth are a function of a given particle’s diameter within that size distribution. Hence there are some questions that need to be addressed: how or should this aspect be taken into account? Is it viable to use an average value for different discrete particles?
1.2 Considerations of Computational TimeGenerally the total computational time for a given simulation is a combination of multiple factors:
In the DEM simulation the time-step, ΔTstep, is defined as the time between each iteration. A simulation is stable only if the timestep employed is lower than a critical timestep ΔTcritical, which is generally defined as a fraction of the natural frequency of an equivalent mass-spring system. For the commercial EDEM software the critical time step is defined as the critical Rayleigh timestep, TR as follows:
| (1) |
Where TR is function of the particles properties such as radius and density, R is the particle radius, ρ the particle density, G the shear modulus and ν the Poisson’s ratio.
A number of techniques (tricks?) exist for reducing the computational cost DEM simulations. One of this is “particle scaling”, also referred to in the literature as particle ‘coarse graining’ or ‘parcel-based, by artificially increasing the particle size (Radl et al., 2011; Shigeto et al., 2011; Hilton and Cleary, 2012). The shear modulus, G, also has a major influence in the critical time step. It is common practice in the literature to arbitrarily reduce the value of the shear modulus as this dramatically increases the critical time-step, thus reducing the simulation time. However, it has been shown that this approach could lead to errors on the simulation results depending on the actual system that is modelled (Lommen et al., 2014). Using an angle of repose test and a penetration test it was shown that by reducing the shear modulus below to 107 Pa or keeping the average normal overlaps below 0.3 % of the particle radius these parameters would not alter the simulation results. On the other hand, it was shown using a bulk compression test that properties such as the bulk stiffness and bulk restitution changed as a result of the stiffness reduction. This example showed that DEM users should always be cautious and verify the impact of varying the input parameters case by case and possibly verify the approach through some real validations.
In the context of industrial applications for DEM it is maybe more prudent to consider is the actual “total modelling time” rather than only the simulation time where the former includes the time required to evaluate the DEM input parameters, the geometry and simulation set-up time, the time for the simulation itself plus the data post-processing and “validation time”. To provide justification and economic value to many typical industrial applications, the total modelling time to set up, solve and analyse a DEM model has to be relatively quick and faster than actually carrying out experiments, unless it is justifiable by the fact that computer simulations can give insight of information at the particle level that would be difficult or impossible to obtain by experiments.
In this paper we assess the validity and time requirement of the calibration approach to model parameter selection, and test the quality of the model by evaluating the accuracy of the predictions for different scenarios against experimental data. Specifically, a shape optimisation for cylindrical α-alumina pellets and a subsequent parameter calibration using the “sand-pile” test were performed. The model was then used to simulate two different geometries of a rotating drum and the results compared to experimental data.
The work in this paper was done using the commercial DEM code (EDEM Version 2.4.3) using a E5645 2.40GHz DUAL-QUAD PC machine with 12GB of RAM using a 2 core licence.
The handling of particulate material in rotating drums is a fundamentally important unit operation in many industries such as pharmaceutical, chemical, minerals, food and agricultural industries, as well as many others. The aim of this work is to model the motion of 3.12 × 2.96 mm alumina cylindrical pellets within two rotating drum systems: a conventional cylindrical rotating drum and a baffled attrition test rig.
As discussed above, the first major step was to determine the input parameters for the simulations and additionally, in the case of a complex shape, what is the best combination of multiple spheres to represent the right pellet geometry in order to capture the its effect during the tablet motion. The approach that has been applied in this work is represented in Fig. 1.
Approach for the DEM model: shape calibration, input parameters calibration and model validation.
The intent was thus to establish an adequate representation of pellet shape (using the multi-sphere representation), then “calibrate” the DEM input parameters and subsequently validate the model—all using a different experimental method:
Each of these activities and their results will now be described along with comments on the time required to complete the task.
The ultimate aim of the present study was to model the motion of cylindrical alumina pellets (3.12 × 2.96 mm and pellets density of about 1800 kg/m3 within two rotating drum systems: a baffled rotating drum and a plain (i.e. no flights) rotating drum. The first aspect of building the model was to assess and optimise the representation of the cylinders using the multi-sphere approach. This was carried out by comparing the model results to experiments for loading the pellets into a horizontally and vertically aligned cylinder. This approach has previously been used for the validation of a semi-stochastic particle loading simulation software (Xu et al., 2008).
Considering the Hertz-Mindlin contact model employed in the EDEM software, if only gravity is considered as an external force then the particle characteristics and contact input parameters that must be defined are:
Pellet shape representations were generated using the multiple sphere approach; the cylindrical pellet template was generated in CAD and approximated with conjoined spheres. For the specific purpose of establishing the best representation of pellet shape a single set of DEM input parameters was arbitrarily assumed. These values are listed in Table 1, which were used in all simulations relating to the optimisation of the pellet shape representation. Note that the shear modulus is set to a low value of 2 × 106 in order to shorten simulation times. The only pellet properties that were actually known were the pellets size and density whereas all the other parameters were chosen according to the authors’ prior DEM experience.
| Parameter | Value |
|---|---|
| Particle (p) density, ρp [kg/m3] | 1800 |
| Particle shear modulus, G [Pa] | 2 × 106 |
| Particle Poisson’s ratio, νp | 0.25 |
| Vessel (v) density, ρv [kg/m3] | 1000 |
| Vessel Poisson’s ratio, νv | 3 × 109 |
| Particle-particle static friction, μs-pp | 0.25 |
| Particle-particle rolling friction, μr-pp | 0.65 |
| Particle-vessel static friction, μs-pv | 0.6 |
| Particle-vessel rolling friction, μr-pv | 0.65 |
| Particle-particle restitution, εpp | 0.01 |
| Particle-vessel restitution, εpv | 0.85 |
An example of the effect poor particle representation on packing of cylindrical pellets into a rotating drum is given in Fig. 2. The experiment was carried out using a horizontal Perspex drum (D × L: 240 × 90 mm) which was filled to 20 % volume fill level. Using the known pellet density and the weight of material this corresponded to about 22,000 pellets in the DEM. Two pellet shapes, Shapes 1 & 2 were generated and the pellets properties such as mass, volume, moments of inertia were calculated using the grey cylindrical template reported in Fig. 2.
Lower packing density (fill level) is obtained if a poor pellet shape representation. Left side: Perspex drum (D × L: 240 × 90 mm) filled with 20 % fill level cylindrical pellets. Right side: DEM results using 22,000 pellets.
Due to the poor representation of the real shape the overall bulk or packing density is evidently wrong (dotted line in Fig. 3). This would have a significant effect on the overall flow and behaviour in a dynamic system. The main reason is probably excessive rounding of the sharp corners, which it can easily be envisaged will result in closer packing of the pellets and thus a higher packing density—as witnessed by the results. Evidently a geometrically more correct representation of the pellet shape is required.
Cylindrical pellet packing into a 52 mm tube using refined pellet shape representation by increasing conjoined spheres number and decrease sizes.
Representing complicated shapes by using a simple sphere can thus lead to errors. Glued spheres clusters do not duplicate the overall non-spherical particle shape exactly, especially the sharp corners. Although algorithms have been developed to approximate particle shapes using glued spheres (Radshaw and O’Sullivan, 2004; Wang et al., 2006), the fine details of the overall non-spherical particle surface remain irregular even with a large number of constituent spheres.
This inaccuracy can be improved by using for example a higher number of spheres with smaller or mixed diameters to achieve a better representation of the square sharp corners and flat edges, Table 2.
| Type | n° (radius-mm) |
|---|---|
| Shape 1 | 12 (0.765) |
| Shape 2 | 24 (0.765) |
| Shape 3 | 24 (0.765) |
| 18 (0.35) | |
| Shape 4 | 48 (0.765) |
| 34 (0.35) |
By more accurately defining the particle details better precision and a more accurate representation of the packing would be achieved. However this considerably increases the computational time since in the multi-sphere (MS) a higher number of particles would be considered into the simulation.
Experiments for comparison were run by loading pellets into a vertical cylindrical Perspex container (diameter 52 mm). Three different fill levels were used: 200, 250 and 300 grams corresponding respectively to 5212, 6515 and 7818 pellets. Equivalent DEM simulations were run for the different shape representations, each with an increasingly complex multi-sphere assembly.
The qualitative output from this is shown in Fig. 3. Again, it is shown that an inadequate particle shape representation can lead to clear errors and in this case artefacts in terms of the packing structure of pellets. The pellet shape is exported from the DEM software using a pure cylindrical basis (centroid and orientation angles).
The poor representation of the particles corners leads to apparent overlaps due to the rounded corners and to pellets penetrating the cylindrical container wall.
The simulation and experimental results are compared quantitatively in Fig. 4. Shape 4 appears to provide a reasonably good approximation for the packing behaviour of the real pellets in terms of bed height since the DEM results fall within the experimental error, whereas the other two shapes were statistically different. Shape 4 also shows little evidence of significant overlap and intersection of the cylindrical wall.
Plot comparison bed height measured experimentally and DEM simulation results for shape 2, 3 and 4 (experimental results = mean ± SD, n = 3).
However, to run just 2 seconds with about 8,000 pellets the time was just over 45 hours on a 2-core machine, 9.5 hours using 8 cores. This does not bode well for future simulations over longer times for actual applications also involving many pellets.
3.2 DEM input parameters calibrationTo determine the DEM input parameters for the cylindrical pellets the typical “calibration” strategy was used to determine their values, as discussed in the introduction, using the “sand-pile” test as the exemplar. For the experiments, see Fig. 5, a Perspex funnel was filled with 300 grams of the pellets and the pellets were dropped on to flat Perspex surface. The angle and the height of the pile were compared to the results of DEM simulations carried out using the “Shape 4” representation from the preceding section.
Parameter calibration: experimental and simulation set-up plus representative results.
| 300 grams (7818 pellets) | Shape 2 | Shape 3 | Shape 4 |
|---|---|---|---|
| 2 cores license | 3.5 h | 12.6 h | 45.1 h |
| 8 cores licence | 9.43 h |
It is expected that varying the various DEM input parameter values would result to a different final results for the discharge of pellets from the container in terms of shape, angle and height of the pile. The successive parameter sets used in simulations for the calibration of the shape 4 pellets are reported in Table 4. This of course considers only a subset of the whole parameter set. The artificially low value for the shear modulus has been maintained and the literature value for Poisson’s ratio has also been accepted. Initially the coefficients of rolling friction were set to zero as in previous work regarding tablet motion in a pan coater it was assumed zero for shaped particles as the rolling resistance was expected to be accounted for by the actual shape representation (Ketterhagen, 2011).
| Run | εpp | μspp | μsrpp | εpw | μspw | μrpw |
|---|---|---|---|---|---|---|
| 1 | 0.4 | 0.625 | 0 | 0.4 | 0.625 | 0 |
| 2 | 0.5 | 0.625 | 0 | 0.45 | 0.625 | 0 |
| 3 | 0.55 | 0.625 | 0 | 0.5 | 0.625 | 0 |
| 4 | 0.65 | 0.625 | 0 | 0.55 | 0.625 | 0.05 |
| 5 | 0.85 | 0.625 | 0 | 0.75 | 0.625 | 0.05 |
| 6 | 0.4 | 0.625 | 0 | 0.4 | 0.625 | 0.05 |
| 7 | 0.4 | 0.625 | 0 | 0.4 | 0.625 | 0.02 |
| 8 | 0.35 | 0.625 | 0.05 | 0.35 | 0.625 | 0.05 |
| 9 | 0.35 | 0.625 | 0.05 | 0.35 | 0.625 | 0.05 |
| 10 | 0.4 | 0.625 | 0.065 | 0.4 | 0.625 | 0.025 |
| 11 | 0.4 | 0.625 | 0.065 | 0.4 | 0.625 | 0.04 |
| 12 | 0.4 | 0.7 | 0.07 | 0.4 | 0.7 | 0.05 |
The results from the simulations with these successive parameter sets are shown in Fig. 6 and the quantitative outputs from each simulation for pile height are compared with the experimental value in Fig. 7.
DEM results for the pellet-pile for the different input parameter sets given in Table 4.
Height for the pile obtained from DEM simulations compared to the one obtained experimentally.
For the early runs the pile could not actually be formed due to excessive rolling of the pellets.
The coefficient of rolling friction actually has a significant impact on the results for this particular particle flow problem as the cylindrical pellets can just roll on the side of the cylindrical shape. It was necessary to introduce a non-zero value for both the particle-particle and particle-surface rolling coefficients. The other parameters were relatively unchanged from their assumed initial values, with only a modest increase in both static friction coefficients.
In total 12 simulation runs were required to achieve a parameter set that yielded a satisfactory representation of the “pellet-pile”. It can be noted that for the last 3 sets of parameters (Simulation 10, 11 and 12) the pile heights were all within the experimental measurement error. This implies that any of the last 3 sets could be assumed as “calibrated” set of input parameters. For the subsequent validation in the next sections set 12 was arbitrarily used. A more extensive systematic statistical and theoretical evaluation should be carried out in order to assess the relationship between the DEM input parameters and the model outputs for a number of several granular flow scenarios. However, this reinforces some of the time scale limitation for industrial applications of DEM.
Each simulation was run for 4 seconds, sufficient time for the pellets to settle into a static state. Each of these simulations required 60 hours on 2 cores, (estimated at 13 hours on 8 cores), thus a total of 720 machine hours plus the time to set up the simulations and do the results post processing. This is very clearly a very time-intensive step. The timestep used in this simulation was 30 % of the critical timestep given by Eqs. 1, and it corresponded to 1e–5 seconds.
The total required (elapsed) time to achieve this will, of course, not fall in linear fashion. The time to set up the simulations and to process the results is likely to be more or less constant and not a function of the power of the computer used to do the simulations. In reality parameterisation of the model has taken several weeks, and this is a significant barrier to the wider application of DEM to industrial problems.
The previous results have established a “calibrated” DEM model: “Shape 4” and input parameter set 12. Results from DEM simulations will now be compared to two pellet flow situations:
The basis for comparison of the experimental results was a simple qualitative and macroscopic comparison of particle distributions in the first case. For the latter, however, a far more quantitative approach was taken and the comparison based on:
The attrition tester is described in an ASME procedure (ASTM D 4058, 2011), for the attrition testing of catalyst particles. It consists of a horizontal rotating drum with D = 254 mm, L = 152 mm and a single baffle with 51 mm height extending the full length of the drum. For the experiments 115 grams of pellets were used and the drum was rotated at 40 and 60 rpm. A high speed video camera (Photron SA4 FASTCAM) was used to capture images of the particle flow through one of the end plates using a frame rate of 1,000 Hz.
For the simulations 3000 pellets (an equivalent to the pellets used in the experiment) were used and the simulation was run for 10 sec. The Pellets were created and allowed to settle from 0–3 sec before the rotation was stated. The rotation was then modelled for 7 seconds. At 40 rpm this corresponds to 4.7 rotations while at 60 rpm this corresponds to 7 rotations. Snapshots from the DEM results after the total 10 second simulations were compared qualitatively with the experimental images in Fig. 8 and 9 for 40 rpm and for 60 rpm respectively.
Comparison pellets flow for attrition test with 3000 pellets between high speed camera and DEM at 40 rpm. 8 baffle positions (represented in yellow).
Comparison pellets flow for attrition test with 3000 pellets between high speed camera and DEM at 60 rpm. 8 baffle positions (represented in yellow).
Both sets of results showed a good qualitative comparison between the general behaviour of the pellets in terms of their locations and the spread of the pellet cluster at all baffle plate angles. Qualitatively we appear to have, therefore, a good and accurate model. However, it is important to note some interesting little discrepancies between experiments and simulations. For example, in Fig. 8 when the baffle is in the vertical position (frame 7), the tablets take off the baffle earlier than the simulations. Instinctively, it could be attributed to the selection for the particle-wall interaction in parameters, perhaps not being high enough. However, this is not valid in Fig. 9, where for the same baffle position the tablets seemed to detach correctly from the baffle. It is also worth noting that in Fig. 8, frame 5, the high speed image shows more tumbling than the DEM results, this might be attributed to particle-particle interaction parameters. Therefore, it could be said that this approach of inferring the set of parameters from the calibration tests employed is not suitable to obtain an optimal calibration. However, it is also important to note that most of the DEM literature would consider this qualitative comparison good enough and it would claim that the model has been validated arguing that a model will never be perfect but good enough to give a useful result.
For the record, each of the simulations required 12 hours machine time on a 2 core computer (estimated at 2.5 hours for 8 cores). In this case the data post-processing is only qualitative and therefore relatively quick—as indeed is the processing of the experimental results. This, of course, is not the case if a quantitative approach to validation is taken.
4.2 Rotating drumThe rotating drum was made of Perspex and of diameter D = 240 mm and length L = 90 mm with flat endplates and no baffles or flights. It was operated at a 20 % volumetric fill level at three speeds: 20, 35 and 50 rpm. High speed videography was used to capture images to allow evaluation of the upper and lower angles of repose (Fig. 10). The digital video images were also used to evaluate Eulerian particle velocity maps at the end plate using a Particle Imaging Velocimetry (PIV) approach to image analysis described later.
Angles of repose, upper and lower, for the Perspex rotating drum 20 % filled.
For the simulations, 22,000 pellets were introduced into the drum, representing also a 20 % fill level. They were allowed to settle for 2 seconds and rotation was then started. The rotation was modelled for 3 seconds corresponding to 2.3 and 5.8 rotations at 20 and 50 rpm respectively. Snapshots for comparison with the experiments were taken at the end of the total 5 seconds simulation, from which the dynamic angles of repose were evaluated, Fig. 11. Each of the simulations required 220 hours machine time on a 2 core computer (estimated at 46 hours for 8 cores).
Angles of repose, upper and lower, from DEM simulation using 22,000 pellets represented by Shape 4 and Input parameters by simulation 12.
Fig. 12 shows the evolution of the average pellet speed though the 3 seconds of simulation time for both simulated drum speeds. This shows that a steady state has ostensibly been achieved after 1.3 and 2 seconds (20 and 50 rpm respectively) and thus the 3 seconds simulation time seems to be adequate.
Average pellets bed velocity when the rotation is started. Showing the bed reaching a steady state quite quickly.
Fig. 13 shows the global quantitative comparison between the experiments and the simulations for the dynamic angles of repose and the bed height.
Comparison angles of repose, upper and lower and bed height from DEM simulation and high speed camera. These values where taken after 3 seconds of drum rotation.
Although general trends for the angles as function of the speed are in good agreement, some discrepancies are apparent. The bed height predicted by the DEM simulations is much higher than the experimental values as is the upper dynamic angle of repose (β). These discrepancies could be related, for example to the value for the static friction value between particle and wall which is probably inadequately captured by the “sandpile” test used for calibration or perhaps the shape representation due to the rounded corners. It could be argued that the parameter calibration should be repeated; but this surely defeats the stated object of the calibration exercise!
It has recently been emphasised that validation of complex models, such as DEM, against a single objective function may be neither adequate nor sufficient (Stitt et al., 2013). Therefore, the model here is compared to detailed particle velocity fields in addition to the global results above. For the DEM results, the velocity is calculated by dividing the side view in a grid 0.5 × 0.5 mm (much smaller than the pellets sizes) and by averaging the velocity for the pellets going through each element of the grid over 1 second (when the bed average velocity is constant in Fig. 12) to generate an Eularian velocity field that can be compared to the Particle Imaging Velocimetry (PIV) results.
PIV is a well-established and widely used technique for velocity imaging of both gases and liquids (Raffel et al., 2007) that can also be applied to particulate flows. It measures whole velocity fields by taking two images shortly after each other and calculating the distance individual particles travelled within this time, tracking individual particles using a cross correlation technique. From the known time difference and the measured displacement the velocity is calculated to generate a Eularian velocity map. More extensive descriptions of the procedure are given in the literature.
Successive digital images at a frame rate of 1,000 Hz for the rotating drum were recorded using the high speed camera were analysed using the PIV software provided by (LaVision, 2014).
Fig. 14a–d, show the comparative velocity fields evaluated from the simulations and the experimental PIV for both drum rotating speed: 20 and 50 rpm. As discussed previously, the main pellet bed features, trends in bed shape and height, are reasonably well captured by the DEM simulations. Moreover, the velocity fields are well captured near the drum wall and the central part of the pellet bed. For both drum speeds the lowest pellet velocities are correctly captured within the central core part of the pellet bed. However, the top bed velocity is much higher from the DEM. In Fig. 14e the distribution for the velocity field is given to show the general good agreement again in terms of general trends but poor validation in capturing the top bed pellet velocity.
a–d) Comparison for the velocity field between DEM simulations and PIV experiments. e) Velocity distribution. The pellets velocity on the top of the bed is over predicted by DEM.
It is interesting to note that comparisons of PEPT and DEM also indicate that DEM systematically over-predicts velocities by approximately 10 % in free flowing systems (Hassanpour et al., 2011; Marigo et al., 2013). The discrepancies in this case do, however, seem to be larger than this. To what extent this depends on the selection of the input parameters (perhaps choice of static/rolling friction between pellets material and drum wall?) is of course not clear. What is clear is that the model calibration carried out using the sand-pile test does not appear adequate for its application to a rotating drum.
4.3 DiscussionDespite the lengthy and careful approach to model parameterisation, in line with literature methods, the derived model shows significant discrepancies when compared to experimental validation data. While it is possible to hypothesise on the sources of these errors, ascribing the blame to one or two given parameters, due to the nature of the model parameterisation it is not simple to assess the origins of these errors.
On a broader front, this raises questions over the efficiency and accuracy of the “calibration” approach to parameter estimation. In this study the “sand-pile” type test was used—as it has been in the literature. Two measured parameters were used to optimise the values of 6 adjustable fitted parameters. This is an inherently ill-posed problem and as such multiple combinations of parameter values will provide an equally good fit. So, while a good representation of the “sand-pile” test was obtained in the model calibration, there is no way of assuring that the fitted parameter values are correct without detailed statistical evaluation of the quality of fit (e.g. 95 % confidence intervals). Is the “sand-pile” test a reasonable basis for parameter evaluation for a rotating drum model application, as also suggested in the literature? The response to this question has, for the moment, to be equivocal. What is clear is that significant improvements are required in methods to assign the values of the model parameters (which all have clear physical definitions) in DEM applications in order to achieve a quantitatively correct model.
DEM is still in its infancy when compared to other modelling tools such as Finite Element Analysis (FEA) and Computational Fluid Dynamics (CFD). However, the development of DEM and its increasing applications are following a similar pattern to those previously seen with CFD. Many limitations and difficulties regarding the application of DEM for real industrial scale applications exist.
With advances in computational power and programming, DEM is becoming more and more accessible to academia and industry alike. At present, on commercially available desktop computers, simulations of up to a million particles can be performed. On large high performance computing clusters, the trajectories of millions of particles can be computed.
The use of this modelling technique is accordingly expanding and it is being exploited in many applications. Despite this, as extensively reported and demonstrated, there are still some limitations that need to be addressed to use DEM as a possible fully predictive rather than indicative tool for modelling real scale industrial systems. Some of these issues are highlighted as follows:
This paper has investigated the particle shape representation for a system involving cylindrical pellets using tube loading and the model calibration for DEM input parameters using the “sand-pile test. The resulting model has then been subjected to a qualitative and quantitative validation against experimental data for two rotating drum systems; one with a single baffle, one unbaffled. The lengthy procedure from getting the input parameters to the validation DEM application has proven quite challenging for these real industrial (albeit rather small scale) applications.
The results indicate that the use of calibration is not wholly successful in allowing quantitative prediction of the applications. It is, however, not possible to interrogate the root causes of the errors in the final application simulations. The paper also highlights the high demands on computer simulation time and human resource necessary to complete this exercise. There is a need to improve in both of these areas if DEM is to gain more traction as a design rather than a research tool.
The implication is that the current methods of model parameterisation are not adequate. This begs the question as to whether model calibration using a dissimilar system to the final application is appropriate. It has been highlighted that the parameterisation by “calibration” is inherently ill posed. This means that in terms of the fitting exercise there will inevitably be questions over the parameter confidence intervals and cross correlation between some of their values. There are few reported attempts to understand the quality of the fit other than by means of simple residuals.
As a DEM modelling community, there is an over-riding need to better establish protocols for the evaluation and estimation of the key particle properties—the input parameters. This paper is intended to highlight and exemplify the shortfalls with current methodologies and to encourage further work in this unattractive but yet absolutely vital area.
Michele Marigo
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/32_2015016_15.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Michele Marigo is a Principal Scientist at Johnson Matthey’s Technology Centre. He obtained an undergraduate degree with Master’s in Mechanical Engineering from University of Padua, Italy and a doctorate in Chemical Engineering EngD from University of Birmingham, UK. Michele’s expertise include materials science, particle engineering, discrete element modelling (DEM) and finite element modelling (FEM).
Hugh Stitt
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/32_2015016_16.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
Hugh Stitt is a Scientific Consultant in Johnson Matthey’s Technology Centre. He holds BEng and PhD from the University of Bradford, UK. He is a Visiting Professor at the University of Birmingham and at Queen’s University Belfast. He is a Fellow of the Institution of Chemical Engineers and a Fellow of the Royal Academy of Engineering. He has 25 years of industrial research experience across a variety of themes related to catalytic reaction engineering and catalyst manufacture with an emphasis on both experimental techniques and modelling with over 100 refereed publications.
