Prediction of flowability of cohesive powder mixtures at high strain rate conditions by discrete element method

Abstract Flow behaviour of powders plays an important role in defining their performance in many industries. In this study, we present our work on prediction of flow characteristics of binary and ternary mixtures under dynamic conditions using Discrete Element Method and simulating the expended mechanical work of a rotating impeller penetrating a packed bed. Three commonly used methods for calculating a mean have been explored to express the mixture Bond number, namely arithmetic, geometric and harmonic mean. This is done by introducing a weighting factor based on the fractional surface area of each component of the mixture. The mixture Bond number is dependent on the interfacial surface energy, particle size and density; and a wide range of Bond number is covered in our study by varying all three. The expended work of impeller shows the clearest trend with the mixture Bond number when it is expressed in terms of the arithmetic mean. Although we only used binary and ternary mixtures (40 different mixtures in total) in this study, the trend should be applicable to multi-component mixtures and therefore useful as a design aide for powder formulation.


Introduction
Flow behaviour of powders plays an important role in performance of various powder processing operations, such as discharge from hoppers, feeding, dosing and pneumatic conveying. Shear cell testing has long been established for designing suitable hoppers for powder discharge [1], where the shear resistance is characterised as a function of consolidation stress or state of packing. Demand for faster production rate and newer and more advanced products has led chemical, pharmaceutical, and food industries to move to more complex powder formulations, making it more difficult to characterise the powder mixtures at early stages of development. Consequently, there is great interest in developing predictive models for flow behaviour of powder mixtures based on the properties of their individual components. A large number of studies have been reported in the literature for prediction of flowability of powders and their mixtures under quasi-static conditions, mainly focusing on effect of particle size distribution [2][3][4], shape [5][6][7][8], moisture [9][10][11] and cohesion [12][13][14] . The work reported by Capece et al. [15] stands out and is highly relevant, as they have developed a predictive model for flowability of adhesive mixtures using quasistatic shear cell testing. They relate the flow function coefficient measured by shear cell testing to powder mixture Bond number. Their approach takes account of the inter-particle forces, density, particle size, and surface roughness and has been validated for a number of mixture systems, such as single component, binary and ternary mixtures. The level of adhesiveness of the powder mixture was assessed by averaging the granular Bond number of each components and their interaction with other components using a harmonic weighted mean. They proposed a weighting factor based on the fractional surface area of each component of the mixture. In their work, almost all tested materials had a similar adhesive characteristic (i.e. the values of surface energy of components were between 35 and 50 mJ/m 2 ). Harmonic averaging method, as used by Capece et al. [15], tends to give an average value towards the smaller numbers that are being averaged. Hence, their methodology was not tested for extreme cases, where for example one component is free-flowing and another highly adhesive and the bulk powder flow is under medium or high strain rate conditions.
In recent years, the attention of the particle technologists and physicists has focussed on gaining a better understanding of the flow behaviour of particulate solids at intermediate and high strain rate conditions. There are many cases where the dynamic flow behaviour is critical for process design and operation, such as in screw conveyors, mixers, fast feeding for tabletting machines, and discharge from storage devices. In contrast to quasi-static test methods, there is no shear cell that can characterise the bulk parameters at high strain rates, i.e. in the dynamic regime. However, over the past few years, several powder rheometers, i.e. Freeman FT4 powder rheometer, Anton Paar Modular Compact Rheometer and the Couette device have been developed to characterise the resistance of bulk powders to motion as a function of speed of shearing (strain rate). For more detail on the dynamics and analysis of cohesive powder flow at high strain rate conditions, the reader is referred to a recent review by Ghadiri et al. [16]. The Freeman FT4 powder rheometer (referred to as FT4 hereinafter) characterises the flow resistance by measuring the amount of mechanical work expended for a rotating impeller penetrating a certain distance into a powder bed. Several studies have been carried out by several researchers looking at the relevance of this measurement method to bulk powder flow behaviour [5,[17][18][19][20][21][22][23][24].
In this study, we expand on the work of by Capece et al. [15] for calculation of the granular Bond number of powder mixtures by using three averaging methods: harmonic, arithmetic, and geometric weighted means. The proposed averaged granular Bond numbers are then used to predict the expended work (flow energy) of various binary and ternary mixtures using Discrete Element Method (DEM) simulations.

Model materials
SiLibeads® Type S glass beads, provided by Sigmund Linder GmbH, Warmensteinach, Germany, are used in this work for single particle characterisation. Their physical and mechanical properties are characterised for use in DEM simulations rather using calibration approach [25]. The use of spheres (representing glass beads) makes numerical simulations fast and reliable as contact laws for particles interactions are well established. A quantity of the beads was first subdivided into smaller samples by a spinning riffler to obtain representative small samples for particle characterisation. The beads were initially sieved into 425-500 μm and 850-1000 μm size classes by mechanical sieving method. Observing them by Scanning Electron Microscopy (SEM) revealed that there were particles in the sample that were far away from being spherical, as shown in Fig. 1a. The spherical glass beads were separated from non-spherical ones by using a tilted pan with a smooth surface, in which the spherical beads easily rolled down the pan and could be separated from the non-spherical ones. They are shown in Fig. 1b. In DEM simulations, the beads are treated as spheres, which are generated using a random size distribution in the range of 850-1000 μm and 425-500 μm. They are also assigned two densities of 2500 and 1000 kgm −3 to investigate the effect of density by simulations only.

Interfacial surface energy
Since our interest here is to predict the flow behaviour of adhesive and non-adhesive powder mixtures, the surfaces of glass beads are chemically modified to make them more adhesive. This is done by silanisation by applying a commercially available silane coating on their surfaces, known as Sigmacote®, supplied by SigmaAldrich®, Missouri, USA, using the approach of Zafar et al. [26]. The 'apparent' interfacial surface energy of glass beads is then measured experimentally by Drop Test Method developed at the University of Leeds [26]. In this method the particles are subjected to a transient tensile force, produced by impact and the detachment force is calculated, from which the interfacial surface energy is inferred using the analysis of Johnson et al. [27]. The measurements were actually done using a smaller size silanised glass beads (100-200 μm), as the large ones all detached in the test within the impact velocity range tested. The beads were dispersed onto three different surfaces: glass, silanised glass and stainless-steel slides that were stuck to aluminium stubs. Following this step, images of the particles were captured by SEM, as shown in Fig. 2a. The particles were then subjected to impact tests at different velocities, ranging from 0.5 to 2.0 m/s, and the contact time was measured by a Photron Mini AX 200 high speed video camera at 70,000

Coefficient of restitution
This was measured by analysis of impact against stainless steel and glass slides using Photron Mini AX200 high speed video camera at 10,000 frames per second. Upon impact of glass beads to the target of interest, the pre-and post-rebound velocities are calculated using ImageJ software and a particle tracking algorithm developed by Sbalzarini and Koumoutsakos [28]. The coefficient of restitution for each interaction is measured for approximately 100 particles and its average values are 0.95 and 0.93 for impacts against glass and stainless steel plates, respectively.

Coefficient of sliding friction
A simplified methodology compared to work reported by Nan et al. [29] is used here for measuring the sliding friction coefficient between particle-particle and particle-wall for use in DEM simulations. The glass beads are deposited and glued onto a flat glass slide in a packed state using super glue as shown in Fig. 3a. Care is taken to ensure the particles top surfaces are not covered with super glue, which could significantly affect the sliding behaviour of the particles against the smooth surface. The prepared slide is placed faced down against another flat glass surface. The plates are then tilted until the top plate starts sliding down as shown in in Fig. 3b. At this point, the corresponding angle of sliding is recorded. Since the surface of the plate is flat and smooth, coefficient of sliding friction can be calculated using the tangent of the recorded angles. This method of measurement is used for both particle-particle and particle-wall friction coefficient. Sliding friction coefficients are measured at least five times for glass-glass and glass-stainless steel surfaces and are summarised in Table 1.

Discrete element method
The particles are modelled as discrete entities and their translational and rotational motions are tracked individually by solving Newton's laws of motion [30,31]: where m i , I i , v i and ω i are the mass, moment of inertia, translational velocity and angular velocity, respectively; F c,i is the contact force, originating from its interaction with neighbouring particles or walls; M c,i is the contact torque, arising from the tangential and normal contact forces; R i is the rotation matrix from the global to the local coordinate system in which the calculation of the rotation expressed by the Eq. (2) is accomplished. In this work, the elastic contact force is described by Hertz-Mindlin contact model [32], and the adhesive interaction is accounted for by JKR theory [ where Γ is the interfacial surface energy; E* is the equivalent Young's modulus; R* is the equivalent radius; a is the contact radius; α is contact overlap. In the unloading process, the normal contact force F n is not zero when the normal overlap α is negative, as further work is required to separate the cohesive contact. To make the computational time practical, Young's modulus in the simulation is 1000 times smaller than the one measured in the experiment, and the surface energy is correspondingly scaled down in the simulation by keeping the Cohesive number of particle constant [33], as shown in Eq. (5) similar to work of Nan et al. [29]. For calculation of mixture Bond number, the actual measured value of interfacial surface energy is used here rather using the modified value of interfacial surface energy.

Analysis of flow of adhesive mixtures by DEM
The flow behaviour of multi-component powder mixtures as a function of strain rate is studied by simulating the expended mechanical work, commonly termed 'basic flow energy (BFE)' as measured by Freeman FT4 rheometer, for several binary and ternary mixtures using Discrete Element Method modelling. EDEM® software package (DEM-Solutions, Edinburgh, UK) is used here. In this study, the 25 mm diameter Freeman FT4 vessel with the 23.5 mm impeller was used, as shown in Fig. 4. The standard test procedure is applied to the powder, whereby the bed is initially conditioned by rotating the impeller clockwise to gently slice the bed and produce a reproducible and low stress packing state. The cell is then split to remove any material above a bed height of 80 mm. Following this step, the blade is driven down into the bed anti-clockwise with a helix angle of 5°and a tip speed of 100 mm/s. The vertical force acting on the base, and the torque acting on the blade are measured at approximately 200 μm increments of vertical displacement. The BFE is calculated using Eq. (6).   where T, R, α, and F are blade torque, radius, helix angle, and force acting on the base of FT4 rheometer, respectively.

Powder mixture bond number
In order to describe the flow behaviour of adhesive powder mixtures at bulk-scale based on single particle properties, it is essential to consider the particle-scale interactions. In the absence of moisture and electrostatic effects, the inter-particle adhesion is dominated by van der Waals (vdW) forces [6,34]. There are a number of adhesion force models in the literature [27,35,36]. In this study, Johnson, Kendall and Roberts (JKR) model of elasto-adhesive contact force [27] is used: where R ⁎ and Γ are reduced radius and interfacial surface energy of the two particles in contact. The effect of adhesion on particle behaviour under gravitational field is commonly described by granular Bond number [37], When adhesive force is larger than the gravitational force (Bo g > 1), the two particles tend to stick to each other unless separated by additional external forces. For granular Bond numbers smaller than unity, the two particles could easily detach from each other due to gravity and hence they are considered non-adhesive. Bond number has been used successfully for describing bulk-scale powder properties [2,15,[38][39][40]. In particular, Capece et al. [15] used it to predict flowability of powders under quasi-static conditions. Their work is most relevant to our study here, since they considered multicomponent powder blends and showed that the powder mixture Bond number can satisfactorily describe the flowability of mixtures obtained by quasi-static testing. Here, we extend on their work for prediction of flow rules of powder mixtures under dynamic conditions, i.e. as a function of shear strain rate. In addition, we investigate different methods of averaging the Bond number as Capece et al. [15] only used the harmonic mean. Here, we define the powder mixture Bond number based on the arithmetic, geometric and harmonic means, where w ij is a weighting factor for the interaction between two particle types i and j. Bo ij is the granular Bond number between the two particle types i and j: where the adhesion force (F ad,ij ) is calculated from Eq. (7), characteristic particle weight (W ij ) is calculated using a harmonic mean approach based on the weight of two particles (W i and W j ), Given the fact that the van der Waals forces act on particles in close proximity and the interactions depend on the coordination number, the weighting factors in Eqs. (9) to (11) are based on the fractional surface area (f SA,i ) of each component (i and j), as shown in Eq. (14), The weighing factor can also be calculated based on the coordination number; however, in practice quantification of coordination number is very challenging. For a binary powder mixture, the mixture granular Bond number based on arithmetic, Eq. (4), geometric, Eq. (5), and harmonic, Eq. (6), averaging methods can be written as, where subscripts '1' and '2' refer to first and second component of the mixture, respectively.

Case studies
Four types of particles (i.e. non-adhesive, poorly adhesive, adhesive, and highly adhesive) are considered as mixture components to form a range of binary (see Table 2) and ternary mixtures (TMS1: Ternary mixture of non-adhesive component from BMS1 and adhesive components from BMS2 and BMS3). of the mixtures refers to binary, ternary, mixture, same size ratio, and different size ratio of the components, respectively. It is also worth mentioning in DEM simulations, the mixtures are generated such that all the components are generated randomly in the special domain. For binary mixture systems of BMS1 to BMS3, having a size ratio of unity, seven mixing ratios based on number density are considered (i.e. 0, 10, 25, 50, 75, 90, and 100% of the more adhesive component). For BMD1 to BMD6 cases, having a size ratio of 2, only three mixing ratios of 25, 50, and 75% of the more adhesive component are considered.
In the case of BMD1 to BMD3 the larger size component of the mixture is adhesive, whilst BMD4 to BMD6 the smaller size component is adhesive.
In the case of ternary mixtures (size ratio of unity for all components) three components are considered: highly adhesive, adhesive, and nonadhesive. Six mixing ratios of the components are considered as indicated in Fig. 11. For indicating which component is adhesive, reference is made to interaction Bond numbers Bo ij for each system as shown in Table 2.

Results and discussions
The approach adopted here is to link flow behaviour of powder mixtures to the single particle properties of their components (i.e. density, size, interfacial surface energy). Here the flow behaviour of powder mixtures is simulated as a function of the shear strain rate using the Freeman Technology FT4 rheometer. Initially, a range of DEM simulations are run for a single component powder bed to calculate the BFE for each component separately. When particles are made more adhesive, it is expected the value of the BFE to increase, due to the fact that adhesive particles do not flow well. However, the calculated flow energies, shown in Fig. 5a, suggest that more work is actually required for the impeller to penetrate into the particle bed for the free-flowing component of BMS1 case than the adhesive component of BMS3 case. Looking at the properties of individual components, this contradiction can be explained by the density difference between the two cases, where the free-flowing component of BMS1 and adhesive component of BMS3 cases have a density of 2500 and 1000 kg/m 3 , respectively. The latter low density is artificially assigned for numerical simulations to check the influence of density on expended work. This suggests for comparison of flowability of adhesive powders by FT4, the BFE should be normalised by the mass of the powder in the vessel. Here, we refer to the normalised BFE as 'Specific Downward Flow Energy (SDFE)', where a consistent trend can be observed as shown in Fig. 5b. As expected, more adhesive particles (larger granular Bond number) have larger specific flow energies as compared to those particles with smaller granular Bond number. Therefore, for analysis of flowability of the mixtures, the SDFE is considered in this study.
The three binary mixture systems (i.e. BMS1, BMS2, and BMS3) are considered to represent the extreme cases for the range of granular Bond numbers used in the simulation. The BMS1 is the base model system, where both components in the mixture have a density of 2500 kg/m3 and a random particle size distribution ranging from 850 to 1000 μm. Based on the values of interfacial surface energy chosen, the granular Bond numbers for the adhesive and free flowing components in the base model are 2.9 and 0.01, respectively. In the case of BMS2, the level of interfacial surface energy of adhesive component is increased by a factor of ten compared to the base model, leading to granular Bond numbers of 29.3 and 0.1 for the adhesive and non-adhesive components. For the BMS3 case the same interfacial surface energy as in the case of BMS2 is used, but the density of both components is reduced from 2500 kg/m 3 to 1000 kg/m 3 to increase the granular Bond number of adhesive component to 73.3, according to Eq. (7). As it can be seen from Fig. 6, in the case of BMS1, increasing the concentration of the adhesive component does not change the SDFE. On the other hand, in the case of BMS2 and BMS3, the SDFE does indeed increase with increasing the concentration of the adhesive components and the extent of increase is much higher in the case of BMS3 than BMS2. This is as expected since the granular Bond number of the BMS3 adhesive component is much larger than BMS2. For each data point in Fig. 6, the mixture Bond number is calculated based on Eqs. (6) to (8), and the results are shown in Fig. 7. It is evident that the harmonic mixture Bond number, shown in Fig. 7a, does not provide a good unification of SDFE. In contrast, both arithmetic and geometric mixture Bond numbers provide a good unification of the data points. When using the arithmetic Bond number (see Fig. 7b), the SDFE is initially constant around 4.1 mJ/g, corresponding to the non-adhesive mixture system, and starts rising when the mixture Bond number exceeds 10. A similar trend prevails for the geometric mixture Bond number case, as shown in Fig. 7c. The SDFE starts to rise when the mixture Bond number is in the range of 1 to 4. So the transition clearly depends on the definition of the average value. Nevertheless, the existence of a critical range at which the transition takes place is evident, and has practical significance in formulating powder mixtures.
For validation purpose, a series of experimental tests were carried out for different mixing ratios of BMS1 case, as shown in Fig. 8. A good agreement is observed for BFE as function of concentration of adhesive component in the case of BMS1 giving confidence in the conclusions derived from this study.
Further investigation is carried out for establishing the link between the mixture Bond number and SDFE for binary mixtures when the two components have different particle size ratios (i.e. BMD1 to BMD6 cases). This is done by considering a particle size ratio of two and making either larger size (BMD1 to BMD3) or smaller size (BMD4 to BMD6) components adhesive. Three levels of adhesion are considered here similar to BMS1 to BMS3 cases. When the smaller size particles are made adhesive, the ratio of adhesive force to gravitational force becomes larger, resulting in larger values of granular Bond number for the component under consideration. The relationship between the concentration of adhesive components and SDFE for binary powder mixtures with particle size ratio of two is shown in Fig. 9.
When the larger particles in the mixture are adhesive (i.e. BMD1, BMD2, and BMD3), there is a mixing ratio (50-50%) at which the maximum SDFE is expended. The underlying cause of this trend requires further investigation, considering a larger number of mixture ratios. This type of behaviour is also observed by Fulchini et al. [13]  Specific Downward Flow Energy (mJ/g)

Components:
NA: Non-Adhesive A: Adhesive HA: Highly Adhesive Fig. 11. Calculated SDFE for various mixing ratios of non-adhesive, adhesive, and highly adhesive components in a ternary powder mixture.
silanised glass bead surfaces with fine zeolite particles. They found a critical surface coverage of zeolite particles beyond which the flowability of the mixture decreased. The SDFE as function of mixture granular Bond number for the BMD1 to BMD6 is shown in Fig. 10. Similar to the powder mixtures with size ratio of unity, a good unification of the data points is achieved when considering the arithmetic mixture Bond numbers. The SDFE remains constant for the arithmetic mixture Bond numbers up to 10 for all the cases analysed. The unification of the data is poorer for geometric mean of mixture Bond number. Increasing the mixture Bond number beyond these values results in an increase in the SDFE, meaning the flowability of adhesive mixtures decreases with an increase in their Bond number. A similar analysis is also carried out for a ternary powder mixture with size ratio of unity to check the validity of our approach. The calculated SDFE for different mixing ratios of the components in the case ternary mixtures is shown in Fig. 11. It can be seen that the highest SDFE is expended for the largest concentration of highly adhesive component in the mixture and the lowest for the largest ratio of the non-adhesive component. In between, there are a number of mixtures for which the expended SDFE does not vary greatly.
The relationship between the mixture Bond number and SDFE for all mixture systems analysed here is compared in Fig. 12. It can be seen again that the harmonic mixture Bond number does not unify the data points and therefore is not representative of the bulk flow properties. Among the three averaging methods, the arithmetic mixture Bond number seems to give the clearest trend. Hence, it can be concluded that the arithmetic averaging is the most suitable method of averaging the mixture Bond number for prediction of the SDFE as measured by FT4 rheometer. It is worth noting that inter-particle and particle-wall friction coefficients have been kept constant here. A much wider range of single particle properties, including friction and shape needs to be analysed to arrive at a general correlation. However, the focus of this study is on establishing a methodology by which the mixture Bond number can be used to give an indication of flowability of powder mixtures at high strain rate conditions. The results show that there exists a narrow range of mixture Bond numbers below which the mixture behaves as free flowing.
When representing the mixture bulk cohesion by a mean Bond number, a complete random mixture is considered, with contributions from the Bond numbers of the components weighted with respect to surface area of the mixture composition. According to work reported by Berger et al. [42], bulk cohesion and friction are both influenced by strain rate. Therefore, a thorough analysis is needed to address the coupled effect of these two bulk powder properties, structuring of cohesive powders [41], and underlying reasons why arithmetic averaging gives the best unification.

Conclusions
A wide range of DEM simulations have been carried out in this study for relating the flowability of adhesive powder mixtures to the properties of particles constituting the components under dynamic conditions. Initially binary powder mixtures were considered and the size and density of all components were kept constant (i.e. BMS1 and BMS2 cases) whilst the level of interfacial surface energy was increased. Later, the size and interfacial surface energy of the components were kept constant and the density was varied (i.e. BMS3). For each of these cases, a particle size ratio of 2 was also considered (i.e. BMD1 to BMD6 cases). Finally, for the proof of concept, a ternary powder mixture was also considered.
The DEM results indicate that the flowability of powder mixtures can be linked to the level of adhesion of individual components in the mixture. The granular Bond number of individual components is used to calculate an average value for binary and ternary powder mixtures. It is shown that their flowability is described well by the mixture Bond number by considering a weighting factor based on fractional surface area of individual components in the mixture. Three averaging methods were investigated to calculate the powder mixture Bond number (i.e. harmonic, arithmetic, and geometric). In contrast to the previous work reported in the literature [11,12] indicating that harmonic averaging of the Bond number for the powder mixture can be used for calculating a representative powder mixture Bond number for describing flowability, our study shows that among the three averaging methods, the arithmetic mixture Bond number gives the best unification of data for a very wide range of Bond numbers. The dependence of flowability of powder mixtures on the properties of individual particles constituting the components, such as interfacial surface energy, size, and density, is well described in this way. The methodology developed here is helpful for describing the flow behaviour of powder mixtures, when characterising powder formulations.