Discrete element modelling of granular materials incorporating realistic particle shapes

This paper proposes an approach to generate realistic particle shapes considering the major plane of orientation of particles in discrete element modelling (DEM). The particle generation framework includes capturing high-quality scanning electron microscope (SEM) images, followed by image processing and generation of clumps using a commonly used multi-sphere (MS) approach in particle flow code (PFC3D). A set of experimental direct shear tests (DST) and subsequent DEM simulations were performed by incorporating realistic particle shapes. The simulation results show a good agreement with those obtained in the laboratory. In addition, the normal stress showed a significant effect on the structural anisotropy of the granular materials.


Introduction
Morphological features of granular materials have attracted many researchers in the field (e.g., [1,2]); [3,4].It was evident from the experimental investigations that the mechanical behaviour of granular materials including compressibility, shear strength, dilation, and crushability, was highly influenced by the morphological features (e.g., [5][6][7][8][9]).However, performing laboratory tests is time-consuming and may not be feasible to relate the morphological features of realistic particles.
For continuum-based numerical methods such as finite-element or finite difference methods, the physics occurring at the grain scale is difficult to capture because of its discrete nature and may not be possible to explore microscale insights.On the other hand, the discrete element method (DEM) pioneered by Cundall and Strack [10] has progressed rapidly over the decades and is now capable to capture microscale response of granular materials (i.e., fabric orientations, force chain networks and associated displacement vectors) considering realistic particle shapes.Using DEM, several researchers have considered different traditional element tests, such as triaxial (e.g., [1,4,[11][12][13][14] and direct shear tests (DST) (e.g., [15][16][17][18][19] to explore the microscopic responses of granular materials.In general, to incorporate the particle shape effect into DEM, one can adopt either of the two approaches that are commonly used.The earlier one is to introduce a rheology-type rolling resistance contact law between the particles, while the latter one is to incorporate irregular particle shapes (i.e., using multi-sphere (MS) approach in particle flow code (PFC3D)) for representing the constitutive relationships [20][21][22].When compared with the rolling resistance contact model, a more robust method is to incorporate realistic particle shapes, but how effectively the particle geometric shapes can be modelled using the various approaches is still a big concern.
In the recent past, the micromechanical behaviour of granular materials using direct shear test (DST) was explored by several investigators (e.g., [15-17, 19, 23-27]).However, very limited studies have focused on generation of realistic particle geometries to compare the effectiveness of the approach because of the difficulties in characterizing the complex particle shapes and utilization of highly sophisticated instrumentation in the laboratory [19,24].Therefore, based on the previous attempts made by several researchers, the present investigation is aimed at proposing a simple approach for the generation of realistic shapes of sand particles considering the major plane of orientations and incorporating them into a commercial DEM package i.e., PFC3D.This approach involves capturing the high-quality scanning electron microscopy (SEM) images of the sand particles (i.e., the major plane of orientation), image processing, analysis using MATLAB subroutines and finally generation of realistic particle shapes using a built-in clump logic available in PFC3D using MS approach [4].In summary, quantitative comparisons at the macroscopic level and subsequent microscopic evaluations were investigated by incorporating realistic particle shapes.

Properties of sand
Figure 1a shows the particle size distribution (PSD) of sand ranging between the minimum and maximum diameters of 0.075 mm and 5 mm, respectively.The physical properties of sand used are as follows: The specific gravity (G) of sand is 2.60, coefficient of uniformity (C u ) and coefficient of curvatures (C c ) are 4.28 and 1.21, respectively.Based on the Unified Soil Classification System (USCS), the gradation curve is classified as well-graded and it mainly consists of medium sand (i.e., greater than 60%).The mean particle diameter (D 50 ) is noted as 1.2 mm.In addition, X-ray powder diffraction (XRD) analysis results confirmed that the selected sand consists of a major proportion of quartz with minor microcline compositions, as presented in Fig. 1b.

SEM and image processing
The development of digital imaging technology has enabled the growth of advanced 3D imaging approaches along with numerical techniques, such as X-ray computer tomography (μCT) and laser scanners (e.g., [28][29][30][31][32][33]).However, these techniques require highly sophisticated instrumentation and complicated methods.On the other hand, 2D approaches which are precise and robust can be used as an effective method for particle shape characterization and understanding the interaction of contacts using DEM [34].About 10 different representative particle shapes were selected to capture the geometric particle shape of sand.These particles were cleaned thoroughly with alcohol solution before inserting them in a scanning electron microscope (SEM) JSM-6510 (JEOL).The sand particles were mounted on a specialized sample holder with adjustable height to optimize the operation.To reduce thermal damage, the air-dried particles were coated with a conductive thin film by the secondary electron signal beam under a voltage of 15 kV.For representative purposes, the magnification was fixed at 30 and subsequently yielded the image with pixels size of 1024 × 819.To obtain the particle shape/ geometry information, image processing techniques (e.g., cropping, thresholding etc.) were applied before analyzing the image using an open-source image processing package called ImageJ [35].An example of a processed micrograph of the typical sand particle is shown in Fig. 2a.A subroutine using MATLAB was written and executed to trace the boundary of the particle.The illustration of the obtained particle outline is shown in Fig. 2b.The outline of the particle closed boundary extracted from the binary image was opened up at a constant interval of 0.1 degrees in the anti-clockwise direction, as shown in Fig. 2c.To fit multiple circles to a series of data points along the periphery of the particle, an approach proposed by Gander et al. [36] was utilized and the representations are provided in Fig. 3.A detailed procedure on how to find the corner and non-corner regions, particle outline and best-fitted circles can be found in Zheng and Hryciw [37].

Generation of clumps
There has been growing interest in the packing structure of non-spherical (irregular) particles using DEM.Researchers, such as Thakur et al. [38], Indraratna et al. [24], Coetzee [23], Zhou et al. [39] and Gong et al. [17], have made attempts to model the complex particle shapes of granular materials using clump logic i.e., an approach for generating irregular particle shapes by joining and overlapping number of individual spheres of various sizes, and by assigning the corresponding radii and coordinates [40].An additional subroutine was written using MATLAB to build the clump templates.The information obtained from the subroutine includes radii and coordinates of the fitted circles (see Fig. 3) in the Cartesian coordinate system were then exported to PFC3D.A directory of the generated realistic particles (clumps) used in the present investigation is shown in Fig. 4.

DEM simulation of the direct shear test
The direct shear box with dimensions of 100 mm (length) × 100 mm (width) × 50 mm (height), split into two halves, modelled using rigid walls in PFC3D.The simulated DST specimen has two independent boxes namely 'box top' and 'box bottom' with corresponding planes covering both the ends (top and bottom).In addition, two more flanges were added (left and right) on either side of the shear box to overcome the particle expulsion during shear.Figure 5 shows the model configuration of the clumped assembly.To nullify the boundary effects, the ratio of the testing chamber dimension to the maximum particle size ratio is maintained between 7 and 8 [24,41].Following this criterion, a total of 4,989 discrete particles (clumps) with sizes ranging from 13 to 0.3 mm were randomly distributed within the domain, which satisfies the requirements of a representative elemental volume (REV) in DEM [42].The basic micromechanical parameters used are presented in Table 1.A simple linear contact law (i.e., linearly elastic in both normal and tangential directions) was implemented between the particles as adopted by many prior investigations (see [38,43,44]).Once clumps are generated, the assembly was cycled to reach equilibrium by arranging particles to form contact with each other.A built-in servo-control mechanism was activated to maintain the desired stress levels during one-dimensional compression under different normal stresses.While shearing the specimen, the 'box top' was permitted to move towards a rightward direction at a very low rate of shearing velocity (i.e., 5.0 × 10 −2 mm/s) to achieve quasi-static condi- tions.An auto timestep scaling mode is enabled to minimize the computational cost associated with DEM simulations [16,18].

Results and discussion
The experimental and subsequent numerical validations were done using the clumps generated using the proposed method.The stress-strain responses were compared at three different normal stresses i.e., 75, 179, and 318 kPa.The global shear strain ( γ s ) is defined as the ratio of displacement of the 'box top' to the initial height of the sample while the vertical strain is considered as the vertical displacement of the 'box top' to the original height of the sample.One may refer to Cui and O'Sullivan [16] and Kodicherla et al. [18] for more details on how to calculate stress and strains in DEM.

Evolution of the macroscopic response
Figure 6 shows the comparison of stress-strain behaviour of granular materials at three different normal stresses.For all the specimens, irrespective of normal stress, stressdependent and strain-hardening behaviour is observed.The volume change in terms of vertical strain exhibits the initial contraction followed by volumetric dilation.The higher the normal stress, the higher the peak stress is, and subsequently smaller dilation is Normal to shear stiffness ratio, k * (≡ k n /k s ) 4/3 Local damping constant 0.7 noticed.It is found that the DEM trends of the curves are generally in agreement with the experimental findings (except for the strain variation for 179 kPa) and minor variations are observed, which may be attributed to the rigidity of the loading plates and ignoring the particle breakage criteria during simulations.These observations are generally in agreement with Indraratna et al. [24], a study on three-dimensional DEM simulations of fresh and fouled railway ballast subjected to direct shear testing.
Figure 7 shows the best-fit lines between the applied normal stress and peak shear stress.The peak angle of shearing resistance, φ′ , is calculated as where τ is the peak shear stress and σ n is the applied normal stress.Both the experimen- tal and numerical simulation results of the peak angle of shearing resistance are found to be 31.29  , respectively.As compared with the experimental peak angle of shearing resistance, the numerical simulation results are overestimated.

Distribution of forces and displacement vectors
Figure 8 shows the distribution of normal contact forces along the XZ plane at a normal stress of 179 kPa.In PFC3D, the inter-particle forces are represented by solid lines with their thickness and colour proportional to the force magnitudes [40].The lines connecting the particles can form a force chain.At the initial state (i.e., γ s = 0% ), the normal contact forces are uniformly distributed and transmitted vertically from the top to bottom of the shear box, while shearing advances, the normal contact forces intensified in the diagonal direction (see Fig. 8b).This can be attributed to the reduction of contact numbers of each particle associated with an increase in dilation and the corresponding drop in the shear strength.
Figure 9 shows the displacement vectors for the applied normal stress of 179 kPa.As observed in Fig. 9a, at small strain (i.e., γ s = 1%), particles in the 'box top' moved  horizontally and particles in the 'box bottom' tended to move in the downward direction causing a contraction in the assembly.On the other hand, at higher global shear strain, (i.e., γ s = 18%), particles in the 'box top' tend to move up or dilate (see Fig. 9b).These microstructural features in DEM clearly illustrate the likely mechanism of volumetric changes during the shearing process, which may not be possible to capture during laboratory experimentation.

Average coordination number
The average coordination number (Z), is a microscopic scalar parameter that can provide information on the packing structure, material fabric and stability of the granular assembly, which is defined as [18,45]: where N c is the number of contacts, and N p is the number of particles in the assembly.It should be mentioned that the only contacts between the clumps are considered for the evolution of Z. Considering the friction between the particles, assuming no sliding anywhere, the number of degrees of freedom of a particle is 6, which means that the total number of degrees of freedom in the system is 6N p for a three-dimensional assembly.The number of constraints at each contact is 3, thus the total number of constraints in the system is 3N C .If the assembly is said to be statically determinate, then Z is equal to 4 (i.e., 3N C = 6N p ).According to Gong [46], if the value of Z is greater than 4, it can be treated as redundant, otherwise, it is structurally unstable.Figure 10 shows the variation of Z against global shear strain under different normal stresses.It is observed that the Z increases with increasing normal stress.However, for a given normal stress, it doesn't mean that a higher Z value can lead to higher shear stress, for independent of applied normal stresses above 4 of Z indicates structurally stable assembly at any stage during the entire simulation.
Figure 11.shows the probability distributions of Z for different normal stresses.The curves are well-fitted by Gaussian distribution functions.It should be mentioned that as normal stress increases, the distribution moves toward the right corresponding to an increasing Z.

Evolution of fabric anisotropy
The contact network-based geometric anisotropy can be defined as [47][48][49][50]: where n k is the unit normal contact vector of k with i, j = 1, 2, 3 .N c is the number of contacts in the assembly.The overall contact network can be subdivided into strong and weak force subnetworks.If the value of contact normal force is greater than that of the average contact normal force F , it can be considered a strong subnetwork, otherwise, it is a weak subnetwork [13,14,51].The fabric tensors related to subnetworks can be given as: where N s c and N w c are numbers of contacts in the strong and weak subnetworks, respectively.The principal values and principal directions of fabric tensors, φ 1 , φ 2 and φ 3 of the fabrics can be achieved by solving the eigenvalues and eigenvectors respectively.The deviator fabric φ d (= φ 1 − φ 3 ) , used as the difference between major and minor eigenval- ues of φ i , which is suggested by Thornton [52].
Figure 12 presents the evolution of deviator fabrics of overall, strong and weak force subnetworks under different normal stress conditions.It can be noticed that the overall and subnetworks are developed significantly and the degree of anisotropy increases with increasing shear strain.Also, it is found that the increase in normal stress can limit the development of fabric anisotropy, thereby increasing shear strength.Furthermore, it can be attributed that the normal stress has a significant effect on structural anisotropy and higher normal stress lead to a nominal magnitude of the deviator fabric, which is in line with the triaxial behaviour of granular materials [53,54].

Conclusions
An approach has been introduced for generating realistic particle shapes considering the major plane of orientation.An attempt has been made to see the effectiveness of the proposed method and constitutive relationships were compared.Based on the analysis of results and discussions, the following key findings are drawn: • A commonly used MS approach has been employed in PFC3D to model realistic particle shapes considering the major plane of orientations of the particle.• Irrespective of the applied normal stresses, stress-dependent and strain-hardening response is observed.Moreover, the higher the normal stress the higher the peak stress is, and subsequently smaller dilation is noticed.• Although the simulation results for shear stress are higher than those in the laboratory DST results, they are found to agree reasonably well with each other.• The normal stress has a significant effect on the structural anisotropy.In addition, both normal stress and deviator fabric are found to be inversely proportional to each other.Future research needs to focus on the development of sophisticated three-dimensional particle generation approaches to incorporate realistic particle geometry of granular materials and establish the correlations among the geometric shape parameters (i.e., aspect ratio, sphericity, convexity and overall regularity).In addition, the particle breakage/crushing effects should also be considered when the normal stresses reach a cutoff level for a given granular material.

Fig. 1 Fig. 2 Fig. 3
Fig. 1 Characteristics of sand used: a particle size distribution; b XRD analysis

Fig. 6
Fig. 6 Macroscopic comparison between experimental and DEM results: a shear stress against global shear strain; b vertical strain against global shear strain

Fig. 7
Fig. 7 Best-fit lines between normal stress and peak shear stress

Fig. 8
Fig. 8 Distribution of contact normal force networks: a global shear strain of γ s = 1%; b global shear strain of γ s = 18%

( 1 )Fig. 9
Fig. 9 Displacement vectors of particles: a at a global shear strain of γ s = 1%; b at a global shear strain of γ s = 18%

Fig. 10
Fig. 10 Evolution of average coordination number