Investigation of the Sand Porosity via Oedometric Testing

Investigation of the porosity of Klaipeda sand by odometric test is presented. The Klaipeda sand is typical Baltic seashore sand consisting of grains the average diameter of which varies from 1.18 mm to 0.3 mm. Variation of the porosity during odometric compression of the whole mixture of the sand and of three separate fractions were investigated experimentally. Porosity was characterized by the maximal (initial) and minimal (after compression) values of the void ratio. It was proved experimentally that porosity of the sand mixture is practically predicted by the coarse-grained fraction with grain diameters ranging between 1.18 mm and 0.6 mm. The role of microstructure in the densification mechanism is explained by employing the discrete element method simulations. The spherical particles and commercial EDEM code were used for modeling purposes. Discrete element method simulations confirmed generally the macroscopic experimental results and yielded additional data on microscopic behavior. The non-smooth deformation behavior was observed during detailed numerical time-history analysis. The detected instabilities are explained by rearrangement of the sand grains.


Introduction
A proper understanding of the mechanical behaviour of soils is of the major importance in the construction engineering, including construction of bridges and roads. In major cases, the soil (i.e. a basement for engineering structures) is treated as a continuum. However, the non-conventional behaviour, i.e. the transient state of the soil varies between the liquid and deformable solid states. Rigid particles, voids between them that could be fi lled by water and/or gas are the main contributors, conditioning the state of stress and deformability properties of this complex mixture in whole.
Th e macroscopic behaviour of soil as homogeneous compressible continuum is described in terms of the internal state variables, namely, stresses, internal pressures, strains, displacements, velocities, and etc. Th e relationship between them is governed by the various models which are ranging from the simplifi ed engineering approaches up to the sophisticated elastic-plastic-viscous theories. Description of them may be found in numerous textbooks, e.g. in Terzaghi et al. (1996), Ortigao (1995), Mitchell (1993), Fredlund and Rahardjo (1993).
Porosity (a parameter describing the void contribution to soil volume) variability is of a major importance when describing its deformability in stress ranges compatible with conventional design.
Verifi cation of the soil models requires experimental evidence. In major cases, the physical and mechanical properties of soils are site-specifi c (Juknevičiūtė, Laurinavičius 2008;Ždankus, Stelmokaitis 2008). Th erefore, usually the requested information on the soil properties is obtained for each construction site. Th is procedure is relatively expensive in both fi nancial and time resources. Th erefore, any approach (e.g. numerical simulation) minimizing the extent of fi eld experiments for identifying soil mechanical properties is actually expected.
Recently the behaviour of soils has been extensively investigated via analytical, experimental and numerical approaches. It is worth to note that standard well defi ned testing procedures (Head 1986) have been developed for the experimental evaluation of the soil behaviour. Th e oedometric, the triaxial and the direct shear tests are at present the most common tests for determining the macroscopic soil parameters in laboratory (Amšiejus et al. 2009;Amšiejus, Dirgelienė 2007;Lade, Prabucki 1995;Lade, Wasif 1988;Peric, Su 2005;Vervečkaitė et al. 2007). Th e oedometric test is acknowledged to be the most wide-ly employed general method for evaluating the soil compressibility by measuring porosity changes vs. pressure.
Despite the huge eff ort, adequacy of theoretical models of the soils to reality is still problematic. Soils in general and sands in particular are of the discrete nature and present heterogeneous compositions of grains. Th e continuum approach used in a frame of the classical mechanics does not provide insight into variation of porosity occurring by the local motion at the scale of individual grain.
Th e experimental investigation of particle's interaction is quite costly, therefore, application of numerical simulations are frequently involved to enhance macroscopic physical experiments. Th e discrete (distinct) element method (DEM) became recognized as a tool for simulating the particulate matter aft er the publication of the work by Cundall and Strack (1979). Th e DEM concept off ers a unique approach capturing the various particle shapes and physical models by a discrete set of the quantities. Th e fundamentals of DEM and the important details may be found in books and review papers of Allen and Tildesley (1987), Pöschel and Schwager (2004), Džiugys and Peters (2001), Kruggel-Emden et al. (2007), Zhu et al. (2007), Zhu et al. (2008).
Discrete concept and available commercial DEM codes like the DEM Solutions, 2009. EDEM v2.2, Edinburgh, UK, DEM Solutions L;Itasca, 2003, PFC 3D User's manual, v3.0, Minneapolis, Minnesota, USA, Itasca Consulting Group Inc. and numerous non-commercial research codes made DEM as an attractive research tool.
At the same time, computational capabilities limiting a number of particles remains the main disadvantages of the DEM technique. Simulation technique and soft ware issues are discussed by Raji and Favier 2004, Balevičius et al. 2006, Kačianauskas et al. 2010.
Th is tool is effi cient when simulating micro-scale response and explaining the physical nature of sand as composite discrete material.
Some applications of the DEM to modelling of dry sands may be emphasised. Th ey comprise behaviour of sand in technological processes such as fi lling and compaction of sand by vibration (Rojek et al. 2005) or simulation of various tests. In many cases, three-dimensional sample in a form of rectangular parallelepiped was examined. Th e eff ect of the use of fl at boundaries during one-dimensional compression was considered by Marketos and Bolton (2010). A cuboidal sample to simulate biaxial compression with rigid walls was studied by Yan et al. (2009). Similar approach in numerical simulations of the triaxial compression was employed by Belheine et al. (2009). DEM simulations of standard oedometric tests with cylindrical samples were presented by Oquendo et al. (2009).
Th e paper presents experimental investigation of the Klaipėda sand by oedometric test supported by the DEM simulations. It is organised as follows. An investigation concept and the basic data are given in section 2. Experimental set-up and testing results are given in section 3. Th e DEM methodology with particular emphasis on de-scription of motion and inter-particle forces is presented in section 4. Th e simulation results and the discussion are given in section 5; conclusions are presented in section 6.

Porosity problem and basic data
Th e state analysis of soils requested for engineering application is much more diffi cult when compared to structural analysis. Th e deformation of the traditional construction material such as steel or concrete occurs without significant volume change because of relative material continuity. However, sand as the natural granular material possesses the considerable compressibility and is subjected to densifi cation during deformation with a continuous change of volume. Th erefore, investigation of porosity is the important step related to fi nal evaluation of the stressstrain state of sand and its physical constants. It could be reminded, that densifi cation of the sand layers predefi nes the settlements of the structures and evaluation of the porosity is compulsory in foundation design (Murthy 2002;Terzaghi et al. 1996).
Moreover, all the natural geological materials especially sands are characterised by diversity of the grains sizes and shapes, conditioning the actual physical and mechanical properties. Since the microstructure of soils is geologically site-specifi c, therefore, their properties should be tested for each general construction site.
Th e Klaipėda sand as typical Baltic see-shore sand is under investigation. Its mineralogical composition consists basically of dominating ingredients, namely, of ~ 85% silica and of ~ 6% sunstone with remaining contribution of carbonate, mica and some other minerals.
Composition of the Klaipėda sand was evaluated via standard granulometric testing according to Geotechnical Investigation and Testing. Identifi cation and Classifi cation of Soil. Part 2: Principles for a Classifi cation. When considering the sand microstructure, it was found that the grains of sand are of a rounded shape, the grain surface is of abrade smooth character.
Average diameter d particles vary in the range from 1.18 mm to 0.3 mm. Size distribution of sand grains is presented in Fig. 1. Th e sand properties due to above mentioned standard are: uniformity coeffi cient C U = 1.47, coeffi cient of curvature C C = 0.93. Th e porosity in current investigation is characterised by a void ratio e obtained as the ratio of voids volume to the volume of the grains V gr , m 3 , i.e.: Finally, the void ratio is (1) where -the volume occupied by soil as solid material which is measured in oedometric device, m 3 .
The current study is aimed not only to quantify max/min porosity of the dry virgin Klaipėda sand under uniaxial compression but also to explain the role of microstucture in the densifi cation mechanism on the basis oedometric testing and employing the DEM simulations.

Experimental set-up
Th e uniaxial oedometric compression tests have been carried out in standard manner following (Atkinson et al. 1997).
An oedometric testing device consists of a metallic cylinder being fi xed to a rigid base (Fig. 2). Th e internal volume of the cylinder is of height H = 35 mm and diameter D = 70 mm. It responds to the volume of the sample of the soil. All cylinder walls are assumed to be rigid. Th e friction between the sample and the device walls is not eliminated.

Fig. 2. Oedometric device
For the experimental procedures the sand was dried up and test samples have been prepared by the fi lling. Th erewith, the grains were slowly poured into the cylinder of device. Th e top surface was fl attened by removing horizontally the grains above the surface. Th en, an initial vertical pressure (loading) was applied on the top wall by the rigid stamp. Th e loading was realized via the step-bystep procedure therewith each loading increment was imposed by a specifi ed loading rate and kept constant aft er some time to relax a sample. Th e incremental loading was continued until the required fi nal loading level has been reached.
Th e axial deformations were measured by the controlling the height decrease of the sample height Δh during the test.
Evaluation of void ratio e described via Eq (1) was done by having employed the following assumptions: the volume of sand grains remains constant  V gr = const, i.e. their volume changes due to contact deformation are negligibly small; the particles are homogeneuous and have the con- stant density of the grains ρ; the macroscopically observed volume changes  of the specimen are specifi ed by reduction of pores because of particle's rearrangement. As a consequence of the fi rst two assumptions, the volume of grains V gr = const could be calculated from the weight of the tested material.
Th e initial void ratio e 0 is characterised by its max value directly obtained by expression (1) assuming that entire volume of the cylinder is occupied by granular solid, i.e. V sol = V cyl .
At each load level i defi ned by pressure increment Δp, the porosity e i is related to the measured height h change Δh i of the specimen. It is calculated as follows (2) Oedometric modulus E 0 (explicitly employed in geotechnical engineering) can be obtained in the same manner by , ( where Δσ zi = Δp i and Δe i -specimen axial stress and void ratio change, respectively.

Experimental results
Th ree samples denoted hereaft er as Sample 1, Sample 2 and Sample 3 were prepared from three segregated fractions of the sand grains, namely, the coarse-grained fraction consisting of relatively large particles, varying in the range between 1.18 mm and 0.6 mm; the medium-grained fraction consisted of particles varying in the range between 0.6 mm and 0.425 mm and the fi ne-grained fraction consisted of relatively fi ne particles varying in the range small grains varying in the range between 0.425 mm and 0.3 mm were extracted and tested. Sample 4, the composition of all three fractions of equal parts, i.e. consisting of particles varying in the range within 1.18 mm and 0.3 mm, was also prepared and investigated in addition. Each of the above mentioned samples were subjected by the identical loading history up to a max pressure magnitude p = 400 kPa (Fig. 3). Th e load was increased via the four loading steps specifi ed by the equivalent pressure increment Δp i = 100 kPa, ensuring the axial strain rate of 0.11 mm/min. Th e relaxation at each loading increment lasted 60 s. Th e measured decrease of the height Δh of solid material, e.g. for Sample 1, during the testing is shown in Fig. 4.
For defi nition of the densifi cation nature the graphs presenting variations of the void ratio e and oedometric modulus E 0 for various loading levels can be given. Th e above mentioned parameters obtained by expressions (2) and (3) on the basis of experimental results for Samples 1 and 3 given in Figs 5 and 6 for illustration.
Th e graphical summary of porosity results is illus trated in Fig. 7. Th e three-dimensional histogram represents the max (initial) and min (the end of test) values of the void ratio of segregated fractions (Samples 1-3) and that of mixture (Sample 4) of sand.
Results show that max porosity (e 0m = 0.844) as well as min porosity (e 04 = 0.714) of the sand mixture (Sample 4) are practically predicted by the coarse-grained fraction (Sample 1).

Th eoretical background
Th e DEM methodology considered in this paper is aimed at simulating the dynamic behaviour of dry non-cohesive frictional visco-elastic particles. Sand sample is presented as a system of the fi nite number of particles which is synonymous to discrete elements. Each particle is treated as deformable body with the given geometry and material properties. As the particles move, they impact each other and undergo the contact deformations. Th e particle shape is restricted to spheres. Th e DEM methodology employed in below simulations is based on the original proposals of Cundall and Strack (1979) while an evaluation of some force components involves later developments of Tsuji et al. (1992), Raji and Favier (2004). Th e latter was implemented into the commercial code EDEM developed by DEM Solutions Ltd. (DEM Solutions 2009).
Th e DEM proceeds separately each particle denoted hereaft er by a subscript i. Th e motion of the particle as a rigid body in the global Cartesian coordinates is described using a framework of classical mechanics. Th e translation- . al behaviour of arbitrary particle i is characterized by the global parameters: positions X i , velocities and accelerations of the mass centre; as well as the resultant force vector F i acting on the particle. Th e motion of particle i in time t obeys the Newton's second law and is formulated for the mass centre of the particle as .
( 4) Th e equations describing the rotational motion depend on the shape of particle but for spheres they may be considered in the same way as translation. Th e rotation is governed by three independent rotational degrees of freedom θ i (t) presenting Euler angles, angular velocities w(t) and accelerations (t). Finally, rotation equations related to global coordinates of particle are: , where m i -the mass, kg; I i -the moment of inertia, of particle i, m 4 . Vectors F i and T i present the resultant force and torque acting on the particle i. Th e gravity force and all contact forces between contacting neighbour particles will be taken into account. Methodology of calculating contact forces and torques in (4-5) depends on the particle's size, shape and mechanical properties as well as on the constitutive model of the particle interaction.
Contact between the particles i and j will be considered to illustrate DEM methodology. Geometry of particles is defi ned by the radii R i and R j while elasticity properties by elasticity moduli and Poisson's ratios E i and E j , ν i and ν j , respectively. Interaction between colliding bodies is defi ned by the coeffi cient of restitution e end coeffi cient of friction μ.
Th e force and displacement vectors at the particle contact point of particles i and j can be separated into the normal and tangential components to the contact surface which are denoted hereaft er by the superscripts n and t, respectively. Th e normal direction of the contact surface is defi ned by the unit vector n ij extending through the centre of the contact area C ij . Th e normal contact direction n ij for the spherical particles always coincides with the line connecting the particle centres. Th e unit vector t ij of the tangential contact direction is perpendicular to n ij .
Th e particle's deformation due to collision is assumed to be approximated by the overlap of the particles (Fig. 8). Th e depth of the overlap h ij is considered as normal displacement component. It is provided that overlap is signifi cantly smaller than the particle size, h ij << R i , R j . Th e size of the overlap h ij is defi ned by considering the distance between the centres of the spheres . (6) During a contact, the particles move relatively to each other with velocity along the distance in the tangential direction, termed hereaft er as tangential displacement. It is obtained by integrating velocity as .
( 7) Th erewith, partial slip in the case, if the tangential force exceeds the limit defi ned by static friction is captured.
Decompose force vector into a normal and a tangential components by . (8) Th en both normal and tangential components present a combination of the elasticity, damping and separation or sliding eff ects. An inter-particle contact is defi ned by stiff ness coeffi cients and , by the damping coeffi cients and and the friction coeffi cient μ, respectively.
Th e normal elasticity force follows Hertz model while damping force is defi ned by the Tsuji model (Tsuji et al. 1992). Explicitly, the normal force is expressed as , where the normal stiff ness parameter related to the eff ective Young's modulus E eff given by .  sliding. Th e static force comprises elastic and viscous ingredients identical to those of normal force (9), while the dynamic force confi rms the friction force expressed by the Coulomb' low. Th e tangential force can't exceed the friction limit, the real force responds to the min of the two above forces: . Explicitly . (10) Formally, the right-hand terms in Eq (5) present a composition of two independent torques. Th e contact torque (Fig. 8) with particle j is defi ned as vector product as originally proposed by Cundall and Strack (1979), namely: . (11) Th e second term in Eq (5) presents rolling resistance torque which is introduced to account the nonspherical character of real particles. Th e latter may be interpreted as a frictional resistance torque defi ned in local plane: , . (13) where -the normal contact force; -the min radius of the contacting spheres; -a rolling friction coeffi cient (Iwashita, Oda 2000).

Simulation and analysis
Th e numerical experiments were performed to investigate the behaviour of the soil specimen. Simulation capability is limited by a number of particles and lack of inter-particle data, therefore the several simplifi cations are introduced and applied for the current DEM simulations performed by EDEM 2.2.1 Academic code. DEM solutions, 2009. Simulation of the full-scale sample would be enormously time consuming. Consequently, reduction of the problem size was done by reducing the size of computational domain. Th e cylinder of oedometric device was reduced seven times while retaining the ratio between height H and diameter D, i.e.
. Th e reduced cylinder is defi ned by H = 5 mm and D = 10 mm. Th is simplifi cation allows to retain the original scale of grains, while aspects ratio of cylinder and grain diameters is .
Th e data on particle physical properties of the investigated sand grains are defi ned indirectly. Diff erent values of the density and the elastic constants for various modifi cations of silicon grains and may be found in available internet data sources and references.
Viscous damping eff ect is evaluated by a coeffi cient of restitution. Th e average value of the normal grain-wall coeffi cient of restitution was obtained by the grainbed collision experiments performed in the wind tunnel (Wang 2008) was taken for current simulations. Th is value was also employed for characterisation of the particleparticle interaction, thus the normal restitution coefficient magnitude is c n = 0.5. Th is introduced assumption is mainly conditioned by the experimental diffi culties arising in measurements of the grain-grain restitution. Th e tangential coeffi cient of restitution e t is assumed to be fraction of the normal coeffi cient Th e role of friction of granular materials was reviewed by Zhu et al. (2007). Its infl uence on densifi cation is emphasised here. It is obvious, that because of scattering of the experimental results (Horváth et al. 1996), only a rough approximation of the friction of sand is available. In many cases (Belheine et al. 2009;Oquendo et al. 2009;Rojek et al. 2005) various values of the friction coeffi cient between 0 and 1 are employed for numerical simulations. Here, the friction coeffi cient between the particle and wall (microscopic scale) is chosen to be the same, i.e.  =  w =  gr , while a relatively low value  = 0.2 is assumed regarding the smoothness of particles. Th e data values applied in the DEM simulations are summarized in Table 1. Numerical simulation is naturally divided into two stages -generation of the initial state of particles, or packing stage and the compression stage.
Packing problem was comprehensively discussed in review paper of Zhu et al. (2007), packing of spheres in cylinder was discussed by Mueller (2005). Th e generation of cylindrical packing considered hereaft er is subjected to specifi c requirements. Adequacy of the generated sample to the experimental one is retained by the specifi ed value of the porosity e 0 . From this the fi xed volume of sand material in terms of the volume of grains V gr is recalculated.
Another issue is the design of particles. In a case of the fractioned sample, diameters of particle are restricted by the specifi ed min and max values d max and d min , respectively.
Th e initial state of the sample, i.e. location of particles in the cylinder, is obtained by the DEM, simulating the particle's deposition. Particle's diameters and placing above the cylinder are generated randomly. Later they fall free due to the gravity force forming a packing structure. Th e process of fi lling was controlled by considering the evolution of the volume (or mass) particle and interrupted when it reaches specifi ed values V gr or m gr , respectively. Finally, the total number N of particles is obtained.
As it was explained above, the dynamic state of particles at the time t is obtained by the numerical integration of the equations of motion (4-5). Th e explicit integration scheme with a constant time step t = 0.1 μs is applied. Th e magnitude of t is chosen to be a smaller than the Rayleight time step calculated by EDEM. Th e total time of fi lling t f = 0.18 s is continued to stabilise of the particles motion aft er fi lling.
Two samples of the specimens composed of the coarse-sized and medium-sized grains were generated numerically. Th e characteristic data values for both specimens are summarized in Table 2. Th e initial states of particles are shown in Fig. 9.
A simulation of compressing the specimens was performed in the second stage. It was realized in a slightly different manner comparing with the real experiment. Firstly, the compressive loading p(t) was imposed by vertical moving the upper rigid wall of the cylinder with proportional increase of its vertical displacement u t (t). Secondly, the loading history was simplifi ed by restricting to one stage load increment i.e. max load pressure 400 kPa was reached in shorter time without relaxation. An integration of the equations of motion (4-5) was performed with a constant time step, t = 10 -7 s.
Th e loading histories of both of the samples in terms of the top wall pressure p(t) are given in Fig. 10. It could be observed that in Sample 1 with the coarsegrain particles, the max pressure p max = p(t L1 ) = 400 kPa was reached aft er t L1 = 1.08 s at the lower displacement u max1 (t L1 ) = 0.076 mm ( %), while in Sample 2 with the medium-sized particles, the max pressure p max = p(t L2 ) = 400 kPa was reached aft er t L2 = 0.8 s at the lower displacement u max2 (t L2 ) = 0.2 mm. Compression results are presented in Fig. 11 in terms of the relationship between the porosity and pressure e(p). Void ratio magnitudes are recalculated by expression (2) as e(u t (t)), while pressures are taken from the DEM simulations. Here, DEM results are also compared to those being obtained experimentally (bold curve).

Simulation results and discussion
A numerically simulated macroscopic behaviour in terms of pressure-displacement relationship contains descending intervals exhibiting instability of material. Th e above character is hardly to explain from the continuum point of view but we recognize it on the results of microscopic analysis of the particles motion.
Let us consider the central vertical section of Sample 1 (Fig. 12). Th e fragment of fi ve red coloured particles is selected for analysis purposes. Th e central particle is denoted by capital letter C, while the four remainder particles by L (left ), T (top), R (right) and B (bottom).
Porosity was characterised by the max (initial) and min (aft er compression) values of the void ratio. It was proved experimentally that porosity of the sand mixture with grain diameters ranging between 1.18 mm and 0.3 mm is characterised by max (e 0 = 0.844) and min (e 04 = 0.714) values of the void ratio which are smaller compared to those of separate fractions. Moreover, these values are practically predicted by the coarse-grained fraction with grain diameters ranging between 1.18 mm and 0.3 mm.
DEM simulations confi rmed generally the macroscopic experimental results and yielded additional data on microscopic behaviour. Th e non-smooth deformation behaviour was observed during detailed time-history analysis. Th e detected instabilities are explained by failure of microstructure and rearrangement of the sand grains.
Identifi cation of the infl uence of other factors such as non-sphericity of particles, deformation rate, friction, etc., requires, however, the further research. A motion of all selected particles during compression is characterized by velocity histories v C , v L , v T , v R and u B (Fig. 13a). Particular velocity v i (t) stands for displacement magnitudes . Th e relative motions of neighbour particles i with respect to the central particle C are characterized by diff erences of the velocity magnitudes shown in (Fig. 13b). Comparing time instants of the pressure drop in (Fig. 10) it is easy to fi nd that velocity-displacement jumps indicate the bifurcation from one dynamic equilibrium to the other.

Conclusions
Variation of the porosity of the Klaipėda sand during oedometric compression was investigated experimentally and by the Discrete Element simulations. Behaviour of separate (segregated) fractions and the whole mixture of the sand were treated independently.