Numerical modelling of powder caking at REV scale by using DEM

This work deals with numerical simulation of powder caking process caused by capillary condensation phenomenon. Caking consists in unwanted agglomeration of powder particles. This process is often irreversible and not easy to predict. To reproduce mechanism involved by caking phenomenon we have used the Discrete Elements Method (DEM). In the present work, we mainly focus on the role of capillary condensation and subsequent liquid bridge formation within a granular medium exposed to fluctuations of ambient relative humidity. Such bridges cause an attractive force between particles, leading to the formation of a cake with intrinsic physicochemical and mechanical properties. By considering a Representative Elementary Volume (REV), the DEM is then performed by means of a MULTICOR-3D software tacking into account the properties of the cake (degree of saturation) in order to establish relationships between the microscopic parameters and the macroscopic behaviour (tensile strength).


Introduction
Powder solids are regularly faced to the problem of caking during their periods of storage and/or transport.The agglomeration of particles that compose a powder depends not only on their water content, the temperature and the applied pressure, but also on the interactions between the solid substance and water molecules present in the atmosphere, i.e., on relative humidity (R h ) at which the product is stored [1].It is well known that ambient humidity plays an important role in most stages leading to caking: capillary condensation of water at contact points between particles, subsequent dissolution of a solid and formation of saturated solution eventually followed by precipitation of solid during evaporation of water.Almost all types of powder solids are subjected to caking mechanism, including crystalline or amorphous materials, soluble or insoluble, fine or coarse, etc. Caking is involved by attractive forces acting between particles under environmental conditions over time.As mentioned above, caking arises due to interparticle forces depending on powder state and physical aspects.For example, van der Waals forces in the case of dry caking, capillary forces for wet caking, material bridge formation for solid caking and more occasionally, electrostatic, magnetic or interlocking forces are also involved.According to many theoretical studies and experimental works carried out on this issue, capillary forces are at least 1 or 2 orders of magnitude higher than van der Waals forces.In conclusion whatever the main mechanism of caking is, the presence of water either in vapour, liquid, adsorbed or absorbed state could largely accentuate the process.In fact, capillary condensation is the first e-mail: mohamed.guessasma@u-picardie.frstep of caking for both amorphous and crystalline powders [2][3][4][5][6].This phenomenon considerably affects the amount of water which appears in a bulk media at a given R h .The water will then act as a plasticizer in case of amorphous materials or as a solvent for crystalline powders [7].In the present work, we are interested in simulating capillary force effects due to caking mechanism by using DEMbased approach.Indeed, in view of the nature of the powder, the DEM seems to be the best way for modelling a large number of particles taking into account physical and mechanical aspects.The state of art on DEM modelling of discrete media with capillary interactions reveals many works dealing with numerous applications fields such geomechanics and process engineering [8,9].With DEM, the second law of Newton is solved for each particle and the contact interactions between particles are computed by specific constitutive laws.This allows to characterise the behaviour of discrete media at a macroscopic scale starting from interactions at the binary contact scale.By considering a REV, the DEM is then performed using a numerical tool to study the influence of physical parameters on the mechanical behaviour, for example saturation degree and tensile strength.All numerical simulations discussed in this contribution have been carried out with MULTICOR-3D software [10].

Theoretical modelling
At relatively higher water activities, a w > 90%, water molecules become mobile enough to form liquid bridges, subsequently, the adsorbed molecules on the solid surface form a liquidlike a meniscus [2,3].Geometric as-Figure 1.Schematic representation of a pendular liquid [14] sumptions and analytic models to characterise the capillary bond are given below.

Fundamental equations
Usually the liquid interface is described by two principal radii, ρ 1 and ρ 2 , that define the capillary radius r c 1.
Kelvin model is commonly used to describe the capillary condensation phenomenon at high a w .Under mechanical equilibrium condition at the liquid interface, the Kelvin equation is given by 2.
with γ the vapour-liquid surface tension, R the universal gas constant, T the temperature and V l the molar volume of liquid.
Figure 1 gives a schematic representation of liquid meniscus in a pendular bond state when two identical particles (R 1 = R 2 = R) at a distant of 2h are joined by a liquid bridge.Considering the capillary torque profile [11,12], the first and second principal radii of curvature are given as follows 3.
where ψ and θ are the filling and average liquid-solid contact angles, respectively.
The volume of liquid V l is obtained by integrating the function expressing the arc > ABC (figure 1) with respect the variable x.The dimensionless form of liquid volume, V * = 3V l 4πR 3 , is commonly considered [14].As evoked in introduction, the fluctuations of ambient R h led to the formation of liquid bridge within a powder, which causes an attractive force between particles.According to Laplace theory for capillary forces, the total attractive force [13] between the two particles including both the capillary and the surface tension components is given by equation 4.

Parametric studies on a binary agglomerate
When a binary agglomerate is subjected to humidification the capillary condensation happens between two limiting values of the filling angle ψ.The variation of principal radii, ρ 1 and ρ 2 , as a function of ψ, in the case of a perfectly wetting system (θ = 0), leads to extremal values of ψ given by the intersection of curves.As shown on figure 2, the filling angle ψ, for 2h * = 0 (h * = h R ) and θ = 0, reaches extremal values, ψ min and ψ max at 0 • and 53 • , respectively.Beyond these values, the equivalent capillary radius r c becomes negative which means that the capillary condensation is brought to an end.Besides, if the relative gap 2h * > 0.3, capillary condensation phenomenon stops [14].
The capillary force acting between two particles is very sensitive to the water activity, a w , due to the liquid amount of the capillary bridge.The influence of the liquid volume is shown through the evolution of the capillary force as function of the separating gap for different volume of water (figure 3).Obviously, the attractive force increases when the relative gap between particles decreases.Furthermore, the attractive force varies slightly when the amount of liquid reaches its maximum (5% of particle volume) and exponentially when the amount of liquid becomes lower.The next study deals with the variation of capillary force with respect to the filling angle for all considered relative gap 2h * (figure 4).The volume of liquid is kept constant.As shown on the figures a lower value of filling angle involves a capillary force close to zero, except when the separating gap 2h * is nil.In such case, the dimensionless attractive force reaches a symbolic value, namely F * = π.

DEM modelling
In this section we are interested on modelling of capillary forces acting into a Representative Elementary Volume  (REV) composed of a set of particles subjected to a given water activity.The numerical modelling is made with Discrete Element Method (DEM) approach by using the explicit formulation [15], which is also called smoothed formulation.

Contact law
The spherical particles are rigid, a local contact deformation is considered and the interactions (F n , F t ) are governed by analogies with damped springs (K n,t , C n,t ) in normal and tangential directions n, t .Therefore, by using discrete approach, the contact forces between particles, without capillary effect, are described with a contact model depending on elastic force displacement law, Coulomb friction and viscous damping coefficient 5.
δ n,t , v n,t : relative displacements and velocities In the presence of water vapour condensation, the particles in the REV are subjected to attractive normal forces due to the formation of capillary bounds.The normal component of the contact force is, then, composed of two parts;  repulsive one and attractive one 6.The viscous effect of liquid is neglected in the modelling.
On figure 5-a, representing a snapshot from MULTICOR-3D software, are shown four binary contacts for different separating distances 2h * (0.02, 0.08, 0.2, 0.3).At each contact, the capillary bridge shape is drawn with respect their respective to geometric parameters, namely the principal radii ρ 1 , ρ 2 and the filling angle ψ.

Discrete simulations of capillary effects at a REV scale
A REV of 10 000 monosized particles with radius R = 2.285 μm and density ρ v = 1500 kg/m 3 is generated with MULTICOR-3D software by using a radius expansion algorithm [16] until the porosity of packing reaches a desired value.In the present case, the porosity is fixed to 50%.Saturation degree, υ, of the REV is of the order of 0.64% and considered as a reference value for simulations with multiple saturation degree (υ varies in the interval [0.64% − 12.83%]).The amount of liquid is spread equally at each contact and kept constant during the simulations (figure 5-b).It should be recalled that only the pairs of particles with a relative gap 2h * ≤ 0.3 are subjected to capillary forces.Thereafter, the attractive forces are computed numerically using an iterative algorithm.
The influence of saturation degree is investigated by considering tensile tests for several values of υ.The tensile strength of the REV is obtained numerically and compared to analytical one, according to Rumpf's model [13], for a random packing of monosized spherical particles and an even distribution of liquid.Regarding the tensile test, the REV is placed between two walls, at the lower and upper sides 6.Before starting the tensile test, the REV should first be in an equilibrium state.In fact, when the saturation degree is affected to the REV the particles start to move slightly as a result of the capillary forces.Once the REV reaches its equilibrium (minimum of kinetic energy) the tensile test is started.The lower side of the REV is fixed during the test while the upper one moves vertically with a constant velocity (10 −4 m.s −1 ).The particles of both lower and upper sides are subjected to the same kinematic conditions as the fixed and mobile walls.On figure 7 are plotted the tensile curves for six values of the saturation degree.The tensile stress, σ, is estimated from the resultant of capillary force acting on the cross section at the middle of the REV.The relative displacement, , is computed with respect to the particle radius.As expected, the higher the saturation degree of the REV, the higher the tensile strength is.This is due to both increasing of liquid amount at binary contact scale and coordination number of particles.Comparison made with analytical model exhibits discrepancies about 30% with numerical predictions.This is mainly due to the statistical parameters of the Rumpf's model, namely the coordination number and the mean capillary force.Nevertheless, the DEM remains more suitable to model mechanical behaviour of powder agglomerate starting from binary contact scale.

Figure 4 .
Figure 4. Dimensionless capillary force F * as function of filling angle ψ

Figure 7 .
Figure 7. Tensile curves for several saturation degrees