A mesoscopic model for compression of granular materials

The design of protective structures often requires numerical modeling of shock-wave propagation in the surrounding soils. Properties of the soil such as grain-grading and water-fraction may vary spatially around a structure and among different sites. To better understand how these properties affect wave propagation we study how the meso-structure of soils affects their equation of state (EOS). In this work we present a meso-mechanical model for granular materials based on a simple representation of the grains as solid spheres. Grain-grading is prescribed, and a packing algorithm is used to obtain periodic grain morphologies of tightly packed randomly distributed spheres. The model is calibrated by using experimental data of sand compaction and sound-speed measurements from the literature. We study the effects of graingrading and show that the pressures at low strains exhibit high sensitivity to the level of connectivity between grains. At high strains, the EOS of the bulk material of the grains dominates the behavior of the EOS of the granular material.


Introduction
The study of the mechanical behavior in granular materials has a large importance in various fields.Each field requires accurate modeling of different aspects of the complex mechanics.Issues in powder mechanics in silos and soil strength in civil engineering require the understanding of strength and stability of soils under low and intermediate pressures.The study of shock-wave propagation for the defense industry [1] requires models for the response of different soils under compression at a wide range of pressures.In the work presented here we focus on the volumetric response of granular materials to pressure, i.e., their EOS.
Experimental studies of the dynamic compression of soils are complex and limited.Dynamic measurements in the split Hopkinson bar test and measurements of spherical-shock propagation present one-dimensional cases of the strain state [2][3].Thus, neither retrieve the material behavior under different spatial loading scenarios, nor enable the separation of material response into a deviatoric and a hydrostatic part.Static measurements of the mechanical response of soils are regularly performed in tri-axial cell devices.These experiments enable different spatial loads but not at the full range of pressures required for shock dynamics.In addition, there is an open discussion in the literature whether the material response of soils is rate-dependent [4][5].For that reason, the applicability of static measurements for the modeling of dynamic phenomenon is controversial.Of special interest for our work is static tri-axial experiments done by Laine and Sandvik on sand from the shores of Sjöbo, Norway [6].In their experiments the EOS was isolated by applying a hydrostatic load, and strains at pressures up to 100 MPa were measured.The experimental results were later implemented in the commonly used hydrodynamic code AUTODYN [7].We use these results to calibrate our model.
In this work we focus on capturing the full range of compaction including large deformations of the grains up to full closure of the pores in the material and subsequent compression of the bulk material (the material from which the grains are formed).The aim is to find the dependence of the EOS on different material parameters such as: grain-size distribution, porosity, and properties of the bulk material.
Preliminary work to ours was done by Larcher and Gebbeken to compute the full compaction range by a mesoscopic model [8].Unfortunately, convergence of the model was not achieved, and no conclusions regarding the effect of the grain meso-structure or material properties on the EOS were derived.

The model
The mesoscopic model presented here predicts the EOS of granular materials by utilizing finite element simulations (FEM) of random grain-morphologies.In order to create the morphology of a granular material, we apply a random packing algorithm and create different realizations of random grain morphologies.Then, a material model is defined for the bulk material of the grains.The morphologies are loaded and unloaded by displacement boundary conditions, and the pressure is measured.This section provides the details of the implementation our model.

Morphology generation
A granular material is generally composed of grains of different shapes, sizes and composition.Under compression, the grains are packed and deform.For simplicity we degenerate the shape of the grains into solid spheres.To form the morphology, the spheres are packed by using a random sphere-packing algorithm [9][10].The initial condition for the packing algorithm is a periodic ensemble of non-overlapping spheres created by a Matérn hard-core process (see Fig. 1).Then, each sphere is given a random initial velocity vector, and a constant growth-rate.On each computational step of the algorithm, the collision times between all spheres are calculated and the spheres are repositioned to their locations at the nearest collision time.The post-collision velocities are then computed and the computational step is repeated until the ensemble locks.The resulting morphology is then cut to obtain a cubic representative volume element (RVE) to be used in the FEM simulation (see Fig. 1).
Different sphere-size distributions may yield different porosities.The porosity for a random packing of equal diameter spheres is ϕ 0 = 0.38, and interestingly corresponds to the porosity of SJOBO sand.In correspondence we chose an equal diameter spheremorphology in our simulations.

Boundary conditions
The RVE is compressed via applying displacement boundary conditions on all its faces.Since the morphology is periodic, an exact boundary condition to use would be a periodic one.For simplicity we choose to use an approximate boundary condition of prescribed velocities on the faces: where x i , v i , l, v 0 are the components of the location, components of velocity, RVE edge length, and boundary velocity respectively.The conditions are implemented using moving rigid walls which compress the RVE.The contact forces between the grains and the rigid walls are then measured and averaged by area to obtain the pressure.
Two main issues may arise from the choice of boundary conditions: (1) dynamic effects may occur due to the finite time is takes the elastic wave to propagate throughout the RVE, and (2) the approximate boundary conditions may give rise to errors in the pressure.These issues are addressed in convergence tests in the results section.

Material model
We use Quartz as the bulk material from which the sand is composed.Since the elastic behavior of Quartz is almost isotropic, we assume isotropic linear elasticity with respect to the true-strain measure.The averaged bulk and shear moduli used in the simulation are G = 46.9GPa, K = 37.7 GPa [11].The bulk density is ρ Quartz = 2650 kg/m 3 and corresponds to initial sand density of The granular material undergoes irreversible deformation in the form of pore-closure.This irreversibility is caused by fracture and micro-plasticity in the grains.In order to obtain a first approximation of this irreversibility we use a simple continuum plasticitymodel.Some energy is stored elastically and released upon unloading.A simple linear-hardening plasticity model is used where σ, σ y , β, and ε p are the stress, the yield stress, the hardening coefficient, and the plastic strain measure respectively.

Simulation parameters
Simulations are performed in the commercial explicit finite-element code LS-DYNA.The geometry is meshed with 4-node constant-strain tetrahedron elements, and the ratio of element-size to RVE length is taken as 1/40.Sensitivity of the results to the stiffness of the contact algorithm between the grains was tested and shown to be low.

Results
A representative tri-axial compression simulation of a morphology of 80 grains is presented in Fig. 2. The initial contact area between the spherical grains is theoretically zero and the forces between grains cause high local deformation (as classically studied in the Hertz contact between two spheres).Most of the mechanical work done on the RVE is transformed into irreversible plastic deformation and some into elastic energy.As compression increases, the grains deform, the contact area between them increases, and the pore-space decreases.Finally all the pore-space is filled and the spheres reach full contact.At that point, the response of the granular material degenerates into the bulk material response.Any additional work on the RVE goes only into elastic energy of the bulk-material.Thus, the process of compression can be viewed as a combination of two processes: an irreversible pore-closure process, and a reversible elastic deformation of the bulk material.

Convergence tests
Convergence tests were performed to verify the level of sensitivity of the simulation to the approximate boundary conditions and the dynamic wave-propagation in the RVE.
To examine the effects of the boundary conditions, two simulations describing the uni-axial compression of RVEs of a different size were performed.Morphologies from ensembles of 80 and 640 grains were studied with an arbitrary constant yield stress Quartz (σ y = 200 MPa).As the effect of the approximate boundary condition increases with the surface-area to volume ratio, it is expected that any significant error originating from them will cause a deviation between the results of the two simulations.Fig. 3 presents the pressure vs. density curve of the simulations.It is clear that the sensitivity to the ensemble size is very low.Thus, an ensemble of 80 grains is sufficient.A small initial ramp in the pressure is visible in Fig. 3.This is a result of strength effects as the loading is not hydrostatic.Strength effects will be studied in a future work.Simulation results at different strain-rates are presented in Fig. 4. The smooth increase in pressure at the beginning of the loading shows that there is sufficient time for the elastic waves which form at the moving boundaries to reverberate and reach equilibrium in the grains.In addition, the curves at different rates coalesce.Hence, a strain-rate of 100 µs -1 is sufficiently converged, and is used in all other simulations presented here.

Calibration of the grain material model and the Loading EOS
Two parameters of the bulk material require calibration, both define the plastic dissipation: σ y , β.In Fig. 5 we present the pressure-density curves for several parameter sets.The low yield-stress and the significant hardening required for calibration (σ y = 40 MPa , β = 400 GPa) may insinuate that the spherical shape of the grains gives rise to an underestimation of the increase in contact area between the grains.Alternatively, using a fracture model for the grains instead of plasticity might deliver better calibration parameters.Since calibration is performed on a single realization of the mesoscopic structure, it is important to validate that the variation of the calibrated EOS among different realizations is low.The sensitivity to the randomness of the grain morphology is examined on 9 realizations of grain ensembles containing 80 grains.The pressuredensity curves of all simulations along with the deviation from the average EOS are presented in Fig. 6.The variation in pressure is up to 5%.Hence, the accuracy of the calibrated EOS can be roughly taken as 5%.

The unloading EOS
The unloading portion of the EOS is computed in simulations of a loading stage up to a prescribed density, followed by an unloading stage down to zero pressure.The unloading portion is obtained in simulations by reversing the velocities on the boundaries.This yields different unloading curves for different peak previous densities (ρ max ).Generally, the EOS can be described by Notice that according to this model, a soil that passed a loading-unloading sequence will, when loaded again, follow the unloading branch of the EOS until it reaches the peak previous density, and then continue on the loading branch.Fig. 7 presents the unloading curves from 8 different densities.The simulation shows that unlike the loading process which exhibits both accumulation of elastic energy and the plastic dissipation, the unloading process is purely elastic.The initial slope represents the drop in forces between grains.As unloading progresses, some grains disconnect, and as a result, the overall stiffness decreases down to zero-stiffness -the point at which the grains are not in contact.( Since the model is calibrated from the loading curve, the speed of sound is a useful parameter for validating the model.SJOBO experimental data and prediction of the herein model for the volumetric speed of sound at the beginning of an unloading phase are presented in Fig. 8.Note that the five lower density points of the SJOBO data were calculated from wave propagation measurements.Since the experimental apparatus did not enable pressures above 100 MPa (corresponding to a density of 2200 kg/m 3 ), data points of higher density were extrapolated linearly.The two upper data points were taken as the sound of speed of Peaks Pike Granite [12] (assuming full pore closure) and show a clear discontinuity which suggests that the extrapolation is inaccurate.Our Mesoscopic model presents an increase in the sound of speed up to the Quartz speed of sounds at the point of full pore-closure.The reasonable fit between the speeds of sound predicted by the model and the SJOBO measurements supports the accuracy of the model and shows that the EOS can be reasonably predicted even with simplified physics and grain morphology.

Extrapolation of the model to soils with higher initial density
The density of the soil depends on the tightness of packing and the grain size-distribution.Other properties of the soil may also vary due to gain size-distributions.For example, an ideal one-sized spherical grain ensemble in a tight random packing reaches a volumetric packing of 62.5%.A random tight-packing of spheres with a bimodal distribution may theoretically reach up to 85.9% [10].One useful application of the model is to study how the density of the soil affects the EOS.In order to investigate this effect we choose a packing algorithm which assures a tight packing with a relatively low number of spheres.The initial grain ensemble is of tightly packed equally-sized spheres.Then, an unoccupied point in the RVE space is chosen randomly and a sphere is placed in it with initial random velocity.The sphere is then increased in size while colliding with adjacent spheres up to the point of locking.Then, if the sphere fulfils a certain size criterion, the process is repeated until the desired density is reached.Fig. 9 presents the progress of the packing as a function of the number of grains in the RVE.Notice that reaching very low porosities require a large number of grains, i.e., grain distributions with very large deviations in sizes.Simulations results for the compression of low porosity (high-density) sand with a porosity of ϕ 0 = 0.3 are presented in Fig. 10 along with the standard-density sand.The low-porosity sand presents a higher stiffness, initially threefold the stiffness of the regular sand.As expected, the stiffness increases asymptotically up to the stiffness of Quartz.Note that this result is applicable to sand that has a low porosity due to its grain-grading.Reduced porosity may also occur due to past loads as analyzed in section 3.3 in this work.To ease the comparison between the two sands, a translation along the abscissa of the low-porosity curve down to the initial density of the regular sand is plotted in a dashed green line.

Fig. 1 .
Fig. 1.Right -the initial conditions for the packing algorithm; left: the periodic packed-sphere morphology.

Fig. 4 .
Fig. 4. Pressure-density curves at different strain-rates.The EOS of Quartz is represented by a Purple line.

EPJFig. 5 .
Fig. 5. Calibration of the EOS: pressure-density curves at different bulk-material parameters.In black markerexperimental measurements of SJOBO sand.The black dotted line represents the extrapolation used for the sand model in the AUTODYN hydrodynamic code.

Fig. 6 .
Fig. 6.Pressure-density curves for different realizations of the mesoscopic structure with calibrated parameters.

Fig. 7 .
Fig. 7. Pressure-density curves for unloading from different densities.The volumetric speed of sound c in the soil is related to density-pressure curve of the unloading by c 2 = dP/dρ.(5)

Fig. 8 .
Fig. 8. Volumetric speed of sound c at the peak previous density ρ max .

Fig. 9 .
Fig. 9. Packing fraction as a function of the number of grains.

Fig. 10 .
Fig. 10.Pressure-density curves of low and regular porosity sands.To ease the comparison between the two sands, a translation along the abscissa of the low-porosity curve down to the initial density of the regular sand is plotted in a dashed green line.