Influence of dry cohesion on the micro-and macro-mechanical properties of dense polydisperse powders & grains

Understanding how cohesive granular materials behave is of interest for many industrial applications, such as pharmaceutical or food and civil engineering. Models of the behaviour of granular materials on the microscopic scale can be used to obtain macroscopic continuum relations by a micro-macro transition approach. The Discrete Element Method (DEM) is used to inspect the influence of cohesion on the micro and macro behaviour of granular assemblies by using an elasto-plastic cohesive contact model. Interestingly, we observe that frictional samples prepared with different cohesion values show a significant difference in pressure and coordination number in the jammed regime; the differences become more pronounced when packings are closer to the jamming density, i.e. the lowest density where the system is mechanically stable. Furthermore, we observe that cohesion has an influence on the jamming density for frictional samples, but there is no influence on the jamming density for frictionless samples.


Introduction
Granular materials behave differently from usual solids and fluids.One way to understand the macroscopic particle behaviour is using the Discrete Element Method (DEM) [1][2][3].Models based on this method rely on the contact force models used and solve a coupled system of equations of motion for all the interacting particles [4,5].Even though millions of particles can be simulated, this number (for fine powders) is generally too small to regard the systems modelled as truly macroscopic.Therefore, a transition from the micro to the macro level should be established.The microscopic properties can be used to derive macroscopic constitutive relations.These relations can be used to describe the particle behaviour on the large scale application/process level [6].The testing and characterization of dry, non-sticky powders are well established.The main challenge comes when the powders are cohesive and less flowable like those relevant in food industry [7].Research has already been done on cohesive granular materials (see refs [8,9]), however the influence of cohesion on granular packings is still poorly understood.There are two cases where cohesion becomes important: When particles become very small the cohesive forces become larger than the other forces on each particle, as is the case for dry fine powders [5,10].Not only the size of the particles contributes to the influence of cohesive attractive forces, but also a liquid between the particles, as is the case for wet granular materials [11][12][13].The research presented here will focus on dry cohesion and DEM is used to study granular packings made of polydisperse e-mail: r.g.h.kievitsbosch@student.utwente.nlcohesive particles.The question arises how does the presence of attractive forces affect macroscopic properties of the packings?So far, only a few attempts have been made to answer this question.Gilabert et al. [14] focussed on a two-dimensional packing made of particles with shortrange interactions (cohesive powders) under weak compaction.Yang et al. [15] studied the effect of cohesion on force structures in a static granular packing by changing particle size.Singh et al. [16] studied the effect of friction and cohesion on anisotropy in granular materials under quasi-static shear.The goal is to understand the influence of the microscopic parameters on the macroscopic properties of the packings.Knowing the influence of cohesion on particulate systems will advance of development of new constitutive models to predict the macroscopic material behaviour, to be used to model real life applications and to understand and optimize processes.

Simulation details
We use the Discrete Element Method (DEM) to understand the behaviour of granular systems.In the model, we relate the force interacting between the particles to the overlap δ that the particles have with each other.DEM solves Newton's equations of motion for all forces f i = f n n + f t t acting on particle i for the translational and rotational degrees of freedom.

Contact model
For the sake of simplicity, the linear visco-elastic normal contact force model can be used, but for cohesion the model in Sect.2.2 is applied.It involves a linear repulsive and a linear dissipative force: f n = kδ + γ 0 δ with k as spring stiffness, δ = (a i + a j ) − (r i − r j ) • n > 0 as particle overlap, n = n i j = (r i − r j )/|r i − r j | as normal unit vector, γ 0 as viscous damping coefficient and δ the relative velocity in normal direction v n = −v i j • n = δ.For simulations regarding frictional particles, we used the Coulomb friction law as explained in ref [5], with tangential force f t and tangential direction vector(s) t.
An artificial damping force f b is introduced to reduce dynamic effects and shorten relaxation times: This force acts not on contacts but directly on particles, proportional to their velocity v i .
The integration time-step Δt MD used for simulations needs to be much smaller than the contact duration t c to make sure that the integration of the equations of motion is stable.Note that in extreme cases of an overdamped spring, t c can become extremely, artificially large, i.e. dissipation γ should be neither too weak nor too strong.

Adhesive, elasto-plastic contact model
In this paper, a linear variant of the hysteric spring model is used (see refs [4,5,17,18]).This model is a simpler version of nonlinear-hysteric force laws [17][18][19][20].Different spring constants are assumed for loading and unloading.The (hysteric) force can be written as: with k 1 ≤ k 2 .k 2 is variable, depending on the previous deformation amplitude, see refs.[5,11]: There are different phases during the particle contact (see figure 1).First there is an initial loading; The force increases linearly with the overlap δ.The unloading starts after the overlap reached its maximum value δ max .This point defines the maximum force between two particles.During unloading, the force decreases linearly with a slope k 2 , to zero, at overlap δ 0 = (1 − k 1 /k 2 )δ max , which resembles the plastic contact deformation.Reloading will result in an increasing force with the same slope k 2 .When unloading takes place below δ 0 , then attractive forces will occur until the minimum cohesive force branch −k c δ min is reached.An overview of the parameters used in the DEM simulations can be seen in table 1.

Preparing samples
The sample preparation procedure consists of several steps.The first step of the preparation process is to compress isotropically the loose packing that was created randomly (with small overlaps) with a volume fraction of ν = 0.3 up to a volume fraction of ν = 0.5.The system  is then relaxed at a constant volume fraction of ν = 0.5 and the particles are allowed to dissipate their kinetic energy and reach a zero-pressure, relaxed configuration before jamming.Isotropic compression is then applied to reach a volume fraction of ν = 0.82 and, in the last phase, the compressed packing is decompressed isotropically until a volume fraction of ν = 0.5 is reached again [21].The compression and decompression are performed by applying a constant strain rate to each particle.We use a simulation time of 4000 [μs] per phase for the whole procedure, which results in a strain rate of ε = 3.8 • 10 −5 , for loading and unloading.By following the above procedure frictional and frictionless samples were studied for varying cohesion, from 0 ≤ k c /k ≤ 20.It is worth to mention that some particles have overlaps during the initial random generation.These overlaps form small clusters due to high cohesive force between particles which will last during sample evolution when k c /k > 1.This causes inhomogeneity for samples with extremely high cohesion.

Loading
Here, we present the general definition of the investigated quantities.Pressure and coordination number were calculated during sample compression to observe the influence of cohesion.We use the dimensionless pressure which is Powders & Grains 2017 calculated by using the average normal stress: P * = P/k * , where P = (σ xx + σ yy + σ zz )/3 and k * = k/(2 a ).The coordination number is defined as the average number of contacts per particle (C = M/N, where M is the total number of contacts and N is the total number particles).Particles with zero number of contacts and particles having a too small number of contacts, so called rattlers, were excluded, because they do not contribute to the mechanical stability of the packing [22].So the coordination number becomes: C * 4 = M 4 /N 4 , with M 4 = total number of contacts of particles with at least 4 contacts and N 4 = number of particles with at least 4 contacts.Note that C * 4 and C * 3 show similar behaviour for the results presented later, for that reason C * 4 is just shown here.The jamming density is the next quantity, it can be obtained from the performed simulations using the pressure and coordination number data.It is defined as the transition of the material from a fluid-like state to a solid-like state.When the constituent particles of such packings are so crowded that they are in mechanically stable contacts with one another, the whole packing undergoes a sudden dynamical arrest, where the static pressure becomes non-zero.show the fractions of particles with a certain number of contacts during compression and decompression.Looking at figure 2(a), there is only a small number of particles having more than three contacts at low volume fraction.It is clear that by increasing the volume fraction, the number of particles with a high number of contacts increases (C ≥ 4) and particles with a small number of contacts decreases (C < 4).Having a closer look, we can see that there is a regime of volume fractions, around the so-called jamming density, where the strongest changes occur in the number of contacts.This transition regime (around ν = 0.56) displays in particular a change in slope for C = 4.

Unloading
Figure 2(b) shows the particle fractions during the sample decompression.As expected, the number of particles possessing a high number of contacts is greater at high volume fraction.By moving along the decompression path, it can be seen that there is not a big change in the number of particles with high number of contacts (especially for C = 4, 5, 6) and very few particles with contacts C = 0 and 1; this proves the forming of agglomerates during decompression due to the cohesive forces.Particles sticking together will stay attached to each other during decompression due to the high cohesive forces among them.Note that the compression path is used for further results, because decompression path is inhomogeneous, but this effect and the jamming transition of cohesive packings will be studied elsewhere.In the case where the samples are prepared without friction, we see no significant changes in the pressure, throughout the whole regime (figure 3(a)) and in the coordination number at high volume fraction (figure 4(a)), for different cohesion values.For the lower densities, below jamming, a systematic increase in C * 4 can be seen with respect to k c .On the contrary, we observe a change in  ).This is less pronounced for highly dense samples.Cohesive forces are relatively stronger at low volume fractions, so that the change of coordination number can be observed for both frictionless and frictional samples.Since, there is an extra force (tangential) in frictional samples, it causes particles to slide and rotate, causing rearrangements and a smaller C * 4 .The effect of cohesive forces is thus already important at loose packings (close to the jamming density) since particles have more space to attract, repel and restructure.Far from the jamming density, there is not enough space for particles to establish new cohesive forces since they are strongly compressed and their normal force is dominated by elastic loading.

Conclusion
In this study the influence of cohesion on the macroand microscopic properties, pressure, coordination number and jamming density, was observed.Frictionless and frictional samples were prepared with different interparticle cohesion.Sample preparation plays a key role to obtain a homogeneous medium; for more cohesive packings we observed stronger formation of clusters and agglomerates during sample preparation in particular during unloading.During compression, the influence of cohesion is more pronounced for samples with friction due to tangential forces that cause rotations and rearrangements.Moreover, cohesion plays a more significant role at low volume fractions, since samples have not been compacted yet and translational and rotational movement of particles can activate cohesive forces.The jamming density is a crucial state variable at which the transition from fluid to solid-like behaviour of configurations occurs.We observed that variation of cohesion does not change the jamming point of frictionless samples, but frictional samples get affected by cohesion.Increasing cohesion leads to a jamming transition at considerably lower densities.Future studies will focus on the effect of cohesive forces on the jamming density.Some packings with strong cohesive forces formed some tiny clusters from the beginning due to cohesion.We are currently quantifying how such clusters in the packings affect the results by creating cohesive packings without the early formation of clusters.

Figure 1 .
Figure 1.Schematic graph of the piece-wise linear, hysteric model.The adhesive force-displacement for normal collision.The non-contact forces ( f 0 ) are kept equal to zero in this study and also the line for negative δ is neglected in this paper

7 Figure 2 .
Figure 2. Fraction of particles possessing number of contacts from 0 to 6 or more, for the frictional sample with k c /k = 20, during the sample (a) compression and (b) decompression.Where the arrow indicates the jamming density around ν = 0.56

Figures 3 and 4
Figures 3 and 4 depict the evolution of pressure and coordination number, respectively, for frictionless (figure 3(a) and 4(a)) and frictional (figure 3(b) and 4(b)) samples.All data is taken during isotropic compression with a few cases during decompression.

Figure 3 .
Figure 3. (a) Frictionless (μ = 0) and (b) frictional (μ = 0.5) simulation results of the dimensionless pressure P * plotted against volume fraction ν for different cohesion values k c , during compression (solid lines with points) and decompression (dashed lines), where the insert shows a zoom into the jamming regime.The arrows in the zoom indicate the range of φ j ∈ [0.56, 0.58] for μ = 0.5 and φ j ≈ 0.65 for μ = 0.The blue arrows in the plot indicate the compression and decompression paths for the largest and smallest k c

Figure 4 .
Figure 4. (a) Frictionless and (b) frictional coordination number C * plotted against volume fraction ν for different cohesion values k c during compression

Table 1 .
The microscopic contact model parameters values