CJRM (Figure 1A) contains two sets of joints (joint sets 1 and 2), which cut the rock mass into columns. The mechanical behaviors of such rock mass are determined by the joints and rock columns (Wang and Huang, 2014; Chang et al., 2019; Do and Wu, 2020). In the enhanced model, the assumed failures of CJRM include three categories: 1) the shear and tensile failures on the weak planes of joints, 2) shear and tensile failures of the rock columns, and 3) bending failure along the rock column axes. In addition, the assumed failure sequences between the weak planes and rock columns are such that 4) the weak planes fail before the rock columns, and 5) the weak plane at the greater distance from the stress state point to the corresponding strength envelope fails first if both weak planes fail in the same iteration process.

The reasons for making such assumptions are as follows: firstly, in the cross-section plane of the rock columns (Figure 1B), the rock mass is blocky so that the strength and deformability are mainly controlled by the joints (Singh et al., 2002; Singh and Rao, 2005). Assumption 1) ensures that material failures in such plane occur only on the weak planes of joints and that plastic deformation in the cross-section plane is mainly controlled by sliding deformation on the weak planes and tensile deformation perpendicular thereto. Secondly, in the axial direction of rock columns, the failures of the rock mass mainly occur in the rock columns. The failure modes include shear and tensile failures of the rock columns and bending failure (Adhikary and Guo, 2002; Neff et al., 2008) along the rock columns, as per assumptions 2) and 3), respectively. Thirdly, studies (Jiang et al., 2006; Hatzor et al., 2015) have shown that jointed rock masses usually fail on the joints first, followed by failure of the intact rocks, so assumption 4) guarantees a reasonable sequence of failure events between joints and rock columns.

Similar to the principle of U-J model, the enhanced model is a completely equivalent continuum model (Mühlhaus, 1993; Adhikary and Dyskin, 1998; Sitharam et al., 2001) that the influence of the joints is smeared into the continuum description of the rock mass, so the distributions (such as the spacing and length) of the joints, except the orientation, are not necessarily defined explicitly (Agharazi et al., 2012; Wang and Huang, 2014; Zhou et al., 2021). In the enhanced model, the local coordinate system is redefined based on the orientations of the two sets of joints the better to describe the mechanical properties of CJRM.

One of the rock columns is taken as an example to illustrate the method for establishing the local coordinate system (Figures 1B,C). First, the axial direction of the rock column is defined as the z ′ -axis; second, the direction perpendicular to z ′ -axis and in the plane of joint set 1 is defined as the x ′ -axis, and the direction perpendicular to z ′ -axis and in the plane of joint set 2 is defined as the y ′ -axis. Since the orientations of the two joint sets are not necessarily perpendicular, the angle between the x ′ -axis and y ′ -axis is between 0 and 180° on the cross-section of the rock column (Figure 1D). The directions of the local coordinates x ′ -, y ′ -, and z ′ -axes are determined as described below.

Vectors J 1 and J 2 are defined as direction vectors normal to the weak planes of joint sets 1 and 2, respectively, and vectors x ′ , y ′ , and z ′ are defined as the positive direction vectors of the local coordinate x ′ -, y ′ -, and z ′ -axes, respectively. Based on the perpendicularity relationships between the five vectors, the directions of the local coordinate axes are given by z ′ : J 1 ⋅ z ′ = 0 J 2 ⋅ z ′ = 0 | z ′ | = 1 , x ′ : J 1 ⋅ x ′ = 0 z ′ ⋅ x ′ = 0 | x ′ | = 1 , y ′ : J 2 ⋅ y ′ = 0 z ′ ⋅ y ′ = 0 | y ′ | = 1 , (1) where J 1 and J 2 are calculated from the dip angles ( D i p 1 and D i p 2 ) and dip directions ( D D 1 and D D 2 ) of joint sets 1 and 2, respectively. J 1 = ( sin D D 1 , cos D D 1 , cot D i p 1 ) , J 2 = ( sin D D 2 , cos D D 2 , cot D i p 2 ) . (2)

Matrix [ C ] is defined as the rotation tensor between the local and the global coordinate systems, and then the transformation of stress components between the two coordinate systems (Itasca Consulting Group, 2011) is given by [ σ ] ′ = [ C ] T [ σ ] [ C ] [ σ ] = [ C ] [ σ ] ′ [ C ] T , (3) where [ C ] is [ C ] = [ cos ( x ′ , x ) cos ( x ′ , y ) cos ( x ′ , z ) cos ( y ′ , x ) cos ( y ′ , y ) cos ( y ′ , z ) cos ( z ′ , x ) cos ( z ′ , y ) cos ( z ′ , z ) ] . (4)

Weak planes of joints in the enhanced model are ubiquitous in the rock mass based on the theory of the equivalent continuum model (Mühlhaus, 1993; Adhikary and Dyskin, 1998; Sitharam et al., 2001). The stress and strain components at any point in the rock mass can be expressed as shown in Figure 2 without consideration of the spacing and length of the joints in the local coordinate system. On the weak planes of joint set 1, shear stress τ 1 ′ and shear strain γ 1 ′ have the following forms: τ 1 ′ = σ 2 ′ 1 ′ 2 + σ 3 ′ 1 ′ 2 , γ 1 ′ = ϵ 2 ′ 1 ′ 2 + ϵ 3 ′ 1 ′ 2 . (5)

The equivalent forms are readily obtained for joint set 2: τ 2 ′ = σ 1 ′ 2 ′ 2 + σ 3 ′ 2 ′ 2 , γ 2 ′ = ϵ 1 ′ 2 ′ 2 + ϵ 3 ′ 2 ′ 2 . (6)

On the cross-sections of the rock columns, the similar forms are τ 3 ′ = σ 1 ′ 3 ′ 2 + σ 2 ′ 3 ′ 2 , γ 3 ′ = ϵ 1 ′ 3 ′ 2 + ϵ 2 ′ 3 ′ 2 . (7)

In what follows, the failures of rock mass are detected by the six components of the generalized stress vector: σ 1 ′ 1 ′ , σ 2 ′ 2 ′ , σ 3 ′ 3 ′ , τ 1 ′ , τ 2 ′ , and τ 3 ′ whose components of the corresponding generalized strain vector are ϵ 1 ′ 1 ′ , ϵ 2 ′ 2 ′ , ϵ 3 ′ 3 ′ , γ 1 ′ , γ 2 ′ , and γ 3 ′ .

The failure criterion for the weak planes of joints used in the enhanced model is a composite M–C criterion with a tension cut-off expressed in terms of σ i ′ i ′ and τ i ′ (where i = 1 ,   2 ,   3 ), as illustrated in Figure 3.

When shear failures occur on the weak planes of joints, the stress states calculated using the initial elastic estimate are located in Domain 2 (Itasca Consulting Group, 2011). The corresponding failure criteria meet the following conditions:

For joint set 1, f J 1 s = τ 1 ′ + σ 1 ′ 1 ′ tan ϕ J 1 − c J 1 > 0 , h J 1 = τ − τ J 1 P − a J 1 P ( σ 1 ′ 1 ′ − σ J 1 t ) > 0. (8)

For joint set 2, f J 2 s = τ 2 ′ + σ 2 ′ 2 ′ tan ϕ J 2 − c J 2 > 0 , h J 2 = τ − τ J 2 P − a J 2 P ( σ 2 ′ 2 ′ − σ J 2 t ) > 0 , (9) and τ J 1 P = c J 1 − σ J 1 t tan ϕ J 1 , a J 1 P = 1 + tan ϕ J 1 2 − tan ϕ J 1 , τ J 2 P = c J 2 − σ J 2 t tan ϕ J 2 , a J 2 P = 1 + tan ϕ J 2 2 − tan ϕ J 2 , (10) where c J 1 , c J 2 , ϕ J 1 , ϕ J 2 , σ J 1 t , and σ J 2 t are the cohesion, angle of internal friction, and tensile strength of the weak planes of joints, respectively. For weak planes with non-zero internal friction angles, the maximum values of the tensile strengths are given by σ J 1 t ≤ c J 1 tan ϕ J 1 σ J 2 t ≤ c J 2 tan ϕ J 2 . (11)

The potential functions g J 1 s and g J 2 s corresponding to a non-associated flow rule have the following forms:

For joint set 1, g J 1 s = τ + σ 1 ′ 1 ′ t a n ψ J 1 . (12)

For joint set 2, g J 2 s = τ + σ 2 ′ 2 ′ t a n ψ J 2 , (13) where ψ J 1 and ψ J 2 are the angles of dilation of the weak planes.

When tensile failures occur on the weak planes of joints, the stress states calculated from the elastic estimate are located in Domain 3. The corresponding failure criteria meet the following conditions:

For joint set 1, f J 1 t = σ 1 ′ 1 ′ − σ J 1 t > 0 , h J 1 = τ − τ J 1 P − a J 1 P ( σ 1 ′ 1 ′ − σ J 1 t ) ≤ 0. (14)

For joint set 2, f J 2 t = σ 2 ′ 2 ′ − σ J 2 t > 0 , h J 2 = τ − τ J 2 P − a J 2 P ( σ 2 ′ 2 ′ − σ J 2 t ) ≤ 0. (15)

The potential functions g J 1 t and g J 2 t corresponding to an associated flow rule have the following forms:

For joint set 1, g J 1 t = σ 1 ′ 1 ′ . (16)

For joint set 2, g J 2 t = σ 2 ′ 2 ′ . (17)

When no failures occur on the weak planes of joints, the stress states calculated using the elastic estimate are located in Domain 1 and are considered as the final incremental stress without any correction by way of potential functions. The corresponding failure criteria meet the following conditions:

For joint set 1, f J 1 s = τ 1 ′ + σ 1 ′ 1 ′ tan ϕ J 1 − c J 1 ≤ 0 , f J 1 t = σ 1 ′ 1 ′ − σ J 1 t ≤ 0. (18)

For joint set 2, f J 2 s = τ 2 ′ + σ 2 ′ 2 ′ tan ϕ J 2 − c J 2 ≤ 0 , f J 2 t = σ 2 ′ 2 ′ − σ J 2 t ≤ 0. (19)

In addition, it should be noted that both of the weak planes may fail in the same iteration process under a certain stress state in the rock mass. The choice of potential functions for stress correction is determined by the failure sequence and affects the corrected stress state. In the enhanced model, the failure sequence is judged by the distances, d J i s and d J i t , from stress state points A and B to the corresponding envelopes, f J i s and f J i t , as illustrated in Figure 3. The greater the distance (either d J i s or d J i t ), the greater the amount by which the stress calculated by the elastic estimate exceeds the corresponding failure envelope; so, the weak plane with the greater value of d J i s or d J i t is more prone to failure than the other and it fails first. The distance functions have the forms: d J i s = f J i s 1 + t a n 2 ϕ J i , d J i t = f J i t . (20)

The failure mode of the rock columns involve shear and tensile failures of the columns and bending failure along the columns as assumed (Adhikary and Guo, 2002; Neff et al., 2008). When shear and tensile failures occur, a similar procedure to that used in M–C constitutive model (Itasca Consulting Group, 2011) can be used to calculate the incremental stress and strain. This section mainly discusses the bending failure of the rock columns.

Figure 4 illustrates the bending failure of a rock column in the local coordinate system. When a bending moment acts on the rock column, as shown in Figure 4, the rock column is subjected to compressive stress on the upper side and tensile stress on its lower side. With increasing bending moment, the tensile stress on the lower side increases until it exceeds the tensile strength, resulting in a bending failure on the lower side (Agharazi et al., 2012; Sainsbury and Sainsbury, 2017). However, the tensile stress is difficult to calculate accurately, as joints are smeared into the continuum description of the rock mass by an equivalent continuum method and have no boundaries and thickness (Agharazi et al., 2012; Wang and Huang, 2014; Zhou et al., 2021). Thus, the normal stress on the column cross-section σ ′ 3 ′ 3 ′ is used as an approximate replacement for the tensile stress. The bending failure criterion has the form: f b = σ 3 ′ 3 ′ ′ − σ 0 b ≥ 0 , (21) where σ 0 b is the bending strength of rock columns. Generally, the bending strength is less than the tensile strength although it has similar forms of failure criterion in the z ′ -direction. The potential function g b corresponding to an associated flow rule has the form: g b = σ ′ 3 ′ 3 ′ . (22)

Numerical triaxial compression tests were conducted to verify the enhanced model. The samples used in the tests were cylinders with a diameter of 1 m and a height of 2 m. The compressive strength, plastic zones, and deformation of the samples were analyzed: contrasting simulations were conducted using the enhanced model and U-J model. The material parameters used in these simulations are listed in Table 1. The confining pressure used in all such simulated tests is 1 MPa.

The difference of the triaxial compression test samples between use of U-J model and the enhanced model mainly lies in the consideration of the number of joint sets and local coordinate system. For the sample simulated using U-J model, only one set of joints is considered at a time (Figure 5A); the three local coordinate axes are mutually orthogonal and meet the right-handed spiral criterion. By contrast, two sets of joints are considered at a time for the sample simulated using the enhanced model (Figure 5B); the included angle of the local coordinate axes x ′ and y ′ is determined by the occurrence of joints, which ranges from 90° to 270°, and the coordinate axis z ′ is perpendicular to the plane formed by axes x ′ and y ′ .

During the numerical simulations, the compression stress and confining pressure acting on the samples are initially 0.1 MPa and gradually increased with an increment of 0.1 MPa. The confining pressure stops increasing when it reaches the target value of 1 MPa and then stays the same. The compression stress keeps increasing until the calculation no longer converges. The total compression stress of the final converged solution is taken as the compressive strength of the samples (Neff et al., 2008; Chang et al., 2019).

The case study involves the underground excavation simulation of a pumped storage power station in China. The underground caverns are located in a medium-coarse-grained granite (Figure 7A). Two conjugate sets of joints, J1 and J2, occur therein (Figure 7B). The axial directions of the powerhouse and transformer chamber are parallel and lie along the north-south direction (Figure 7C). The in situ stresses in the rock mass are such that 1) the horizontal major principal stress σ 1 ranges from 12 to 18 MPa with the east-west direction; 2) the horizontal minor principal stress σ 3 ranges from 7 to 11 MPa with the north-south direction; and 3) the vertical principal stress ranges from 6.6 to 9.3 MPa. The dimensions and excavation sequence of the underground caverns are shown in Figure 7D. The main support measures include 1) mortar anchors with lengths of 6 or 9 m at 1.5 m × 1.5 m grid spacings, 2) prestressed anchor cables with a length of 25 m at 4.5 m × 4.5 m spacings, and 3) steel fibreconcrete liners with a thickness of 150 mm sprayed onto the surface of the cavern walls. The mechanical properties of the rock mass are listed in Table 2 and those of the support structures are listed in Table 3.

The main engineering problems encountered in the underground caverns during the excavation are as follows: 1) with the downward excavation of the powerhouse and transformer chamber, the stress field of the surrounding rock mass evolves, resulting in joint slippage and cracking of the liners (Figure 8A); 2) after the excavation of IV step is finished, the total deformation of the surrounding rock mass with a maximum value of 100 mm is much greater than the designed value of 50–60 mm. The monitored total deformations on some main cross-sections of the powerhouse are shown in Figure 8B. Compared with other underground engineering with similar lithology, in situ stress, size of excavation, and other conditions (Zhu et al., 2008; Zhu et al., 2010), the most prominent geological feature of the present engineering is that the conjugate joints in the rock mass are relatively well-developed. Therefore, it can be inferred that one of the possible factors inducing the aforementioned engineering problems is the strength reduction of the rock mass due to the conjugate joints, making it necessary to use the enhanced model to study the plastic zones and deformation characteristics of the surrounding rock mass.

In the excavation simulation of the underground caverns, the excavation sequence is I through VI, as shown in Figure 7D. The anchors and cables are simulated by cable elements and the liners are simulated by shell elements (Itasca Consulting Group, 2011). The parameter settings of the anchors, cables, and liners are consistent with those listed in Table 3. The enhanced model, as well as M–C model and U-J model used as contrasting models, are each used to calculate the distribution of plastic zones and deformation of the surrounding rock mass. The rock mass properties used in the simulation are those listed in Table 2.

M–C model and U-J model are selected as contrasting models for the enhanced model due to both the similarities and differences between them. The similarities lie in the failure criteria and the potential functions for joints and intact rocks, i.e., they are all derived from M–C yield criterion. The differences lie in the consideration of the numbers of joint sets and the enhancements in the local coordinate system, failure modes, and failure sequences of the enhanced model. In terms of the numbers of joint sets, neither joint set J1 nor J2 can be considered in M–C model, either J1 or J2 can be considered in U-J model, but both J1 and J2 can be considered in the enhanced model. Thus, by comparing the results calculated using the three models, the effects of joints on strength and deformation behaviors of the surrounding rock mass and the responses of the enhanced model in underground excavation simulations can be determined.

The distributions of plastic zones calculated by the three models are presented in Figure 9: when U-J model is used to simulate either joint set J1 or J2, the distributions of plastic zones are mainly concentrated on the upper left and lower right corners of the powerhouse (when considering J1) or on the lower left and upper right corners (when considering J2). By contrast, when the enhanced model is used to simulate both joint sets J1 and J2, the plastic zones are distributed all around the powerhouse. The notable difference in the distributions of plastic zones indicates that the plastic zones are significantly affected by the joints and it is necessary to fully consider the effects of the two sets of joints in the simulation of excavation-induced failure in CJRM.

The volumes of plastic zones calculated using the enhanced model are significantly larger than those calculated using the U-J model (Figure 9): the difference in the volumes is quantified by comparing the normalized volumes calculated using the three models (Figures 10A,B). Compared with the plastic zone volumes calculated using M–C model considering no joints, those calculated using U-J model considering either joint set J1 or J2 are almost 1.8 times greater, and those calculated using the enhanced model considering both joint sets J1 and J2 are 2.5 times greater (Figure 10A). As the excavation proceeds, the growth rates of plastic zone volumes calculated using each of the three models become significantly different (Figure 10B). From the first excavation step to the final excavation step, the plastic zone volumes increase 1–2.6 times when using M–C model, 1.8–5 times when using U–J model, and as many as 2.2–7 times when using the enhanced model.

The proportional occurrence of the failure types of the rock mass is different when using different models. As shown in Figures 10C,D, when M–C model is used, the failure of the rock mass is mainly tensile, and the proportions of shear failure and tensile failure are 20 and 80%, respectively. On the contrary, when the U-J model or the enhanced model is used, the rock mass is mainly subjected to shear failure. The proportions of shear failure and tensile failure calculated using U-J model are 60 and 40%, respectively, while the proportions calculated using the enhanced model are 90 and 10%, respectively.

The displacement contours (Figure 11) calculated using each of the three models are similar in terms of the deformation mode: the deformation of the surrounding rock mass develops towards the opening and the largest deformation occurs on the sidewalls of the powerhouse upon completion of the excavation. The deformation characteristics are consistent with those reported elsewhere (Jia and Tang, 2008; Hatzor et al., 2015), indicating the correctness of the enhanced model in terms of the deformation response. However, there are some differences among the deformation magnitudes calculated using the three models. Compared with the displacement magnitude contour calculated using M–C model (Figure 11A), those calculated using U-J model (Figures 11B,C) are increased on the sidewalls (especially where the arrows indicate in the figures) and the magnitude calculated using the enhanced model (Figure 11D) is increased on almost all of the sidewalls (especially where also found to have been increased when using U-J model, as the arrows indicate in the figure).

For further study, the magnitudes of the total displacements calculated using the three models at each excavation step are compared (Figure 12). The total displacements calculated using the three models contain both of the elastic and plastic components. The elastic components are almost the same since they are calculated under the same stress condition, but the plastic components are different since they are derived from plastic flows of different plastic zones. As shown in Figure 12, the displacement magnitudes of the monitored points calculated using the enhanced model are generally similar to, but slightly larger than, those calculated using M–C model and U-J model. This indicates that the total displacements are mainly controlled by the elastic components and the plastic components which are derived from plastic flows of plastic zones entailing only minor corrections to the total displacements.