Discrete element simulation of wet granular materials: plastic compression

789 Discrete element simulation of wet granular materials: plastic compression IV International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2015 E. O˜nate, M. Bischoff, D.R.J. Owen, P. Wriggers & T. Zohdi (Eds) DISCRETE ELEMENT SIMULATION OF WET GRANULAR MATERIALS: PLASTIC COMPRESSION V.-D. THAN1,2, J.-N. ROUX1, A.-M. TANG1AND J.-M. PEREIRA1 1Universit´e Paris-Est, Laboratoire Navier (UMR 8205), CNRS, ENPC, IFSTTAR, F77455 Marne-la-Vall´ee Cedex 2, France email: [email protected]c.fr 2The University of Danang, College of Technology, Department of Civil Engineering, 48 Cao-Thang St., Hai-Chau Dist., Da Nang, Vietnam Key words: Plastic compression, Capillary force, Consolidation, DEM Abstract. We use Discrete Element Method (DEM) simulations in three dimensions (3D) 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. Previous results obtained in 2D [8] are generalized to systems that may be observed in the laboratory. We study the influence of the initial assembling process and of various micromechanical parameters on the plastic compression curves, from very loose states assembled at low P∗to maximally compressed ones in which capillary cohesion is negligible at large P∗. We also show how the plastic response along those compression curves is influenced by rolling resistance in contacts. 1 INTRODUCTION Many natural and industrial processes involve granular materials, in which grains are bonded, due to a variety of physical effects near intergranular contacts: capillary bridges, solid precipitation, artificial solid bridges. The water menisci joining solid particles in wet granular soils play a key role in the overall behavior studied in geotechnical engineering. Capillary cohesion bestows to these materials specific mechanical features that do not exist with dry grains, such as the ability to form stable structures at very low density, and a strong sensitivity to stress intensity as well as to stress direction. Published studies of bonded granular materials investigate the mechanical properties of cohesive soils (clays and silts) [18, 16], wet beads [6], cohesive powders [7, 8], loessic soils [4, 10], cemented sands [17, 11, 12, 2], or wet sands [15, 3]. The contact networks are hardly accessible to experiments, and intergranular forces are also inaccessible to measurements. Generally, the behavior of materials under a growing external load (oedometric or isotropic tests) is characterized by a compression curve showing an irreversible density increase with pressure. 1 790 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira “Discrete element”simulations, as introduced 35 years ago by Cundall & Strack [5], have become a valuable and efficient tool to investigate the microscopic mechanisms and to classify mechanical properties of granular systems. Cohesive granular materials, especially in loose states, have less frequently been investigated by numerical simulation than cohesionless ones. Inspired by the two-dimensional model for cohesive powders of Gilabert et al. [7, 8], a 3D model 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 behaviour (plastic compression curve) and for the evolution of microstructural and micromechanical parameters. 2 MODEL MATERIAL 2.1 Intergranular forces We consider assemblies of spherical beads of diameter aand mass m, with the elastic properties of glass: Young modulus E= 70 GPa, Poisson ratio ν=0.3. They interact in their contacts by the Hertz law, relating the elastic normal force FE Nto the normal deflection has FE N=E 3(1 −ν2)√ah3/2.(1) This corresponds to a force-dependent normal stiffness KN=dF H N dh =E 2(1−ν2)√ah1/2. The tangential force model, combining elasticity and friction, is described by [1]. Tangential stiffness KT(relating increments of tangential elastic force and of tangential elastic displacement is proportional to KN,KT=2−2ν 2−νKN. Alternatively, a linear contact elasticity might be implemented, with h-independent stiffness coefficients KN,KT. The Coulomb condition ||FT|| ≤ µFE N, involving the sole repulsive elastic part of the normal force (Fig. 1(b)), is enforced with friction coefficient µ=0.3. The attractive capillary force is only present if grains have been in contact and the meniscus has not broken since. A meniscus of volume Vmremains until separation distance reaches its rupture value, Dr≃V1/3 m[9]. We use the Maugis approximation [14] for capillary forces: Fcap N=−F01−1 1+ 4Vm πaD2.(2) Ddenotes here the distance between the surfaces of the particles joined by the liquid meniscus. F0=πaγ cos θis the maximum tensile force, involving surface tension γof the water-air interface (7.27 ×10−2N/m at 20oC), and θis the wetting angle, with θ=0 for a perfectly wetting liquid. Most calculations here are carried out with Vm/(a3)= 10−3, and thus the force range extends to Dr=a/10 (corresponding to a very small capillary force). This ratio of meniscus volume Vmto a3is one of the dimensionless 2 791 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira (a) (b) Figure 1: (a) Static normal force, FN=FE N+Fcap, versus deflection hor distance D=−h. (b) Coulomb cone limiting the value of the tangential force. control parameters in the present study, determining the number of liquid bonds and the range of the corresponding capillary attraction. Our work is limited to the pendular state of isolated menisci [15]. The total static normal force combines the elastic normal force (Hertzian force) and the capillary force, Ftot N=FE N+Fcap N, as shown in Fig. 1(a). A viscous force is added in contacts, opposing the normal relative velocity of grains as in [1], in order to damp the vibrations about equilibrium states – with 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∗[7], with V∗=F0Dr m.(3) The influence of rolling resistance (RR) at contacts, as investigated in 2D in [8], 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 constant KRwhich 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 KP(chosen equal to KR), relating similarly the pivoting moment to the pivoting (i.e., along the normal direction) relative rotation angle; (iii) a rolling friction coefficient µRwith the dimension of a length, setting the maximum norm of the rolling moment ||ΓR|| to µRFE 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 µPFE N. In most simulations with RR, we set µR/a=0.02, KR/(KNa2)=2.5×10−5while KT/KN=1 and KN/a = 4 GPa. 3 792 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira 2.2 Stress control and equilibrium Numerical samples comprise N=4000 equal-sized spherical beads of diameter aand mass m, with no gravity, within a cuboidal cell with periodic boundary conditions. The cell edges are parallel to the coordinate axes (xα)α=1,3. Their lengths (Lα)α=1,3vary simultaneously with the grain positions until a mechanical equilibrium state is achieved for all particles with externally imposed values (Σα)1≤α≤3of principal stresses: σαα =1 Ω i mivα ivα i+ i<j F(α) ij r(α) ij .(4) Here, Ω = L1L2L3is the sample volume, r(α) ij is coordinate αof vector rij joining the centers of neighbor beads iand jand F(α) ij the corresponding force coordinate. Velocities viof 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 [1]. A typical intergranular force value F1=max(F0,Pa2) sets the tolerance levels for individual grain equilibrium, where Pis 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−4F1; (ii) the total moment on each sphere is lower than 10−4F1a; (iii) the difference between imposed and measured stresses is less than 10−4F1/a; and (iv) the kinetic energy per grain is less than 10−7F1a. 2.3 Dimensionless control parameters While calculations are carried out with glass beads of diameter a=0.11 mm, assuming perfect wetting (θ= 0), it is convenient to express dimensionless results as functions of dimensionless input parameters, thereby achieving greater generality. Aside from the friction coefficient, the important dimensionless combinations are the reduced pressure, P∗=a2P/F0, comparing the applied pressure to the tensile strength of contacts, and the stiffness parameter κ=[E/P(1 −ν2)]2/3(for Hertzian contacts), or κ=KN/P (for linear elasticity). P∗1 in cohesion-dominated systems, for which attractive forces may stabilize loose structures. Confining forces dominate for P∗1, and attractive forces become negligible. For the chosen value of a,P∗= 1 corresponds to P= 2 kPa. κ, on the other hand, sets the typical scale of deflections hunder confining forces, as h/a ∼κ−1[1]. The simulation parameters are chosen such that P∗reaches large values before κ−1decreases to 10−3. Additional parameters are µR/a =µP/a in systems with RR, as well as the rotation angles for which rolling and pivoting friction thresholds are reached (small enough to be irrelevant in our case). 4 793 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira Table 1: Initial configurations Φ0V0/V ∗Vm/a3 0.30 0.2041 1.00 ×10−3(reference case) 0.30 0.4082 1.00 ×10−3 1.2247 4.0825 12.2474 40.8248 0.30 0.2041 5.00 ×10−4|2.50 ×10−4 1.25 ×10−4|6.25 ×10−5 3.13 ×10−5|1.56 ×10−5 7.80 ×10−6 0.32 0.2041 1.00 ×10−3 0.35 0.2041 1.00 ×10−3 0.40 0.2041 1.00 ×10−3 0.45 0.2041 1.00 ×10−3 3 NUMERICAL PROCEDURES 3.1 Specimen preparation We first use a hard sphere event-driven method to prepare disordered, low density configurations of 4000 grains. All particles are then launched with Gaussian-distributed random velocities with quadratic mean V0. 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. This initial structure depends on one dimensionless parameter, characterizing the agitation intensity in the assembling stage. It is defined as the ratio of V0to the minimum receding velocity V∗introduced in Eq. (3). As in Ref. [7], larger initial agitation levels are observed to produce better connected structures. The initial parameters are listed in Table 1. 3.2 Compaction process The final packing structures is compacted under growing external isotropic pressure. 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 confinementdominated 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 101/4≃ 1.7783, and one waits until the new equilibrium configuration is reached, within the tolerance criteria stated in Sec.2.2. The compression program is pursued until P∗ max, well beyond complete plastic collapse is obtained. Then, the effect of P∗back to minimum value from its highest value is also simulated. 5 794 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira 4 NUMERICAL RESULTS 4.1 Reference case A typical test is run with a low initial solid fraction, Φ0=0.30. The results are shown in the from of the conventional compression curve e−log(P) in Fig. 2. Three regimes (a) (b) 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Voidratio,e Reducedpressure,P* Referencecase y=-0.3576x+0.9727 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Voidratio,e Pressure,P(kPa) P*=1 Referencecase Typicalresultofcohesionlessgrains Figure 2: (a) Compression and decompression curves in reference case, with (b) comparison with cohesionless systems. can be distinguished in the compression curve, for different pressure ranges. A first regime, thereafter called regime I, is observed for low reduced pressures P∗, in which the initial structure still sustains the increasing pressure without rearranging, and void ratio eremains nearly constant. In regime II, roughly corresponding to interval 0.02 ≤P∗≤2, void ratio strongly decreases, and might be described as linearly varying with log P∗. The loose structures formed at low P∗are no longer able to support the increasing confining pressure, they collapse and restructure. Finally, in regime III,egradually approaches some minimum void ratio emin at the end of the loading process. The slight decrease of efor P∗≥100 is similar to the cohesionless result, and due to elastic deflections in a stable contact network. The chosen parameters are such that κremains large enough not to influence the compression process taking place in regime II. Upon decompression, e increases slightly, remaining very close to emin: the compaction is essentially irreversible, unlike in the cohesionless granular assembly of Fig. 2(b), for which loading and unloading branches are not distinguishable. The compression curve is similar to the ones obtained by Gilabert et al. [8]. 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 zcand the coordination number of distant interactions, through menisci joining non-contacting grains, zd. Fig.3 plots zc,zd, and zversus P∗in the compression cycle. Initially, one has zd= 0, as the velocities in the aggregation process, in 6 795 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira (a) (b) 4.00 4.50 5.00 5.50 6.00 6.50 7.00 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Totalcoordinationnumber,z Reducedpressure,P* 0.00 1.00 2.00 3.00 4.00 5.00 6.00 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Coordinationnumbersofzdandzc Reducedpressure,P* zd zc Figure 3: Coordination numbers in reference system subjected to the pressure cycle: z(a); and zcand zd(b), versus P∗. the reference case, are low, and do not allow contact opening. zcand zdremain unchanged in regime I. Both coordination numbers start to increase as the structure collapses and reorganizes in regime II.zexhibits 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∗. 4.2 Effect of initial solid fraction Among the features affected by the assembling process, the competition between compression and aggregation is the most important one. Fig. 4, obtained with standard values V0/V ∗=0.2041, Vm/a3= 10−3, compares the compression curves of specimens with different initial solid fractions Φ0(the red arrow denotes the increase of Φ0). Denser systems are able to support larger pressures before rearranging, whence larger regime Iplateaus (see Fig. 4(a)). However, in regime II, these compression curves start to converge and very nearly coincide in regime III. Coordination numbers zcand zdare and remain almost equal whatever the initial Φ0(see Fig. 4(b)). 4.3 Effect of initial agitation intensity Ratio V0/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. Fig. 5 shows that this parameter mainly affects the beginning of the compression curve. Six different values of V0/V ∗are used, the difference between 7 796 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira (a) (b) zd zc 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 Voidratio,e Reducedpressure,P* Φ =0.30 Φ=0.32 Φ=0.35 Φ=0.40 Φ=0.45 y1=-0.3427x+0.9926 y2=-0.1101x+0.9187 0.00 1.00 2.00 3.00 4.00 5.00 6.00 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Coordinationnumberszdandzc Reducedpressure,P* Φ=0.30 Φ=0.32 Φ=0.35 Φ=0.40 Φ=0.45 Figure 4: (a) Void ratio eand (b) coordination numbers zcand zdversus P∗in compression cycle for different values of Φ0. minimum and maximum values is two-hundredfold. Under low P∗, the wider regime I plateaus correspond to the higher initial agitation velocities. I If initial velocity V0is of the order of V∗, the initial packing structure is strongly perturbed. I n other words, the larger the velocity, the stronger the initial structure (see Figure 5(a)), with more contacts – relatively strong increases in initial coordination numbers zcand zdare observed at low P∗(see Figure 5(b)). Conversely, with low agitation velocities (V0≤V∗), grains stick gently to one another and clusters of aggregated grains are not disturbed. In other words, lower agitation velocities induce more tenuous aggregates. 4.4 Effect of meniscus volume We now report on investigations of the influence of meniscus volume Vm/a3, which is kept low enough for the material (with a saturation hardly exceeding 1%) is maintained in the pendular state. 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 behaviour of the granular specimens under isotropic compression (no apparent change in the compression curve, eversus log P∗). This is corroborated by the fact that very little change in the coordination number of contacts zcis observed. Only the coordination number of distant interactions is notably influenced by such a change, especially in regimes II and III, and in the decompression process, as shown in Fig. 6. As menisci break at larger interparticle distance with larger meniscus volume, more liquid bonds are present. 8 797 V.-D. Than, J.-N. Roux, A.-M. Tang and J.-M. Pereira (a) (b) zd zc 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Voidratio,e Reducedpressure,P* V 0 /V * = 0.2041 0.4082 1.2247 4.0825 12.2474 40.8248 y1=-0.3576x+0.9727 y2=-0.5575x+0.8876 0.00 1.00 2.00 3.00 4.00 5.00 6.00 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Coordinationnumberszdandzc Reducedpressure,P* V0/V*= 0.2041 0.4082 1.2247 4.0825 12.2474 40.8248 Figure 5: (a) Compression and decompression curves and (b) coordination numbers zcand zdfor the reference case with different values of V0/V ∗. zd zc 0.00 1.00 2.00 3.00 4.00 5.00 6.00 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Coordinationnumberszdandzc Reducedpressure,P* Vm/a3= 1.00x10-3 5.00x10-4 2.50x10-4 1.25x10-4 6.25x10-5 3.13x10-5 1.26x10-5 7.80x10-6 Figure 6: Coordination number zc,zd,z=zc+zdfor standard values Φ0=0.30, V0/V ∗=0.2041 and for different values of Vm/a3. 4.5 Effect of rolling resistance Fig. 7 shows the effect of initial velocities V0and rolling/pivoting friction coefficients µR/a =µP/a =0.02 on the initial assembling process and coordination numbers for Φ0=0.30 and Vm/a3= 10−3. For the smallest velocities, the appearance of RR in contacts creates denser systems under low P∗, with smaller coordination numbers. 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 9