Numerical Simulation on Shear Behavior of Double Rough Parallel Joints Under Constant Normal Stiffness Boundary Condition

The shear characteristics of rock joints under constant normal stiffness (CNS) boundary condition are critical to the stability of underground engineering rock mass. Five different rock joint profiles with random morphology were constructed using the independent segmentation method of the Hurst exponent. Based on the continuously yielding joint model, the discrete element calculation models of double rough parallel joints with different joint surface roughness and spacing of joints under CNS boundary condition were established using UDEC software to investigate the shear effect of joints under CNS boundary condition. The results show that under CNS boundary condition, the shear stress, normal displacement, normal stress, and surface resistance index (SRI) all increase with the increase of joint surface roughness coefficient (JRC) for both single- and double-joint specimens, and these of the single-joint are greater than these of the double-joint. With the increase of the joint spacing d, the peak shear stress τ y , and the peak surface resistance index (SRI p ) show a gradually increasing trend, and the influence of d on τ y and SRI p is greater with the increase of JRC, and both τ y and SRI p show a linear increasing trend with the increase of JRC. With the increase of d for the same JRC, both normal displacement and normal stress show a gradually increasing trend, and also, the increased amplitudes gradually increase; with the increase of JRC for the same d, the two parameters also show an increasing trend.


INTRODUCTION
A large number of joints that exist in natural rock masses can play a vital role in controlling the shear and deformability properties of engineering rock masses. A critical failure mode in rock masses is the shearing of rock joints, and thus, the shear characteristics of joints have been an important research topic in the field of rock mechanics (Bahaaddini, 2017). For most underground projects, the shear failure of rock joints is constrained by the surrounding rock, and the normal stress increases with increasing shear dilation displacement. Therefore, the constant normal load (CNL) boundary condition is no longer applicable, and the constant normal stiffness (CNS) boundary condition should be used (Indraratna and Thirukumaran, 2015;Thirukumaran and Indraratna, 2016;Shrivastava and Rao, 2018;Liu R. C. et al., 2020). At present, the study of the shear properties of single joints has been very mature. However, joints in natural rock masses cannot often exist alone, and the interaction between multiple joints has an important effect on the mechanical properties of rock masses.
Conducting a direct shear test is a standard experimental methodology to evaluate the shear behavior of rock joints (Saadat and Taheri, 2020). Shear failure of rock joints under CNL boundary conditions has been widely investigated due to its simplicity and ease of experimentation. Patton (1966) conducted direct shear tests on saw-toothed synthetic rock joints and proposed a bilinear envelope to predict the shear resistance of rock joints under different constant normal stresses.  proposed a fractal model to predict the shear behavior of a large-scale rock joint during shearing under CNL conditions. Oh et al. (2017) numerically studied the shear dilation of a saw-toothed rock joint under CNL condition. Zhou et al. (2019) studied the shear deformation characteristics of marble dentate joints with different dentate heights under different normal stresses and proposed an empirical formula to evaluate the shear dilation effect. These works suggested that the shear-induced dilation depends strongly on the relative normal stress. Compared with the CNL boundary condition, the CNS boundary condition is more suitable for deep-seated joints. GU et al. (2003) found that the joint surface wear was related to the joint surface roughness and normal boundary conditions. Lee et al. (2014) studied the effects of loading conditions and rock properties on the surface resistance index (SRI) and normal deformation characteristics of the joints. Li et al. (2018) established an analytical model to predict the shear behavior of rough rock joints under the CNS condition. Cui et al. (2019) carried out shear tests under CNL and CNS boundary conditions and investigated the effects of boundary conditions, initial normal stress, and joint roughness on the shearing behavior of the artificial joints. Han et al. (2020) experimentally investigated the effects of initial normal stress, normal stiffness, and shear rate on the shear properties of rough fracture surfaces and proposed empirical equations to assess the shear stress, normal stress, and normal displacement under cyclic shear loading. Yin et al. (2020) investigated the effects of initial normal stress and joint surface roughness coefficient (JRC) on shear stress, normal displacement, normal stress, and shear wear characteristics of fracture surfaces. Thirukumaran and Indraratna (2016) and Han et al. (2021) summarized the recent progress on shear characteristics of rock joints under CNS boundary conditions. The studies discussed earlier on CNS boundary conditions are mainly based on laboratory experiments and elucidate the important influence of CNS boundary conditions on the shear properties of the joint. As an alternative to experimental tests, many research scholars used advanced numerical techniques to characterize the shear failure of rock joints. Wang et al. (2019) investigated the effect of loading direction and normal stress on shear anisotropy and shear dilation deformation of jointed coal rocks by PFC software. Indraratna and Haque (2000) examined the shear behavior of soft saw-toothed rock joints under CNS conditions using UDEC software. Timothy (2018) investigated the effect of shear dilatancy on shear properties of rock joints after yielding under CNL and CNS boundary conditions using UDEC software. The numerical simulation discussed earlier mainly uses the discrete element method (DEM). This method provides an important reference to successfully carry out the numerical simulation of shear characteristics of rock joints under CNS boundary conditions. However, their study is limited to a single joint without taking into account the multi-joint interaction. A few studies (Han, 2019;Huang et al., 2021) have focused on the shear behavior of specimens containing two or more joints by laboratory experiments. However, the impacts of JRC and joint spacing on the shear behavior cannot be considered systematically due to the lack of a wide range of specimen tests.
To date, few studies have been reported on the shear properties of multi-rough rock joints under CNS boundary conditions (Liu et al., 2017). Although laboratory testing is the most common approach for investigating the shear mechanism of rock joints, setting up experiments is expensive and time-consuming (Saadat and Taheri, 2020). Experimental investigations require researchers to generate a wide range of specimens using various materials and conduct laboratory experiments using advanced laboratory equipment. Compared with laboratory tests, numerical simulations can consider more factors and have significant advantages in time, cost, and complexity (Li et al., 2021a;Li et al., 2021b). DEM based on discontinuous media is a direct modeling approach of joint shear damage and has been widely and successfully applied in modeling the shear behavior of rock joints. Numerical codes based on DEM include UDEC, 3DEC, PFC2D, PFC3D, etc. This method uses contact relations between discrete units instead of complex intrinsic relationships, and it can effectively simulate the macroscopic mechanical behavior of block or granular structural materials. In addition, due to the great complexity in natural jointed systems in surface morphology and spatial distributions, tests on models with parallel joints are still an effective approach for understanding the shear behavior of jointed rock masses. In this paper, a distinct element code UDEC was used to study the shear characteristics of double rough parallel joints under CNS boundary conditions. The evolutions of shear stress, normal stress, normal displacement, and SRI of parallel joints with different JRCs and joint spacing during the entire shear process were analyzed.

CONSTRUCTION OF RANDOM MORPHOLOGY ROCK JOINT
The JRC was proposed by Barton (1973). JRC is used to describe the surface irregularities, and it ranges from 0 (completely smooth and planar) to 20 (rough and uneven). Generally, JRC is obtained by examining the joint surface using Barton's comb and comparing it with the 10 standard rock profiles proposed by Barton and Choubey (1977). This method is obviously empirical and significantly depends on personal experiences. Therefore, many scholars (Tse and Cruden, 1979;Li and Huang, 2015;Wang Frontiers Tse and Cruden (1979) is widely accepted. Based on this method, the independent segmentation method of Hurst exponent was used to simulate random morphology rock joint, and the specific implementation process is as follows (Zhao et al., 2013).
(1) Suppose the initial joint profile be a straight line; (2) Split the line with a random point P; (3) On both sides of the point P, bisect the line with 0.5 mm equipartition distance; (4) In the vertical direction, the offset of each equipartition point on both sides of the split point with respect to the previous equipartition point is P(x).
where W is the amplitude-dependent parameter; R is a normally distributed random variable with a mean of 0 and a variance of 1; x is the distance from the equipartition point to the point P; H is the Hurst exponent.
(5) Perform JRC calculation of random morphology rock joint. Tse and Cruden (1979) showed that for the Barton standard profile, the relationship between JRC and the root mean square of the first derivative of the profiles Z 2 satisfies Eq. 2, and the correlation coefficient between the two reached 0.9863.
where L is the total length of the joint profile; (x i + 1x i ) is the measurement step width; y i and y i+1 are the coordinates of the y-axis of the ith and i + 1-th joint discrete points; n is the number of discrete points of the joint profile. The joint profile parameters in Table 1 were used. The calculation program was written in MATLAB language to construct joint profiles with five different JRC values. These JRC values were calculated according to Eqs. 2 and 3. Figure 1 shows five random morphology rock joint profiles called J1, J2, J3, J4, and J5. The height of the joint undulation is magnified by a factor of 3 with respect to the joint length. It is found that the JRC of the joints increases with the increase of undulation. The construction of different rock joint profiles provides a basis for studying the shear characteristics of double rough parallel joints.

NUMERICAL MODEL UNDER CONSTANT NORMAL STIFFNESS BOUNDARY CONDITION Continuously Yielding Joint Model
The continuously yielding (CY) model developed by Cundall and Lemos (1990) is expected to simulate the nonlinear behavior of rock joints observed during the direct shear testing. Compared with the conventional Mohr-Coulomb slip model, the essential features of the CY model in shear loading can account for joint shear and normal stiffness dependence of normal stress and nonlinear hardening and softening behavior in the post-peak stage.
The incremental normal stress Δσ n in response to normal loading can be expressed as where K n is joint normal stiffness; Δu n is joint normal displacement increment. The model assumes that the joints are not tensile, and K n is a quantity related to the normal stress.
K n a n σ en n (5) where a n and e n are constants that are experimentally determined.
This model can show irreversible nonlinear behavior from the onset of shearing. The shear stress increment Δτ is expressed as where K s is joint shear stiffness; Δu s is joint shear displacement increment. Factor F is the governing parameter of the tangent modulus, which may be considered as an index of stress path history, and is expressed as  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 819290 where τ is the current shear stress; τ m is failure stress at a given plastic displacement; factor r is intended to restore the elastic stiffness immediately after a load reversal. The initial value of r is 0. r is set to τ/τ m when the direction of the applied shear load is reversed. Therefore, the CY model can represent the hysteresis in cyclic shear loading. τ m is expressed as where φ m is the mobilized friction angle, and it decreases with the increase of shear displacement.
where φ i m is the initial internal friction angle; Δu p s is the plastic displacement increment; φ r is residual internal friction angle; R represents the speed at which the initial friction angle changes to the residual friction angle. The smaller the R, the faster the rate of change.

Numerical Model
The practical implications of CNS boundary conditions are movements of unstable blocks in an underground excavation ( Figure 2). In Figure 2A, the failure mode of the unstable rock block between two parallel joints is shear-slip failure along the joints. The sliding of such block inevitably increases the normal stress, which in turn increases shear resistance. The size of the unstable block is largely controlled by the joint spacing of a particular set of joints. The CNS boundary condition is usually simulated as an elastic spring with normal stiffness k cns , and the value of k cns is externally controlled by the adjacent rock mass across the joint interface (Thirukumaran and Indraratna, 2016). Apart from this boundary effect, the surfaces of natural rock joints in hard rock are always rough, which has an influence on the shear behavior of joints. Based on Figure 2A, the conceptual model of double rough parallel joints subjected to shear stress under CNS conditions is illustrated in Figure 2B. The joint surface is aligned parallel to the shear direction. The compression of the spring indicates joint dilation and results in an increase in normal stress. The normal stress σ n during the shear process varies as a function of initial normal stress σ n0 , spring normal stiffness k cns , and joint dilation δ n . The shear behavior is governed by initial normal stress, joint surface roughness, normal stiffness magnitude, shear-slip displacement, and joint spacing.
To simulate the direct shear test of double rough parallel joints under CNS boundary conditions, the laboratory scale numerical model in UDEC software was established, as shown in Figure 3. The model consists of two blocks where the spring is simulated by the top block to implement the CNS condition, and the rock specimen is   Figure 1 were selected. The rough joints were digitized in AutoCAD software and imported into the numerical model by "crack commands" using the FISH function of UDEC. The two joints with the same JRC values in each numerical model represent one set of joints. They are equally spaced on both sides of the midline in the height direction of the rock block. The CY joint model was used for the joint constitutive model. The normal stiffness corresponding to the spring block is taken as 6 GPa/m, and its stiffness magnitude is mainly determined by the bulk modulus and shear modulus of the spring block. The microparameters of the CY joint model and the macroparameters of the spring block and rock block are shown in Table 2 from the research results of Masivars (2006). The initial normal stress of 3 MPa was applied to the top of the spring block, and the model was allowed to achieve equilibrium. The spring block was fixed in the x-direction on the sides and in both the x-direction and y-direction on the top boundary after the stress was applied. The block located between the spring block and the upper joint was fixed in the x-direction on the vertical sides of the block. The bottom boundary of the lower rock block was fixed in the y-direction. A horizontal velocity of 0.001 m/s was applied on the left wall of the bottom block between the lower joint and the bottom boundary to simulate the shear procedure.

Model Validation
Under the CNS boundary condition, the increment of normal stress is linearly related to the increment of shear dilation displacement, and the ratio is the normal stiffness. To verify that the normal stiffness remains constant during shear, a specimen with JRC = 12.481 and the joint spacing of 5 mm was used as an example to plot the relationship between shear dilation displacement and normal stress, as shown in Figure 4. It can be observed that the normal stress increases with an increase in the shear dilation displacement. The two are well fitted by a linear function with the slope and intercept approximately equal to the set normal stiffness and initial normal stress, respectively. The coefficient of determination R 2 is 0.9998. The results mentioned earlier show that the model can better achieve the CNS boundary condition required for the present study.

Comparison of the Shear Properties of Single and Double Joints
To investigate the influence of CNS condition on the overall shear performance of double joints, a series of numerical CNS direct shear tests of single and double joints were carried out on different rock joint profiles J2, J3, and J4. Three different rock  Haque (1999) used the spring system to model the CNS direct shear boundary and selected the CY joint model in UDEC; the numerical simulation results are in an excellent agreement with the laboratory data to ensure the capability of the numerical method in reproducing the shear mechanical behavior of the physical specimen. Figure 5 gives the comparison of shear stress-displacement relations of two kinds of specimens under CNS boundary conditions. It can be seen that the relation curves are divided into four main stages: the elastic stage, the transitional stage, the post-peak stage, and the shear-hardening stage. At the elastic stage, the shear stress increases linearly with the increase of shear displacement. At the transitional stage, the curve is an arc-shaped line and tends to flatten, indicating that the shear stress increases slowly and the initial shear stress peak τ y appears. At the postpeak stage, the shear stress decreases gradually with shear displacement. At the shear-hardening stage, the shear stress increases slowly with the increase of shear displacement, and the shear-hardening phenomenon is more obvious in a singlejoint specimen than a double-joint. Due to the interaction between the joints, τ y is reduced by 11.43% (JRC = 7.593), 15.49% (JRC = 12.481), and 9.69% (JRC = 15.241) for the double-joint specimens compared with the single-joint, respectively. These findings correspond with the laboratory experiment results of Han (2019). In addition, τ y increases by    83.64 and 87.24%, respectively, for single-and double-joint specimens as JRC values increase from 7.593 to 18.224. Figure 6 shows the comparison of normal-shear displacement relations of the single-and double-joint specimens under CNS boundary conditions. At the initial stage of shear, the shear shrinkage phenomenon occurs in both single-and doublejoint specimens. With the increase of shear displacement, the normal displacement gradually increases from negative to positive, and the phenomenon of shear dilation appears. The greater the joint roughness, the more significant the shear dilation, which is because the increase of JRC leads to the greater undulation angle of joints and the more significant climbing effect. The increase rate in normal displacement decreases with the development of shearing due to the progressive increase in deformation of the joint asperities. Figure 7 shows the variation of normal stress with shear displacement for the single-and double-joint specimens. Compared with the single-joint, the normal displacement and normal stress of the double-joint specimens are smaller. When the shear displacement is 6 mm, the normal displacement of the double-joint specimens decreases by 15.67% (JRC = 7.593), 27.98% (JRC = 12.481), and 21.71% (JRC = 15.241), and the normal stress of the double-joint decreases by 5.73% (JRC = 7.593), 17.04% (JRC = 12.481), and 12.81% (JRC = 15.241).
It is difficult to define a peak point under the CNS condition because test results under certain CNS conditions cannot clearly    show maximum shear stress values. Therefore, Lee et al. (2014) defined the ratio of shear stress to normal stress during shearing of rock joints under CNS boundary conditions as the SRI because it always showed a clear peak point under all conditions. Figure 8 represents the variation of SRI with shear displacement for singleand double-joint specimens. The curve is divided into three phases: pre-peak phase, post-peak phase, and residual phase.
In the pre-peak phase, SRI gradually increases with the increase of shear displacement, and its growth rate gradually decreases. In the post-peak phase, SRI gradually decreases with the increase of shear displacement. In the residual phase, SRI fluctuates in small increments. For both single-and double-joint specimens, the SRI variation curve gradually shifts upward as the JRC increases and its peak surface resistance index (SRI p ) gradually increases. For the same JRC, the SRI of the singlejoint is greater than that of the double-joint, and the shear displacement corresponding to the SRI p of the single-joint is smaller than that of the double-joint. The SRI p values of the single-joint were 3.98% (JRC = 7.593), 12.84% (JRC = 12.481), and 13.61% (JRC = 15.241) lower than those of the double-joint, respectively.

Effect of Joint Spacing on Shear Behavior of Double Joints
To investigate the effect of joint spacing on the shear behavior of double joints, four different spacings of 5, 10, 15, and 20 mm were set for each of the five different joints of J1, J2, J3, J4, and J5, respectively. The variation characteristics of shear stress, normal displacement, normal stress, and SRI of the double-joint were investigated. Figure 9 gives the variation of shear stress with shear displacement for joint J3 with different spacings. τ y gradually increases with the increase of the joint spacing. Figure 10 gives the variation of normal displacement with the shear displacement of J3 at different spacings. The normal displacement gradually increases with the increase of joint spacing, and also, the larger shear displacement leads to the more obvious increase of normal displacement. Figure 11 gives the variation of normal stress with shear displacement at different spacings of J3. The normal stress gradually increases with the increase of the joint spacing. Figure 12 gives the variation of SRI with shear displacement for different spacing of J3, and SRIp gradually increases as the joint spacing increases.
All data under each scenario were recorded to obtain τ y , normal displacement δ v6 , and normal stress σ n6 at the shear displacement of 6 mm and SRI p , as summarized in Table 3.
The effect of joint spacing and JRC on τ y in the double-joint is shown in Figure 13. From Figure 13A, it can be seen that for the same JRC, τ y tends to increase gradually as d increases. As d increases from 5 to 20 mm, τ y increases by 3.45% (JRC = 2.398), 15.51% (JRC = 7.593), 16.04% (JRC = 12.481), 17.69% (JRC = 15.241), and 28.51% (JRC = 18.224), respectively. As the JRC increases, the effect of the joint spacing on τ y becomes greater. As can be seen from Figure 13B, for the same joint spacing, τ y tends to increase linearly as JRC increases, with fitted straight-line slopes of 0.242-0.325 and R 2 = 0.9607-0.9821. Figure 14 shows the variation of the normal displacement δ v6 with the joint spacing d and JRC at the shear displacement of 6 mm. It can be seen that for the same JRC, δ v6 shows a gradual increase with an increase of d.  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 819290 = 18.224). For the same d, δ v6 tends to increase as JRC increases. δ v6 and JRC can be better fitted using a linear function with the slope of 0.0474-0.0586 and R 2 = 0.9742-0.9893. The increased amplitude in δ v6 increases gradually with the increase of d. Figure 15 shows the variation of the normal displacement σ n6 with the joint spacing d and JRC at the shear displacement of 6 mm. It can be seen that for the same JRC, σ n6 shows a gradual increase with an increase of d. When d increases from 5 to 20 mm, the increases of σ n6 are 5.64% (JRC = 2.398), 3.64% (JRC = 7.593), 9.36% (JRC = 12.481), 13.73% (JRC = 15.241), and 14.53% (JRC = 18.224). For the same d, σ n6 tends to increase as JRC increases. σ n6 and JRC can be better fitted using a linear function with a slope of 0.284-0.352 and R 2 = 0.9734-0.9881. The increased amplitude in σ n6 increases gradually with the increase of d.
Han (2019) experimentally investigated the influence of two different joint spacing and JRC values on the mechanical behavior of double-joint specimens during shearing under CNS boundary conditions and found that both joint spacing and JRC have significant effects on peak SRI, and the effect of JRC is greater than that of joint spacing. On this basis, we performed numerical tests of four different joint spacing and five different JRC vales. Figure 16 shows the variation of SRI p with d and JRC. It can be seen that for the same JRC, SRI p tends to increase gradually as d increases. As d increases from 5 to 20 mm, the increases of SRI p are 1.85% (JRC = 2.398), 8.53% (JRC = 7.593), 9.86% (JRC = 12.481), 11.92% (JRC = 15.241), and 16.38% (JRC = 18.224). As JRC increases, the effect of d on SRI p becomes greater. For the same d, SRI p tends to increase with the increase of JRC. SRI p and JRC can be better fitted using a linear function with a slope of 0.0341-0.0437 and R 2 = 0.9884-0.9985. The findings discussed earlier of numerical tests are consistent with the results of Han (2019).

CONCLUSION
(1) Under CNS boundary conditions, shear stress, normal displacement, normal stress, and SRI all increase with the increase of JRC for both single-and double-joint specimens, and these of the single-joint are greater than these of the double-joint.
(2) Under the CNS boundary conditions, for the same JRC, the initial peak shear stress τ y and the SRI p show a gradual increase with the increase of the joint spacing d. When d increases from 5 to 20 mm, the increases of the two parameters are 3.45-28.51% and 1.85-16.38%, respectively, and as the JRC increases, the effect of d on the two parameters is larger. For the same d, the two parameters tend to increase linearly with an increase of JRC. (3) Under the CNS boundary conditions, for the same JRC, the normal displacement and normal stress show a gradual increase with the increase of d. When d increased from 5 to 20 mm, the increase of the two parameters are 8.06-22.97% and 3.64-14.53%, respectively, at the shear displacement of 6 mm. For the same d, with the increase of JRC, the two parameters gradually increase, and the increase amplitudes gradually increase with the increase of d.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material; further inquiries can be directed to the corresponding author.