Development of a new DEM contact model for hygroscopic bulk solids

n dieser Arbeit wird ein Kontaktmodell zur Abbildung von kristallinen Festkörperbrücken in DEMSimulation vorgestellt und bisherigen sogenannten Bond-Modellen gegenübergestellt. Mit dem Modell sollen sämtliche Zustände eines hygroskopischen Schüttgutes in einem Silo mit einer allgemein gültigen Kalibrierung abgebildet werden. Anhand von Parametervariationen wird das entwickelte Simulationsmodell qualitativ verifiziert. Außerdem wird für die Ermittlung der notwendigen Modellparameter ein Ablauf für eine einfache Kalibrierung vorgeschlagen.


INTRODUCTION
Even today, hygroscopic bulk solids such as fertiliser urea and other salts are mainly stored in storage sheds or domes, since their behaviour in silos remains challenging. Disadvantages often associated with dome-based storage include time-consuming dosing and mixing using the 'last-in-first-out' principle, sometimes carried out with wheel loaders. These issues lead to reduced product quality and significant economic losses every year. Computer-assisted analysis of bulk material behaviour by means of the discrete element method (DEM) could remedy this situation and provide the necessary data for designing process-safe silos. However, industry still lacks a suitable contact model for the simulation of hygroscopic bulk solids that can reproduce the unconfined yield strength achieved by solid-bridge formation under various operating conditions.
In [3] it is shown that a DEM contact model for hygroscopic bulk solids must reproduce the following properties in order to cover the entire parameter space of possible operating conditions for a given set of material parameters: • strength as a function of moisture • strength as a function of time • strength as a function of pressure The present article concretises and verifies the approach presented in [3]. Figure 9 illustrates this approach as a simulation scheme for the calculation of crystallisation-induced solid-bridges. In addition, we present a simple calibration procedure for determining the required model parameters.

EXISTING DEM COHESION MODELS FOR HYGROSCOPIC BULK SOLIDS
For the development of contact models for the DEM simulation, cohesive bulk solids can be divided into two groups based on the binding mechanism: reversible and irreversible consolidation.
Reversible and an irreversible consolidation can be distinguished based on the time it takes for the formation of tensile forces acting between two neighbouring particles with respect to the duration of the simulation. In both groups, flow of the bulk solid can be caused only by exceeding the adhesion forces. If a bulk solid with reversible consolidation starts to flow, the formation of the corresponding adhesive forces takes place continuously during the simulation between neighbouring particles. This group includes cohesive powders with particles in the micro-to nanometer range, as well as moist bulk solids with particles up to the millimetre range. Corresponding to the large number of differerent binding mechanisms, a large number of contact models already exist for this group. These models are comparatively easy to implement, since in most cases the size of the adhesive forces is only a function of interparticle distance.
However, in the case of bulk solids with an irreversible binding mechanism, no adhesive force is formed between other particles when the adhesive forces are exceeded. The models developed for this purpose are so-called bond models and are used, for example, for chemical compound bulk solids with sinter bridges and hygroscopic bulk solids with crystalline solid-bridges, such as urea. The implementation in DEM programmes is a bit more complicated: In addition to the size of a possible adhesive force, the possible persistence of an upright connection between neighbouring particles must be determined. It should be noted that the term 'irreversible' applies only to typical simulation times of DEM simulations. Over typical realworld storage timescales, which go far beyond possible simulation timescales, consolidation can certainly occur via renewed formation of adhesive forces after previous breakage. A bond model implementation commonly used in DEM programmes is that for rocks after Potyondy and Cundall [7]. The input parameters required for this model are listed in Table 1.
Although bond models seem to be suitable for hygroscopic bulk solids in general, the implementation according to [7] is hardly applicable to industrial applications for designing a silo for urea, for example. In contrast to the model's workings in its original field of application, the dependence of solid-bridge strength on storage time and moisture content is strong for hygroscopic bulk solids. Therefore, the already complex calibration of the models, in particular the bond radius, must be carried out each time under new operating conditions. Furthermore, the bond radius and thus the strength among all particles is assumed to be the same size, with local particle distances not considered. In the following analysis, therefore, a bond model is presented that takes into account these dependencies of bulk solid pressure, moisture and storage time during the calibration process.

DEVELOPMENT OF A BOND MODEL FOR HYGROSCOPIC BULK SOLIDS
The principle of the new bond model has already been described in [3] (see Figure 9) and is based on the physical processes responsible for the formation of a crystalline solid-bridge. These processes occur in the following rough order: 1. Moisture exchange between the particles and the environment across particle surfaces 2. Capillary formation of liquid-bridges at contact points 3. Dissolution of the crystalline particulate solid 4. Formation of crystalline solid-bridges through complete evaporation or local supersaturation and crystal growth The time scale for processes 3) and 4) is many times greater than that for processes 1) and 2). For detailed information, see [3] or [4].
The bond model takes effect only at a user-selected time t create . At this time there is a random bed with various particle distances a. First, temporary liquidbridge radii are calculated on the basis of the predetermined moisture content described by the liquid saturation , from which bond radii are subsequently determined in accordance with the storage time t storage . The required parameters are summarised in Table 2.

Calculation of the liquid-bridge radii
Depending on the amount of liquid contained in a bulk material, different saturation regimes and different types of liquid-bridges between the particles can be distinguished. The amount of liquid is expressed in terms of the pore volume occupied by the liquid saturation S. flussig flussig In (1), ε describes the porosity of the bulk solid. In the so-called pendular-state, i.e., the range of S min < S < 0.08, liquid-bridges occur exclusively in the form of individual annular bridges (see Figure 1). The minimum moisture at which the bulk solid loses its cohesive properties is stated in [5] as S min = 0.002. In the funicular-state 0.08 < S < 0.24, the bulk solid gradually achieves union of individual bridges into more complex structures. Although the global connection of liquidbridges in both cases is interrupted by pores, liquid can be exchanged along the particle surfaces via channels in the size range of the roughness itself [5]. In the equilibrium state, therefore, capillary pressures are similar to one another across all liquid-bridges due to liquid exchange [10], if the influences of surface tension and gravitational acceleration are neglected. According to Pietsch and Rumpf [8], the radii, the volumes and the capillary pressures of the liquidbridges can be calculated given known particle diameter , particle distance , filling angle and contact angle (see Figure 2). With complete wetting, the contact angle , so the liquid-bridges connect tangentially to the particles.
In order to determine the equilibrium local liquidbridge radii, a target value for the capillary pressure as a function of liquid saturation must be obtained. To this end, a mean filling angle is first calculated on the basis of the empirical equation presented in [14]: and , which can be omitted for a = 0 and δ = 0.
Eq. (2) is in good agreement with Pietsch and Rumpf's [8] general model, in which cannot be represented in explicit form for a monodisperse bulk solid in the range 0.05 ≤ S ≤ 0.3. Although nonpenetration of neighbouring particles (a > 0) is considered, several DEM simulations have shown that the average particle distance a remains less than zero. Therefore the correction factor C a need not be taken into account for the calculation of β .
According to the Young-Laplace equation: with liquid surface tension γ i and the required radii according to [8]: the dimensionless positive capillary pressure F pK can be calculated if R 2 takes a negative sign in (5): The liquid surface tension γ i therefore does not need to be determined.
The local filling angles β and liquid-bridge radii R 2 are calculated as a function of the particle distance for a given dimensionless capillary pressure F pK .
Since an explicit equation for β according to Pietsch and Rumpf [8] is not possible, (2) would seem to be an obvious source for the calculation of the local filling angles. However, Figure 3 shows that there is not sufficient agreement between the analytic equation and the empirical equation using the correction factor C a (3).
As can be seen in Figure 4, on the other hand, the analytic function can be very well approximated by the following cubic function: 3 Eq. (9) is written by transformation as a reduced cubic equation without a quadratic term and then solved for β [°] by applying Cardano's method.

. Curve family f(F pK , , a/d P ). Comparison between [8] and by the cubic function (9) in the validity range of the developed contact model
According to Figure 4, equations (9) -(12) provide a very good approximation for (5) or rather (8) in a range of 0.1 ≤ F pK ≤ 18. The range specified for the dimensionless capillary pressure F pK covers the entire range of validity of the liquid saturation of 0.05 ≤ S ≤ 0.3 for a realistic bulk solid porosity of ε ≈ 0.3-0.5. .
The liquid saturation range is limited by the underlying model of individual annular bridge occurrences. Since the approximation currently contains no dependence on the contact angle δ, this work applies exclusively to δ = 0.
With the local filling angles determined, the liquidbridge radii R 2 are calculated according to (7).

Calculation of bond radii
On the basis of the determined liquid-bridge radii and the storage time to be investigated, the bond radii R B can be calculated using the time parameter t 63 according to Tomas [11]: According to (13), crystallisation growth is determined by an inverse exponential function. Therein t 63 describes the time at which the crystallisation progress reaches 63%. For very long storage timescales, the bond radius corresponds to the liquid-bridge radius. Depending on the solubility L of the particulate solid in the liquid, a more or less porous bond cross-section is formed during crystallisation. This circumstance is taken into account by the calibration of the breaking stress in accordance with Chapter 4. Since the solubility L can be strongly temperature-dependent (see Table 3), a direct consideration of this parameter in (13) is possible: When using (14) instead of (13), the calibrated breaking stress takes on a different value, and temperature dependence is taken into account by the actual solubility.

CALIBRATION OF MODEL PARAMETERS
The particle parameters and the interactions parameters without cohesion in Table 1 can be calibrated with experimental devices such as the Draw Dawn Test [9] or the Shear-Slip-Fall-Test [2]. To carry out a DEM simulation with the presented contact model, the parameters listed in Table 2 must be set. As with the parameters in Table 1, the liquid saturation and the storage time t storage must now be specified instead of directly defining the unknown bond radius. The liquid saturation of the bulk solid can be determined by means of so-called sorption isotherms from the relative humidity RH and the ambient temperature T U . Table 4 shows some values of the sorption isotherm of urea at 25 °C. which can be easily converted into the liquid saturation required as input parameter: The typical progression of sorption isotherms according to Table 4 shows at first a very small amount of adsorbed liquid and a strong increase in moisture from a critical relative humidity CRH.
The bulk porosity ε can be obtained in many cases from the results of a first simulation without bonds. However, if the particle size selected in the simulation does not correspond to the actual grain size, this is not possible and the parameter ε must then be calibrated. This case will not be discussed here. However, a possible procedure for its calibration could be a final comparison of the calculated liquid-bridge volume with the adjusted liquid saturation according to (1).
The physical value of the contact angle must be determined in accordance with the wettability of the bulk solid with the surrounding liquid. Therefore, a calibration is not required.
For the bond range a max , a relatively large value can be selected. Since no unphysical results are allowed while determining the correct solution of the cubic equation, no liquid-bridge will be created unless the capillary pressure at that distance reaches the desired target value, even if the value for a max is large. In this sense, a max only defines the range of possible contact partners. For many cases, a max can be chosen in a range of 0.02 ≤ a max /d p ≤ 0.05.
The time for the creation of the bonds t create is chosen by the user on the basis of the simulation process so that the particle bed has already come to rest and the kinetic energy is close to zero.
Since crystalline solid-bridges consist of the same material as the particles themselves, the same stiffness settings can be used for t he bond stiffnesses k N and k T as for the cohesionless particle contact model. If the elastic modulus is used to quantify particle stiffness, it can be converted into the required stiffness per unit area [N/m 3 ] as follows: Finally, the breaking stresses σ max and τ max and the crystallisation parameter t 63 remain from Table 2 for the calibration. Since studies on single grains (e.g., [12] and [4]) are generally very complex, the uniaxial compression test is used here and it must be carried out twice. For this purpose, the sample material is prepared identically in both cases with respect to liquid saturation and temperature. Subsequently, the samples are precompressed by pressure from above and stored in a climate chamber under defined conditions. In the last step, the walls of the mostly cylindrical sample geometry are removed, and the sample is crushed by loading from above. Based on the maximum applied force and the cross-section of the sample, the compressive strength is determined.
The first sample is stored as long as the crystallisation process is complete, at which point no further increase in compressive strength is to be expected with further storage ( Figure 5). The stress value determined at t 1 serves to calibrate the maximum breaking stresses σ max and τ max .
The heretofore unknown parameter t 63 required for the calibration simulations must be determined first. This is done using the second sample shown in Figure 5, which is stored for only a few hours. The exact storage time t 2 depends on the bulk solid and the added liquid. The goal is the determination of a representative point σ 2 on the solidification curve for the determination of t 63 with:

VERIFICATION AND VALIDATION
To verify the developed contact model, we simulate a uniaxial pressure test, and we then qualitatively assess the compressive strengths achieved at different precompactions, liquid saturations and storage times. For this purpose, particles with a uniform diameter of 6 mm are filled into a virtual sample cylinder. The cylinder has a diameter of 61 mm and a height of 85 mm. Figures 6, 7  and 8 show the results of these simulations.
According to Figure 6, a dependence on compression pressure is evident only at very low pressures. This restriction corresponds to the range over which particles rearrange from a loose bed to a dense packing. Further reduction of interparticle distance is relatively minor in comparison. Under typical silo pressures, this dependence is observed to be rather low, coinciding with the observations of Wahl et al. [12].    Figure 9. Simulation procedure for consideration of pressure-and time-dependent solid-state bridges [3] The trend in compressive strength as a function of liquid saturation and storage time (Figures 7 and 8) is in line with expectations. A detailed validation by comparative measurements is still pending.

CONCLUSION
A DEM contact model for hygroscopic bulk solids was developed. An essential feature is its dependence on storage time and moisture. With a one-time calibration, a parameter set valid for all operating conditions can be found. With the help of the presented calibration method, all required parameters can be determined with a simple test bench and little effort. Verification of the model was based on various parameter variations, yielding promising behaviour.
In the next step, the model for various hygroscopic bulk solids, such as urea, sugar, etc., is calibrated. Subsequently, the simulation results are compared with comprehensive measurement series.
Future steps also include an extension of the model with regard to the following possibilities: • Consideration of a possible particle size distribution • Detection of contact angles δ ≠ 0° • Impact of local humidity and temperature differences Приказан је модел елемената у контакту за симулацију кристалног моста чврстих материјала кориш-