Triaxial compression behavior of 3D printed and natural sands

Different particle properties, such as shape, size, surface roughness, and constituent material stiffness, affect the mechanical behavior of coarse-grained soils. Systematic investigation of the individual effects of these properties requires careful control over other properties, which is a pervasive challenge in investigations with natural soils. The rapid advance of 3D printing technology provides the ability to produce analog particles with independent control over particle size and shape. This study examines the triaxial compression behavior of specimens of 3D printed sand particles and compares it to that of natural sand specimens. Drained and undrained isotropically-consolidated triaxial compression tests were performed on specimens composed of angular and rounded 3D printed and natural sands. The test results indicate that the 3D printed sands exhibit stress-dilatancy behavior that follows well-established flow rules, the angular 3D printed sand mobilizes greater critical state friction angle than that of rounded 3D printed sand, and analogous drained and undrained stress paths can be followed by 3D printed and natural sands with similar initial void ratios if the cell pressure is scaled. The results suggest that some of the fundamental behaviors of soils can be captured with 3D printed soils, and that the interpretation of their mechanical response can be captured with the critical state soil mechanics framework. However, important differences in response arise from the 3D printing process and the smaller stiffness of the printed polymeric material. Artificial sand analogs were 3D printed from X-ray CT scans of sub-rounded and sub-angular natural sands. Triaxial compression tests were performed to characterize the strength and dilatancy behavior as well as critical staste parameters of the 3D printed sands and to compare it to that exhibited by the natural sands.


Introduction and background
The mechanical behavior of coarse-grained soils is governed by the skeletal forces that develop at the particle-particle contacts resulting from applied boundary stresses [41]. The distribution and magnitude of these skeletal forces, along with the resulting normal and shear deformation at the particle contacts, are influenced by the inherent properties of the particles, such as gradation, shape, surface roughness, and mechanical properties of the constituent materials. These particle-scale interactions control the global-scale behavior observed in laboratory tests and field conditions. Previous studies have employed experimental and numerical methods to expand our understanding regarding the effects of different inherent particle properties on the engineering properties of coarse-grained soils, such as friction angle [9,21,24,25,28,48,50,56], shear wave velocity, and small-strain modulus [5,8,12,22,27,35,43,54,58]. However, contradicting trends regarding the effect of particle properties have been reported, likely due to challenges associated with the isolation of individual particle properties for their systematic investigation in natural soils.
3D printing technology has advanced rapidly in the last decade. This technology can be used to generate parts composed of different materials such as polymers, metals, and ceramics. 3D printing has been previously used to generate artificial soil analog particles with independent control over particle shape and gradation [1,18]. Morphological comparison studies on natural sand particles and 3D printed analogs have shown that these analogs can replicate the morphology of natural sand particles successfully [1,29]. Triaxial tests on polymeric 3D printed particles of different shapes and sizes have demonstrated that the analog assemblies exhibit a stress-dilatancy behavior that is similar to that typical of coarse-grained soils [1,18,29] and revealed the dependence of the assembly stiffness on the particle shape [1,4,30]. However, 3D printed particles exhibited a more contractive initial response during shearing compared to natural particles [18], and time-dependent compressibility during consolidation [1]. Comparison of compression test results on natural and 3D printed particles showed that 3D printed particles qualitatively exhibit key aspects of 1D compression behavior of natural soils, such as different compression and recompression indices [2,17]. However, the 3D printed assemblies exhibited greater compressibility than the natural sand assemblies [2]. Bender element test results on 3D printed particles revealed the dependency of shear wave velocity and shear modulus on mean effective stress and void ratio, similar to that observed in natural sands [2]. A characterization study of the contact behavior of 3D printed spheres of different materials indicated that 3D printed particles may be used to validate discrete element method (DEM) simulations against experimental results [26]. The above findings demonstrate the ability of 3D printed particle analogs to emulate the macro-scale behavior of coarse-grained soil. However, the research to date also indicates important differences in the behavior of 3D printed and natural sand assemblies, particularly with regards to the compressibility at medium and large strain regimes.
Although assemblies of natural and 3D printed particles exhibit similarities in their mechanical behavior, further insight requires understanding of the particle-scale contact response. The normal force-deformation behavior of two spheres in contact can be described using Hertz theory, which relates the contact deformation, δ, to an applied normal force, F, as follows: where R is given by: 1∕R = 1∕R 1 + 1∕R 2 , and R 1 and R 2 are the radii of the contacting spheres. The effective Young's modulus, E * , is defined as: 1∕E * = 1 − 2 1 ∕E 1 + 1 − 2 2 ∕E 2 , where ν 1 and ν 2 are the Poisson's ratios of the contacting spheres. As shown, δ is proportional to E* 1/6 ; thus, it follows that particles with (1) δ = 9 16 smaller stiffness will undergo greater deformations under any given F. Figure 1 shows a subset of the results of uniaxial particle-particle compression tests by [2]. The force-displacement responses correspond to equal-sized spheres (3.175 mm diameter) of borosilicate glass and 3D printed polyjet polymer. As shown, the contact force-displacement response of the glass spheres is mostly elastic and follows Hertz solution at loads smaller than 80 N. In contrast, the response of the 3D printed polyjet polymer deviates considerably from Hertz solution at contact deformations smaller than about 0.06 mm due to initial plastic deformation of micro-asperities, as shown in Fig. 1b. This initial plastic response due to the deformation of micro-asperities has been reported by a number of researchers [10,13,15,42], who observed plastic yielding until the contact normal force, F, reached a threshold force, N GT . The threshold force depends on the particle's surface roughness, surface radius at the contact point, and the constituent material's stiffness [15]. Once F exceeds N GT , the force-displacement response follows that predicted by Hertz theory. This is shown in Fig. 1b by a second dotted line that represents the Hertz solution with a constant offset in the initial constant displacement. Although a Hertzian model with an offset appears to provide a reasonable approximation over the range shown, the 3D printed particle exhibits a response that is more complex possibly due to the coupling of elastic and plastic deformations. Table 1 summarizes Young's modulus, E, and N GT values for a natural sand, glass spheres, and polyjet 3D printing polymer spheres. The N GT values are smallest for glass ballotini and borosilicate glass spheres, followed by the polyjet 3D printing polymer spheres, and highest for the Leighton Buzzard sands. Since the normal contact deformation depends on E (Eq. 1), the N GT /E ratio can give insight into the relative magnitude of plastic to elastic deformations that can be expected at a given contact. The glass spheres have the smallest N GT /E ratio, indicating that their response is likely to be dominated by elastic deformations as shown in Fig. 1a. In contrast, the polyjet 3D printing polymer has the highest N GT /E ratio, highlighting the important role of the plastic deformation of micro asperities as shown in Fig. 1b. These differences in particle-level response are likely to be responsible for differences in the response of granular assemblies tested in a laboratory setting.
While a number of researchers have investigated the similarities in mechanical behavior of natural and 3D printed particles, there is need to characterize the shearing behavior across a range of effective stresses, void ratios, and imposed boundary conditions. The current study investigates the triaxial compression behavior of angular and rounded 3D printed sands. In addition, the analysis presented here examines the existence of a critical state for the 3D printed sands in the stress and compression planes. Both drained and undrained behaviors of 3D printed analogs are investigated and compared to those of natural angular and rounded particles.

Materials and methods
Drained and undrained isotropically-consolidated triaxial compression tests were performed on specimens of two types of 3D printed sand and two types of natural sand to evaluate their shear response and stress-dilatancy behavior and make comparisons between the materials. The 3D printed particles were generated from X-ray CT scans of the natural particles with the goal of reproducing their shape and size, leaving the constituent material as the main difference between them. This section provides an overview of the 3D printing technology used in this study, the materials tested, and the experimental procedures used to characterize the triaxial compression behavior.

3D printing technology
3D printing technology has advanced rapidly in last few years, and several 3D printing methods have been developed to generate parts composed of different materials. Different methods and materials have enabled the availability of modern 3D printers in a wide spectrum of precision, capability, and cost. Large-scale and specialized 3D printers are capable of printing highly complex geometries using materials such as polymers, metals, ceramics, concrete, and even clay with high accuracy. Some of those printers can mix different materials on demand to achieve the desired mechanical properties and aesthetics [23,32]. More economic desktop 3D printers are typically constrained to printing polymers that can print layers with thickness as low as 10 μm [33]. Contemporary 3D printing technology affords greater design freedom and production flexibility in both research and industrial applications [46].
A given 3D printing method offers certain advantages and drawbacks. This study uses the polyjet 3D printing technology because it provides relatively inexpensive and fast manufacturing of small, detailed parts. Usually, polyjet printers have two print heads that deposit different resins. The resins are liquid photopolymers and are hardened by an ultraviolet laser. One resin creates the desired object while the other resin acts as a support structure. Once the printing is complete, the support structure is removed from the finished 3D object by water jetting and chemical treatment using a 2% sodium hydroxide solution. The interparticle uniaxial compression test result shown in Fig. 1b was performed on two spheres of polyjet-printed polymer.
The polyjet polymer has a Young's Modulus of 2.4 GPa and a Poisson's ratio of 0.3, as shown in Table 2. As reported in the literature, both the printing direction and printing material can significantly affect the strength and surface properties of polyjet printed objects [11,44]. Specifically, the fracture stress and strain, and tensile toughness are affected by the printing orientation [11,19], which may affect the fracture and breakage of 3D printed objects. However, it should be noted that no significant breakage of particles was observed during any of the tests reported here. On the other hand, the contact deformation response of 3D printed particles is not expected to have an anisotropic behavior based on literature reporting that the printing orientation has no significant effect on the modulus of elasticity, and the ultimate tensile and compressive strengths [19,45,51].  [10], b [42], c [13], d [2], e this study

Natural and 3D printed sand particles
The natural sand particles used in this study were obtained by sieving a well-graded quartz sand. The particles that passed through the #6 (3.36 mm) and were retained by the #8 (2.38 mm) sieves were separated to obtain a poorlygraded soil. This soil consisted of both angular and rounded particles, which were then manually separated to create two separate soil samples: one with angular particles and one with rounded particles (Fig. 2). This approach ensured similarity in the gradation and mineralogy of the two samples, and allowed isolating particle shape as the only difference between them. The 3D printed particles were generated from X-ray CT scans of a set of randomly chosen angular and rounded natural sand particles (90 for the former, 70 for the latter). The particles were printed using an Object Eden 260 V printer from Stratasys with VeroWhitePlus rigid acrylate-based polymer resin with a printing resolution of 30 μm. The X-ray CT scans had a resolution of 10 μm, which were reduced to expedite the 3D printing process. The scans of reduced resolution were used to generate the 3D printed particles. A comparison of natural particle scans, scans with reduced resolution, and scans of 3D printed particles are presented in Fig. 3a-c, respectively. The morphology of a particle can be characterized at different scales: shape (medium scale) and surface texture or roughness (small scale) [7]. This study considers the particle shape to compare the similarity between the natural and 3D printed particles. A particle shape is typically quantified by roundness and sphericity, both of which can be described in a number of ways [7,16,31]. The shape parameters are usually estimated from 2D projection of a particle [7]. Roundness is a measure of the smoothness of the angles or corners of a particle, whereas sphericity is a measure of how closely the particle shape approaches that of a circle (in 2D) or sphere. This study considers Wadell roundness [49], area sphericity, and width-to-length ratio sphericity [31,60] to describe the particle shape. The shape parameters were obtained by image analysis of 2D projections of particles using the code by Zheng and Hryciw [60] and are shown in Fig. 4 and summarized in Table 2. The shape parameters of the 3D printed particles compare well with those of the natural particles, indicating the successful reproduction of particle shape.

Coefficient of friction test
Differences in the friction coefficient can lead to differences in the stiffness and strength behavior of granular assemblies [20,34,39]. The frictional resistance between polyjet 3D printed materials was measured by shearing two equal-sized blocks (63 mm length, 25 mm width, and 19 mm thickness) against each other in a modified direct shear test apparatus. Anisotropy in the frictional resistance due to the direction of material layer deposition was examined by shearing the blocks in the directions perpendicular and parallel to the printing direction (Fig. 5a). Every test used newly printed blocks. Five tests in each direction were performed under normal stresses (σ v ) ranging from 50 to 400 kPa.

Triaxial test
Drained and undrained isotropically-consolidated triaxial compression tests were performed to characterize the shear behavior and stress-dilatancy response of specimens of 3D printed and natural sands. This study used an automatic triaxial testing system with digital data acquisition capabilities. Cell and pore pressure and volume change were controlled using two digital pressure volume controllers. The measured volume changes are used to determine the specimen volumetric strain, ε v . Axial load was measured by an external load cell mounted on the load frame. Axial displacement was measured by an external linear variable differential transducer which is used to determine the specimen axial strain, ε a . Pore pressure transducers were used to measure the specimen pore pressure as well as the triaxial confining pressure. The changes in void ratio, e, of the specimens is computed based on the measured volumetric changes following the equation: Tests were performed on specimens of 70 mm diameter and 150 mm height. The dense specimens were prepared by pouring the particles in the mold in three lifts. After pouring each lift, a tamping rod was used to densify the layer to the target void ratio. The loose specimens were prepared by pouring the particles from a small height into the mold in a single lift without any tamping. Once prepared, a small vacuum was applied to stabilize the specimen. The specimen was then placed in the triaxial cell and the cell was filled with de-aired water. The specimen was saturated by applying back-pressure while maintaining a constant small difference between the cell and the back-pressure. The back-pressure was increased slowly until a B-value of 0.95 was obtained. After saturation, the specimen was consolidated isotropically to the target confining pressure; once the consolidation phase finished, the shearing phase commenced. All the specimens exhibited a bulging failure with no visible shear bands, with the exception of two tests on dense 3D printed rounded sand which exhibited shear bands (tests "6-3DPR-30" and "7-3DPR-30" in Table 3). While a membrane penetration correction was attempted for the volume changes, the effect was insignificant (less than 2%); thus, correction was ignored for all the tests. Table 3 summarizes the 34 triaxial compression tests performed. The testing ID convention is such that ID "x-AB-y" Fig. 3 Comparison of X-ray CT scans of a natural particles, b reduced scans for 3D printing, and c 3D printed particle analogs Fig. 4 Comparison of shape parameters for natural and 3D printed particles. Note that standard deviation in shape parameter values is shown by the error bars and the "after TX" parameters were obtained after performing over 15 triaxial tests on the 3D printed particles corresponds to test number x of particular sand AB at an effective confining pressure ( ′ 3c ) of y kPa. For example, "4-NR-530" corresponds to test number 4 on natural rounded sand at ′ 3c of 530 kPa. Given the greater stiffness of quartz as compared to the 3D printed polymer (Table 2), the tests on natural sand were performed at greater confining pressures than the tests on 3D printed sand. Initially, a series of five drained tests were performed on angular and rounded 3D printed sands ( ′ 3c range of 20-90 kPa) and natural sands ( ′ 3c range of 90-742 kPa) to characterize their triaxial compression behavior. Then, five triaxial tests were performed on angular 3D printed and natural sand specimens with similar initial void ratios with the objective of shearing them along specific drained or undrained stress paths to compare their response. Two similar tests were performed on rounded 3D printed and natural sand specimens subjected to undrained shearing. These tests on natural sands were performed at confining pressures that were about 22 times greater than those used for the 3D printed particles (i.e. 650 kPa for the natural sands and 30 kPa for the 3D printed sands). This ratio was chosen to be somewhat close to the ratio of the Young's moduli of quartz and polyjet 3D printing polymer while minimizing the possibility for particle breakage in the natural sands by maintaining the mean effective stress, p', to values smaller than about 2000 kPa throughout the triaxial tests. It is noted that both natural and 3D printed sands showed no visible particle crushing. Figure 4 presents particle shape parameters for the 3D printed sands obtained after being subjected to 15 triaxial compression tests under cell pressures varying from 20 to 90 kPa. As shown, no significant changes in roundness, area sphericity, and width-tolength ratio were observed, implying that damage to the 3D printed particles was not significant.

Results
This section first discusses the friction coefficient measurements of 3D printed polymer and compares it with published values from natural quartz sand particles to provide insight into differences in the inter-particle frictional interactions. Then, the results of drained triaxial tests performed on angular and rounded 3D printed and natural sands are presented to explore the stress-dilatancy behavior and make comparisons between the different materials. Tests along analogous drained or undrained stress paths are presented to compare the mechanical response of the 3D printed and natural sands, and to provide insight regarding the feasibility of modeling the behavior of coarse-grained soils with 3D printed analogs. Section 4 includes a comparison of the response of 3D printed and natural sands along analogous stress paths, along with an estimation of the critical state lines (CSLs) in e-log p' space.

Friction of polyjet 3D printing resin
The results of frictional resistance tests between equal-sized blocks printed with polyjet 3D printing polymer revealed dependence of the friction coefficient, μ, on the orientation of the printing layers and the magnitude of normal stress, σ v , as shown in Fig. 5b and c. The measured μ was significantly greater (2.9 times on average) for the direction Fig. 5 a Schematic of friction test used to determine printed material friction coefficient with different printed layer orientations, b friction coefficient of polyjet 3D printed blocks for different applied vertical stresses, and c relationship between shear stress and applied vertical stress parallel to layer deposition as compared to the μ values measured in the direction perpendicular to layer deposition. Also, μ decreased as σ v was increased in both shearing directions. This may be related to the plastic deformation of the micro-asperities on the polymer surface. Namely, as the σ v is increased the normal contact force increased, leading to yielding of micro-asperities as shown in Fig. 1b. As plastic strains accumulate at the contact, the micro-asperities flatten, likely resulting in a smaller surface roughness that leads to a decrease in interlocking interactions between the asperities. While in reality the failure envelope is curved, the data can be reasonably fitted with two separate linear envelopes for different σ v ranges where the slopes correspond to the average friction coefficient (Fig. 5c). As shown, the average μ for σ v < 200 kPa is 0.44 in the parallel direction and 0.15 in the perpendicular direction, whereas the average μ for σ v > 200 kPa is 0.35 in the parallel direction and 0.11 in the perpendicular direction. Average μ values for Leighton Buzzard and ASTM 20-30 sand particles range between 0.17 and 0.36 (Table 1). While not reported in the literature, these values are likely not direction-dependent. Thus, the results of the block friction tests reported here Table 3 Summary of triaxial testing program for natural and 3D printed particles ′ 3c , effective confining pressure; e 0 , initial void ratio; e c , void ratio after consolidation; e f , void ratio at the end of test; q peak and q end , peak and end of test deviatoric stress, respectively; p ′ peak and p ′ end , peak and end of test mean effective stress, respectively; NA, natural angular; NR, natural rounded; 3DPA, 3D printed angular; 3DPR, 3D printed rounded; CD, consolidated drained; CU, consolidated undrained suggest that the μ in the parallel direction for 3D printed polymer are greater than those expected between natural sand particles; however, the μ in the perpendicular for 3D printed polymer are smaller than those expected between natural sand particles.

Triaxial compression behavior
Five drained triaxial compression tests were performed on specimens of natural and 3D printed angular and rounded particles to characterize their consolidation and shearing behavior. Figure 6a-d show the isotropic consolidation curves of natural angular, natural rounded, 3D printed angular, and 3D printed rounded particles, respectively. As shown, the decrease in void ratio with increasing mean effective stress is greater for the angular and rounded 3D printed particles (Fig. 6c, d) than for the corresponding natural particles (Fig. 6a, b), highlighting the greater skeleton compressibility of the 3D printed soils. The average slopes of the consolidation curves for natural angular and rounded particles are 0.032 and 0.036, respectively, while the corresponding values for 3D printed angular and rounded particles are 0.116 and 0.138, respectively. The greater compressibility of 3D printed sand is likely due to the initial plastic yielding of micro-asperities of the rougher inter-particle contacts, as shown in Fig. 1b and discussed by Ahmed and Martinez [2]. The results presented herein also highlight that the natural and 3D printed rounded particles exhibit a slightly greater compressibility than the natural and 3D printed angular particles, respectively. The consolidation results from the natural sands indicate that the curves that start at different void ratios remain roughly parallel to each other in the range of mean effective stress considered. This is particularly evident in Fig. 6b and is in agreement with Altuhafi and Coop [3] and Yamamuro et al. [57] for mean effective stresses at which no significant particle crushing takes place. In contrast, the consolidation results from the 3D printed soils indicate that the curves tend to converge; this is particularly evident in Fig. 6c. This difference in behavior between the natural and 3D printed sands could be explained by the limiting compression curve (LLC) concept. Namely, the compression curves for natural cohesionless soils with different initial densities tend to converge to a unique curve at high mean effective stress levels, referred to as the LCC [14,36]. At low stress levels, the volume changes of soils are due to elastic compression of the soil skeleton and particle rearrangement, while the LCC response is controlled by particle damage (i.e. asperity breakage and particle crushing) [37]. The damage of particles is affected by particle shape, gradation and mineralogy. In this study, no crushing was observed in the specimens of natural sand and 3D printed sand. However, damage of the 3D printed particles during consolidation likely took place in the form of plastic yielding of the micro-asperities, which Fig. 6 Isotropic consolidation curves for a natural angular sand, b natural rounded sand, c 3D printed angular sand, and d 3D printed rounded sand could be analogous to the breakage of asperities in natural sands. For the 3D printed sands, the results suggest that a mean effective stress of about 50 kPa causes this convergence of the compression curves.
The deviatoric stress and volumetric change responses for the natural and 3D printed sands during drained shearing correspond to those expected for coarse-grained soils. As shown in Fig. 7 for the angular sands, greater deviatoric stresses were developed at higher confining pressure for both natural (Fig. 7a) and 3D printed (Fig. 7e) sands. The q-ε a curves for the natural angular sand tests performed at confining pressures of 90, 200, and 318 kPa show a slight Fig. 7 Drained triaxial test results on a-d natural angular sand and e-h 3D printed angular sand peak followed by strain softening (Fig. 7a) accompanied by dilative volumetric strains (Fig. 7b). On the other hand, the q-ε a curves for the tests at confining pressures of 530 and 742 kPa exhibit strain hardening along with overall contractive volumetric strains (Fig. 7a, b). The q-ε a curves for the tests on angular 3D printed sand exhibit strain hardening and an overall contractive behavior (Fig. 7e, f). However, the specimens exhibited slight dilation at axial strains greater than about 15%.
The points at the end of shearing for the natural and 3D printed sand specimens were taken as the critical state points to estimate the critical state line in the q-p' plane ( Fig. 7c,  g). The tests on natural angular sand yielded CS points in the q-p' plane that can be fitted with a straight line passing through the origin with a slope, M, of 1.35, corresponding to a critical state friction angle, ′ cs , of 33.4°. In contrast, the CS points of the tests on 3D printed sand cannot be fitted with a straight line in the q-p' plane. While the CSL may be curved, two straight lines are used here to define average friction angles for two ′ 3c ranges. A straight line fitted to the tests with a confining pressure smaller than 50 kPa has a M of 0.936 ( ′ cs of 23.9°), while a line fitted to the tests with confining pressures greater than 50 kPa has a slope of 0.792 ( ′ cs of 20.6°). The reason for the decrease in M and ′ cs with increasing ′ 3c could be due to the decrease in inter-particle friction coefficient of the 3D printed polymer as the normal stress is increased, as shown in Fig. 5b and c. The results also indicate that the M and ′ cs values obtained for the 3D printed angular sand are smaller than those for natural angular sand. One possible reason for this is the smaller interparticle friction coefficient of the 3D printed polymer in the direction perpendicular to layer deposition (between 0.12 and 0.19) compared to that of natural sands (between 0.17 and 0.36) ( Table 1). Another reason could be plastic deformation of the particles' asperities, which would decrease the interlocking between particles. However, further research is required to determine the reason for the differences in M and ′ cs values. Note that, post-test assessment of particle shapes shows no statistically significant changes in particle shape parameters (Fig. 4). The stress paths in the e-log p' plane in Fig. 7d and h show the dilation of the natural sand observed during tests at confining pressures of 90, 200, and 318 kPa and the contraction observed during tests at greater confining pressures. In contract, the stress paths show the overall contraction during all the tests on 3D printed sands. Due to this continued contraction, the 3D printed sand specimens reach lower void ratios (0.581-0.683) than the natural sand specimens (0.627-0.750). It is noted that the difference in final void ratios is likely due to a combination of the greater compressibility of the 3D printed skeleton as well as the differences in volume change tendencies dictated by the specific void ratio and mean effective stress of each specimen.
Similar to the angular natural and 3D printed sands, the drained shear behavior of rounded natural and 3D printed sands agrees with that expected for coarse-grained soils (Fig. 8). The q-ε a plots indicate that both natural and 3D printed particles develop higher shear stress at higher confining pressure (Fig. 8a, e). All the q-ε a curves for the natural rounded sand exhibit a peak followed by slight strain softening (Fig. 8a) accompanied by dilative volumetric strains (Fig. 8b). The tests on 3D printed rounded sand exhibit strain hardening (Fig. 8e) accompanied by an overall contractive behavior. However, in a similar way as the 3D printed angular sand, slight dilative volumetric strains take place at axial strains greater than about 15% (Fig. 8f). The trends shown in the q-p' plane ( Fig. 8c, g) are similar to those described for the angular sand. Namely, a straight line passing through the origin can be fitted for the natural rounded sand in the q-p' plane with an M of 1.29 ( ′ cs = 32.1°), which is slightly smaller than the M for the natural angular sand (M = 1.35, ′ cs = 33.4°). This difference reflects the greater roundness and sphericity of the rounded particles (Fig. 4). The CS points from the tests on 3D printed rounded sand can be fitted with two different straight lines in the q-p' plane (Fig. 8g). The fitted M slope for the tests at ′ 3c < 50 kPa is 0.887 ( ′ cs = 22.7°) while the slope for the tests at ′ 3c > 50 kPa is 0.742 ( ′ cs = 19.3°). Both M values obtained for 3D printed rounded particles are lower than those obtained for 3D printed angular particles due to the greater roundness and sphericity of the former. The stress paths in the e-log p' plane shown in Fig. 8d and h indicate similar trends as described for the angular sands. Namely, they show net dilation for the natural sand and contraction for the 3D printed sand as well as smaller void ratios for the 3D printed sand at the end of the tests.
The drained triaxial compression results indicate that the 3D printed sands exhibit stress-dilatancy behavior in agreement with that of natural coarse-grained soils. The stress-dilatancy (R-D) relationships for the tests on the natural and 3D printed angular and rounded sands are shown in Fig. 9, where R = � 1 ∕ � 3 is the stress ratio, ′ 1 is the major principal effective stress and ′ 3 is the minor principal effective stress, and D = 1 − dε v /dε a is the dilatancy. The test data is compared to Rowe's flow rule [38], which is expressed as: Figure 9a and b show that the stress-dilatancy data for natural angular and rounded sands, respectively, closely

Comparison of the behavior of 3D printed and natural sands along analogous stress paths
Triaxial test pairs were performed on natural and 3D printed sands at confining pressures of 650 and 30 kPa, respectively. The tests pairs were performed with void ratios at the beginning of shearing (e c ) that were within 0.02 from each other. This section discusses the similarities in triaxial response for drained and undrained tests on the angular sands and undrained test pairs on the rounded sands.

Drained behavior of natural and 3D printed angular sands
The effect of the void ratio at the beginning of shearing on the natural and 3D printed angular particles was investigated by conducting test pairs on specimens with larger and smaller void ratios. The results of drained tests on natural angular and 3D printed angular sands with e c values of 0.705 ± 0.01 and 0.645 ± 0.01 are shown in Fig. 10a-h. As shown in the q-ε a curves (Fig. 10a, e), the specimens with smaller e c yielded greater deviatoric stresses than the specimens with greater e c for both natural and 3D printed sands. The natural angular sand specimen with smaller e c exhibited an initially contractive response followed by strong dilation (Fig. 10b); the continued dilation at the end of the test indicates that the specimen did not reach a constant-volume state (i.e. critical state). The specimen with the greater e c exhibited overall contractive volumetric strains throughout the entire test. The 3D printed angular sand specimens with both larger and smaller e c exhibit an initial contractile response up to an axial strain of about 15% followed by some dilation (Fig. 10f). The specimen with smaller e c exhibited a smaller amount of initial contraction followed by greater rate of dilation, as compared to the specimen with the larger e c . The initial contraction observed in all the tests on 3D printed sand specimens may be indicative of skeleton compression due to the increase in p' during the initial 15% of axial strains. This suggests an agreement with the consolidation results shown in Fig. 6c and d and with 1D compression results presented by Ahmed and Martinez [2], which indicate greater compressibility of the 3D printed sand skeleton in comparison with the natural sand skeleton. In the q-p' plane, the stress paths for the natural angular sand specimens appear to evolve towards the CSL (Fig. 10c) even though it appears that neither test reached critical state. The stress paths of the tests on 3D printed angular sand also appear to evolve towards the CSL (Fig. 10g), with the test with an initially greater e c converging to the CSL and the test with the smaller e c still above the CSL due to the larger rate of dilation. The stress paths in the e-log p' plane are qualitatively similar between the tests on natural and 3D printed natural sands. Namely, both tests with larger e c show net contraction and both tests with smaller e c show initial contraction follow by dilation (Fig. 10d, h).

Undrained behavior of natural and 3D printed angular and rounded sands
The undrained triaxial compression behavior of natural and 3D printed angular particles with different void ratios at the beginning of shearing was investigated by performing test pairs at e c of 0.705 ± 0.01, 0.650 ± 0.01, and 0.590 ± 0.01 (Fig. 11). The q-ε a curves (Fig. 11a, e) indicate that the specimens of both natural and 3D printed particles with an e c of 0.590 ± 0.01 mobilized the greatest deviatoric stresses while the specimens with an e c of 0.705 ± 0.01 mobilized the smallest deviatoric stresses. The natural angular sand specimen with the highest e c exhibited positive Δu that increased to about 310 kPa and then decreased to about 210 kPa, while the specimen with an e c of 0.649 exhibit initial positive Δu and decreased to slightly negative values by the end of shearing (Fig. 11b). The densest natural sand specimen generated negative Δu that reached a value of about -330 kPa at the end of shearing. The trends exhibited by the angular 3D printed sand specimens follow a similar trend where the specimen with the greatest e c (0.714) generated the greatest magnitude of Δu, followed by the specimen with an e c of 0.654 and then by the specimen with the lowest e c of 0.591 (Fig. 11f). The Δu of the 3D printed Fig. 11 Undrained behavior of a-d natural angular sand and e-h 3D printed angular sand sand specimens during the tests increased rapidly at axial strains smaller than 5%, after which the magnitudes tended to decrease slowly throughout the end of the tests. These trends are in agreement with the rapid rate of contraction observed in the drained tests in Figs. 7f and 10f at small axial strains followed by the slight dilation at greater axial strains.
In the q-p' plane, the end-of-test points for the natural angular sand specimens converged towards the CSL (Fig. 11c). The stress paths of the angular 3D printed specimens evolve towards the CSL; however, the points are above Fig. 12 Undrained behavior of a-d natural rounded sand and e-h 3D printed rounded sand the CSL at the end of the test, indicating that critical state has not been reached (Fig. 11g). In the e-log p' plane, the end-of-test points for the natural particles show the increase in p' which is greatest for the specimen with the smallest e c (Fig. 11h). The stress paths of the 3D printed sand specimens follow similar paths; however, the increase in mean stress is considerably smaller than for the tests on natural angular sand.
The results of the undrained triaxial tests on natural and 3D printed rounded sands show similar trends as described for the angular sands (Fig. 12). Namely, the specimens with smaller e c mobilize greater deviatoric stresses (Fig. 12a, e) due to the greater magnitude of negative excess pore pressures (Fig. 12b, f) than the specimens with greater e c . All the q-ε a curves exhibit a hardening response, which is characteristic of dilative undrained sand. The tendency of the 3D printed sands to contract at axial strains smaller than about 5-10% is evident in the results; however, this tendency is stronger for the specimen with the greater e c . In the q-p' plane, the end-of-test points for the natural sand specimens converge toward the CSL (Fig. 12c) whereas they are above the CSL for the 3D printed sands, indicating that critical state has not been reached (Fig. 12g). In the e-log p' plane, all the stress paths of the tests on both natural and 3D printed sands show a net increase in the p' (Fig. 12d, h).
The change in mean effective stress with axial strain during undrained shearing is typically associated with the volumetric deformations of the specimen. Constant p' during undrained shearing indicates volumetric deformation of the soil skeleton itself while skeleton shear deformations in the form of dilation or contraction due to particle rearrangement have not yet taken place [55]. Figure 13 shows the p'-ε a plots of the undrained tests on natural and 3D printed sands. As shown, p' for both angular and rounded natural sands (Fig. 13a, b) begins to change at very small ε a (smaller than 0.5%), indicating that skeleton shear deformation starts at a very small ε a . However, p' remains relatively constant for ε a up to about 3% for the 3D printed angular sand (Fig. 13c) and about 7% for the loose 3D printed rounded sand (Fig. 13d), indicative of the initial volumetric contraction of the skeleton observed up to a larger ε a compared to the natural sands. This suggests that skeletal shear deformation for the 3D printed sands initiates after greater initial skeleton volumetric contraction and at a greater ε a compared to the natural sands. This observation for the 3D printed sands is in agreement with the initial contraction up to a larger ε a observed in drained tests (Figs. 7f, 8f, 10f).

Estimation of critical state lines in e-log p' space
The critical state lines in e-log p' space for both the natural and 3D printed sands are approximated by best fitting the critical state points obtained from both drained and undrained tests. The CS points for CD and CU tests were Fig. 13 Mean effective stress evolution during undrained triaxial tests on a natural angular sand, b natural rounded sand, c 3D printed angular sand, and d 3D printed rounded sand obtained by extrapolating the end-of-test results following methods described in Zhang et al. [59] and Torres-Cruz and Santamarina [47], respectively. Examples of the extrapolation procedures are presented in "Appendix A", which consist of extrapolating the dilatancy to a value of 1.0 in drained tests and the rate of pore pressure change to a value of zero in undrained tests.
The end-of-test points in the e-log p' plane for natural angular and rounded sands are shown in Fig. 14a and b, separated into specimens that were deemed to reach critical state and those that did not. A curved form of the critical state line (CSL) is considered here to better capture the shape across a wide range of p' values. Hence, the CSL is presented by a power function [52] as: where e f is the critical state void ratio, e Γ is the critical state void ratio at p' = 0 kPa, p a is the atmospheric pressure (≈100 kPa), and λ and ξ are material constants describing the material compressibility and the non-linearity of the CSL, respectively. The CSL equations obtained for both natural angular and rounded sands are provided in Fig. 14a and b, respectively, showing a greater e Γ intercept for the angular sand, a slightly greater λ for the rounded sand, and similar ξ exponents. Similar to the natural sand, the CS points in e-log p' plane obtained for the 3D printed sand can be approximated by Eq. 3 as shown in Fig. 14c and d for the angular and rounded sands, respectively. Comparison of the equations for angular and rounded 3D printed sands indicates similar trends as the natural sand. Namely, the e Γ intercept is greater for the angular sand and λ is slightly greater for the rounded sand. However, the ξ exponent is slightly smaller for the rounded sand. Comparison of the CSL equations for the natural and 3D printed sands reveals important trends. The e Γ intercepts are relatively close in value albeit slightly greater for the natural sands: 0.785-0.746 for natural sands and 0.745-0.725 for 3D printed sands. The λ values, which are related to the material compressibility, are considerably greater for the 3D printed sands: 0.0215-0.024 for natural sands and 0.148-0.165 for 3D printed sands. These differences agree with the greater compressibility during isotropic and 1D compression shown in Fig. 6a-d and by Ahmed and Martinez [2]. Finally, the ξ exponents were smaller for the 3D printed sands: 0.78 for natural sands and 0.70-0.65 for 3D printed sands. This indicates a greater non-linearity of the CSL of the 3D printed soils. The estimation of the CSLs for the 3D printed soils suggests that their behavior can be captured within the critical state framework in the same way as for natural soils. It can be envisioned that the behavior of natural soils could be modeled with 3D printed analogs using critical state concepts. For example, tests could be performed at a combination of void ratio and initial mean effective stress such that the state at the beginning of shearing for the natural and 3D printed sand specimens have the same state parameter with respect to their corresponding CSLs (as described by Been and Jefferies [6]).

Discussion on the modeling of soil behavior with 3D printed particle analogs
One of the greatest potential benefits offered by the 3D printed sands is the ability to systematically control different particle properties, such as the shape, size, and constituent material. Additionally, 3D printed particles may also enhance validation procedures for discrete element modeling simulations against experimental data by ensuring the use of particles with similar morphology in both experimental and numerical investigations [26]. The results presented here indicate that polyjet 3D printing technology can be used to successfully reproduce the shape of natural sand particles (Figs. 3, 4). It can thus be envisioned how synthetic particles could be generated with methods such as spherical harmonics [53] and 3D printing to systematically investigate differences in strength and stiffness of granular assemblies due to changes in particle shape or particle size distribution alone. The triaxial compression tests indicate that the 3D printed sands exhibit many of the fundamental behaviors that characterize sands and gravels, suggesting that they can be used to model the behavior of natural soils. Namely, the stress-dilatancy behavior conformed to established flow rules (Fig. 9 There are important differences in the response of 3D printed sands with respect to that of natural sands. The polyjet 3D printed particles employed in this study exhibit greater plastic deformation at inter-particle contacts compared to that observed between glass or sand particles due to the greater surface roughness of the particles resulting from the layer deposition process (Fig. 1). In addition, greater elastic deformation is also experienced at the contact due to the lower Young's modulus of the 3D printed polymer. These are likely the reasons for the greater compressibility of the 3D printed sands (Fig. 6,  14), along with the greater skeleton compressive deformation observed at small axial strains during triaxial testing (Figs. 7, 8, 13). The layer deposition process also results in direction-dependency of the friction coefficient of the 3D printed material (Fig. 5), which is likely not the case for natural soil particles. The smaller friction coefficient in the direction perpendicular to layer deposition may be responsible for the smaller friction angles of the 3D printed soils (Figs. 7 8). Understanding of these differences is required for assessing how closely the macro-scale behavior of natural soils can be modeled or reproduced by the 3D printed soil analogs.
The micro-scale behavior of 3D printed particles is determined by the surface created by the printing process and by the mechanical properties of the printed material. Multiple types of 3D printing technology exist (e.g. stereolithography, fused deposition modeling, selective laser sintering), and new technologies are being developed rapidly. Because each technology is capable of printing different materials and each one employs a different manufacturing process, the possible effects on the response of individual particles (i.e. contact normal stiffness, friction coefficient) and granular assemblies should be evaluated and understood if the behavior of soils and other granular materials is to be modeled in a quantitative manner.

Conclusions
This paper investigates the feasibility of using 3D printing technology to generate analog particles to model the triaxial compression behavior of coarse-grained soils. A total of 34 drained and undrained triaxial compression tests on specimens of natural and 3D printed angular and rounded particles were performed. The main findings of this study are summarized as follows: • The polyjet 3D printing technology can accurately reproduce the shape and size of natural coarse sand particles. However, the surface texture of the polyjet 3D printed particles is dependent on the printing layer direction, which results in different surface roughness that affects the inter-particle frictional coefficient. The friction coefficient of the polyjet polymer was also observed to be dependent on the magnitude of applied normal stress. • The 3D printed sands are more compressible compared to the natural sands, which is likely due to lower Young's Modulus of the polymer and to plastic yielding of microasperities on the surface of the 3D printed material. During shearing, this greater compressibility likely causes larger initial volumetric contraction (or positive excess pore pressure generation) followed by skeletal shear deformation that begins at a greater axial strain compared to that of natural sand. • The drained test results show that the 3D printed sands follow Rowe's flow rule, demonstrating that they can replicate the stress-dilatancy behavior observed in natural sands. Also, the angular 3D printed sand mobilizes greater critical state friction angle than that of rounded 3D printed sand, in agreement with the results on the natural sands. • The results from test pairs performed at similar void ratios but different confining pressures (30 kPa for 3D printed sand, 650 kPa for the natural sand) indicate that analogous drained and undrained stress paths are followed by the test pairs in both the q-p' and e-log p' planes. This included contractive volumetric changes and generation of positive excess pore pressures in loose specimens, and dilative volumetric changes and generation of negative excess pore pressures in dense specimens. • The critical state line in the q-p' plane is curved for the polyjet 3D printed sands, with a slope that decreases as p' is increased. The critical state line in the e-log p' space can be described with a power-law function. Owing to the greater compressibility of the 3D printed sands, their critical state void ratios are smaller than those for the natural sands.