Modelling of Rock Joints Interface under Cyclic Loading

The most important effects related to monotonic and cyclic response of contact interfaces of brittle materials are analyzed in this paper. The class of the materials is quite wide and includes rock joints, artificial and natural joints, existing cracks in rocks, brittle materials, masonry and other cementitious joints. Rock masses contain discontinuities, such as faults, joints, cracks and other plane of weakness.

The mechanical behavior of rock masses depends not only on the properties of intact structure but also on the properties of the interfaces. The mechanical behavior of such contact interfaces (shear strength), shear and normal deformation dominates the mechanical behavior of in situ rock masses under different loading condition.

Frictional properties of joints are usually investigated by performing shear test under constant normal load CNL^[11] (Fig. 1) and constant normal stiffness CNS condition.

Figure 1

Two types of interfaces can be considered: first as a natural joint, where the surface morphology is described by the joint roughness coefficient JRC^[6,35] (Fig. 2a), and second, where the interface is artificial sawcut material joint and the surface morphology is represented by periodic piecewise linear or curvilinear shape^[11,34] (Fig. 2b, c). The contact interface issue was under investigation both in laboratory tests [10], [16], [19] and from the modeling point of view [2], [7], [12], [17], [29]. Most of the researchers focus on monotonic loads [11], [25], [26], but cyclic loads, where degradation of the contact surface occur [9], [27] can also be found.

Figure 2

In several studies, the surface roughness is considered as a composition of asperities of two different length scale: primary asperities (largest one) and secondary (or subordinated) asperities (Fig. 3). On each contact surface, there are some kind of surface asperities, which can be seen in Fig. 2. Each rock produces distinctive structure during cracking. The resulting asperities correlate with the length of the crack. In the literature, there are many works dealing with the influence of the shape of these asperities on the behavior of the cyclic shear process.^[1,26] For the analyzed issues, the most important influence on tangential stresses tendencies and changes in the contact layer height have primary asperities.^[4,18]

Figure 3

The basic model used in the analysis of soil deformation processes (digging, compaction) is a model of perfectly plastic material that meets the Coulomb or Drucker-Prager condition of plasticity and associated or non-associated flow rules. In this model, the hardening or softening of the material is omitted and unlimited deformability in the state of plastic flow is assumed. Coulomb’s linear condition gives a good estimate for dense materials. However, for loose and susceptible materials, when there is simultaneous compaction and loosening during deformation processes, nonlinear models have found widespread use. In this type of processes, models with a closed yield surface are used (i.e., ellipse). Surfaces with a closed yield criterion, which can describe the effects of compaction and loosening, and deformation of the soil for complex load paths, have been used by many researchers, that is, [3], [5], [24], [25], [27], [28] and [30]. To describe the state of deformation, both associated and non-associated flow rules were used, as well as various concepts of density and deviatoric hardening. To describe the yield surface, elliptical surfaces, super ellipse or irregularly shaped surfaces were used.

In the presented simulation, the elliptical model defined in the space of contact stresses (σ_n, τ_n) was used to describe the material in contact layers. Models with elliptical yield surface were successfully used to describe the behavior of the shear band deformation, which is formed in the rocks [9], [20] and [21]. The evolution of the state of stress and strain in connection with the effects of wear and damage is presented in the paper [19].

Typical dilation and stress response in cyclic shear test is shown in Fig. 4. The shearing displacement is usually accompanied by the dilation of the joint. The character of dilation response depends on the joint type (primary asperity shape). With the change of the sliding direction, the dilation decreases until the shear displacement changes its sign and then increases again. For repeated cyclic shearing at very low pressures, the dilation response is reversible (Fig. 4a), while under high normal pressure, the dilation decreases with increasing number of cycles (Fig. 4b). This effect is related to wear and crushing of asperities and granular fraction formation on the contact surfaces.

Figure 4

The granular layer consists of free particles with different shapes and sizes. In cyclic shear test, the effect of compaction associated with the shear direction changing is observed (Fig. 4c). To describe such a behavior, the phenomenological 2D model of a granular layer is proposed (3D is not considered). This model is based on elliptic yield surface with density and deviatoric hardening described by the rotation of yield surface.

In this work it is assumed that the primary asperities are responsible for interfacial dilation, while the secondary asperities govern the frictional conditions at the inclined contact of primary asperities (Fig. 3).

The rock joint interface is treated as the contact zone with properties different than the rock block with shape described by asperities. During shearing, this zone adopts the properties of granular material. The granular layer with thickness h is sheared as the result of tangential displacement u_t with an amplitude u0 under normal stress load σ_n (Fig. 5).

Figure 5

Material under cyclic shear test is in a certain state of stress and strain in accordance with (1). (1)σ=[σnτn],ε=[εnγn]\sigma = \left[ {\matrix{{\sigma _n } \cr {\tau _n } \cr } } \right],\varepsilon = \left[ {\matrix{{\varepsilon _n } \cr {\gamma _n } } } \right] where strain rate is equal to ε˙n=−u˙nh,γ˙n=u˙th{\dot \varepsilon _n} = - {{{{\dot u}_n}} \over h},{\dot \gamma _n} = {{\dot ut} \over h} .

In the presented work, σ_n, ɛ_n are treated as positive and displacement u_n – as negative; so, increase in layer height has a positive value.

To describe the strain of a layer, it was assumed that the strain rates components are decomposed as elastic and plastic: (2)ε˙=ε˙e+ε˙p\dot \varepsilon = {\dot \varepsilon ^e} + {\dot \varepsilon ^p} or (3)[ε˙nγ˙n]=[ε˙neγ˙ne]+[ε˙npγ˙np]\left[ {\matrix{ {\dot \varepsilon _n } \cr {\dot \gamma _n } } } \right] = \left[ {\matrix{ {\dot \varepsilon _n^e } \cr {\dot \gamma _n^e } } } \right] + \left[ {\matrix{ {\dot \varepsilon _n^p } \cr {\dot \gamma _n^p } } } \right] where ε˙p{\dot \varepsilon ^p} and ε˙e{\dot \varepsilon ^e} are elastic and plastic strain rate component, ε˙ne{\dot \varepsilon_n^e} , γ˙ne{\dot \gamma_n ^e} are normal and shear elastic strain rate components and ε˙np{\dot \varepsilon_n^p} , γ˙np{\dot \gamma_n ^p} are plastic components.

Also, it is assumed that constitutive relations for elastic deformation can be expressed in incremental form: (4)[ε˙neγ˙ne]=Ce[σ˙nτ˙n]=[1/Kd001/Gd][σ˙nτ˙n]\left[ {\matrix{{\dot \varepsilon _n^e } \cr {\dot \gamma _n^e }} } \right] = C^e \left[ {\matrix{{\dot \sigma _n } \cr {\dot \tau _n }} } \right] = \left[ {\matrix{{1/K_d } & 0 \cr 0 & {1/G_d }} } \right]\left[ {\matrix{{\dot \sigma _n } \cr {\dot \tau _n }} } \right] where K_d, G_d are the normal and shear moduli.

Plastic strain rate is described by the associated flow rule (5)ε˙p=λ˙∂F∂σ, λ˙≥0, F≤0, λ˙F=0{\dot \varepsilon ^p} = \dot \lambda {{\partial F} \over {\partial \sigma }},\dot \lambda \ge 0,\,\,\,F \le 0,\,\,\,\dot \lambda F = 0 where F is given by (6).

In this study, an elliptical yield function was assumed, which in the space of contact stresses σ_n, τ_n, can be formulated as: (6)F(σn,τn)=(σn−σ0)2a2+(τn)2b2−1=0F\left( {{\sigma _n},{\tau _n}} \right) = {{\left( {{\sigma _n} - {\sigma _0}} \right)^2} \over {{a^2}}} + {{{{\left( {{\tau _n}} \right)}^2}} \over {{b^2}}} - 1 = 0 where a, b are the sizes of semi-major and semi-minor axes of an ellipse, and σ₀ is the position of the center of an ellipse on the axis of normal stresses.

It was also assumed that the size of the ellipse is related to the location of its center σ₀ in accordance to: (7)a=kσ0a = k{\sigma _0} where k is a material constant.

Along with changes in the density of the material, the size and position of the yield surface changes according to the following formula: (8)σ0=C1[(ρmax−ρminρmax−ρ)C2−1]{\sigma _0} = {C_1}\left[ {{{\left( {{{{\rho _{\max }} - {\rho _{\min }}} \over {{\rho _{\max }} - \rho }}} \right)}^{{C_2}}} - 1} \right] where C₁, C₂ are material constants and ρ_max, ρ_min mean the highest and lowest density of the material.

It was assumed that the critical state line for which plastic flow occurs without volume changes is given by: (9)τ0=A1⌊1−e−A2σ0⌋=b\tau _0 = A_1 \left\lfloor {1 - e^{ - A_2 \sigma _0 } } \right\rfloor = b where A₁, A₂ are material constants and σ₀ is the stress determining the state of compaction.

The relation given above also defines the size of the semi-minor axis.

Examples of elliptical yield surfaces and the critical state line (csl) are shown in Fig. 6a. The figure shows three yield surfaces corresponding to material with densities ρ1 <ρ2 <ρ3. The critical state line separates stress states for which the material is softened by loosening (positive dilation angles) and hardened by compaction (negative dilation angles).

Figure 6

Fig. 6b shows the changes in the position of the center of ellipse σ₀ and the size of the semi-major and semi-minor axes (a, b) as the function of density ρ.

Figures 7, 10, 12, 15 present simulation results for each stage of model building. All cyclic shear tests of the contact layer were performed with an amplitude of u₀ = ±5 mm and normal load σ_n = 2, 10 and 20 MPa. Assumed material parameters are: Young modulus E = 23.13 GPa, Poisson’s ratio υ =0.23, initial density ρ_ini=2.3·10³kg/m³, A₁=200 MPa, A₂=2.887·10⁻²MPa⁻¹, ρ_max = 3·10³kg/m³, ρ_min = 2·10³kg/m³, C₁ = 10 MPa, C₂ = 2.

Figure 7

Fig. 7a shows the changes in tangential stress value as a function of tangential displacement. For low values of normal stress (σ₁) in the initial load cycle, there is a sudden increase in tangential stress and then decrease and stabilization at a much lower level limited by formula (9). As the normal load increases, this effect does not disappear. At high contact pressures (for σ₃), the tangential stress gradually increases rather than decreases. Fig. 7b shows changes in the material dilation layer thickness. For low stresses (σ₁), the material loosens, while for high pressures (σ₃) the compaction of the material starts immediately.

This simple elliptical model is not able to describe the complex nature of rock joint interface, in particular, when changing the load direction. The response stabilizes just after the first load cycle.

During the change of the tangential movement direction, the grain configuration changes at the contact layer. This effect is shown in Fig. 8.

Figure 8

The concept of the description of the configurational rearrangement was presented in the work [31] assuming the description of the contact layer by the Coulomb condition, the non-associated flow rule and rotation of yield surface according to sliding direction. In this model, the idea of changing the position of the yield surface is shown in Fig. 9.

Figure 9

The rearrangement of particles is unequivocal with the change in the location of the yield surface.

The change of the failure surface position is related to the introduction of the rotation angle θ (Fig. 9). With the increase of the u_t displacement, the surface position is changed by the angle θ to the maximum value ± θ_max depending on the deformation direction ut described by the parameter s(11). Rate of θ is described by (10). (10)θ˙=η˙kt(sθmax−θ)\dot \theta = \dot \eta {k_t}\left( {s{\theta _{\max }} - \theta } \right) (11)s=γ˙np|γ˙np|=±1s = {{\dot \gamma _n^p} \over {\left| {\dot \gamma _n^p} \right|}} = \pm 1 (12)η˙=|γ˙np|\dot \eta = \left| {\dot \gamma _n^p} \right| Where k_t [-] is rotation parameter and γ˙n\dot \gamma _n is shear strain rate.

In the local coordinate system, the state of stress and strain is described by (13). (13)σ′=[σn′τn′], ε˙′=[ε˙n′γ˙n′]\sigma^\prime = \left[ {\matrix{{\sigma _n^\prime } \cr {\tau _n^\prime } } } \right],\;\;\;\dot \varepsilon^\prime = \left[ {\matrix{{\dot \varepsilon _n^\prime } \cr {\dot \gamma _n^\prime } } } \right] The relationships between global and local (rotated) coordinate systems are given by (14). As shown before, in the rotated ellipse system also, plastic strains are described by the associated flow rule (15). (14)σ′=Qσ,σ=QTσ′ε˙′=Qε˙,ε˙=QTε˙′\matrix{{\sigma ' = Q\sigma ,} \hfill & {\sigma = {Q^T}\sigma '} \hfill \cr {\dot \varepsilon ' = Q\dot \varepsilon ,} \hfill & {\dot \varepsilon = {Q^T}\dot \varepsilon '}} (15)ε˙′p=λ˙∂F′∂σ′ ε˙p=QTε˙′p{{\dot \varepsilon^\prime}^p} = \dot \lambda {{\partial F'} \over {\partial \sigma '}}\,\,\,\,\,{\dot \varepsilon ^p} = {Q^T}{\dot \varepsilon '^p} where Q is the rotation matrix connected with angle θ.

The adaptation of the yield surface rotation in cyclic shear test is shown in Fig. 10. Similar as before, the simulations were performed for three different values of normal stress at a constant initial density and θ_max = 25°.

Figure 10

The changes in tangent stress are similar to simulations for a model that does not include rotation of the failure surface. Considering yield surface position changes by adding the angle of rotation θ, results in qualitative changes in the dilation of the layer at load reversal. When the direction of shear changes, the initial compaction of the material occurs, followed by the relaxation in the further stage to the level dependent on the state of normal stress σ_n. In subsequent load cycles, the response is the same.

In the model discussed here, proposed shape of the primary asperity, similar as in [31], is given by formula: (16)ua=f(ut)=δAtanh(ut2tan2ψ0+g02−g0δA){u_a} = f({u_t}) = {\delta _A}\tanh \left( {{{\sqrt {u_t^2{{\tan }^2}{\psi _0} + g_0^2 - {g_0}} } \over {{\delta _A}}}} \right) where δ_A asperity height, g₀ and ψ₀ are shape parameters.

The influence of asperity shape on the parameter g₀ is shown in Fig. 11. Higher the value of the parameter g₀, the more gentle are the asperity slopes.

Figure 11

Simulations in Fig. 12 show the results with respect to the shape of asperity together with the elliptic yield condition and the change in θ angle value of the failure surface position. The total rotation angle θ_max of the elliptic yield condition (17) is the sum of the angle θ and the angle ψ_as resulting from the shape of the primary asperities (18). (17)θmax=θ+ψas{\theta _{\max }} = \theta + {\psi _{as}} (18)ψas=duas(ut)dut{\psi _{as}} = {{d{u_{as}}\left( {{u_t}} \right)} \over {d{u_t}}} The results show specific changes in the tangential stress state (Fig. 12a) and their step change when passing through the starting point for which u_t = 0. The changes in the dilation of the contact layer (Fig. 12b) are the sum of changes in the asperity height and changes in the behavior of the newly created third body granular layer. In addition, Fig. 12c shows the asperity shape assumed for calculations, which is symmetrical with respect to “0” assuming model parameters: δ_a = 2 mm, g₀ = 0.03 mm, ψ₀ = 10°.

Figure 12

During the cyclic shearing of rocks, formation of the granulate fraction due to abrasive wear is observed. Wear occurs both by changing the shape of the asperities and by friction wear of the contact layer.

Figure 13 shows the generation of third body granular layer due to the shear process.

Figure 13

The progressive abrasive wear considerably contributes to the formation and development of this structure. Therefore, a certain form of friction wear was introduced into the model: (19)w˙=|u˙t||τn|kw(hmax−h)\dot w = \left| {{{\dot u}_t}} \right|\left| {{\tau _n}} \right|{k_w}\left( {{h_{\max }} - h} \right) where k_w is a wear parameter and h_max is maximal thickness of the granular layer.

Changes in layer height are also related to the process of wear (attrition) of primary asperities. The wear function is formulated as: (20)g=g(κ)=gf⌊eκκ0−1⌋g = g(\kappa ) = g_f \left\lfloor {e^{{\kappa \mathord{\left/ {\vphantom {\kappa {\kappa _0 }}} \right. } {\kappa _0 }}} - 1} \right\rfloor where κ is accumulated frictional work increment given by: (21)κ˙=|τn||u˙t|\dot \kappa = \left| {{\tau _n}} \right|\left| {{{\dot u}_t}} \right| Asperity wear is described by: (22)was=δAtanh(g2+g02−g0δA){w_{as}} = {\delta _A}\tanh \left( {{{\sqrt {{g^2} + g_0^2} - {g_0}} \over {{\delta _A}}}} \right) Using the above relationship, the function describing the shape of asperity with respect to their wear was modified as: (23)uas=f(ut)=δA[tanh(ut2tan2ψ0+g2+g02−g0δA)−tanh(g2+g02−g0δA)]{u_{as}} = f\left( {{u_t}} \right) = {\delta _A}\left[ {\tanh \left( {{{\sqrt {u_t^2{{\tan }^2}{\psi _0} + {g^2} + g_0^2} - {g_0}} \over {{\delta _A}}}} \right) - \tanh \left( {{{\sqrt {{g^2} + g_0^2} - {g_0}} \over {{\delta _A}}}} \right)} \right] Figure 14 shows a gradual change in the shape of asperities as a function of the progressive wear described by (20). The higher value of wear function g as the load cycles increase, the greater degradation of the asperity shape.

Figure 14

Figure 15 shows simulation results for the material subjected to cyclic shear described in the model with respect to assumptions presented so far (for k_w = 0.1·10⁻⁴mm²/J, h_ini = 5 mm, h_max = 7 mm, δ_a = 2 mm, g₀ = 0.03 mm, ψ₀=10°, g_f = 4 mm, κ₀ = 3 J/mm²(20)). Variations of shear stress, dilatancy and rotation angle are clearly visible. Similar to the non-wear plots (Fig. 12), sudden changes in the state of shear stress when passing through the point “0” can be observed (Fig. 15a). These step changes are smaller as the number of load cycles increases. Considering the variations of dilation (Fig. 15b), it can be seen that the contact layer height is gradually reduced. As a result of progressive wear, the plot is not symmetrical with respect to the point “0”, so that there is no single material response to the change of load direction.

Figure 15

Figure 15c shows the failure surface rotation angle θ changes. Also, in this case, due to wear, this parameter value changes are seen as well.

In addition, Fig. 16 shows graphs that complement the response of the material for applied normal stress σ_n = 10 MPa. Figure. 16a contains the contact layer height increase from the assumed initial (5 mm) to the final value, that corresponds to the end of the process. Fig. 16b shows changes in the shape of asperities as the number of load cycles increases. This shape is subject to changes due to progressive wear. On the other hand, Fig. 16c shows the characteristic of changes in the behavior of the third body granular layer.

Figure 16

It was decided to carry out a comparison of the proposed model with the experimental data (Fig. 17). The results of research presented by Kamonphet et al. in the paper [10] were selected as the comparative data (Fig. 17b and Fig. 18b). In the experiment, the rock sample (Tak granite) was placed in the apparatus for conducting shear tests, loaded with a normal stress σ_n = 0.5 MPa and σ_n = 4 MPa and forced to move with the displacement amplitude u₀ = ±5 mm. Test was performed on natural joint. Simulations were carried out with the assumption of uniform, symmetrical with respect to “0”, shape of the asperity. The following model parameters were assumed: h_ini = 0.1 mm, ρ_ini = 2.1·103 kg/m³, ρ_max = 2.71·10³kg/m³, ρ_min = 2·10³kg/m³, A₁ = 80 MPa, A₂ = 0.01 MPa⁻¹, C₁ = 1.1 MPa, C₂ = 3, k_t = 1, k_w = 0 mm²/J, h_max = 10 mm, δ_a = 0.1 mm, g₀ = 0.0 mm, ψ₀ = 12.5°, g_f = 1 mm, θ_max = 15°, κ₀ = 0.25 J/mm². These parameters were chosen so as to best fit the real changes.

Figure 17

Figure 18

When comparing the characteristics with each other, good qualitative agreement obtained by simulation and experiment can be noticed.

In case of tangential stress, a step change when passing through the point “0” in the first load cycle can be seen. When the number of cycles increases, these steps are gradually flattened and the changes in contact layer height are getting smaller. The adoption of assumptions regarding the abrasive wear of both the asperities and the creation of the third body granular layer were the most appropriate. This can be observed comparing the asymmetrical (relative to “0”) dilation diagrams obtained in the simulation and the experiment.

Differences in the quantitative agreement occurred, especially in the first shear stage. This is due to the assumption of uniform, symmetrical with respect to “0”, shape of the asperity. As can be seen in the above figures, the model is able to describe the nature of changes in both tangential stress and dilation with wear, but not the first peak, which may result from fitting together two parts of the tested rock or from sudden destruction in the first stage of surface asperities. Apart from the first peak, the rest of the shear process is modeled correctly.

For cyclic shearing under high normal pressure, the dilation decreases with increasing number of cycles. This effect is due to wear and degradation of asperities and formation of granular fraction on the contact surfaces. In the cyclic test, the effect of the compaction when changing the shear direction is observed.

A large number of laboratory tests are available regarding the influence of the shape of the surface asperities’ profile on the cyclic shear process, whereas insufficiency is based on numerical modeling. Most of the studies on the shear behavior of rock joints are concerning the pre-peak and peak shear strength in monotonic loading. Presented model provides the ability to study the post-peak behavior of rock joints in cyclic test. However, further studies are required to investigate the correct mapping of real changes, especially in the first shear stage. This will be achieved by applying the damage function.

This paper analyzes the deformation of brittle materials on contact surfaces subjected to cyclic loads. It was shown, that some of the effects observed in experiments may result from more than one mechanism. Simple conditions were proposed providing description of these mechanisms.

A comprehensive model containing: –

elliptic yield condition with density hardening/softening with associated flow rule,

This analytical model should be treated as a preliminary, covering the effects of degradation of surface asperities.