The excavation in layered rocks is an issue for a number of geoengineering applications; these kinds of rocks all exhibit transverse isotropic features due to the process of metamorphic differentiation. This paper focuses on providing two methods, i.e., the 3D printing technology and the combined finite-discrete element method, to simulate the anisotropic characteristics of layered rocks. The results showed that both the 3D-printed samples and the FDEM numerical models are considered as a good match, and both revealed that as the inclined angle increased, the UCS of the sample first decreased and then increased, showing a U-shaped pattern. The results of this paper served as a reference to the promotion of the 3D printing technology and the combined finite-discrete element method in the geotechnical engineering field and laboratory test research.

1. Introduction

Shale formations widely exist in nature. These kinds of rocks, including schists, sandstones, shales, and basalts, all exhibit features of transverse isotropic, which means that these materials all show the same phenomenon along a certain plane. During the diagenetic process of compaction with mineral particles, many joints, fractures, and tectonic structures will be formed; these structures have great influences on the strength and failure characteristics of shale formations and, hence showing great anisotropy [1, 2]. To study the mechanical behavior of shale formations is important for a number of geoengineering applications, including the drilling industry and the excavation of tunnels and the underground cavities, as well as the shale gas mining [3].

Laboratory experiment has been widely used to study the anisotropy of layered rocks, e.g., Attewell and Sandford [4] on shale; Brown et al. [5] on slate; Nassri et al. [6] on schist; and Al-Harthi [7] on sand stones. These studies all show a similar phenomenon, that is, the maximum strength can be observed at angles α = 0° or α = 90°, whereas the minimum strength can be observed at the angle of α between 20° and 30°. Here, α is introduced as the angle between bedding plane and the loading direction applied on the sample. However, the size of the laboratory test sample is small (φ50 mm × 100 mm), and it is difficult to obtain such a small size of natural rock sample due to the inevitable damage induced by the low strength characteristic of soft rock and the high-stress release disturbance [8]. In addition, the heterogeneous characteristic of rock makes it extremely difficult to obtain repeatable experimental samples. To avoid the difficulty of sampling natural layered rocks, some researchers have used rock-like samples produced by portland cement, quartz sand, gypsum, and water to investigate mechanical behaviors of layered rocks [9], but experimental accuracy control is still an inevitable problem.

The 3D printing technology, enjoying rapid development in recent years, can be used to produce various 3D entities with complex structures and unique shapes. The repeated prints in the same batch have high consistency in accuracy, detail, size, performance, etc. A new way for the physical simulation of rock test has been created by this technology since no two samples are identical in nature, and it can provide authoritative and effective comparative verification for the parameter calibration of numerical simulations.

Appearing in the mid-1990s, 3D printing technology has received close attention from all walks of life. This technology has been widely used in many fields, including biomedicine [10], aerospace [11], automobile manufacturing [12], and electronic components [13], and preliminary exploration has been launched in the field of geotechnical engineering for the past few years. 3D printing can be divided into fused deposition molding (FDM), stereolithography appearance (SLA), and powder-ink binders (PIB) in terms of the formation methods [14]. Different materials are used in various printing technologies. PLA plastic is the material used in FDM printing technology. The samples produced by it have obvious elastic-plastic characteristics which are not suitable for simulating rocks [14, 15]. Photosensitive resin was used in the SLA technology. Zhou and Zhu [16], and Zhu et al. [17] compared and analyzed several different 3D printing technologies. They found that the samples printed by SLA technology had better brittleness and were best suited to simulate hard rock. Gypsum powder was used in PIB technology, which was shaped by binder. Research has shown that its mechanical properties are more similar to those of sandstone [18, 19]. However, there are still two major problems. Firstly, with low intensity, the strength of the printed sample is even lower than that of the weakest rock in nature; secondly, with strong ductility, the stress-strain curve has no obvious elastic brittleness characteristics.

Therefore, when applying the PIB 3D printing technology to the study of rock properties, the first work needed is to increase its strength and brittleness as much as possible. Through changing the thickness of each layer of powder and the saturation of binder, Vaezi and Chua [20] studied the surface accuracy and tensile strength of the samples. Research carried out by Lu et al. [21] showed that the samples printed with finer powder particle size had higher strength and surface accuracy. The conclusion that the interval of each layer of powder laid in the printing process had certain impacts on the strength of the sample can be seen in the study of Farzadi et al. [22].

Besides laboratory testing, numerical simulation method has also been widely utilized to investigate the failure mechanism of rocks [23, 24]. These numerical approaches are typically classified either as continuum- or discontinuum-based methods. Continuum-based methods include nonlinear finite element method (FEM), Lagrangian finite difference method (FDM), and boundary element method (BEM). Discontinuum-based methods include universal distinct element code (UDEC), particle flow code (PFC), and discontinuous deformation analysis (DDA) method [25]. British scholars, Munjiza et al. [26, 27], combined the advantages of finite element method and discrete element method and proposed a new way to simulate multiple interacting deformable solids, which is called the combined finite-discrete element method (FDEM). Theoretically, FDEM can be a better way to simulate the gradual change process of materials from continuous to discontinuous under load conditions, which has been widely used in the field of rock mechanics in recent years.

In this paper, laboratory tests and numerical tests were performed using the 3D printing technology and the combined finite-discrete element method. By combining these two methods, the anisotropic mechanical behaviors of natural shale rock formations were studied. The structure of this paper is as follows: In Section 2, we introduce the equipment which is needed in the laboratory test. In Section 3, an optimized printing scheme is proposed to improve the strength and brittleness of the sample by enhancing the binder saturation level. On this basis, by changing the direction of 3D printing, the anisotropic characteristics of natural shale rock formations are simulated. In Section 4, the combined finite-discrete element method is used to simulate the mechanical characteristics of shale formations under uniaxial compression test, and the 3D printing test results are used to calibrate the required parameters for further research.

2. Equipment and Methodology

2.1. 3DP Equipment and Methodology

Produced by the 3D system company, the equipment used in the test is the ProJet 860 printer (Figure 1), which is based on a digital model and uses VisiJet PXL gypsum powder to construct the model. The forming process can be divided into the following steps [28, 29].

The first step is to build the model by 3D modeling software (ProE 5.0) and save it as a STL format; after reading the established STL format files by the printer, the 3D model will be decomposed into cross-sections with an interval of 0.1 mm along the z-axis. Then the printer will determine the binder spraying area of each cross section. After warming up the 3D printer, a certain thickness of powder will be prepaved on the build bed to provide a smoother foundation for sample printing. The specific molding process is shown in Figure 2.

Firstly, the roller is rolled counterclockwise to lay an extremely thin layer (0.1 mm thickness) of powder on the build bed. The excessive powder will be collected by the overflow tank for reusing. Then the build bed descends 0.1 mm and gets ready for the next layer. After this, the print head will spray binder on the corresponding area of each layer. The entire volume of the sample will be constructed by cycling these three steps.

When the sample printing is finished, the printer runs a drying circle for about 2 hours. Once the sample dries, it needs to be cleaned by removing the excess powder. After that, the strength of the sample can be further enhanced by submerging it into the infiltrant (ColorBond™) which is produced by 3D Systems company. However, it was hard to infiltrate fully, so the infiltrant was not used in the research.

2.2. RMT Equipment and Methodology

The RMT-301 rock and concrete mechanical testing system (Figure 3), developed by Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, was used as the loading equipment. The instrument contains two sets of force sensors with different measuring ranges. A small force sensor with a measuring range of 100 kN was used during the test due to the low strength of the 3D-printed sample. And four displacement sensors with a measuring range of 5 mm were used to test both the vertical and horizontal deformation. During testing, the loading mode was initiated through displacement control with a ratio of 0.002 mm/s.

3. Anisotropy in 3D-Printed Samples

3.1. Effect of Binder Saturation Level on the Strength of 3D-Printed Samples

To detect the strength of printed sample, we performed a UCS test which is widely used in rock mechanics. The samples were several cylinders with a diameter of 50 mm and a height of 100 mm. Before testing, the cylinders were naturally dried in a room with temperature of 28°C and air humidity of 60%. This process lasted for two weeks until the weight of the samples became stable.

After the uniaxial compression test, the failure pattern of the specimens printed by the default printing method of the system can be seen in Figure 4. As shown in the picture, it was found that the skin of the specimens was peeled off and the central part was in the form of incomplete bond powder. This was caused by the way of spraying the binder. In Figure 5(a), it can be seen that when spraying the binder on each layer of a cross section, the default spraying method is to spray the shell more densely, while the center part is only partially sprayed. In fact, the interior of the printed specimens is not completely bonded, even though its surface seems very firm. Hence, the default printing mode of the system is not applicable to the simulation test of rock. What is needed for the test is a more homogeneous sample. Therefore, the amount of binder inside the printed sample must be increased.

The 3D printing software of the system can adjust the amount of binder inside the sample by changing the parameter called “binder saturation level.” 100% of binder saturation level is the default parameter of the system, which represents the internal binder volume of the sample accounting for 11.9% of the total volume. In addition, we recorded the situations when using 125%, 150%, 175%, and 200% of the binder saturation. 15 samples were printed to test each case thrice. When the binder saturation level increased to 200%, the internal binder volume of the sample accounted for 23.8% of its total volume. Visually, the binder has been spread evenly, as shown in Figure 5(b).

When the binder saturation level increased from 100% to 200%, the UCS of the samples increased obviously and the ductility decreased according to the uniaxial compression test results (Figure 6). The UCS was 5.13, 6.11, 7.53, 9.62, and 9.99 MPa, respectively, which exhibited an increasing trend. When the binder saturation level was 175%, the UCS of the sample was close to the maximum value. It was nearly double the sample printed with the default printing method and was only 3.7% lower than the strength of the sample printed with 200% binder saturation. The strain-stress curve of the 3D-printed samples under the uniaxial compression test was similar to those of rock-like materials, as they underwent the stages of elastic deformation, crack development, strain softening, and residual deformation [30]. Considering the cost of the binder ($350 per liter), 175% of the saturation level will be used for further research. The UCS of 3D-printed sample was close to 10 MPa at that time, which has been improved compared with 6.7 MPa in the study by Jiang et al. [31] and 7.93 MPa in the study by Jiang Quan et al. [14].

3.2. Effect of Printing Direction on the Strength of 3D-Printed Samples

Previous studies have shown that the layering orientation has great influence on the strength of shale samples, since the natural shale is laminated [32]. At the same time, PIB 3D-printed samples are also stacked layer by layer. During the process of printing, the printing rate was approximately 2–4 layers/min; therefore, the sprayed binder took 15–30s to be concreted. As a result, the obvious directionality can be seen in the final printed sample. This forming process is similar to that of natural shale formations. By changing the printing direction of the 3D-printed sample, the physical and mechanical properties of natural shale formations can be well simulated.

Figure 7 shows the definition of inclined angle α, which is introduced as the angle between bedding plane and the loading direction. In order to study the influence of printing direction on the anisotropy characteristics of 3D-printed specimens, the inclined angle α was divided into four equal parts on average from 0° to 90°, and five different α angles, namely, 0°, 22.5°, 45°, 67.5°, and 90°, were tested. According to the research results in section 3.1, 175% of binder saturation was selected in this part, and after printing finished, the samples were naturally dried at room temperature of 20°C for two weeks. To reduce the uncertainty of the obtained results, two samples were printed in each direction.

Figure 8 shows the stress-strain curve of 3D-printed samples. As shown, the compressive deformation properties of the 3D-printed samples were similar to those of rocks when α = 0° and α = 22.5°, and more ductile deformation was shown at other angles according to the tests. The maximum failure strength of 3D-printed sample is at α = 0°, and the minimum value is at α = 22.5°; this type of anisotropy is observed in the research of layered rocks [6].

4. Numerical Simulation

4.1. Model Description

The combined finite-discrete element method (FDEM), proposed by A. Munjiza et al. in 1995 [26], has been widely used in geotechnical engineering. In order to study the mechanical properties of layered rocks, A. Lisjak et al. proposed two different modeling methods: discrete approach [33] and smeared approach [34]. Compare with the samples simulated by the discrete approach, the failure patterns of the samples simulated by smeared approach are closer to that of the natural layered rocks which will serve as the theoretical basis for the numerical simulation in this paper.

For the purpose of comparison with the laboratory test, the inclined angles α were set as 0°, 22.5°, 45°, 67.5°, and 90°, respectively. The sample at α = 67.5° was taken as an example (Figure 9). The calculation model was 50  100 mm rectangular specimen according to the real test size. Uniaxial loading platens at both ends of the sample were moving in opposite directions with a constant velocity of 0.05 m/s. Although the loading speed is much higher than 0.002 m/s in the actual laboratory test, the research carried out by Lisjak et al. [34] shows that as long as the loading rate is less than 0.25 m/s, the simulation results tend to be constant.

In order to capture the delamination process along bedding planes, the model was divided into several blocks, as depicted in Figure 9(a). Moreover, to reduce the constraint imposed by the mesh topology on the model behavior, each small block was meshed with an unstructured Delaunay triangulation scheme by GMESH software. Taking an enlarged view of the model, as Figure 9(b) shows, the resultant mesh is a combination of a random triangulation for the intralayer material together with edges preferably aligned along the layering direction for the bedding planes. The definition of t is the layer thickness, and this value can only affect the potential least-energy failure path but does not seem to interrelate with the values of UCS when there are 8–10 meshes in one diameter region [34, 35]. Therefore, the layer thickness was kept constant (t = 10 mm) in the simulated test.

The definition of θ is the angle between crack element and the bedding plane, as shown in Figure 9(c). It is assumed that the values of the cohesive strength parameters fθ (GI, GII, c, ft) conform to a linear variation law according to the bedding plane [34]. The maximum and minimum values of these parameters are at θ = 90° and θ = 0°, respectively; when θ is between 0° and 90°, the linear relationship conforms to the formula as follows:where f and f90° are the values of the cohesive strength parameters at θ = 0° and θ = 90°, respectively. The friction angle, ϕ, is kept constant. In addition, five independent elastic constants are required to describe the elastic deformation of the triangular elements: two orthogonal elastic moduli Ex and Ey; two Poisson’s ratios and ; and the shear modulus G'.

4.2. Mechanical Test Modeling and Parameter Calibration of Layered Materials in FDEM

In FDEM, the requirement of simulation accuracy and calculation amount should be considered comprehensively to select a reasonable mesh size. The study of Liu and Deng [35] showed that the mesh size of the uniaxial compression standard sample with a diameter of 50 mm should not be greater than 1.8 mm, so in this study, the specimens were meshed with the grid of 1 mm average element edge length. Correspondingly, the time step size was set as 3 × 10−9 s.

The basic input parameters are based on the laboratory tests. However, similar to the traditional discrete element modeling method, there are some microscopic parameters defined in the FDEM (i.e., various penalty values, fracture energy) which cannot be measured through laboratory tests, and these parameters have important effects on the macroscopic mechanical performance of rock models. The FDEM input parameter calibration procedure proposed by Tatone et al. [36, 37] was used to calibrate these parameters. The final values of mechanical parameters are reported in Table 1. After the parameters are selected, the preprocessing software Y-GUI [38], developed by Grasselli’s Geomechanics Group of Canada, was used to assign the corresponding properties to the rock sample, load plate unit, and node.

4.3. Analysis of Simulation Results

Figure 10 shows the failure patterns of both simulated samples and 3D-printed samples during UCS tests. The simulated fracture patterns are considered as a good match with the laboratory results and all show obvious anisotropy.

During the process of uniaxial compression test, the generation of tensile crack and shear crack was different for simulated specimens with different inclined angles. When the inclined angle α = 0°, the orientation of bedding plane was parallel to the loading direction. A vertical splitting crack parallel to the bedding plane was generated and was intersected with the shear failure plane. When the inclined angle α = 22.5° and α = 45°, most of the cracks at failure were shear cracks, which are basically produced from the ends of the specimen and propagate along the bedding plane. When the inclined angle α = 67.5°, the sample no longer broke along the bedding plane, the split tensile failure dominated the failure process of the sample, and the cracks penetrated through the bedding plane. When the inclined angle α = 90°, the orientation of bedding plane was perpendicular to the loading direction, and the sample produced failure patterns of uniaxial compression test similar to that of the rock material [39]. The sample broke along the angle of , but, due to the influence of the bedding plane, the shear failure had some slip along the horizontal bedding.

Figure 11 illustrates that the response of all samples showed a steep drop after reaching the peak strength values, which can be characterized by elastic-brittle failure mode. The peak values of the simulated samples are 11.23 MPa, 4.91 MPa, 8.26 MPa, 8.60 MPa, and 9.43 MPa when the inclined angles are 0°, 22.5°, 45°, 67.5°, and 90°, respectively. These values are considered as a good match with the laboratory results (Figure 12). They showed high consistency and both revealed that the bedding plane direction had a great influence on the strength of the samples. As the inclined angle increased, the UCS of the simulated sample first decreased from 11.23 MPa (α = 0°) to the minimum value of 4.91 MPa (α = 22.5°) and then increased to reach 9.43 MPa (α = 90°), showing a U-shaped pattern. The influence of inclined angle with respect to the sample anisotropy was similar to the previous research results of shale formations.

The failure of rock mass is a change process from continuous to discontinuous rock mass, which is hard to be characterized by single continuum mechanics theory or discontinuous medium theory. However, the theory of continuum mechanics can well simulate the process before rock mass failure, while the failure mechanism can be more accurately reflected by the method of discontinuous medium mechanics [40]. The FDEM combines the advantages of the above two methods to establish a full scale analysis from microscopic fracture mechanism to macroscopic failure process, which is undoubtedly a better choice.

5. Conclusion

In this paper, laboratory tests and numerical tests were performed using the 3D printing technology and the combined finite-discrete element method. In order to systematically study the anisotropic mechanical behaviors of natural shale rock formations, we changed the printing direction of 3D-printed samples, as well as the inclined angle of FDEM simulated samples. Compared with the layered rock materials, the following conclusions are obtained:(1)It is necessary to increase the internal binder spraying of the sample, since the sample printed by the system default printing mode is not suitable for simulating rock materials. The test results show that when the binder saturation level is 175%, the UCS of the sample is close to the maximum value, reaching 9.62 MPa, and is only 3.7% lower than that printed with 200% binder saturation level. At this time, 175% binder saturation level was considered as the optimal value to maximize the UCS of the sample, and it is of little significance to increase it any more.(2)The 3D printing layering direction has a significant influence on the uniaxial compressive strength of the sample. When the printing inclined angle α is changed from 0° to 90°, the uniaxial compressive strength of the sample first decreases and then increases. When α is 0°, the UCS of the sample is the biggest, while it is the lowest when α is between 20° and 30°. When α is 45°, 67.5°, and 90°, the stress-strain curve is similar to that of ductile materials. When α is 0° and 22.5°, the ductility is greatly decreased.(3)A method based on the selection of a set of parameters reflecting the anisotropy of materials was introduced in the FDEM in order to simulate the physical and mechanical properties of layered rock. Due to the repeatability of 3D printing technology, the calibration of these microscopic parameters by 3D printing experiment is more accurate than the traditional calibration method of laboratory test, which is more suitable for subsequent large-scale simulation.

This paper initially explored the application of 3D printing technology in geotechnical engineering and made some breakthroughs in the form of the strength and stress-strain curve of printed samples. However, due to the limitation of printing technology and materials, the strength of 3D-printed samples is still much lower than that of natural rocks. Therefore, 3D printing technology can only be used as a qualitative analysis method at present. With the development of technology, the settlement of these limitations will be achieved and the 3D printing technology will be widely used in the field of geotechnical engineering.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.