Modelling Pasty Material Behaviour Using the Discrete Element Method

Mixtures of a fine-grained dry solid and a liquid, with a higher solid content in the mixture, show paste-like behaviour. In many technical processes, pasty materials are handled in large quantities. Pasty Materials show the same characteristics as Bingham Plastics, behaving like a rigid body but flowing like a viscous fluid under a certain stress level. This behaviour is due to attracting forces between the particles, resulting from the capillary pressure and the surface tension of the liquid, which forms individual capillary bridges or capillary bridge clusters between the solid particles. The behaviour of granular material can be represented in a discrete element method (DEM) simulation. The calibration of the simulation parameters is achieved by comparing laboratory tests, which reflect a typical material characteristic, with the results of calibration simulations. In this project, several DEM contact models, describing attractive forces between particles depending on the distance between them, were analysed and assessed based on their ability to display the pasty material behaviour of a fine-grained solid water mixture with a high water content by comparing the simulation results to a slump test. The most promising contact model was then optimised to enable a minimal computing time for the simulation of bigger technical processes. Many existing contact models also consider attractive forces between particles (e.g. JKR cohesion) but are based on different physical effects. For this reason, the contact models assessed in the course of this project are, in general, based on the capillary effect.


Introduction
Pasty materials have a wide range of application in technical processes. From assembly pastes, grease-like lubricants with high solid content and good thermal stability, to the mixing process of tailings, a waste product after separating the valuable fraction from the uneconomic fraction of an ore, with crushed waste rock for better waste disposal [1]. The most time and cost-efficient way of optimising a technical process is to simulate it using computer modelling methods suitable for the process.
In DEM simulations, granular media is represented with particles that interact with each other according to predefined contact models. Discrete element simulations have a wide range of application ranging from process engineering (e.g. various mixing processes or chemical processes) to mining engineering, in which, for example, the conveying process of bulk goods (e.g. including breaking behaviour of material) is simulated. The individual particles are approximated by spheres in many DEM simulations in order to combine a short computational time with a very accurate result. The interacting forces between individual particles or between particles and system components, like belts and chutes, are calculated by evaluating the overlap of the particle with other particles or geometries of system components. The interaction forces resulting from this overlap in accordance with the defined contact models are evaluated in discrete time intervals of the simulation, so-called time steps.
Besides DEM simulations, Smoothed Particle Hydrodynamic (SPH) simulations also show promising potential for the simulation of pasty materials. In SPH simulations, the flow of a continuum medium is represented by discrete elements. In comparison to the DEM, the behaviour of an individual element is influenced by not only adjacent elements but the influence of more distant ones is taken into account with a decreasing effect, following a Kernel function [2]. Although this simulation method is most commonly used to simulate pure liquids, it may also be calibrated to simulate highly viscous materials to replicate the behaviour of slurries. However, paste-like media can form certain stable states, which are not possible for pure liquids or pasty media with high water content (highly viscous but still liquid-like). As such behaviour would be relevant for paste-like materials, but cannot be replicated with SPH particles, which do not stop their movements (even after a long period of time) to form said steady states, this method is not suitable for this application and will consequently not be considered any further in the course of the project.
Capillary bridges form between particles of a granular medium when water is added. Due to the surface tension of the liquid and the capillary pressure, attractive forces act between particles [3,4]. The capillary pressure, also known as Laplace pressure, is the pressure difference between the surrounding air and the liquid pressure. As the distance between the particles increases, this attractive force, the capillary force, continues to decrease until it abruptly drops to 0 at a certain distance, the so-called rupture distance [5]. The higher the water content, the more capillary bridges connect to form liquid clusters that eventually fill out all the pores between the solid particles with some entrapped air bubbles [5]. If this state (capillary state) is reached, the fluid pressure is still lower than the surrounding air pressure [6], and the addition of further liquid leads to an increase of the liquid pressure, and therefore to a rapid decrease of the cohesive force [5]. The final state, in which the liquid pressure rises to the surrounding air pressure, is called the slurry state. The dependence of the cohesive force between the particles on the saturation content and the capillary states is shown in Fig. 1.
The capillary force decreases with an increase of the distance between the particles between which the capillary bridges or clusters exist until it drops abruptly to 0 at the rupture distance. This is shown in the example of the capillary bridge theory of Lambert in Fig. 2. The exact course of the function of the capillary force depends on the saturation state but always decreases over the distance [5,8] The simulation of pasty materials in DEM is still quite unexplored. Some capillary based contact models are already included in certain DEM software. Although those contact models are based on the interaction between particles connected by capillary bridges or capillary clusters, measured in laboratory tests, they do not represent the bulk behaviour of a large amount of pasty material in a suitable way, even though the correct physical parameters, e.g. the surface tension and density of the liquid are used. This is due to the fact that these contact models describe the behaviour of isolated capillary bridges [3,10] or small capillary clusters [5,6].
In this project, a method was developed to simulate the bulk behaviour of pasty materials with a water content low enough to still form specific steady states under a certain stress level.  Example of the capillary force over the inter-particle distance depending on the amount of liquid within the capillary bridge (V), cf. Lampert [9] Test In order to calibrate the parameters for particle-particle interactions and the interactions between particles and system components in the simulation, calibration simulations are carried out in which a laboratory test is replicated. This laboratory test is selected based on its ability to depict the characteristic behaviour of the material. Since a standardised test, a so-called slump test, already exists to check the desired viscosity of concrete on construction sites, which reflects the mixing ratio of cement and water, a slightly modified experiment is used to calibrate the simulation. Although the slump test is initially used in civil engineering and has already been simulated on several occasions using the DEM [11], it is used in this project as a relatively simple and fast test for the study on contact models, with the intention to apply that contatact model on industrialsized simulations in the field of conveying technology. Furthermore, such pasty material characteristics are typically not of a common kind for conventional conveying technology applications. The experiment is divided into two steps that are visualised in Fig. 3.
In the first step of the experiment, shown in Fig. 3a, the truncated cone-shaped container is raised at a fixed speed, given in Fig. 3a as v , from the base plate after filling it with the pasty material. Both the container and the base plate are standardised, according to DIN 12350-5.
After the material has flown out of the container and comes to rest, the remaining height of the slumped material, dimensioned in Fig. 3a as h 1 is measured. Due to the particle size in the simulation, the front of the slumped material is not as flat as in reality. This is why the remaining height is measured and not the diameter, as suggested in the norm. The result of the laboratory test is then compared to the result of the calibration simulation. If the dimensions of the slumps deviate above a tolerable level, the calibration simulation is run again with altered parameters. The influences of the simulation parameters on the simulation results are discussed further on. A comparison of the results is shown in Fig. 4 (comparing the slump shapes of the laboratory tests and of the simulation).
In the second step of the experiment, shown in Fig. 3b, the base plate on which the entire test is carried out is lifted as far as possible and then dropped. The base plate, standardised according to DIN 12350-5, is equipped with hinges leading to the plate being only liftable from one side, up to a maximum angle indicated as in Fig. 3b, which results in an asymmetric slump of the material upon impact as seen in Fig. 5. This step is carried out only one time and not several times as suggested in the norm so that the contour of the slumped material remains intact. After the material has come to rest after the impact, the remaining height of the slumped material and the asymmetric slump of the plateau, indicated as h 2 and in Fig. 3b, and not the diameter, is measured again. With this step, the impact behaviour of the material is determined, which is particularly interesting for future simulations on a larger scale, e.g. chute simulations, which are used for designing material handling conveying systems. As in the first step, the results of the laboratory test are once again being compared to the results of the calibration simulation. A comparison of the results is shown in Fig. 5.
In order to reproduce this experiment as precisely as possible in a calibration simulation, the behaviour of the outflowing material from the container is analysed. It was found that the material moves forward in a "rolling" motion like a breaking wave. This behaviour is observed by focusing on small air bubbles on the surface of the outflowing material, as can be seen in Fig. 6. The bubbles on top of the spreading flank move at a higher speed than the material in contact with the base plate since they migrate from the top of the flank over the vertex until they disappear under the front of the flank. It is therefore concluded that the material Fig. 4 Comparison of the results from the first step of the slump test (Fig. 3a), between laboratory (top) and DEM simulation (bottom) results

Fig. 5
Comparison of the results from the second step of the slump test (Fig. 3b), between laboratory (top) and DEM simulation (bottom) results

Fig. 6
Analysation of the outflow during the lifting process of the bucket; focusing on the path of an air bubble on the surface of the material as a reference for the laboratory setup (top), in comparison to the motion path of a tracked particle in the simulation (bottom) in direct contact with the base plate moves at a slower pace due to adhesive forces between the pasty material and the base plate.

Contact Model in DEM
The DEM software ThreeParticle/CAE is used for the simulations [12]. The particles with which the pasty material is simulated are approximated as spheres, with a physical radius and a collision offset which combined result in the contact radius, in order to save computational time. Therefore the particles can overlap each other in two distinctive ways. If the contact radii but not the physical radii of two particles overlap, a negative overlap is detected by the software, as illustrated in Fig. 7a. A positive overlap is only detected if the physical radii of two particles overlap, as seen in Fig. 7b.
Negative forces act attracting between two particles, and positive forces act repelling.

Particle Forces in DEM Simulations
In DEM simulations, the displacements of particles between discrete timesteps are calculated from the forces acting on said particles. Those forces can result from a contact of the particle ( F Contact ) with another particle or with a system component, e.g. the wall of a chute, or are general forces F General from an additional source, like gravitation or a force field; see Eq. (1).
Contact forces (resulting from an interaction) are further divided into master forces ( F Master ) and slave forces ( F Slave ), as according to modern DEM software; explained with Eq. (2).
The master contact force is calculated from the master contact model and represents the primary force in normal and tangential direction acting between particles in contact (see Eq. (3); cf. [13]), whereas the slave force is calculated by an additional slave contact model to add additional forces, which allow the consideration of cohesive, adhesive or in this case capillary effects (as can be seen in Eq. (4)). In this project, an additional cohesive force at a positive particle overlap is added to the capillary force at a negative overlap for reasons explained further on.

Capillary Forces at a Negative Overlap
Capillary forces are only calculated if a negative overlap is detected between two particles. In this case, the collision offset represents the amount of liquid on the surface of a particle which leads to the formation of a capillary bridge or a liquid cluster when two or more particles get into contact.
The capillary force between two interacting particles is calculated with a mathematical function that should be as easy as possible to save a lot of calibration time. It has to be noted that the mathematical function of this contact model does not depend on values measured in the laboratory, such as the surface tension of the liquid or the rupture distance of the capillary bridges, as it is the case with conventional capillary bridge contact models [8]. Although the course of the capillary force over the negative overlap used in this contact model resembles the force curve of a real capillary bridge, it is still an approximation in order to reproduce the bulk behaviour of the pasty material as accurate as possible, to enable an easy calibration of the simulation parameters and to avoid a discretisation effect at simulations with an economic computation time, which will be discussed later in detail. The considered mathematical functions are shown in Fig. 8.
The mathematical function of the capillary force between particles with a negative overlap is chosen as a polynomial function of the second degree and is computed according to Eq. (5): with 0 ≥ d ≥ d c , and where d is the overlap between interacting particles, d c is the rupture distance, and a , b , and c are coefficients responsible for the force-over-distance characteristics, e.g. calibrated by means of the slump test. This decision is based on several calibration simulations carried out with the different models seen in Fig. 8. The chosen function has distinguished itself from the other ones by delivering exact simulation results, being easy to calibrate and being little susceptible to discretisation errors while still resembling the course of a capillary force function based on physical parameters. For better applicability in technical processes, the polynomial function can also be made dependent on the water content of the material. That way, the drying process of the material can also be taken into account in the simulation. The parameters of the function, and therefore the exact course of it, are calibrated in the calibration simulations.

Capillary Forces at a Positive Overlap
If a positive overlap is detected, the contact force between two particles is calculated with the standard Hertz-Mindlin master contact model. Normal and tangential spring-back forces are calculated depending on the overlap as well as damping forces depending on the relative motion of the particles.
The spring-back force of the Hertz-Mindlin contact model results in the particles being repelled back in the event of a collision with an object or other particles, which does not correspond to the behaviour of a pasty material. As it is observed during the second step of the slump test, the  material spreads out laterally instead of being repelled back vertically. In order to counteract this behaviour in the simulation, the Hertz-Mindlin contact model is combined with a constant cohesion slave contact model, which opposes the spring-back force in the normal direction, as seen in Fig. 9.

Discretisation Effect in DEM
The force between two interacting bodies in a DEM simulation is calculated in discrete timesteps and not continuously. In most simulations, the ideal time step corresponds Fig. 11 Visualisation of influences of the time step size on the discretisation effect (analysed in a test simulation); comparing a polynomial model (a) to a polynomial model with a timestep half the size (b), and further to a polynomial model with the same timestep as in (a), but with twice the rupture distance (c) to approximately 20% of the Rayleigh time step, which is dependent on the material data and geometry of the particles. The larger the time step, the shorter the simulation time since the interacting forces between every interacting particle are calculated not as often during the same amount of simulated time. However, one cannot blindly rely on the result of this approximated time step because even though the approximation with the Rayleigh time step proves to be correct in most cases, it can also lead to the following error.
In the case of two particles moving at a high relative speed towards each other, with the contact model between those particles being valid in only a very small range or the contact force changing abruptly over a short distance, a big timestep may lead to an incorrect calculation of the interacting course of the contact force between those particles over distance and time. The more often the capillary force is calculated from the contact model, and the smaller the difference between the values of the calculated force between time steps, the more accurate the representation of the contact model. Depending on the application, the impact of this effect needs to be evaluated. In a chute simulation, for example, particles very rarely exceed relative velocities of 0.5 m/s in relation to their direct neighbours. A trivial experiment with two particles moving away from each other at a relative velocity of v = 0.5 m/s and the timestep being 4e-5 is shown in Fig. 10. While the capillary forces calculated in the discrete time steps of this experiment (visualised as 'x' points in Fig. 10) with a common capillary bridge model do not represent the actual course of the force, the contact model proposed in this project represents the course of the actual force much more accurately.
A comparison of Fig. 10a and b show that the contact model based on a polynomial of the second degree is more resistant to those errors than a common capillary bridge model.
To minimise the errors caused by this effect as much as possible, the coefficients of the mathematical function of the contact model can be altered. This is visualised in Fig. 11, where the same trivial experiment, as described in Fig. 10, is carried out with several versions of the contact model proposed in this project. The influence of the time step size on the contact model with the same mathematical function can be seen by comparing Fig. 11a and b. A smaller time step would minimise this effect, as seen in Fig. 11b, but in industrial-sized simulations, a contact model applicable with a bigger timestep is definitely preferred due to the reduced simulation time. By increasing the scope of the contact model, meaning an increase of the rupture distance as seen in Fig. 11c, requiring a bigger collision offset (which is not typical for general capillary bridges in physical terms), the contact model becomes even more resistant to this effect and thus more suitable for applications where higher relative velocities between particles occur. A good approach for a contact model that is applicable for as many cases as possible in combination with a short computational time is based on a function with a medium rupture distance and a rather flat course of the capillary force over the overlap.
The described aspects must further be taken into account when choosing a contact model for a specific application in order to combine the biggest possible time step with the shortest possible computational time per time step to improve efficiency, especially regarding industrialsized applications.

Parameter Tuning
During the calibration simulation, the simulation parameters are varied, and their influence on the calibration simulation, regarding, e.g. the height of the slump, the discharge behaviour and the impact behaviour, is observed until the results of the simulation match the test carried out in the laboratory.
In the course of this project, the determined parameters for the particle interactions in the DEM simulation, e.g. friction coefficient, cohesion factor etc., are in a value-range typical for DEM simulations. For comparative purposes, some general material properties that are used in the simulations are shown in Table 1. Also, a typical particle size distribution (PSD) is used in the simulation, where the radii of the slurry particles range from 2.5 mm to 3 mm.

Maximum Capillary Force
Increasing the maximum capillary force increases the cohesion of the particles, resulting in a smaller diameter and a larger remaining height of the slumped material. If the maximum force is too low, the plateau of the slumped material from the first step of the laboratory test no longer appears flat in the simulation since the edges of the plateau slide further down due to a lack of cohesion.

Static Friction Coefficient
The sliding properties within the material can be influenced by changing the static coefficient of friction between the particles. The higher the value of the coefficient, the greater the friction forces between the particles moving past each other, thus increasing the viscosity of the simulated material. By increasing the value of the coefficient, the diameter of the slumped material also decreases and the remaining height increases. In order to achieve the exact height of the slumped material measured in the laboratory test, an optimal compromise between the static coefficient of friction and the capillary force must be found.

Rolling Friction Coefficient
The rolling friction coefficient has a significant impact on the simulation result. Although the ability of particles to rotate is switched off in some simulations, it is essential for the rotation to be enabled for a correct simulation result of pasty materials. Otherwise, a particle accumulation could only move translationally in space. An important behaviour of the pasty material, observed in the first step of the slump test, is the spreading of the flank during the outflowing process from the container. By choosing a rolling friction coefficient that is too low, particle layers lose the ability to flow over each other since individual particles lack the grip to climb out of pits from the particle layer below them. As a result, the upper particle layers cannot spread at a higher speed than the particles in the lower layers or the layer in contact with the base plate and therefore preventing the outflowing process described in more detail earlier. The outspreading particle front, therefore, consists of the same particles throughout the whole first step of the simulation, leading to the formation of cracks with increasing diameter of the front, since the gaps between the particles cannot be filled by inflowing particles from upper layers.

Restitution Coefficient
The particles restitution coefficient is set to 0 in the simulation since pasty materials do not rebound after colliding with a wall. In the second step of the simulation, a restitution coefficient not equal to 0 would result in the slumped material first being compressed after the collision with the base plate, due to their inertia and then rebounding into the air due to the restitution of the particles, which does not correspond to the real material behaviour.

Rupture Distance
The rupture distance not only significantly influences the calibration simulation but also all subsequent simulations, as discussed earlier, and must, therefore, be chosen wisely. An increase in the rupture distance leads to particles attracting one another over larger distances. This results in thicker flanks of the slumped material in calibration simulations with an increased rupture distance. In order to obtain the same results in the calibration simulations with an increased rupture distance, the maximum capillary force must be increased and the cohesion factor reduced.

Cohesion Factor
The cohesion factor of the constant cohesion model, taken into account at a positive overlap of the particles, has a major influence on the second part of the slump test. The greater the cohesion factor, the less the material is repelled back after the impact, which corresponds to the actual material behaviour. However, it should be noted that an increase of the cohesion factor also has a significant impact on the first part of the test, as this also results in an increase in the remaining height of the slumped material and a decrease of the diameter. Accordingly, an optimal compromise between cohesion factor, maximum capillary force and static friction coefficient must be found.

Conclusion
Pasty materials are used in a wide range of technical processes, which can often best be optimised through simulations. To be able to simulate pasty material behaviour with a DEM program, a suitable contact model has to be implemented in the software and methods to calibrate the simulation parameters are required. Several contact models following simple mathematical functions based on the capillary bridge effect were implemented into a DEM software, and the model best suited for the simulation of pasty materials was further developed. To evaluate the contact models, a standardised material test, adapted to reproduce a characteristic behaviour of a pasty material, was simulated. To ensure accurate results of industrial-sized process simulations in combination with the biggest possible time step, the discretisation effect was taken into account. In the end, the influence of the parameters from the implemented contact model and the standard Hertz-Mindlin model on the simulation results were documented and analysed to enable a simple and straight-forward process for future applications.

Future Work
By mixing pasty materials with granular media, e.g. when mixing the tailings after the mineral extraction process with the waste rock from the mining process to achieve better long-term storage options [1] for these waste products, the properties of the pasty medium change during the mixing process. Such aspects can further be taken into account by including/implementing additional variables into the contact model.