Sensitivity analysis of the dynamics of fine and ultrafine particles using DEM

In this paper, we focus on particle–particle and particle–wall interactions considering tiny particle dimensions and the processes and phenomena arising from the contact dynamics. This has the important implication that the discrete element method can be used for large-scale computations as well as for tiny particles, i.e. particles with fine and ultrafine dimensions. Particular attention is paid to the granular cohesion dynamics where the particles interact prior to their physical stick. We investigate the sensitivity of the interactions, i.e. we assess how particle size distribution, frictional forms of particle–particle and particle–wall collisions and Van der Waals or liquid cohesive forces shape the particle motions. Through computations, we show how neglecting the above features influences computations of particle positions and particle linear and angular velocities over time.


Introduction
Discrete element method (DEM) [1] is one of the most popular research tools that models the dynamics of granular flows. Unlike the continuous approach [2], DEM distinguishes between trajectories and rotations over time for individual particles, where linear and angular velocities are recorded versus time. The main advantage of this approach is that it more accurately reflects the real dynamics of granular flows. However, its main disadvantage involves high computational time values, especially for large-scale computations where a massive number of particles are involved. The problem of high computational time becomes greater and is exacerbated when we need to consider the behaviour of fine or ultrafine particles. Commercial codes, for example, Altair EDEM [3], ROCKY-DEM Ansys [4], reduce high computations by using specific computer hardware such as multi-Graphical Processor Units as well as parallel computing. Another useful tool, MFiX [5], like the previous ones, also involves many hardware/software tricks that reduce computational time. The above codes are of great interest to the scientific [6] and engineering [7] communities. However, the codes listed above do not pay particular attention to tiny particles where their interactions involve not only collision-repulsive forces [8] but also cohesive forces [9] where the interaction begins before the physical contact. As the particle size decreases, we observe qualitatively different effects, i.e. the phenomenon of cohesion dominates during the interaction of the parti-cles. Such particle cohesion is of interest in many engineering applications, for example: in 3D jet printings [10], particulate matter in polluted air [11], fluidised bed [12,13], toner particles in laser printers and photocopiers [14] and a growing number of medical treatments [15,16]. Therefore, neglecting particle cohesion becomes a shortcoming in the fundamental approach for mathematical modelling of the granular dynamics of fine and ultrafine particle interactions.
In this study, particle-particle and particle-wall interactions play the dominant role, and many simplifications given in the modelling process could not reflect the real dynamics of particle motions. These improper dynamics are highly visible in particles having very low dimensions-so-called fine or ultrafine particles, and for high particle concentrations in space where multi-particle collisions occur. According to our previous experience [17,18], we show how simplifications in the mathematical modelling have a strong influence on the dynamic's particle motions, i.e. drastic changes of particle trajectories and particle velocities. Here, we refer to these simplifications or/and extensions as sensitivity analysis.

Mathematical illustration
According to own investigations [17], we assume population of particles where n indicates the total number of particles which one considers in computer simulations. The individual motion of the centre mass of a particle "k" is described by the following system of ordinary differential equations: for particles moving individually, i.e. without particleparticle or particle-wall interactions and considering particle-particle or/and particle-wall collisions. The above expressions describe the motion for the particle k, where: m k -mass of the particle, J k -mass moment of inertia, x k -position of the centre of mass, ω k -angular velocity, F l -force acting on the particle, M l -torque acting on the particle, P rep j (k) -repulsive/friction forces resulting from a collision, P attr j (k) -cohesive/friction forces resulting from short-range interactions, M rep j (k) -torques resulting from repulsive/friction forces, M attr j (k) -torques resulting from the action of cohesion/friction forces. More explanations can be found in our own work [17]. Figure 1 illustrates the concept of the interaction of two particles. We understand the interaction between particles in a broader sense, i.e. particles can collide as well as interact with each other due to cohesion. This is an extension of the concept of the "soft sphere method" [19] with interactions resulting from the interaction of short-range forces and occurring before the collision itself. Figure 1a shows the radius of the r k particle and the virtual radius of the r attr k particle resulting from the forces of attraction/cohesion. Figure 1b explains the local particle collisions as particle overlapping ζ j (k) as well as local coordinates. We will form a cohesion and contact forces system in these local coordinates.
In the existing literature, the rebound force model of spring-damping [20] is assumed. Given the lack of multi-particle contacts [21], we assume fractional rebound force model [17], where its component in the direction of ζ j (k) is defined by local coordinates ζ j (k) denotes left-sided fractional derivation defined in the Caputo sense [22], t be j (k) and t en j (k) represent times where the collision begins and ends, α j (k) is the fractional order of the fractional derivative, ζ j (k) is the particle' overlap as shown by Fig. 1b, C j (k) and K j (k) reflect the material properties of two colliding bodies, i.e. damping and spring coefficients.
Regarding the cohesive forces, we use Rump's formula [23] for the Van der Waals force P attr−vW where ζ vW j (k) is the overlap of particle cohesion as shown by Fig. 1, C 5 is the smallest distance between contacting particles and according to [24]  is the mean radius of the colliding particles (see Fig. 1),Ĥ j (k) = H k H j (k) is the average value of the Hamaker' constants for the two interacting particles andR j (k) = R k R j (k) represents the average roughness value for the two interacting particles.
Hamaker constants for selected granular materials can be found In [25]. Another cohesive force arises from the surrounding moisture of the particle's surface. According to [26], we use the following expression P attr-LB where ζ LB j (k) is the overlap of particle cohesion as shown is the mean radius of colliding particles (see Fig. 1), ς j (k) denotes surface tension and the quantities r I j (k) , r I I j (k) are functions of the wetting angle ψ j (k) and are related by the following expressions where To find the wetting angle ψ j (k) , we need to solve the following equation is the quotient of the volume V LB j (k) of the moisture layer surrounding the particles to the volume of the particlesV j (k) = V k + V j (k) . Before applying Eq. (5), the nonlinear expression (8) must be solved to obtain wetting angle ψ j (k) . As noted in [26], the wetting angle is ψ j (k) ∈ 0, Π 2 . However, the frictional force comprehensively described in the author's work [17] includes the normal force, in which all forces acting along the normal direction are taken into account, including the cohesive forces described above. This increases the friction force.
All other forces and torques derived from these forces, including those involving various forms of friction, are presented in the author's study [17] and omitted in this paper.
The systems of Eqs. (1) and (2) are strictly nonlinear and cannot be solved except by numerical methods. The Runge-Kutta-Fehlberg method of 4-5 order [27] was used to solve these equations. From the definition of this method, it follows that it is always convergent and has an accuracy of O (Δt) 5 . Based on [27], one can estimate the calculation time step as a function where err is the user-assumed precision of the computation. The linked cell method [28,29] was used to determine the start time of t be j (k) and the contact duration t en j (k) − t be j (k) . According to Allen [28], this method is simple, discrepant and the most effective. In jumping from expression (1) to (2) and vice versa, the leapfrog method [28] was used.
In conclusion, it is stated that the proposed mathematical notation (1)-(9) together with our own precise description [17] allows us to carry out computer simulations of the flow dynamics of granular material, considering the phenomena of cohesion and friction.

Results
We rely heavily on previous experience when assessing the uncertainty of the calculations. Our previous studies [17] show, by comparing four different repulsive force models, that the fractional repulsive force model given by Eq. (3) is a reliable model. This model reflects the viscoelastic properties for the contacting bodies, the contact start and end times, as well as the coefficient α j (k). This coefficient is responsible for the degree of conversion of the impact energy into rebound and dissipation energies. Looking at Eq. (3), it can be seen that for α j (k) = 0 all impact energy is converted into rebound energy, and formula (3) takes the form . This is well known in the literature as a linear-spring model. On the other hand, for α j (k) = 1 all the impact energy dissipates and formula . This is also well known in the literature as a linear-damping model. It seems to be that α j (k) is like the restitution coefficient, which is frequently used in the hard sphere event-driven molecular dynamics [30]. As noted in my own study [17], α j (k) is called the conversion degree, which has more significant meaning for the roughness of contacting particle surfaces. In this paper, we generally rely on previous experience and do not make a local scale comparison with different repulsive force models and friction force forms, which can be found in [17]. Comparison on a global scale requires many assessments of particle motion dynamics and an assessment of the effects of self-scaling. This is a very difficult task, and in this work, we will limit ourselves to the dynamics of particle outflow from a container. In this regard, we stretch our previous considerations given in work [18], where the outflow dynamics of 3000 particles were considered, to the outflow dynamics of 4000 particles from the container. In this case, we refer to comparing the experiment results with computer simulations for invariant material data and increasing the quantity by only 1000 grains. Figure 2 shows such a comparison, i.e. the input and output effects. Experimental studies were carried out according to the following parameters given by Fig. 2a. We consider two cases of the particle size distribution that reflect the height of the particle bed in good agreement with experimental data. In the second case, we achieved greater compatibility of the computer simulation with the experimental data, i.e. the time of emptying the grains from the tank was 1.08 s for the computer simulation and 1.16 s ± 0.04 s for the experiment. In conclusion, we note that it is sufficient to record two observations to check the DEM model on a global scale: the height of the particle bed and the time it takes for the particles to flow out of the container.
Initial comparisons of computer simulations with the results of the experiment allow a smooth transition to the sensitivity analysis. Thus, we start all over again. Figure 3 shows the influence of the conversion degree α value (3) on the dynamics of 3 particle collisions in 1D. In the strong repulsion regime, (α < 0.3), the particle positions are observed to be very far apart. On the other hand, in the weak repulsion regime (α > 0.5), the particle positions are not far apart. Extending our analysis, we show how weak and strong repulsion affect the dynamics of particle collisions. Figure 4 presents a similar situation to the previous one but in 2D. Here, we observe that the trajectories of the particles are almost linear in a strong repulsion regime. However, the particle trajectories change as the α conversion ratio approaches unity. Thus, the particle trajectories are curve-linear in a weak repulsion regime. Consequently, formula (3) of the repulsive force makes it possible to control the dissipation energy during the collision of single particles as well as during the collision of multi-particles with each other. The nonlinear trajectories of the grains extend their contact time. On a macro scale, this manifests itself in the formation of grain clusters, which affects the dynamics of the movement of the entire population of grains.   Extending our case study, we focus on assessing the effect of friction on particle trajectories. All forms of friction (sliding, rolling, torsional) are included in [17]. Here, we limit our considerations to the sliding/rolling of particle-particle or particle-wall friction. Figure 5 shows two scenarios where inter-particle friction is turned on and off. Here, we consider the 500 particles that lie at the bottom of the container forming a layer of particles. The particle dimensions are as noted in case II in Fig. 2. The crosses on the particles indicate the positions of individual points on the particle surface and show the particle rotation resulting from the friction force. We deliver air from the bottom of the container at an average velocity of 1.5 m/s and then record images of the positions of individual particles in a 3D space. However, Fig. 5 shows images of the particle positions flattened down to 2D. The left image in Fig. 5 shows particle positions where we ignore interparticle friction. It can be observed that the crosses on the particles do not change their position, i.e. the particles do not rotate. On the other hand, the right image shows the particle positions in which we consider the inter-particle friction. It can be observed that the crosses on the particles change their position, i.e. the particles rotate. Without taking friction into account, the dynamics of particle motion appear to be a splash type dynamics, i.e. the particles rebound more strongly off each other. However, the particle dynamics are more consistent and compact when inter-particle friction is taken into account. In the images presented above, one may observe the formation of entirely different structures. Taking the friction process into account during the collision of particles converts the kinetic energy of the collision into the kinetic energy of translational and rotational forms of motion. This well-known and obvious fact significantly changes the image of particle movement, which, considering the collisions of many particles with each other, has a significant impact on the formation of particle clusters and other forms of colliding particles.
The next simulations show the interaction dynamics between the particles, considering the Van der Waals cohesion at the nanometric scale. Figure 6 presents 4 scenarios of the formation of particle clusters depending on the initial velocities of these particles and the energy conversion degree α value. It can be observed that the values of the initial particle velocities have a significant effect on the formation and spattering of particle clusters only for the strong replication regime α = 0.09. However, we do not notice a significant influence of the initial particle velocities on the distribution of particle clusters for the weak reflections regime, α = 0.85, which are enhanced by the Van der Waals cohesion force (4). This analysis has a significant impact on 3D printing [10].
Another example presents computer simulations of particle elutriation dynamics under the influence of mean air velocity. Figure 7 shows 3 scenarios depending on the type or absence of cohesive force. Here, three different behaviours of the same population of particles acting under the same conditions can be seen.
The image on the left shows floating particles where the mutual cohesion between the particles is switched off. In this case, the particles are fully raised and do not form any particle clusters. The middle image shows the particle float, considering only the cohesion in the form of the Van der Waals force (4). We notice that the particles float more slowly and form small clus- ters. However, their lifting force is so great that they all move up the container. The picture on the right shows the behaviour of the particles, taking in to account the Van der Waals cohesion force (4), as well as the capillary cohesion force (5). As can be seen, almost all the particles stay at the bottom of the container forming a sticky cluster, i.e. a coherent fluidised layer. The air only filters the resulting layer of agglomerated grains. The individual agglomerates are detached from this fluidised bed and lifted into the top of the container.
By extending our considerations and considering the scenario presented in Fig. 2, we can simulate emptying the container of nano-particles, where we can also turn off cohesion or gradually turn it on by increasing the value of the Hamaker constant as we introduce capillary cohesion. Figure 8 presents 5 scenarios. In this case, we can control the cohesion value, which can lead to a scenario of total outflow of particles from the container (left image-no cohesion) or a slow outflow of particles from the container (middle images-increasing the Van der Waals cohesion level) or to a complete blockage of the outflow of particles from the container (image on the right-capillary cohesion).
The latest analysis assesses the effect of capillary cohesion on the Brazil nut effect [31]. Figure 9 shows 4 images of a vibrating container with a population of fine particles and one coarse particle. The particle data is the same as in Fig. 5 The large particle size is 10 mm. The container was vibrated at 100 Hz but with different amplitudes. The first two images on the left show the behaviour of the particles for the container vibration amplitude A = 3 · 10 −4 m. The following two images show the behaviour of the particles with an increased container vibration amplitude A = 4 · 10 −4 m. Here, we can see that the non-cohesive particles more easily push the large particle into the top of the container, which perfectly illustrates the Brazil nut effect [31]. However, considering the capillary cohesion between the particles significantly reduces the Brazil nut effect, i.e. it is more difficult for the large particle to reach the top of the container.

Conclusions
The use of the discrete element method (DEM) in computer simulations of granular matter dynamics provides valuable information about the behaviour of particles Fig. 9 Dynamics of the Brazil nut effect considering a cohesion-less system and capillary cohesion for the 5 · 10 −1 s observation and significantly strengthens our knowledge in this field. However, before using the DEM, it should always be verified both on a local scale, e.g. by comparing various forms of rebound strength, and on a global scale, e.g. by the dynamics of emptying the particles from the container. In this study, we used the fractional repulsive force model to illustrate strong and weak repulsion regimes of colliding particles. Our sensitivity analysis made it possible to turn the friction force on or off, which significantly distorted the positions and velocities of the particles. In addition, we analysed the sensitivity of granular matter on a nanometric scale, which allowed us to present the effect of mutual cohesion between particles on the behaviour of the entire particle population. Depending on the type of cohesion (Van der Waals cohesion or capillary cohesion), various particle motion scenarios were observed, e.g. the formation or disintegration of particle agglomerates. Depending on the size scale considered for the adopted grain population, we suggest that a sensitivity analysis be performed each time to assess the qualitative and quantitative influence of processes and phenomena important from our point of view on the behaviour of granular matter. More advanced sensitivity analysis can be found in [32], where an iterative Bayesian filtering approach was used for automated calibration. This study is limited to the local or small-scale behaviour of granular matter. The properties of particles in the nanoand microscale can encompass countless regimes of granular matter to illustrate the properties of bulk material. Note the significant influence of the nano/micro scale on the actual state of the bulk material. However, this topic goes beyond the purpose and scope of this study. Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest The author declares that he has no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix
Let us complete some definitions that fully explain the model. According to Fig. 1, we calculate the virtual pre-contacts between particles k and j (k) as where · denotes a norm calculated from the relative coordinate x j (k) − x k . If there is cohesion between particles then ζ attr j (k) ≥ 0. We also calculate the particle collision as a virtual overlap If there is collision between particles then ζ j (k) ≥ 0. We assign global coordinates (x 1 , x 2 , z 3 ) and local ones ξ j (k) , η j (k) , ζ j (k) The index k is the number of the particle whose centre-of-mass movement is described by Eqs. (1) and (2). Thus, the set (k, j (k)), where j (k) is variable, denotes a pair of contacting particles. Let us define the unit normal vector as that connects the particles' centres of mass. On a tangent plane, we define tangential unit vectors as and where ||·|| x 1 x 2 denotes the norm which is calculated in the tangent plane. Given the above expressions, the unit vector matrix for the two interacting particles can be built as Assuming spherical forms of contacting particles, we obtain a point determining the temporary mass centre of the overlapping particles (10) and (11) as (16) for particles in the pre-contact regime, i.e. when Eq. (10) is less than or equal to zero and for colliding particles, i.e. when Eq. (11) is less than or equal to zero. We expand the local coordinate system ξ j (k) , η j (k) , ζ j (k) for the above points as illustrated by Fig. 1b. For ζ attr j (k) = 0, we obtain the begin of pre-contact and Eq. (16) simplifies to Also, for ζ j (k) = 0 we obtain the start of the collision and Eq. (17) simplifies to Next we define the relative particle velocities. At the point x C−attr j (k) , we have w attr j (k) = w lin−attr where w lin−attr is the relative translational velocity of the contacting particles (sometimes called the sliding velocity in the contacting regime) and w rot−attr is the relative rotational velocity of the contacting particles in the contacting regime, where ω attr k = e attr j (k) · ω k and ω attr j (k) = e attr j (k) · ω j (k) are angular velocities defined in the local coordinates in the contacting regime and a attr k , a attr j (k) are branch vectors connecting the centres of contacting particles x k , x j (k) with the temporary point x C−attr j (k) given by formula (18). We define such vectors in local coordinates as follows a attr and a attr j (k) = 0, 0, x C−attr In the collision regime, i.e. ζ j (k) ≤ 0 we define once again the relative particle velocities at the point is the relative translational velocity of the colliding particles (sometimes called the sliding velocity in the colliding regime) and is the relative rotational velocity of the colliding particles in the colliding regime, where ω k = e j (k) · ω k and ω j (k) = e j (k) · ω j (k) are angular velocities defined in the local coordinates in the contacting regime and a k , a j (k) are branch vectors connecting the centres of colliding particles x k , x j (k) with the temporary point x C j (k) given by formula (19). Such vectors we define in the local coordinates as follows and Next we define a set of long-range interacting forces F l , for l = 1, ..., 3 as: -Gravity, for l = 1 where m k is the particle mass, g = 9.81 m/s 2 denotes the gravitational constant, -Hydraulic resistance force, for l = 2 where V k = π d 2 k /4 is the particle volume, p represents fluid pressure -Drag force, for l = 3 where β denotes drag coefficient and U g is the gas velocity.
Simplifying the mathematical illustration, we assume all torques issued from the long-range interaction forces are zero, i.e. M l = 0 for l = 1, ..., 3. A set of repulsive forces P rep j (k) involves all particleparticle collisions, but the force components include all interactions acting in the normal direction and the tangent plane. Thus, the symbolic notation P rep j (k) has the following meaning where · t denotes a norm calculated only in the tangent plane ξ j (k) , η j (k) and P sta where the symbolic notation e T j (k) indicates transposition of the matrix of unitary vectors e j (k) and several forces are defined in the local coordinates ξ j (k) , η j (k) , ζ j (k) as are components of static friction force, T sli ξ j (k) , T sli η j (k) represent components of the coupled sliding-torsion friction force and T rol ξ j (k) , T rol η j (k) denote components of the rolling friction, whereas P rep ζ j (k) represents the component of repulsive force along the ζ j (k) direction and was defined previously by Eq. (3). The static friction force is defined in the local coordinates as where μ sta j (k) denotes the static friction coefficient determined between a pair of colliding particles x k , x j (k) . Following to [33], we propose a mathematical illustration for a joint sliding-spinning motion of two-colliding particles in the form of sliding-torsion friction force where μ is the joint-friction coefficient and μ sta j (k) , μ kin j (k) are static and kinetic coefficients of friction, C 1 is an empirical constant, · t represents the norm of a vector calculated in the local-tangent plane. The function F λ j (k) issues from spinning motions of a pair of colliding particles. Farkas et al [33] define this function as where E I , E I I are elliptic functions of the first and second kind and For F ( ) = 0, we have torsion friction without sliding friction but F λ j (k) → ∞ = 1 we obtain sliding friction without torsion friction. With regard to the rolling friction force P rol j (k) , we have considered the two-dimensional form in the paper [34] and extend it to the three-dimensional one as where N j (k) is the normal reaction being perpendicular to the local-tangent plane and has the following form N j (k) = 0, 0, P rep and s j (k) represents the rolling friction arm and is expressed by Next we consider a set of attractive/cohesive forces P attr j (k) . Taking into account the local coordinates ξ j (k) , η j (k) , ζ j (k) and using the transformation matrix e attr j (k) given by Eq. (15), we have P attr j (k) = e T j (k) · P attr j (k) and the cohesive forces defined in the local coordinates look as follows P attr j (k) = P attr−vW j (k) The above formula includes the vectored form of two types of cohesive forces, where P attr−vW j (k) = 0, 0, P attr−vW is the Van der Waals cohesive force where its component P attr−vW is expressed by Eq. (4) and is the capillary force where its component P attr−L B is expressed by Eq. (5). Here, we finish mathematical illustration of all forces acting on contacting/colliding particles.
Next we define torques resulting from the forces. Thus, we have where M sli j (k) is the coupled torsion-sliding torque and M rol j (k) represents the coupled torsion-rolling torque. Using transformation between global and local coordinates, we obtain The sliding torque M sli j (k) is defined in the local coordinates as The rolling torque M rol j (k) is determined by the following formula The torsion torque M tor j (k) = 0, 0, M tor T is responsible for particle spins. Its component M tor ζ j (k) is defined in the work [33] as where the function T λ j (k) is defined as