Three-Dimensional Combined Finite-Discrete Element Modeling of Shear Fracture Process in Direct Shearing of Rough Concrete–Rock Joints

Abstract

A three-dimensional combined finite-discrete element element method (FDEM), parallelized by a general-purpose graphic-processing-unit (GPGPU), was applied to identify the fracture process of rough concrete–rock joints under direct shearing. The development process of shear resistance under the complex interaction between the rough concrete–rock joint surfaces, i.e., asperity dilatation, sliding, and degradation, was numerically simulated in terms of various asperity roughness under constant normal confinement. It was found that joint roughness significantly affects the development of overall joint shear resistance. The main mechanism for the joint shear resistance was identified as asperity sliding in the case of smoother joint roughness and asperity degradation in the case of rougher joint asperity. Moreover, it was established that the bulk internal friction angle increased with asperity angle increments in the Mohr–Coulomb criterion, and these results follow Patton’s theoretical model. Finally, the friction coefficient in FDEM appears to be an important parameter for simulating the direct shear test because the friction coefficient affects the bulk shear strength as well as the bulk internal friction angle. In addition, the friction coefficient of the rock–concrete joints contributes to the variation of the internal friction angle at the smooth joint than the rough joint.

1. Introduction

In a ground-anchored suspension bridge, the anchorage performs the function of supporting the tension force acting on the main cable. As one of the construction methods of anchorage, gravity-type anchorage (GTA) (see Figure 1) is widely used due to a variety of advantages such as simple mechanical configuration, ease of construction, and its vast applications [1]. The GTA resists the cable load via the weight of the anchorage structure itself, and this depends on the frictional resistance between the anchorages (concrete) and the foundation (rock). For designing the GTA, many guidelines are recommended to assess the stability of a GTA using the gravity method [2,3,4,5,6], and it is coupled with the Mohr–Coulomb (M-C) criterion, which is defined by several parameters such as friction angle and cohesion. However, some guidelines assume the M-C parameters to be a conservative constant value (e.g., $c=0$ and $tan\varphi = 0.6$) [5,6] in the vicinity of the joint between the anchorage and the foundation. The rough geometry of the joint is also neglected in the gravity method. Such oversimplifications related to joint properties may lead to a conservative stability analysis, and it may require a huge amount of earth and rock excavation and concrete consumption, which may result in high construction costs and environmental damage in GTA. Therefore, it is crucial to extensively investigate the frictional resistance and failure behaviors around the joint structure between concrete and rock surfaces (i.e., concrete–rock joint) with rough geometry under shearing.

The problem of the failure of the concrete–rock (C-R) joint under shearing has been considered as a special case of the failure of rock mass containing rough joint surfaces (i.e., rock–rock joint) under shearing. The problem of rock–rock (R-R) joints subjected to shearing has been investigated extensively, and several models for this problem have been proposed [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. These studies have examined the major factors that affect the behavior of the R-R joint, considering the rock mass properties around the rock joint, joint surface roughness, and boundary conditions applicable to joints, and some useful solutions have been proposed.

Based on the research and knowledge regarding the shear behavior of the R-R joints, several investigations of the interaction of C-R joint surfaces have been performed by researchers through direct shear tests on the C-R joints under constant normal loading (CNL) and constant normal stiffness (CNS) conditions. Thomas et al. [11] performed a direct shear test on regular triangular C-R joint surfaces under several CNS conditions. As an extension of their study, Johnson et al. [12] proposed analytical solutions to predict the development of the shear resistance in rock-socketed concrete piles, which can be idealized as a C-R joint problem under CNS conditions. Kodikara and Johnston [27] investigated the shear behavior of C-R joints with irregular triangle asperities, and they compared the obtained results with the case of regular triangular asperities. They revealed that the major difference between irregular and regular asperities is in the fact that an irregular joint presents more ductile behavior than a regular joint does. Similarly, Gu et al. [28] carried out a direct shear test with fractal profiles under several CNL conditions, and the results were compared with a case of regular triangular asperity. The response of the joint with fractal profiles showed more ductile behavior and presented a significantly lower peak overall strength than the case with regular asperity, because each single asperity in the fractal profile was progressively fractured rather than simultaneously fractured. Saiang et al. [29] investigated the shear behavior of a bonded shotcrete–rock joint. They found that for the development of the overall shear strength of the joint, the cohesive strength of the joint played an important role when the applied normal stress confinement was fairly low (0.2~0.5 MPa), while the friction was more important under higher normal stress confinement.

Based on experimental observation and theoretical analysis, several constitutive models have also been developed. Seidel and Haberfield [30] developed a constitutive model to predict the mechanism of the development of overall shear strength of the C-R joint. This constitutive model simulated important physical processes including asperity sliding, asperity shearing, and post-peak behavior of the concrete–rock joint. However, Harberfield and Johnston [31] have pointed out the limitations of existing empirical and analytical approaches. The obtained results posted significant dependencies on specific experimental site and test conditions for the existing empirical models, while the requirements of complex input parameters and oversimplifications showed difficulty in terms of the applicability of existing analytical approaches. As a result, they proposed a new theoretical model that can describe the main mechanism and requires joint roughness, basic rock properties, and boundary conditions for the cases of R-R and C-R joints. However, these constitutive models still could not have investigated the process of asperity degradation and shear fracturing during a direct shear test, indicating that the empirical-based and theoretical-based approaches suffer from significant limitations in their applications.

The viability of continuum and discontinuum-based numerical methods to numerically simulate the shear behavior of R-R and C-R joints has advanced alongside computational geomechanics technology. Tian et al. [32] used a three-dimensional (3D) FEM to develop a cohesive interface model of cemented C-R joints and simulated a direct shear test with a smooth joint surface. Using the discrete element method (DEM), Bahaaddini et al. [33] simulated the shear behavior of rock joints using the bonded particle model (BPM), which is available in the two-dimensional (2D) particle flow code (PFC2D). They modeled the shear behavior of a saw-tooth triangular bonded and unbonded R-R joint and considered the joint roughness coefficient (JRC) profiles [7] as they investigated the rock damage process during a direct shear test. Ultimately, they successfully simulate the propagation of tensile cracks and the asperity crushing. Using the BPM of the 3D particle flow code (PFC3D), Park and Song [34] modeled a series of direct shear tests with rock–rock joints and investigated the effect of the geometrical topology and the micro-properties of the joints on shear behavior. Although the DEM with the BPM successfully modeled the fracture asperity and gouges, this method required a more tedious calibration process to acquire the unknown micro-parameters until the response of the macro-scale model accorded with the laboratory test results rather than the FEM application. Many of the computational resources also accompany simulating the crack propagation of the rock and concrete than the FEM.

Recently, the combined finite-discrete element method (FDEM), originally proposed by Munjiza [35,36], has attracted significant attention in the field of computational solid mechanics. Since then, various FDEM codes have been developed in several research groups across the world [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54]. The FDEM can model the three important aspects of solid fracture, i.e., the continuous deformation in the intact regime, the transition from continuum to discontinuum, and the very complex contact mechanics between solid surfaces, including newly created macroscopic fracture surfaces. For this reason, the FDEM is expected to model the three-dimensionally complex asperity degradation and shear resistance process, and thus, this study employed the FDEM.

FDEM has been applied to the investigation of the fracture process in the direct shearing of the R-R joints. For example, Karami and Stead [53] utilized 2D FDEM implemented in ELFEN software to investigate the process of the R-R joint degradation during direct shear tests with various joint roughness under different normal stress. Similarly, Liu et al. [46] applied the 2D FDEM implemented in the Y-HFDEM IDE code to simulate the direct shearing process of rock with rough R-R joints under various dynamic loading and investigated the effect of different patterns of asperity roughness. They successfully simulated the development process of overall shear resistance including asperity sliding, asperity degradation, gouge formation, and gouge grinding. Furthermore, Tatone and Grasselli [54] implemented the simulation of a direct shear test using the 2D FDEM implemented in Y-Geo code for a 2D profiled rock–rock joint in which the detailed joint topology was extracted from the 3D discontinuity surface data obtained from the micro-CT scanning system, and the simulation results of the direct shear test presented a good agreement with laboratory observations over the first 2–3 mm of shear displacement. According to these application examples, the FDEM has great capabilities for detailed investigations of complex fracture processes in the shear test of R-R joints. However, all the previous research using FDEM was limited to 2D problems and R-R joint problems. In addition, in the experience of the author, the heterogeneous fracturing was shown in a laboratory direct shear test of the concrete–rock joint specimen with uniform and rough tooth-saw-shaped asperity. Therefore, to achieve a better understanding of the development of resistance mechanisms in the C-R joints under shearing, 3D numerical simulation of the fracturing process of concrete and rock around the C-R joints is necessary, because the concrete and rock are fractured involving complex 3D processes. However, no previous research has been reported for the application of 3D FDEM to the investigation of the direct shear test of the rough C-R joint under shearing.

In relation to the work summarized above, this paper aims to apply the 3D FDEM to investigate the detailed development of resistance mechanisms in a direct shear test of a rough C-R joint under shearing including asperity dilatation, asperity sliding, and asperity degradation behavior. The organization of this paper is as follows: Section 2 briefly introduced the 3D FDEM implemented in the Y-HFDEM integrated development environment (IDE) whose computational performance has been greatly improved by the introduction of parallel computation utilizing general-purpose-graphic-processing-unit (GPGPU). In Section 3, calibration of the input parameters for the 3D FDEM simulations in Section 4 was first conducted against the experimental results of the target rock and concrete through the uniaxial compressive strength (UCS) and Brazilian tensile strength (BTS) tests. Using calibrated input parameters, a direct shear test of the rough C-R joint was investigated in detail using 3D FDEM. Finally, the authors discuss the effect of the roughness of asperity as well as the effect of friction coefficient on friction resistance during direct shear development.

2. GPGPU-Parallelized 3D Y-HFDEM IDE

In this study, an FDEM implemented in 3D Y-HFDEM IDE code was applied to simulate the fracturing of a direct shear test of a rough concrete–rock joint under shearing. Originally, the code was developed using the object-oriented programming tool (i.e., Visual C++) [44,45,46], and it was based on non-parallelized open-source Y-2D/3D libraries with OpenGL [36,37,38]. Using this code, the various rock engineering problems have been successfully solved by simulating rock fracturing [44,45,46]. However, due to the non-parallelized nature of the code, its previous applications were limited to small-scale 2D problems using relatively rough meshes. Recently, to overcome these limitations, a parallel programming scheme using GPGPU controlled by compute unified device architecture (CUDA) C/C++ was incorporated to develop the GPGPU-parallelized Y-HFDEM IDE code by the authors [47,48,49,50,51] for 2D and 3D FDEM modeling. Since the full details of the Y-HFDEM IDE code can be found in the literature [49,50], we focus in this section on the essential parts of the 3D FDEM in order to model the fracture process in the direct shearing of the rough concrete–rock joint. In the remaining sections of this paper, the 3D FDEM implemented in the GPGPU-parallelized Y-HFDEM 3D IDE code is simply called the 3D FDEM.

In the 3D FDEM, the continuous, i.e., elastic, deformations of rock and concrete are modeled using continuum solid elements (i.e., 4-noded tetrahedral elements, TET4s) [50].

The fracturing in the rock and concrete are modeled by an intrinsic cohesive zone model (ICZM) with the concept of smeared crack [55] representing opening and sliding cracks as well as mixed-mode I–II fracturing. In this model, 6-node initially zero-thickness cohesive elements (CE6) are inserted into all the boundaries of the TET4s at the beginning of the FDEM simulation to model the fracture of the rock and concrete (Figure 2). The opening and sliding of the faces in each CE6 are governed by tensile and shear softening curves as a function of the crack opening ($w_o$) and sliding displacements ($w_s$) (Figure 3). In the tensile and shear softening curves, the normal and shear cohesive tractions (i.e., ${\mathsf{\sigma }}_{\mathrm{coh}}$ and ${\mathsf{\tau }}_{\mathrm{coh}}$, respectively) acting on each side of CE6 are calculated using Equations (1) and (2), respectively:

$${\sigma }_{coh}=\left\{\begin{array}{ll}\frac{2w_o}{w_o^{overlap}}{T}_s& if{w}_o<0\\ \left[\frac{2w_o}{w_o^p}-{\left(\frac{w_o}{w_o^p}\right)}^2\right]f(D)T_s& if0\le w_o\le w_o^p\\ f(D)T_s& if{w}_o^p<0\end{array}\right.$$

(1)

$${\tau }_{coh}=\left\{\begin{array}{ll}\left[\frac{2\left|w_s\right|}{w_s^p}-{\left(\frac{\left|w_s\right|}{w_s^p}\right)}^2\right](-{\sigma }_{coh}\tan (\varphi )+f(D)c)& if0\le \left|w_s\right|\le w_s^p\\ -{\sigma }_{coh}\tan (\varphi )+f(D)c& if{w}_s^p<\left|w_s\right|\end{array}\right.$$

(2)

where $w_o^p$ and $w_s^p$ are the artificial elastic limits of $w_o$ and $w_s$, respectively. $w_o^{overlap}$ is the representative overlap controlling the artificial elastic behavior of CE6s when o is negative, i.e., overlapping. $T_s$, c, and $\varphi$ are the microscopic tensile strength, cohesion, and internal friction angle, respectively, of the CE6. The $w_o^p$, $w_s^p$, and $w_o^{overlap}$ are defined as $w_o^p=2h{T}_s/P_{open}$, $w_s^p=2hc/P_{tan}$, and $w_o^{overlap}=2h{T}_s/P_{overlap}$, respectively, where $h$ is the element length, and $P_{open}$, $P_{tan}$ and $P_{overlap}$ are the artificial penalty terms of the CE6 for opening, sliding, and overlapping, respectively. Equation (2) corresponds to the Mohr–Coulomb (MC) shear strength model with tension cut-off. The function $f\left(D\right)$ is the characteristic function representing the tensile and shear softening curves (Figure 3), and it changes according to a damage value $D$ of CE6. The damage value $D$ is used to consider not only the opening and sliding of the faces of CE6 but also mixed-mode fracturing using the equation below [55]:

$$\begin{gathered} D=Minimum\left(1,\sqrt{{\left(\frac{w_o-w_o^p}{w_o^t}\right)}^2+{\left(\frac{\left|w_s\right|-w_s^p}{w_s^t}\right)}^2}\right) \\ if{w}_o>w_o^{p}or\left|w_s\right|>w_s^p,otherwise0 \end{gathered}$$

(3)

where $w_o^t$ and $w_s^t$ are the critical values of $w_o$ and $w_s$, respectively, in which CE6 breaks and develops as macroscopic fracture surfaces. Meanwhile, the $w_o^t$ and $w_s^t$ in Equation (3) satisfy the mode I and mode II fracture energies $G_{fI}$ and $G_{fII}$, respectively, which are expressed by Equations (4) and (5):

$$G_{fI}={\displaystyle \underset{w_o^p}{\overset{w_o^t}{\int }}{\sigma }_{coh}(w_o)d{w}_o}$$

(4)

$$G_{fII}+W_{res}={\displaystyle \underset{w_s^p}{\overset{w_s^t}{\int }}{\tau }_{coh}(\left|w_s\right|)d\left|w_s\right|}$$

(5)

where $W_{res}$ is the amount of work per CE6′s area performed by the residual stress term in the MC shear strength model. In the present formulation, mode II and III fracturing modes are not distinct. Mode II and mode III responses to microcracks (i.e., CE6) are assumed to be simply explained by the parameters $G_{fII}$, because the clear definition of the crack tips and crack planes is difficult to define in the 3-dimensionally complex fracture networks in FDEM. This is an aspect that is contrary to ordinary fracture toughness tests.

3. Calibration of the Combined Finite-Discrete Element Method by Modeling UCS and BTS

An appropriate determination of the model parameters is essential for accurate simulation of the shearing process in the C-R joint using FDEM. In this section, the failure processes in the UCS and BTS tests, which are frequently performed for acquiring mechanical properties of rock-like materials (i.e., rocks and concrete) are simulated using the 3D FDEM. The numerical simulation results are compared with experimental data to calibrate input parameters in the 3D FDEM. Figure 4 shows the 3D FDEM models for UCS and BTS tests. The calibration works are performed for a granite specimen and concrete specimen with nominal UCS = 35 MPa (hereafter, 35 MPa-concrete). The UCS and BTS test models are discretized around 1.5 mm of the tetrahedral element where the element size is sufficient to exclude the mesh sensitivity effect for the specimen size [47]. The experimentally obtained physical and mechanical properties of the granite and concrete specimens are summarized in

Table 1. Experimentally obtained physical and mechanical properties of 35 MPa concrete and granite.

Parameter	Unit	Value
		35 MPa Concrete	Granite
Mass density	$\mathrm{kg}/{\mathrm{m}}^3$	2200	2700
Young’s modulus	$\mathrm{GPa}$	30.0	46.0
Poisson’s ratio	-	0.20	0.21
Tensile strength	$\mathrm{MPa}$	2.62	10.0
Cohesion	$\mathrm{MPa}$	7.87	26.0
Internal friction angle	Degrees	35.0	58.1
Uniaxial compressive strength	$\mathrm{MPa}$	35.0	202

. The input elastic parameters for the 3D FDEM models for granite and concrete were taken from the values in Table 1. Considering the literature [47], the normal contact penalty (P_n) and cohesive penalty terms ($P_{open}$, $P_{tan}$, $P_{overlap}$) of the granite, concrete, and steel plates are assumed to be $P_n=10E$, $P_{open}=100E$, $P_{tan}=100E$, and $P_{overlap}=1000E$ where $E$ indicates the Young’s modulus of either rock (E_rock) or concrete (E_conc). The granite model or 35 MPa concrete model was placed between the top and bottom rigid loading plates with a slow velocity of 0.01 $m/s$ [47] along with a critical damping scheme to approximately achieve a quasi-static loading condition [47,48,49,50].

Figure 5 and Figure 6 shows the comparison of the FDEM simulation and experiment during the UCS tests of the granite and the concrete, respectively. Figure 5a and Figure 6a present the plot of axial stress versus axial strain, and Figure 5b and Figure 6b indicate the final fracture pattern between numerical simulation and experiment test. The FDEM simulation result of the UCS test shown in Figure 5 and Figure 6 presents good agreement with the experiment result. Meanwhile, Figure 7 and Figure 8 show the comparison of the FDEM simulation and experiment during the BTS test of granite. In Figure 7a and Figure 8a, the plot of indirect tensile stress versus axial displacement is presented. Whereas in Figure 7b and Figure 8b, the final fracture pattern between numerical simulation and experiment test is on display. The FDEM simulation results of the BTS test shown in Figure 7 and Figure 8 are in agreement with the experiment results. Through the aforementioned calibration procedure of granite and 35 MPa concrete against the UCS and BTS tests, the calibrated input parameters for FDEM simulation of the direct shear test of the rough concrete–rock joint under shearing in the next section are summarized in

Table 2. Calibrated FDEM input parameters for the simulation of fracture behaviors of concrete and rock elements in the direct shear tests. Note that the loading platens are rigid in Section 3 while elastic in Section 4.

Parameter	Unit	Value
		35 MPa Concrete	Granite
Density ($\rho$)	$kg/m^3$	2200	2700
Young’s modulus ($E$)	$GPa$	30.0	46.0
Poisson’s ratio ($v$)	-	0.20	0.21
Microscopic tensile strength ($T_s$)	$MPa$	2.10	20.0
Microscopic cohesion ($c$)	$MPa$	6.70	56.0
Microscopic internal friction angle ($\varphi$)	Degrees	35.0	59.0
Microscopic mode I fracture energy ($G_{If}$)	$N/m$	30.0	60.0
Microscopic mode II fracture energy ($G_{IIf}$)	$N/m$	90.0	120.0
Artificial cohesive penalty ($P_{open}$, $P_{tan}$)	$GPa$	3000	4600
Artificial cohesive penalty ($P_{overlap}$)	$GPa$	30000	46000

and

Table 3. Calibrated FDEM input parameters for simulating contact behavior in the simulation of direct shear tests.

Parameter	Unit	Value
Normal contact penalty between rock and rock surface ($P_{n\_con}^{RR}$)	$GPa$	460.0
Normal contact penalty between concrete and concrete surface ($P_{n\_con}^{CC}$)	$GPa$	300.0
Normal contact penalty between concrete and rock surface ($P_{n\_con}^{CR}$)	$GPa$	300.0
Normal contact penalty between rock and platen surface ($P_{n\_con}^{RP}$)	$GPa$	460.0
Normal contact penalty between concrete and platen surface ($P_{n\_con}^{CP}$)	$GPa$	300.0
Friction coefficient between rock and rock surface (${\mu }_{\mathrm{fric}}^{\mathrm{RR}}$)	-	0.6
Friction coefficient between concrete and concrete surface (${\mu }_{\mathrm{fric}}^{\mathrm{CC}}$)	-	0.6
Friction coefficient between concrete and rock surface (${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$)	-	0.6
Friction coefficient between specimen and platen surface (${\mu }_{\mathrm{fric}}^{\mathrm{RP}}$, ${\mu }_{\mathrm{fric}}^{\mathrm{CP}}$)	-	0.1

.

4. 3D FDEM Modeling of Rough Concrete–Rock Joint Failure in a Direct Shear Test

4.1. Modelling Procedure

In this section, the 3D FDEM was employed to model the asperity dilatation, sliding, and degradation during the direct shear test of a rough C-R joint. As shown in Figure 9a, the 3D FDEM model has been generated and originally adopted as the direct shear test specimen. To compare the numerical simulation result with the direct shear test result, a simulation model was generated as an equivalent dimension on the experimental specimen. To consider the effect of asperity roughness of the C-R joint, we prepared three direct shear test models having the asperity angles of α = 40, 25, and 10 degrees (see Figure 9a). Note that the asperity angle is defined as the angle from the specimen base to the inclined plane. There is no suggested method or standards for the direct shear test of the C-R joint. The size of the direct shear test specimen was determined by considering the International Society for Rock Mechanics and Rock Engineering (ISRM) suggested method for the direct shear test of rock, which is presented in this literature [56]. Based on the ISRM suggested method, the test plane of the rock specimen should be exceeded by 10 times the maximum of the asperity height. In this study, we determined the asperity height between 5 and 25 mm, and the specimen width, length, and height were determined as 125 mm, 300 mm, and 230 mm, respectively. The FDEM model consists of discrete 35 MPa concrete and granite blocks. Both the concrete and granite blocks have a saw-tooth shape at the one end of each block. Along the saw-tooth shaped ends, the two blocks are placed in a way that each end perfectly matches to form the C-R joint. To apply normal and shear loadings to the blocks with the rock–concrete joint, we considered two half-boxes (hereafter called loading platens) partially covering the concrete and granite blocks (Figure 9b).

To provide a better understanding of the mechanical behavior of the C-R joint under self-weight loading, only a direct shearing test under CNL condition is considered here. The CNS conditions employed in most of the previous research [27,28,29] for describing C-R joints were not considered because the main interest of this paper is in the joint between the foundation and gravity anchorage in which the normal stress of the rock joint is expected to remain constant during shear deformation. In such cases, it is recommended to adopt the direct shear test with a CNL [57].

In the 3D FDEM simulation, normal confinement stress σ_n was applied to the top surface of the top loading platen (Figure 9b), while the bottom loading platen was kept in space to model the CNL condition. Then, the FDEM simulation continued until the static stress equilibrium was achieved in the C-R joint specimen. Note that the local damping scheme using a damping coefficient = 0.8 (coupled with mass scaling) [58] was utilized to accelerate the convergence rate to a static stress state. The cases of σ_n = 1.0, 2.0, and 3 MPa were considered. After this, the top loading plate was moved along the shear directions (i.e., horizontal direction in Figure 9b) with a constant velocity v = 0.01 m/s, while the bottom loading plate was fixed in space. Accordingly, shearing the discrete concrete–rock blocks under the CNL condition can be realized along the saw-tooth shape joint. It is worth mentioning that the damping scheme was switched from the local damping scheme (with mass scaling) to the critical damping scheme [36] at this stage to achieve the quasi-static loading condition. Otherwise, unrealistic spurious fracturing occurs due to unrealistic stress pile-up near the loading region. The material properties of concrete and granite blocks are set to the calibrated parameters in the previous section (Table 2 and Table 3).

4.2. Results

First, we present the results for case of the asperity angle α = 40 degrees. Figure 10 presents the shear stress–shear displacement curve obtained during our direct shear test simulation. Note that the shear displacement d_s was defined by the horizontal displacement of the top-loading platen while the bulk shear stress τ_bulk was computed by dividing the horizontal contact force of the top-loading platen as an area (i.e., Area = Specimen width × Specimen length). Figure 11 shows the macroscopic crack propagation (i.e., CE6s with D =1 in Equation (3)) by white color (first column in Figure 11) and the maximum principal stress distribution (PS1) and minimum principal stress distribution (PS3) at different stages of the simulation. The tension stress is positive in PS1 and PS3. At the onset of the direct shear test, the constant σ_n was initially applied in the vicinity of the C-R joint as the boundary condition (i.e., CNL condition). When d_s starts to increase, the concrete asperity sliding along the inclined C-R asperity surface and simultaneously surface friction is applied to the contact area of C-R joint surfaces, as shown in Figure 11a-I. In this stage, the tensile stress is induced near the contact asperity, as illustrated in Figure 11b-I. As these tensile stresses reached the tensile strength of CE6 in the vicinity of the asperity, the tensile crack was first initiated and propagated in a perpendicular direction to the contact area (Figure 11a-II). As shown as Figure 11a-IV, with the minor increase in d_s, the existing tensile crack was continuously propagated in a perpendicular direction to the contact area, while a new tensile crack was created by the tensile stress reaching the strength in the vicinity of the asperity at different locations. The accumulation of damage induced by the continuous tensile and shear crack propagation significantly reduced the ability of asperity resistance by the further increase in the d_s. As the tensile crack extends into the bases of the asperity (Figure 11a-V), the asperity is completely separated from their bases. Meanwhile, during the single asperity fracturing, other asperities were also fractured by the same fracturing process (see Figure 11III–V). When all of the asperities are completely separated from their bases as shown in Figure 11b-VI (marked as the yellow line), the overall behavior of the vicinity of the C-R joint behaves as a smooth surface. Figure 10 reflects the above-mentioned shear failure behavior of the C-R joint in terms of the bulk shear stress–shear displacement curve. Initially, the bulk shear stress increased linearly with the shear displacement (Figure 10-I), which was caused by dilation along the normal direction of the C-R joint. The dilatation was induced by sliding the C-R joint along the inclined surface and the dilatation effect of the friction force between the C-R joint by increasing normal forces to the friction area. When the tensile crack is initiated, the τ_bulk is slightly decreased with the d_s increment (Figure 10-II). In this stage, the PS3 distribution changed along the tensile crack because the single asperity behaved as sliding along the crack surface rather than dilatation as the d_s increases. However, the τ_bulk still increased due to local normal stress increment in the other asperities (see Figure 10II–III), until the τ_bulk approached the shear resistance limit, τ_lim (Figure 10-III). In this stage, the fluctuation of τ_bulk is presented because during the single asperity fracture, the other asperities are in various fracturing processes (i.e., crack initiation and propagation). After further d_s, the accumulation of damage induced by continuous tensile crack propagation significantly reduced the ability of asperity resistance. This stage corresponds to the transition from the post-peak reduction of τ_bulk to residual stress regime (Figure 10-IV–V). As the tensile crack extended into the bases of the asperity (Figure 10-VI), the asperity was completely separated from their bases and behaved as a smooth surface passing through the base of the concrete asperity. In this condition, the effect of the roughness of the concrete–rock joint became almost negligible, and the main factor for shear resistance was changed to sliding, which was subjected to normal loading. In Figure 10, before the shear softening, the numerical simulation presented good agreement with the experiment result in the τ_lim and the d_lim. However, in the regime of residual shear resistance, the shear resistance behavior was different from the laboratory direct shear test result. In the regime of the residual resistance, the numerical simulation presented as τ_bulk remained almost in constant residual shear stress, but the experiment showed a continuous decrease after the peak τ_bulk. Comparing to the resultant fracture pattern between the simulation and experiment result as shown in Figure 12, the smooth surface was generated in both cases. However, additional damages (i.e., cracks and fragmentation) were found in the experiments, as shown in Figure 12b.

5. Discussion

5.1. Direct Shearing of Concrete–Rock Joint with Different Asperity Roughness

The direct shear test of concrete–rock joints with various roughness and regular triangular asperity was simulated using the 3D FDEM to investigate the effect of joint roughness on the shear failure process. Figure 13 and Figure 14 compare the three concrete–rock joint models with the asperity angles of α = 10 and 25 degrees in terms of the distributions of the maximum principal stress (PS1) and the macroscopic cracks. The case of α = 40 degrees was already presented in Figure 11-VI. Similarly to Figure 10, the τ_bulk−d_s curves for the three asperity angles α = 10, 25, and 40 degrees in terms of σ_n = 2.0 MPa are compared in Figure 15. As already explained using Figure 11-VI for the case of rougher asperity (α = 40 degrees), the main fracture mechanism in overall shear resistance was joint asperity degradation, which was induced by tensile fracturing within the intact asperity. As can be seen from Figure 13 and Figure 14 for the relatively smoother asperity (α = 10 and 25 degrees), sliding was the main mechanism in the overall shear resistance. Since the smoother single asperity was more free from the bending moment, which more easily induced the tensile crack than the rough single asperity, the only factor of the shear resistance in the smooth asperity was the friction between the top (concrete) and bottom (rock) asperity. Figure 15a indicates that the relatively rougher asperity (α = 40 degrees) shows clear shear softening behavior (transition from peak bulk shear stress to reduction) due to the joint asperity degradation. With the decrease of the asperity angle α, the shear softening behavior due to the joint asperity degradation became less obvious and in the case of α = 10 degrees only showed constant τ_bulk, i.e., a residual stress state immediately after the initial elastic behavior. Compared to the development of shear resistance between the simulation and experiment result in terms of the τ_bulk−d_s curves as shown in Figure 15a,b, the overall development of shear resistance presents similar behaviors at each case of the α. However, the fluctuation of the τ_bulk is observed around the peak bulk shear stress in the numerical simulation, while the τ_bulk shows the stress-drop after the peak bulk shear stress in the experiment.

To interpret the obtained results more clearly, the FDEM simulation results for various values of α (α = 10, 25 and 40 degrees) and σ_n (i.e., σ_n = 1.0, 2.0 and 3.0 MPa) were fitted using a linear fitting line based on the M-C shear failure criteria as shown in Figure 16. To avoid confusion caused by notation between the microscopic cohesion and internal friction angle, which is summarized in Table 2, and the parameters of the M-C shear failure criteria, the M-C parameter is presented as the bulk cohesion (c_bulk) and bulk internal friction angle ($\varphi$_bulk). In addition, the c_bulk and $\varphi$_bulk reflect the overall shear resistance in the vicinity of the C-R joint, unlike the microscopic cohesion and internal friction angle represented in a single material. In this figure, the shear strength for each case was computed as the peak value of each of the τ_bulk−d_s curves for different conditions. In this M-C shear failure criteria, the bulk cohesion (c_bulk) is considered as zero because of an unbounded condition between the concrete and rock surface. In addition, the bulk internal friction angle ($\varphi$_bulk) was determined by calculating the slope of each linear fitting line. The figure clearly shows that the increase in normal stress increments based on the M-C theory and the internal friction angle ($\varphi$_bulk) also increased with asperity angle increments. These simulation results follow Patton’s theoretical model [10], which represents the effect of joint roughness as the increment of the internal friction angle, $∆\varphi$, (e.g., ${\tau }_{bulk}=tan\left(\varphi +∆\varphi \right)$) when the asperity angle, i.e., joint roughness, increases. Ultimately, such results validate the applicability of the 3D FDEM to this class of problem involving joint surface interaction.

5.2. Effect of Fiction Coefficient on Direct Shearing of Concrete–Rock Joint

The effect of the friction coefficient on the direct shearing of the rough C-R joint was examined by performing a direct shear test on the C-R joint with a different friction coefficient. In an FDEM simulation, the friction coefficient has been applied for calculating shear forces on the contacting macro crack surface in a single material and contacting surface between different materials. The detailed description for the calculation of the contact and shear force can be found in the literature [49,50]. For calculating the shear forces on contacting surfaces, several FDEM numerical simulations [39,47,50], including this study, have applied the constant friction coefficient between rock and rock as 0.5–0.6 when a macroscopic crack surface is newly generated. In addition, the constant friction coefficient between the rock and loading platen is determined as 0.1. Using these values, they acquire quite similar results in essential quasi-static test simulations such as UCS and BTS. However, the friction coefficient between different materials (e.g., C-R joint) should carefully be considered, especially in the direct shear test. The friction force, varied by friction coefficient, may present a different fracture pattern and bulk shear strength in the vicinity of the C-R joint. Moreover, the friction coefficient may affect the internal friction angle ($\varphi$_bulk) due to how the basic friction angle [6,10], which was acquired by direct testing a flat joint, determines the internal friction angle in the rock joint as a similar parameter with a friction coefficient in the FDEM. In this sense, the inspection of the effect of the friction coefficient on the C-R joint should carefully be considered.

In this study, the friction coefficient, expressed as ${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$ in Table 3, varied between 0.5, 0.6, and 0.7. Under different friction coefficients, the simulation of a direct shear test was performed on the C-R joint with α = 10 and 40 degrees of saw-tooth asperity under various values of σ_n (i.e., σ_n = 1.0, 2.0, and 3.0 MPa).

As can be seen in the Figure 17, the τ_bulk is presented as a function of σ_n in terms of different ${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$ (i.e., ${\mathsf{\mu }}_{fric}^{CR}$= 0.5, 0.6, and 0.7) and α (i.e., α = 10 and 40 degrees). Figure 17 shows that the bulk shear strength increases with ${\mathsf{\mu }}_{fric}^{CR}$ increments in both α = 10 and 40 degrees. In Figure 17a,b, the $\varphi$_bulk may be more dependent on the ${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$ in the α = 10 than 40 degrees because the $\varphi$_bulk was almost constant with the α = 40 degrees, while the $\varphi$_bulk increases with the friction coefficient increment in the α = 10 degrees. Sliding is the main mechanism in the smooth joint (i.e., α = 10) and the ${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$ directly increased the friction as well as the bulk shear strength at the concrete–rock joint. The asperity degradation occurs in the rough joint (i.e., α = 40) in the vicinity of the concrete asperity base. The friction applied at the macroscopic crack in the concrete part and $\varphi$_bulk may be independent with the variation of ${\mu }_{\mathrm{fric}}^{\mathrm{CR}}$.

6. Conclusions

In this study, a GPGPU-parallelized FDEM was applied to model the shear failure process of a rough concrete–rock joint in direct shear tests under constant normal loading conditions. Before simulating the direct shear test, the combined finite-discrete element method was calibrated by simulating the concrete and rock fracture process in the UCS and BTS tests and was compared with the experimental results. Then, the direct shear test of the concrete–rock joint was modeled using the combined finite-discrete element method to investigate the shear resistance process of a concrete–rock joint. We focused on the asperity dilatation, sliding, and degradation in terms of various asperity roughness under different constant normal loading (CNL) conditions. The effect of the joint roughness and friction coefficient between the concrete and rock interface on the behavior of the shear resistance was also discussed. In this study, the following conclusions can be drawn:

(1)

The GPGPU-parallelized FDEM reproduced the asperity dilatation, sliding, and degradation on the direct shear test of the concrete–rock joint with rough saw-tooth triangular asperity. Compared with laboratory observation, the simulation result presented a good agreement on the bulk shear stress–horizontal displacement curve before the regime of residual shear resistance. However, the laboratory test results presented a continuous decrease of bulk shear stress unlike the numerical simulation result, which showed constant residual stress after shear softening. In this study, we concluded that the additional damages found in the experiment results induced a continuous stress drop and presented different behaviors in the numerical simulation in the regime of the constant residual stress statement.

(2)

The roughness of asperities in the concrete–rock joint affects the overall joint shear resistance, where the main fracture mechanism is asperity sliding in the smooth asperity and asperity degradation in the rough asperity. The smoother single asperity is freer from the bending moment, which was easy to induce in the tensile crack in concrete asperity base, than the rough single asperity. The only factor of the shear resistance in the smooth asperity is the friction between the concrete and rock joint interface. The effect of the increment of asperity roughness affects the shear strength increment and causes the asperity degradation, which presents the shear softening behavior after the peak bulk shear stress, while the smooth asperity roughness presents the residual stress statement immediately after the initial elastic behavior. In addition, when the direct shear test results are presented as Mohr–Coulomb criterion, the bulk internal friction angle was clearly increased with an increasing asperity angle. In this study, the effect of asperity roughness on the bulk internal friction angle follows the Patton’s theoretical model.

(3)

The friction coefficient is an important parameter for simulating the direct shearing of the concrete–rock joint. Especially, the friction coefficient was more effective on the internal friction angle in terms of the smooth joint than the rough joint. The asperity degradation was the main mechanism in the rough joint, the concrete asperity base was fractured, and the friction was applied at the macroscopic crack surface in the concrete asperity base. Thus, the internal friction angle in the vicinity of the rough concrete–rock joint may be independent of the friction coefficient between the rock and concrete. However, in the smooth joint, the main mechanism of shear resistance was sliding without fracturing, and the different friction coefficient was directly reflected to the friction force as well as the internal friction angle.