Influence of stress-path on pore size distribution in granular materials

. Pore size distribution is an important feature of granular materials in the context of filtration and erosion in soil hydraulic structures. Present study focuses on the evolution characteristics of pore size distribution for numerically simulated granular assemblies while subjected to various compression boundary constrain, namely, conventional drained triaxial compression, one-dimensional or oedometric compression and isotropic compression. We consider the effects initial packing of the granular assembly, loose or dense state. A simplified algorithm based on Delaunay tessellation is used for the estimation of pore size distribution for the deforming granular assemblies at various stress states. The analyses show that, the evolution of pore size is predominantly governed by the current porosity of the granular assembly while the stress path or loading process has minimal influence. Further it has also been observed that pore volume distribution reaches towards a critical distribution at the critical porosity during shear enhanced loading process irrespective of the deformation mechanism either compaction or dilation.


Introduction
Pore size distribution (psd) is one of the key microstructural components of geomaterials that influences various flow related geotechnical problems such as erosion, piping, infiltration as well as filter design [1]. Even the nature of soil moisture interaction, which is crucial for unsaturated soils, depends on the distribution of pore spaces [2].
Estimation of psd in a granular geomaterials requires intense geometrical understanding of the amorphous pores [3,4]. Triangulation or tessellation (either Voronoi or Delaunay) approach is widely used for the determination of psd in granular materials. Chan and Ng [5] used tetrahedral tessellation to estimate the pore size distribution (psd) and the pore network. Extensive image analysis technique was employed by Sweeney & Martin [6] to capture the psd of a simulated particulate assembly involving discrete element method (DEM). Reboul et al. [7] introduced a novel pore merging criteria among the neighbouring pores and improved the psd estimation process. Alternatively, a maximum ball algorithm was introduced by Silin and Patzak [8] and further improved by Dong and Blunt [9] to quantify the pore network. Recently Gao et al. [10] and Vincens et al. [11] performed rigorous analysis such as identifying the specific geometry of pore spaces, and the connectivity using a tube network of pore spaces for the realistic characterization of psd.
Despite a large number of noteworthy contributions in the area of psd analysis, only few have attempted to study the effect of external loading on the evolution of psd(s). Kuhn & Bagi [12] identified the association of pore volume evolution with the deformation mechanism of granular materials using void cell distortion technique. Sánchez et al. [14] argued that the variation in pore size or constriction size distribution is a result of a microstructural alteration under external mechanical loadings. Some studies on the pore distribution by Kang et al. [3] and Shire et al. [15] demonstrated the relationship between pore orientation and principal stress anisotropy during the shearing of granular materials. In the present study, we numerically simulate various loading conditions in DEM and simultaneously track the psd evolution. The aim is to capture the characteristics of psd at various stages of deformation in a granular material, precisely to understand the effects of stress path.

Numerical simulations of compressive loading
We employ Particle Flow Code in Three-Dimensions (PFC-3D, ver. 5, 2015) that implements DEM for the present numerical analysis. The code generates cylindrical (76 mm in height and 38 mm in diameter) assembly of rigid spherical particles. Two different packing of particles are constructed, dense specimen with porosity 0.4 consists of 4300 particles and loose specimen with porosity 0.45 consists of 4080 particles. Rigid boundary walls and linear inter-particle contact law is considered for the simulations. Other relevant parameters for DEM simulations are listed in Table 1.  The loading program consists of conventional drained triaxial compression test (TXD), oedometric compression test (ODC), combined oedometric and triaxial compression (ODC+TXD), and isotropic compression test (ISC). TXDs and ODCs are performed on both loose and dense specimens, while ODC+TXD and ISC are performed only on the loose specimen. 500 kPa and 100 kPa initial isotropic confinement is applied prior to the deviatoric compression in TXDs. Fig. 1 shows the volume compression curves for all the three loading programs. The volumetric responses of TXDs follow dilative and contractive nature in the case of dense and loose specimens respectively. In both the TXDs, irrespective of the initial density, a residual porosity is achieved, which is an indication of the critical-state condition. On the other hand, the volume compression curves in ODCs and ISC converge towards the normal compression line for both the loose and the dense samples. 3 Evolution of pore size distribution

Algorithm for pore size analysis
In order to quantify the pore sizes within the particle assembly generated from DEM, a Delaunay tessellation based triangulation is adopted (see the typical sample cross section after triangulation in Fig. 2a). In this process, every particle center virtually connects with three other neighbouring particle centers forming a tetrahedron (Fig. 2c). The void space is the volume of tetrahedron after subtracting the solid part occupied by those four particles (Fig 2d). Further, we employ a merging technique, which initially proposed by Reboul et al. [7] to identify the overlapped void spaces within the tetrahedrons. The algorithm identifies the overlapped pores by measuring the central distance between the equivalent void spheres. If the sum of the radii of any two such void spheres is greater than the central distance, those two void spaces are merged into a single void space (Fig. 2e). Thereafter, the search for overlapped pores is further carried out to the next neighbouring Delaunay cell and so on. However, unlike Gao et al. [10], we do not distinguish between the various levels of pore overlap to identify the discrete pore size and constriction size. Fig. 3(a-b) represents the initial and final psd(s) for triaxial compression tests carried out for two different initial confinements. The evolving psd(s) indicate irrespective of the initial particle packing density of the specimens, near-identical probability density distribution is attained at the critical state i.e. at porosity 0.426 for 100 kPa confinement and at 0.42 porosity for 500 kPa confinement. In the case of dense specimens, the sample dilates during shearing, eventually the small pores diminish, and more number of larger pores evolve and hence modal pore density reduces. Alternatively, in the case of the loose specimen, disintegration of big pores during contraction generates new small pores and modal pore density slightly increases. Pore radius (mm) Fig. 3. PSD evolution during triaxial drained compression at the initial confinement of (a) 100 kPa; (b) 500 kPa.

Evolution of pore size distribution in ODC
Oedometric compression results in a contractive pore volume change for both loose and dense sample (see Fig.  4). However, unlike TXDs, here the evolution of psd(s) demonstrates significant alteration, in terms of the shape of the curve, modal values, and the peak. Such distinct shift in the modal values (~0.25 mm) in ODCs indicates complete alteration in pore distribution and not just a local pore rearrangement as previously observed by Reboul et al. (2008) for randomly packed spheres in the absence of external loading. The climbing peak value in the psd(s) is the evidence of increasing smaller pores because of pore compression and large pore disintegration. The porosity reduction in ODC also decreases the relative dispersion in psd. Additionally, multimodal nature in the psd(s) at the low porosity (0.3 to 0.24) condition is also noticed. ISC and dense ODC type loading show qualitatively similar trend in the psd evolution as the loose ODCs, thus not presented here.

Stress path dependency on psd
Stress path dependency on the characteristics of psd(s) is analysed by comparing the pore distributions evolved through various loading programs at identical porosities. However, such choice may lead to large variation in stress-states. In Fig. 3 it has already been demonstrated that at the critical state condition both the loose and dense samples converge to an unique psd during TXDs. Similar to the TXDs, the pore distributions obtained from ODCs are compared in Fig. 5a for loose and dense samples at porosity 0.35. Despite, completely different stress conditions in these cases, the psd(s) are matching.
Interestingly one can notice the differences in the pds(s) obtained from TXD and ODC given in Fig. 5b at porosity ~0.42. Despite the similarity in the overall shape of the psd(s), they are not as identical as we observed at the critical state in TXD. Similarly, the comparison given in Fig. 5c between the psd(s) from oedometric compression and isotropic compression show dissimilar characteristic at porosity 0.4. These observations indicate psd(s) emerged from different loading process are not same even at identical porosity and stress state. Similar stress paths possibly leads to similar microscopic stress orientation within granular assemblies. Earlier studies [3,16,15] already established that principal stress direction governs pore orientation. Hence, it is likely that the relative magnitude of the principal stresses dictates the pore sizes as well. However, the influence of the stress path on the variation of the pore distribution is less as compared to the porosity variation.

Conclusions
In the present study, we demonstrate that porosity is the major factor that dictates the characteristics of pore size distribution in deforming granular materials irrespective of its initial packing. In addition, minor influence of stress path on the pore size is also noticed. Numerical simulations using DEM show that, a unique pore distribution is obtained at the critical-state during drained triaxial compression. Similarly, oedometric compression on loose and the dense specimens converses to identical pore distribution at the same porosity level. However, variation in stress path leads different pore distribution even at the same porosity level.