Investigation into the isotropic compression of wet granular soils using discrete element method

. In this paper, the authors aim to present a 3D model using the Discrete Element Method (DEM) to study the quasistatic response of very loose assemblies of frictional spherical grains to an isotropic compression in the presence of a small amount of an interstitial liquid, which gives rise to capillary menisci and attractive forces. A reduced pressure P * is defined, comparing the confining pressure to the tensile strength of contacts. The effects of initial assembling process and various micromechanical parameters on the plastic compression curves are first studied. The influences of rolling resistance at contacts and the size polydispersity along those compression curves are also mentioned.


Introduction
Bonded granular materials are present in many natural and industrial processes, in which this bonding may have different physical origins, and proceed for instance from capillary effects, from solid precipitations, or from artificial solid bridges at the intergranular contacts. In geotechnical engineering, for the case of granular soils, the presence of water menisci between soil particles plays a key role in the overall behavior. Hence, the cohesive materials exhibit some specific features that do not appear in non-cohesive materials, such as the ability to form stable structures at very low density, and they are very sensitive to stress intensity as well as the external force. So far, there are a lot of studies in the macromicromechanical behavior of bonded granular materials, such as: cohesive powders [1,2], wet beads [3], loessic soils [4,5], wet sand [6,7], etc. The contact networks are hardly accessible to experiments, and intergranular forces are also inaccessible to measurements. Generally, the behavior of materials under the growing of external load (oedometric or isotropic tests) is characterized by the compression curve which is the irreversible compaction under increasing pressure.
The first numerical simulation method of the "discrete element" type was introduced 35 years ago by [8]. Since then, the DEM has become a valuable and efficient tool to investigate the microscopic mechanisms and to classify mechanical properties of granular systems. Yet, cohesive granular materials, especially in loose/very loose states, have less frequently been investigated by numerical simulation than cohesionless ones. Based on the 2D model for cohesive powders of Gilabert et al. [1,2], a 3D model [9] is developed to investigate the mechanical behavior of a model bonded granular soils (glass beads with capillary bonds in the pendular state). Results are presented for simulations of isotropic compression cycles, both for the macroscopic behavior (plastic compression curve) and the evolution of the arrangement of microstructure upon growing of pressure.

Force model
The elastic and frictional forces are jointly implemented in contacts as in [10]. The simplified Hertz-Mindlin-Deresiewicz force model is used for calculating the contacts between grains, relating the elastic normal force F H N to the normal deflection h, as described in [11,9]. The attractive capillary force is only present if grains have been in contact and the meniscus has not broken since. A meniscus of volume V m remains until separation distance reaches its rupture value, D r ≈ V m 1/3 [12]. The Maugis approximation [13] is chosen for capillary forces: D denotes here the distance between the surfaces of the particles joined by the liquid meniscus. F 0 = πaγcosθ is the maximum tensile force, involving surface tension γ of the water-air interface (7.27x10 -2 N/m at 20 o C), and θ is the wetting angle, with θ = 0 for a perfectly wetting liquid. Most calculations here are carried out with V m /(a 3 ) = 10 -3 , and thus the force range extends to D r = a/10 (corresponding to a very small capillary force). This ratio of meniscus volume V m to a 3 is one of the dimensionless control parameters in the present study, determining the number of liquid bonds and the range of the corresponding capillary attraction. This work is limited to the pendular state of isolated menisci [6]. The total static normal force combines the elastic normal force and the capillary force, F t N = F H N + F cap , as shown in Figure 1-a. A viscous force is added in contacts, opposing the normal relative velocity of grains as in [14], in order to damp the vibrations about equilibrium stateswith no notable influence on quasistatic rheology.
Two particles moving away from each other after a collision will only separate if their receding relative velocity is large enough to overcome the capillary force, i.e. larger than a threshold proportional to V * [1], with The influence of rolling resistance (RR) at contacts, as investigated in 2D in [2], is also studied. For simplicity, we only implement this feature with linear contact elasticity. The existence of RR is related to the particle surface roughness, such that contact regions are larger than the ones deduced from contact elasticity. Four additional parameters are necessary: (i) a rolling spring constan K R which expresses the proportion between relative rotation and rolling moment (in the tangential plane), as long as the rolling friction threshold is not reached; (ii) a pivoting spring constant K P (chosen equal to K R ), relating similarly the pivoting moment to the pivoting (i.e., along the normal direction) relative rotation angle; (iii) a rolling friction coefficient µ R with the dimension of a length, setting the maximum norm of the rolling moment ‖Γ R ‖ to µ R F H N , proportional to the elastic part of the normal force; and (iv) a pivoting friction coefficient µ P (chosen equal to µ R ), requesting similarly the absolute value of the pivoting moment to stay below µ P F H N . In most simulations with RR, the values of µ R /a = 0.02, K R /(K N a 2 ) =10 -4 while K T /K N = 1 and K N /a = 4 GPa.

Stress control and equilibrium
Numerical samples comprise N=4000 monodisperse spherical beads of diameter a and mass m, with no gravity, within a cubical cell with periodic boundary conditions. The cell edges are parallel to the coordinate axes (x α ) α=1,3 . Their lengths (L α ) α=1,3 vary simultaneously with the grain positions until a mechanical equilibrium state is achieved for all particles with externally imposed values (Σ α ) 1≤α≤3 of principal stresses: Here, Ω = L 1 L 2 L 3 is the sample volume, r (α) ij is coordinate α of vector r ij joining the centers of neighbor beads i and j and F (α) ij the corresponding force coordinate. Velocities v i of grain centers comprise, in addition to a periodic field, an affine term corresponding to the global strain rate. Equations of motion for dimensions L α are written in addition to the ordinary equations for the dynamics of a collection of solid objects, and they drive the system towards an equilibrium state in which conditions σ αα = Σ α are satisfied [14].
A typical intergranular force value F 1 = max(F 0 , Pa 2 ) sets the tolerance levels for individual grain equilibrium, where P is the applied pressure. A configuration is deemed equilibrated when the following conditions are simultaneously satisfied: (i) the net force (the total force) on each spherical grain is less than 10 -4 F 1 ; (ii) the total moment on each sphere is lower than 10 -4 F 1 a; (iii) the difference between imposed and measured stresses is less than 10 -4 F 1 a; and (iv) the kinetic energy per grain is less than 10 -7 F 1 a.

Specimen preparation
A hard sphere event-driven method is first used to prepare disordered, low density configurations of 4000 grains. All particles are then launched with Gaussiandistributed random velocities with quadratic mean V 0 . They collide and stick to one another within a cell of constant size, forming larger and larger aggregates. Finally, all grains are connected to one another by cohesive contacts and reach an equilibrium position. All specimens are prepared in the pendular state, just to create liquid bridges at the contact points of grains. Cohesive forces act through the liquid bridges.
This initial structure depends on one dimensionless parameter, characterizing the agitation intensity in the assembling stage. It is defined as the ratio of V 0 to the minimum receding velocity V * introduced in Eq. (2), this ratio is also called initial coordination number and it will be used priorly. As in [1], larger initial agitation levels are observed to produce better connected structures. The initial parameters are listed in Table 1. These values  0 =0.30, V 0 /V * =0.2041 and V m /a 3 =10 -3 will be set for the standard values.

Figure 2.
A cubical specimen before compaction process. Table 1. Initial parameters with initial solid fraction  0 , initial coordination numbers of contacts z 0 c , and of distant interactions In this table, the level of initial agitation in the end of the assembling stage, expressed by ratio V 0 /V * , determines the initial solid fraction.

Compaction process
The final packing structures are compacted under growing external isotropic pressure. The effect of the external pressure should be compared to the adhesive strength F 0 in contacts. This defines a dimensionless reduced pressure P * , P * = Pa 2 /F 0 , in which F 0 is already given in connection with Eq. (1), F 0 = πaγ (assuming a null contact angle), so P * = Pa/(πaγ). In the present study, P * = 1 means P = πγ/a ≈ 2 kPa, with beads of diameter 118µm in water. For P * << 1, initial loose structures are stabilized by adhesion and the relevant force scale is F 0 . For large reduced pressures P * >> 1, externally imposed forces dominate the adhesive strength, the relevant force scale is Pa 2 [9].
The main objective of the present paper is the study of the effect of a gradual compression, starting from cohesion-dominated loose states at small P * , and ending in confinement-dominated denser states at large P * . A stepwise pressure-controlled loading path is applied. In each compression step, pressure P * is multiplied by a constant factor 10 1/4 , and one waits until the new equilibrium configuration is reached, with the criteria stated in Sec. 2.2. The compression program is pursued until P max * (maximum value 10 4 ), well beyond complete plastic collapse is obtained. Then, the effect of P * back to 10 -3 from its highest value is also simulated.

Effect of initial solid fraction
All the tests are run with low packing solid fraction, from  0 =0.25 to  0 =0.45. In this study, the results are priorly shown in the conventional compressiondecompression curves in soil mechanics e-log(P * ), in which e = 1/ -1.  Three regimes can be classified in the compression curve, depending on the range of the pressure P * [9]. A first regime, thereafter called regime I, is observed in the range of low reduced pressures P * . In this regime, the initial structure can sustain the increasing pressure without rearrangement and void ratio e remains nearly constant. The regime II starts when void ratio starts to decrease sharply. Finally, void ratio e decreases slightly in regime III and approaches the minimum void ratio e min at the end of the loading process. Upon decompression, e increases slightly, remaining very close to e min . It can be observed from the decompression that the compaction cycle is an irreversible process. For different initial solid fraction  0 , the denser systems are able to support larger pressures before rearranging, whence larger regime I plateaus. However, in regime II, these compression curves start to converge and very nearly coincide in regime III.
The initial states of structures are strongly influenced by external pressures while the microstructures seem to have no effect on the microstructural arrangement. An important microstructural characteristic, the coordination number z, defined as the average number of interactions per grain, is the sum of the contact coordination number z c and the coordination number of distant interactions, through menisci joining non-contacting grains, z d . Figure 3-b plots z c , z d versus P * in the compression cycle. Initially, one has z d = 0, as the same minimal initial coordination number in the aggregation process, are low, and do not allow contact opening. z c and z d remain unchanged in regime I. Both coordination numbers start to increase as the structure collapses and reorganises in regime II. They exhibit little change in regime III, showing that the increase of the number of contacts is mainly due to the closing of narrow gaps between pairs of grains joined by a meniscus, as the structure is further compressed --a moderate effect, partly reversed upon unloading. Further unloading back to small pressures takes place with nearly constant coordination numbers of both types, reflecting the stability of the dense structure formed at high P * .

Effect of initial coordination number
Ratio V 0 /V * , characterizing the intensity of initial agitation and its ability to break adhesive contacts, strongly influences the initial assembling process and the resulting coordination number. Figure 4-a shows that this parameter mainly affects the beginning of the compression curve. Six different values of V 0 /V * are used, the difference between minimum and maximum values is two-hundredfold. Under low P * , the wider regime I plateaus correspond to the higher initial coordination numbers. If initial velocity V 0 is of the order of V * , the initial packing structure is strongly perturbed. In other words, the larger the initial coordination numbers, the stronger the initial structures.
Moreover, the relatively strong increases in initial coordination number of z c and z d (under low pressure P * ) prove that more contacts are created when large initial coordination numbers are considered (see Figure 4-b). Conversely, with low intensity of initial coordination numbers (V 0 <V * ), grains stick gently to one another and clusters of aggregated grains are not disturbed. In other words, lower initial coordination numbers induce more tenuous aggregates.

Effect of meniscus volume
The investigations of the influence of meniscus volume V m /a 3 are now reported, which is kept low enough for the material (with a saturation hardly exceeding 1%) is maintained in the pendular state. UNSAT Despite the importance of capillary bridges in the stabilization of loose structures, changing the meniscus volume (or liquid content) seems to have no effect on the macroscopic behavior of the granular specimens under isotropic compression (no apparent change in the compression curve, e versus log(P * )). This is corroborated by the fact that very little change in the coordination number of contacts z c is observed. Only the coordination number of distant interactions z d is notably influenced by such a change, especially in regimes II and III, and in the decompression process, as shown in Figure  5. As menisci break at larger interparticle distance with larger meniscus volume, more liquid bonds are present. Figure 6 illustrates the influence of initial coordination numbers V 0 /V * and rolling/pivoting friction coefficients µ R = µ P = 0.02 on the initial assembling process and coordination number for standard values of  0 and V m /a 3 . For the smallest initial coordination number, the appearance of RR at contacts between grains creates weaker system under low P * , with smaller coordination. Unlike the initial loose structures formed without RR, which can sustain a small non-vanishing pressure in regime I, the compression curve of the tenuously connected (with coordination number approaching 2, i.e., virtually no loop) initially formed with RR, and a low coordination number, responds by some plastic collapse at the first pressure increment, as shown in Figure 6-a. Its regime I, if it exists, is confined too much smaller pressures. It reappears on the pressure scale of Figure 6-a in systems assembled with larger initial coordination numbers, which are better coordinated (Figure 6-b), and collapse more abruptly in plastic compression. Interestingly, the systems assembled with the lower coordination numbers first collapse quickly by plastic compression (i.e., the structure is very fragile with RR), as z c is gradually increased meanwhile z d remains equal to zero (no contact opening creating a distant interaction). More complete analyses of the corresponding structural changes are being carried out, and their results will be published in forthcoming paper [15].

Effect of the size polydispersity
In granular materials, polydispersity means "several sources of dispersion", which is conventionally not only summarized by particle size distribution but also by particle shape and density. In this study, the effect of particle size distribution is only mentioned in order to compare with the monodisperse case. The specimen was also prepared in the same procedure and initial parameters, for 4000 frictional spherical grains with a uniform distribution of diameters (d max = 2d min ) by volume fractions. In the present study, a cubical specimen is created with d max = 2, d min = 1, as shown in Figure 7.   Figure 8. A comparison between monodisperse and polydisperse cases. Figure 8 illustrates a comparison of compression and decompression curves between the monodisperse and polydisperse grains. This result confirms that, with the same initial parameters, the effect of the particle size distribution is negligible for the monodisperse grain, at least in the macroscopic behavior. More detailed analyses in the microstructural changes will be mentioned in forthcoming paper [15].

Conclusions and Perspectives
These results show that the simple model of wet granular materials exhibit the familiar plastic compression properties of cohesive soils or powders, with a welldefined limit of gently assembled (minimally coordinated) structures, and low initial densities, such that the behavior does not depend on dynamical model ingredients. A small level of rolling and pivoting friction is enough to produce strong differences in the plastic behavior, with the collapse under growing confining pressure becoming notably more gradual. Moreover, the change of particle size seems to be no effect on the macroscopic behavior. The research program developed here can be extended in several directions. For example, the quasistatic behavior of a 3D model accounting for the effect of rolling resistance under oedometer and triaxial tests, in pendular state, is a straightforward extension of the present work. A comparison with an experimental investigation will be mentioned for macro-microscopic behavior. The role played by solid bridges instead of capillary bonds could also be investigated.