3D Morphology Distribution Characteristics and Discrete Element Simulation of Sand-Gravel Mixtures

Sand-gravel mixtures are typical binary materials, exhibiting highly heterogeneous, discontinuous, and significant structural effects. The contact state between sand and gravel particles has a significant influence on the mechanical properties of the mixtures. This article focused on the complex internal structure and its mesostructural behavior of the mixtures, and a systematic statistical analysis was carried out to study the shape, size, and angularity of the coarse particles. The three-dimensional (3D) shapes of coarse aggregates were approximated to be hexahedron, pentahedron, and tetrahedron. An indicator called angularity and surface texture (AT) index was developed to characterize the combined effect of the coarse aggregate angularity and surface texture. Based on the screening testing and digital image processing, the particle size and AT index of aggregates were extracted, and their means, standard deviations, and statistical distributions were studied. An algorithm for generating 3D aggregates was developed based on the statistical results of the coarse aggregate 3D morphology. The coarse aggregate generating code was written using the fish language in PFC3D. The numerical model was then applied to conduct three typical monotonic or cyclic triaxial test simulations. Retrospective simulation of the laboratory tests using the proposed model showed good agreement, and the reliability of the model is effectively verified. The results interpreted well the mechanism of particle motion and the distribution of interparticle contact force during shearing from mesoscale of the mixtures, which can give better understanding and modeling of the nonlinear behavior of the sand-gravel mixtures.


Introduction
Hydraulic fracturing is a promising method for increasing the production of oil and gas wells. Hydraulic fracturing is currently the main form of natural gas extraction. It requires a large amount of water mixed with chemicals to be poured into the shale layer for hydraulic fracturing to release natural gas. At present, most of the world's natural gas and oil wells are drilled by hydraulic fracturing [1][2][3][4]. Due to the characteristics of deep underground drilling operation of hydraulic fracturing, it often crosses undesirable engineering geology such as a sandy gravel mixture layer [5,6]. It is of great significance to guide the construction of engineering to establish a set of fast and effective methods to analyze the mechanical behavior of sandy gravel mixtures.
Sand-gravel mixtures are a special engineering geological material between the soil and the fractured rock; the sandgravel mixture structure presents the mesostructure between sand particles and gravel particles [7]. The contact state and particle morphology of the mixtures have a significant contribution to their mechanical properties [8][9][10][11][12]. If the heterogeneous and discontinuous structural characteristics of the sandy gravel mixtures are ignored and the particle morphology and particle distribution are not considered, the authenticity of the mesostructure characteristic simulation of the material will be reduced, and it is difficult to describe the microscopic mechanical behavior between the particles [13]. Kennedy and Lin [14] constructed polygons to approximate the shape of particles and found that the number of sides of the polygon is proportional to the circumference and inversely proportional to the length of the sides. This relationship is linear in double logarithmic coordinates, and the slope indicates the irregularity of the particles. Bowman et al. [15] use the Fourier transform and analysis theory to describe the particle morphology of sand. According to the research of Sukumaran and Ashmawy [16], the surface texture reflects the fine morphology of the particle surface and can be used to describe the roughness of the particle surface and the angular microscale. Cho et al. [17] described the particle morphology based on three aspects: sphericity, edge angle, and roughness. Ferellec and Mcdowell [18] replaced the average radius of curvature of the corners with the average radius of all inscribed circles on the projection contour of the particles, which improved the reproducibility of this parameter.
With the development of computer technology and the emergence of the discrete element method, numerical simulation has gradually become an important method to analyze the particle shape; the elastoplastic finite element method (FEM) was introduced in the initial stage. However, due to the discontinuity of sand, FEM based on continuous medium theory cannot reflect the influence of material parameters such as particle size and particle morphology on the simulation results [19,20]. Therefore, the discrete element method (DEM) based on discontinuous medium mechanics has been widely used in the numerical simulation of granular media [21]. This method is low in cost and easy to operate and is very suitable for numerical simulation of discrete media such as sand and gravel.
Mcdowell and Harireche [22] found that adjusting the percentage of pellets empirically reflects the size effect of individual pellet strength. For the particle cluster model, when the particle interaction force exceeds the connection strength, the small particles are detached, which can simulate local or whole fracture. Cheng et al. [23] conducted a simulation of the particle flow of crushable soil, in which the crushed particles used a cluster model of particles formed by the bonding of several small particles, and the randomness of single particle strength was reflected by randomly deducting a certain percentage of small balls. Deluzarche and Cambou [24] used the discrete element code PFC2D and considered breakable clusters of 2D balls. The different parameters were determined from experimental data obtained from laboratory tests performed on rock blocks. The model was validated by comparing the results of the simulation of shearing tests with actual triaxial tests on rockfill material published in the literature. Alaei and Mahboubi [25] focused on numerical modelling of rockfill material with the discrete element method (DEM). Because the DEM models the interaction of separate elements, it was capable of modelling discrete structures of granular materials and particle breakage. Alonso et al. [26] established a distinct element method model developed in which grains were characterized in the clump model. The resulting particle shape approaches real geometries and allows a reasonable breakage evolution, and the particle breakage criterion involves the subcritical propagation of fissures in the grain.
In this study, taking the sand-gravel mixtures distributed at the project site of Binjiang Subway Station of Metro Line 3 in Chengdu as an example, using discrete element numerical simulation and digital image processing technology, the 3D form and distribution of gravel particles are systematically analyzed. Based on the PFC3D particle flow software, a 3D discrete element simulation method for sand-gravel mixtures was proposed. The triaxial numerical test and actual triaxial tests were compared to verify the accuracy of the method and discuss the mesomechanical behavior of the sand-gravel mixtures.

Basic Physical Properties of Sand-Gravel Mixtures
The sand-gravel mixture samples were taken from Binjiang Subway Station of Metro Line 3, Chengdu. The soil layers in this site from top to bottom are Quaternary Holocene (Q h ) artificial fill, alluvial silt, sand and pebble (medium sand lenticel sandwiched in some sections), Quaternary Pleistocene (Q p ) alluvial pebbles, and Upper Cretaceous (K 2 ) Guankou Formation mudstone. Three groups of sand-gravel mixtures at different depths are used for the particle-size analysis test, and the gradation curve is shown in Figure 1. According to the Unified Soil Classification System (ASTM D2487 [27]), the sand-gravel mixtures belong to poorly graded gravel (GP) with the minimum or maximum dry density of 1.437 or 1.740 g/cm 3 , the specific gravity is 2.68, mean particle size d 50 = 6:31 mm, coefficient of uniformity C u = 34:5, and coefficient of curvature C c = 0:37.

The 3D Morphological Characteristics of Sand-Gravel Mixtures
3.1. Test Materials and Test Methods. Gravel particle morphology is an important parameter of mesophysical and mechanical properties of the sand-gravel mixtures. 110 particles with the diameter range from 50 to 20 mm, 117 particles with the diameter from 20 to 10 mm, and 128 particles with the diameter from 10 to 5 mm were selected as the research objects. The 3D shape of gravel particles was captured by digital image processing technology. In order to distinguish the aggregate from the image background more accurately, the background was set to white, and the gravel particles were blackened with ink, as shown in Figure 2.

3D Morphological Evaluation Parameters of Particles.
The important size parameters of gravel particles include the length of the long axis, the length of the central axis, and the length of the short axis. In order to measure the size of the particles accurately, first, the particles should be neatly arranged on the white paper and a high-definition image was taken. Then, the particles were tiled in situ at a different angle and photographed again. The detailed steps are as follows: 2 Geofluids (1) Select gravel particle samples of all kinds of particle sizes (2) Soak gravel particle samples with a certain proportion of ink and put them into the oven for full drying (3) Arrange gravel particle samples on the white paper neatly, and use a digital camera to take highdefinition sample images (4) Replace the particles in situ at a different angle and take a picture again (5) Use graphics software to sort out the particle images and extract 3D morphological parameters. The specific parameter determination method will be described in detail below The size characteristics of the gravel particles are reflected by the aspect ratio, the short axis length, the central axis length, and the long axis length of particles. To be specific, the long axis length L is the relatively large diameter in the 3D dimension of the particles, the central axis length W is the maximum size of the vertical elevation with the long axis in the three-dimensional image of the particles, and the short axis length t is the minimum diameter in the 3D diameter, which means t is the longest diameter perpendicular to both the long axis and the central axis. The details are exhibited in Figure 3. The size ratio between the long axis and the short axis (aspect ratio, L/t) is selected to reflect the shape characteristics of the particles.
In the particle image, the angular characteristics of gravel particles are represented by the smoothness of the contour lines. As shown in Figure 4, on the contour line of the particle image, the perimeter of the contour line increases with the increase of the number of fluctuations. Based on this, a two-dimensional angular coefficient (AC 2D ) is proposed to characterize the angularity of particles. The angular coefficient AC 2D is defined as the ratio of the difference in the surface contour circumference (SCC) and convex contour circumference (CCC) of the particles to the convex contour circumference (CCC) in the image, as shown in in which SCC and CCC denote all the boundary lines of the particles themselves and the boundary line with high flatness that is obtained by ignoring the concave parts on the surface contour line, and the curve of the concave part of the particle surface is directly replaced by a straight line as the convex contour circumference. In order to more realistically reflect the angular coefficient of a single particle, AC 2D of each side image is calculated, and the three-dimensional angular coefficient (AC 3D ) of the particle is obtained by area weighted average, as shown in where n is the number of particles, A i is the area of the image i, and AC i is the angular coefficient of the image i. Figure 5 shows the frequency histograms and normal distribution P-P graph of L, W, and t for particle size gravel of 20~10 mm. The horizontal axis and vertical axis of the frequency histogram of the 3D morphological parameters with different particle sizes represent the value interval and value frequency of the statistic, respectively. The P-P graph is a statistical graph that detects the degree of conformity of distribution types. If the P-P graph is displayed as a straight line from the origin of the vertical axis to the upper right, it means that the variable values follow a normal distribution [28]. In order to obtain the statistical distribution characteristics of L/t and AC, the P-P graphs of Laplace, Logistic, lognormal, Pareto, Student t, Weibull, and uniform distribution of L/t and AC are explored, respectively. Figure 6 shows that L/t and AC satisfy the lognormal distribution.

Results and Analysis.
Using the same method, the probability statistics of the remaining two groups of gravel particles are performed, and the statistical parameter results are shown in Table 1   able boundary conditions. At the same time, the effect of internal particles is ignored during calculation, so it can be regarded as a rigid whole for calculation. Zhang et al. [29] used clump super particles to simulate coarse aggregates based on PFC3D and proposed an algorithm for randomly generating the 3D morphology of coarse aggregate particles. The results showed that clump super particles can effectively make up the difference between traditional spherical particles and real materials. In order to improve the calculation accuracy and calculation efficiency, the author referred to reference [30,31], focusing on the simulation of the 3D shape of gravel particles (greater than 5 mm), using the equal-size particles with a diameter of 5 mm to simulate the sand particles, which reached the balance of accuracy and efficiency.
In order to better reflect the morphological characteristics of gravel particles, polyhedron (hexahedron, pentahedron, or tetrahedron) and angular coefficient were adopted to describe the 3D form of gravel particles. Hexahedron, pentahedron, or tetrahedron was used to reflect the main shape of particles. Meanwhile, the angular coefficient was introduced to describe the local morphological characteristics (smoothness of particle contour). As shown in Figure 7(a), the 3D shape of gravel particles was simulated approximately through hexahedron, pentahedron, or tetrahedron in the model. The 3D size of gravel particles was determined by the length of perpendicular edges AB, AD, and AE; meanwhile, the corner angle characteristic of particles was controlled by changing the angle of angle ADH, angle ABF, and angle ABC. Through the change of the above six parameters, the particles in different sizes and sharp edges could be generated. According to the statistical results, L, W, and t of gravel particles with different particle sizes satisfied normal distribution, and L/t and AC satisfied logarithmic normal distribution. The mean value and standard deviation are shown in Table 2. According to the normal distribution, the model randomly selected the size parameters and angular parameters of particles to generate gravel grains with particle sizes of 50~20 mm, 20~10 mm, and 10~5 mm, as shown in Figure 7(b).

Sand Particle Simulation.
The model simulated the sand particles with equal particle size with a minimum particle size of 5 mm. The specific generation process is as follows: (1) Calculate the number N of spheres to be generated according to volume V, soil porosity n, and sphere radius R (2) Use command "gen" in PFC3D to generate N spheres with a radius of 0:2r in the specimen, and assign the parameters of density and stiffness to the spheres. When the "gen" command is used to generate the spheres, the spheres are randomly generated within the sample and do not overlap each other, which ensures the uniformity of the sample (3) Use command "mult" to enlarge the radius of the sphere by 5 times. There is a certain amount of overlap between the spheres

Geofluids
(4) Use command "cycle" to run the program to eliminate the overlap between the particles, so as to obtain a uniformly distributed sand-gravel soil

Triaxial Mesoscopic Simulation Test and Analysis.
After the sample is generated, there will be overlap between the particles. The model will judge whether the sand particles and gravel particles overlap one by one. If the two coincide, delete the former. Figure 8 shows the discrete element samples of the triaxial test of the sand-gravel mixtures obtained by the above method. In order to simulate the loading process of the triaxial shear test, a square wall simulation loading indenter was generated at the upper and lower ends of the sample; a cylindrical wall used to simulate the confining pressure was generated on the side of the sample; during the test simulation, the wall at the upper and lower ends was controlled to apply an axial load to the sample at a certain speed, and the change in the radius of the side wall was controlled so that the surrounding pressure around the sample remained unchanged; the radial speed of the side wall was automatically controlled using "numerical servomechanism." The method can realize the application of confining pressure and the loading of the test sample. Through repeated adjustments, the values of the mesoscopic parameters of sand-gravel mixtures are basically consistent with the test results. In this article, a set of ideal mesoscopic parameters of sand-gravel mixture models were calibrated, as shown in Table 2. Figure 9 shows the drained shear test results and numerical simulation results of saturated sandy gravel mixtures when the initial effective confining pressure is 300 kPa. In the initial stage of loading, the deviatoric stress increases rapidly with the increase of the axial strain. As the axial strain further increases, the deviatoric stress slowly increases until the peak deviatoric stress is reached, and the volume     6 Geofluids strain decreases with the increase of the axial strain; the soil appears to shrink. After that, the stress-strain relationship curve has obvious softening characteristics. As the axial strain increases, the deviatoric stress decreases. At this time, the axial strain-volume strain curve shows obvious shear expansion characteristics. It can be seen from Figure 9 that the numerical test results are in good agreement with the drained shear test results. Figure 10 shows the undrained shear test results and numerical simulation results of saturated sandy gravel soil mixtures when the initial effective confining pressure is 300 kPa. It can be seen from the figure that when ε < 10%, the difference between the numerical test results and undrained shear test results is very small. When ε > 10%, the numerical test results is different from the undrained shear test results, and the difference increases with the increase of axial strain. This is because the sample is not completely stable when it is in large deformation, and the numerical simulation results cannot fully reflect the mechanical properties of the soil. It is worth noting that the effective stress path of the numerical test and the undrained shear test is basically the same. Figure 11 shows the undrained cyclic triaxial test results and numerical simulation results of saturated sandy gravel mixtures when the initial effective confining pressure is 100 kPa. The pore pressure ratio R u of gravel-sand mixtures presents a "fast-steady" growth pattern, which is basically consistent with the typical test results of gravel-sand mixtures. The stress-strain curve of the triaxial test and the numerical test can maintain good consistency before 7 Geofluids liquefaction, and when the sample is liquefied, the sample no longer has shear strength. Meanwhile, it is difficult for the numerical test to accurately and effectively apply axial stress and reflect the true state of the sample, resulting in a certain difference in the stress-strain curve after liquefaction.
Because the particle flow sample is composed of discrete particles, similar to the actual coarse-grained soil, the interaction of the particles in the PFC3D numerical model can better reflect the mesomechanical properties of the actual sandy gravel mixtures. In the undrained shear numerical test, PFC3D provides the display of the force between particles. From the force between particles, the mesoscopic working mechanism can be analyzed. In order to study the particle movement law of the sandy gravel mixture sample, Figure 12 shows the particle displacement vector diagram of the numer-ical sample under shear failure when the initial effective confining pressure is 100 kPa. When shear failure occurs, the directional movement of the particles is obvious, the upper particles move downward and the lower particles move in the opposite direction, and obvious shear bands appear, indicating that structural damage occurs inside the sample. Figures 13(a) and 13(b) show the distribution of sand particle contact force and gravel particle contact force, respectively. The soil composed of sand and gravel with obvious differences in particle size has a skeleton structure jointly undertaken by sand and gravel. However, in the case of this test, the distribution of the contact force between the sand particles is scattered, and the contact force between the gravel particles can form a complete contact force, so the gravel particles are mainly stressed particles.

Geofluids
The mechanical constitutive model of sand-gravel mixture established based on the 3D morphological distribution characteristics presents good applicability. The research results will provide an important scientific basis for the analysis of the stability and seismic performance of hydraulic fracturing drilling through crossing the sandy gravel mixture layer.

Conclusions
This article comprehensively used digital image processing technology and discrete element numerical simulation to study the mesoscopic structural characteristics and the PFC3D particle flow simulation method of the sandgravel mixtures taken from the Chengdu metro area in China. According to the statistical characteristics of the 3D form and distribution of the particles, a 3D discrete element simulation method for sandy gravel mixtures is proposed. The discrete element method was used to analyze the particle movement law of the sand-gravel mixtures during the shear test process. The main conclusions are as follows: (1) The long axis length L, the middle axis length W, and the short axis length t of the gravel particles in accordance with the normal distribution and the aspect ratio L/t and the angular coefficient AC satisfy the lognormal distribution. The PFC3D particle flow model can well simulate the 3D morphological distribution characteristics of gravel particles to form uniformly distributed sand-gravel mixtures (2) The deviatoric stress-axial strain curves obtained from the simulated drained shear test, undrained shear test, and undrained cyclic triaxial test are in good agreement with the test results, which verifies the applicability of the model (3) From the particle contact force diagram between the particles, it is concluded that the gravel particles in the sand-gravel mixtures are forced particles during the drained shear process

Data Availability
The data are generated from experiments and can be available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
(a) Sand particle contact force diagram (b) Gravel particle contact force diagram