DEM modelling of a jointed rock beam with emphasis on interface properties

This paper compares analyses performed via the distinct element method (DEM), employing rigid blocks and compliant joints, with results using finite-difference software (FLAC) obtained previously by other researchers. The paper then examines the capability of the rigid-block DEM at modelling joints realistically, with emphasis on the moment transfer between blocks. The line of thrust from this analysis was found to fit well with the well-established uniform catenary curve and the parabola, which has been used extensively in the rock engineering literature. This is an important verification exercise that is still lacking in the literature, especially for the rigid-block DEM. Finally, a comparison is made between the DEM and experimental work carried out previously by other researchers. The previously reported laboratory data were reinterpreted to derive more accurate contact laws in both normal and shear directions. A strain hardening or continuously yielding model was adopted in the latter. The calibratio...


INTRODUCTION
To predict deformations within a jointed rock mass, it is necessary to adopt accurate contact models for the rock joints.However, in the literature there is still a lack of validated approaches for jointed rock masses with more than one joint.The merits of adopting an accurate contact model and a reliable contact detection algorithm in the distinct element method (DEM) are demonstrated by simulating a well-instrumented laboratory-scale experiment on a singlelayer jointed beam (called a voussoir beam) carried out by Talesnick et al. (2007) in a centrifuge.The approach illustrated here can be applied to several engineering problems in rock mechanics such as tunnelling (Boon et al., 2014a) and rock slope stability (Boon et al., 2014b).
The influence of jointing has been studied by adopting an equivalent rock mass elastic modulus (Diederichs & Kaiser, 1999a) or through semi-analytical solutions based on empirical databases developed from numerical simulations (Nomikos et al., 2002;Tsesarsky, 2012).To account for the influence of jointing more directly, several investigators have used discontinuum analysis (i.e. the DEM or discontinuous deformation analysis (DDA) (Hatzor & Benary, 1998;Tsesarsky & Hatzor, 2006;Tsesarsky & Talesnick, 2007;Alejano et al., 2008;Bakun-Mazor et al., 2009;Barla et al., 2010;Hatzor et al., 2010) but, until now, these studies have been limited to simple contact models between rock joints.Talesnick et al. (2007) measured the horizontal thrust and the vertical deflection of a voussoir beam subject to increasing levels of gravitational acceleration.Then, Tsesarsky & Talesnick (2007) ran numerical simulations of the beam using the finite-difference method (FDM) and DDA.They found that, for both FDM and DDA, the horizontal thrusts and inter-block slip were underpredicted and the block rotations were uniform compared with the experimental measurements in which rotation was greatest at the abutment blocks and least at the midspan blocks.Furthermore, deflections were underpredicted (the smallest measured deflection profile among the three tests carried out by Talesnick et al. (2007) was used as a benchmark case by Tsesarsky & Talesnick (2007)).
This paper presents DEM simulations that allowed better estimates of the deflections for these centrifuge tests to be obtained.The open-source DEM academic code YADE (Kozicki & Donzé, 2008) was employed in this study, together with a new contact detection algorithm recently proposed by Boon et al. (2012) based on linear programming.Second-order conic programming could also be used (Boon et al., 2013).The robustness of using the DEM to model problems involving moment transfers between nondeformable blocks with compliant contacts is also investigated.

NUMERICAL SETUP
The centrifuge model of Talesnick et al. (2007) consisted of six equal gypsum blocks flanked by two fixed larger abutment blocks.Two opposing sides of each block were fitted with standard abrasive paper.Two prism dimensions were used, 46 mm × 46 mm × 46 mm and 46 mm × 46 mm × 23 mm, herein referred to as the full-block and half-block beam respectively.The properties of the blocks are listed in Table 1.An initial horizontal force of approximately 60 N was applied to the blocks at 1g. Gravity was gradually increased up to 40g for the full-block beam and 90g for the half-block beam.
The same loading sequence was applied in the current DEM simulations.The experimentally derived contact laws reported by Tsesarsky & Talesnick (2007) were also used for the sake of comparison with their numerical results.Concerning the stress-dependent stiffness in the normal direction, a contact model with stiffness being a linear function of the exchanged normal stresses best fits the experimental data.This model is expressed as (Tsesarsky & Talesnick, 2007) where σ n is the contact normal stress, u n is the contact normal displacement, and k n and k s are the contact normal and shear stiffness, respectively.

NEW CONTACT MODEL FROM REINTERPRETATION OF THE EXPERIMENTAL MEASUREMENTS OF TALESNICK ET AL. (2007)
Expressing the contact model as a stress-displacement function is convenient for DEM calculations since the overlap distance between blocks is the input for the contact model.Therefore, integrating equation (1) yields To improve the predictions, the experimental data of Talesnick et al. (2007) and Talesnick (2007) were reinterpreted to derive more accurate normal and shear contact models.The normal stress-displacement relationship for the joints can be derived directly from the experimental stressdisplacement data (see Fig. 1(a)) of uniaxial compression tests carried out by Talesnick et al. (2007) on a three-block column.This approach is more direct than working out the stiffnesses from the slope of the curve as done by Tsesarsky & Talesnick (2007) (equation ( 1)).
In a three-block column, there are two rock joint contacts plus one equivalent rock joint contact, counting the two block interfaces at the two ends of the column in contact with the abutments.The deformations of the intact material were calculated assuming linear elastic behaviour for the intact blocks and a Young's modulus of 5600 MPa (see Table 1).The deformations due to the joints were then deduced by subtracting the deformations of the intact material from the total measured deformations.As shown in Fig. 1(b), the stiffness relationship proposed by Tsesarsky & Talesnick (2007) in equation ( 3) is less accurate.As a first estimate, the normal stiffness for the linear model in equation ( 3) is halved so that From Fig. 1(b), it emerges that the stress-displacement response given by equation ( 4) matches the experimental measurements well.
From the data reported by Talesnick (2007) it is apparent that the joint shear stiffness measured in the centrifuge decreases with shear displacement (see Fig. 2).However, shear stiffness degradation is not captured by the equation proposed by Tsesarsky & Talesnick (2007) (equation ( 2)), which reflects the initial (high) shear stiffness only and therefore leads to an underestimation of the shear displacements.The instantaneous shear stiffness can be modelled as where k s,initial is the initial shear stiffness, u s is the shear displacement and u peak is the displacement at which the shear stress reaches its peak.u peak (intercepts of the lines with the x-axis in Fig. 2) can be assumed as while the initial shear stiffness k s,initial can be modelled as (obtained by fitting a curve for the data points in Fig. 3 to match the y-intercept of Fig. 2) For high normal stresses, the initial shear stiffness predicted using equation ( 7) is lower than that from the linear model (equation ( 2)) proposed by Tsesarsky & Talesnick (2007) (Fig. 3).The model predictions for the instantaneous shear stiffness juxtaposed against experimental measurements are shown in Fig. 2.
In the case of complex stress paths, it is possible that equation ( 5) may give zero stiffness before yielding occurs;  for example, if the normal stress increases as the joint shears, the joint may experience u s ≥ u peak although the shear stress is still less than µσ n where µ is the joint friction coefficient.This behaviour is not physically sound.Rather than modelling stiffness degradation, it is numerically more robust to model this behaviour as strain hardening or continuous yielding, so that where μ 0 is the maximum friction coefficient (i.e.40°).In this shear model, the friction coefficient increases from zero with the shear displacement and the elastic shear stiffness is simply k s,initial .The stress-update algorithm for the DEM is illustrated in the Appendix.Figure 4 shows that equation ( 8) predicts the laboratory measurements well.

RESULTS
Comparison with finite-difference software (FLAC) simulations (Tsesarsky & Talesnick, 2007) The same values of joint stiffness as those adopted in Tsesarsky & Talesnick's FDM (FLAC (Itasca, 2000)) simulations (see Table 2) were employed in the DEM simulations.In two of the three simulations, the stiffness was assumed constant; in the third simulation, a linearly varying stiffness was adopted according to the contact law of equations ( 2) and (3).The horizontal thrusts obtained from the DEM simulations for the full-block beam together with the values obtained by the FDM analyses of Tsesarsky & Talesnick (2007) are shown in Fig. 5.A comparison of the midspan deflections is shown in Table 2 and one of the obtained deflection profiles is plotted in Fig. 6.Differences in block rotations (i.e.greater for the blocks at the abutment than at the midspan) are more pronounced in the DEM simulations.The figures show that both the horizontal thrusts and the deflection profiles obtained by the DEM simulations compare well with the results obtained by the FDM where a dense mesh (25 × 25 for each block) was used.
Since the DEM analyses are in agreement with the FDM simulations of Tsesarsky & Talesnick (2007), obviously the DEM simulations cannot provide a closer fit to the experimental measurements than the FDM.The main discrepancies concern the magnitude of the thrust, which is about 30% lower than the experimental measurements, and the displacements at both midspans and abutments, which are underpredicted as in the FDM simulations of Tsesarsky & Talesnick (2007).Assuming that the problem does not lie with the capabilities of the numerical codes, a strong possibility is that the contact models of equations ( 1) and ( 2) are not sufficiently accurate.

Verification of the line of thrust for DEM simulations with rigid blocks and deformable joints
In the study of masonry structures (Timoshenko, 1983;Heyman, 1997), the line of thrust under the structure's selfweight is known to assume the catenary shape of Hooke's inverted chain (refer to Lamb (1916) for the equations).On the other hand, in the rock engineering literature, the line of thrust of a voussoir beam is assumed to be of parabolic shape (Evans, 1941;Beer & Meek, 1982;Brady & Brown, 1993;Sofianos, 1996;Sofianos & Kapenis, 1998;Diederichs & Kaiser, 1999b).The contact points between blocks were thus fitted using both a uniform catenary curve and a parabola for gravitational accelerations of 10g and 40g to extrapolate the line of thrust for the beams.Figure 7 shows that both the catenary curve and the parabola fit the contact points equally well, with very little difference between them.Also, the DEM results compare well with the experimental measurements deduced from strain gauges.
Diederichs & Kaiser (1999b) considered the line of thrust obtained by UDEC (Itasca, 1996) for non-deformable abutments to be unrealistic because they expected, according to their analytical solutions, the contact force between the beam and abutment to be not located at the lower corner.This was subsequently experimentally confirmed by the centrifuge tests of Talesnick et al. (2007).This issue is probably due to the contact detection algorithm employed in UDEC and was not observed in the current analyses reported in this paper.

Results from the new contact model
The parameters and models employed are summarised in Table 3. Figure 5 shows that the more accurate contact models have negligible influence on the magnitude of horizontal thrust for the full-block beam.The discrepancies between the load cell and strain gauge measurements are still not resolved (Talesnick et al., 2007).However, using the proposed more accurate contact model in the normal direction (equation ( 4)) while keeping the shear stiffness model (equation ( 2)) of Tsesarsky & Talesnick (2007), the predicted deflections of the DEM simulations agree better with the experimental measurements (see Fig. 8).
Using the more refined shear contact model of equation ( 8), the calculated shear displacements at the abutment (Fig. 9) are slightly larger than those shown in Fig. 8 simulated using the simpler model of equation ( 2).Overall, for the full-block model, this resulted in an overprediction of approximately 28% at the midspan, averaged among the three sets of tests at 40g (Fig. 9), in comparison with the contact models of Tsesarsky & Talesnick (2007) (equations ( 1) and ( 2)) which resulted in an underprediction of approximately 40% on average (see Table 3 for a summary of the comparison).
The displacements at the abutment were underpredicted by the DEM simulations (Fig. 9).Poor predictions of shear displacements by the advanced shear model (equations ( 7) and ( 8)) are largely due to the lack of experimental shear data for model calibration at the (higher) normal stresses under which the experiments were conducted.In fact, the shear tests for stiffness calibration were carried out up to a normal stress of 0•023 MPa, in contrast to centrifuge tests for which the approximated normal stresses were up to 0•5 MPa for the full-block beam and 1•0 MPa for the half-block beam.

CONCLUSIONS
This paper has presented a rigorous DEM exercise that showcases the importance of using accurate contact laws to predict actual displacements resulting from rock block interactions.From our results, it emerges that discrepancies between DEM and FLAC predictions are small when the same contact constitutive model is used (see Fig. 6) and that the DEM predictions of the experimental behaviour improve significantly when more accurate contact models are used.The calibration of the models carried out in this paper shows the benefit of using customised approaches to derive specific contact models for the DEM.These approaches can be applied in other rock engineering applications as well.
Furthermore, it has been shown that the DEM with the contact detection algorithm of Boon et al. (2012) is capable of correctly calculating the line of thrust for jointed singlelayer rock beams subject to gravitational acceleration.

APPENDIX: STRESS-UPDATE ALGORITHM (SHEAR DIRECTION)
The total shear displacement increment du consists of elastic and plastic components, du e and du p , such that from which the increment of plastic displacement can be expressed as  Talesnick et al. (2007).The forces estimated from strain gauges were obtained from test 180-E (Talesnick, 2006).The magnitude (417 N derived from strain measurements) stated in figure 20 of Talesnick et al. (2007) seems to be inconsistent with the strain plots in figure 20, where the value is closer to approximately 349 N.Among the three DEM simulations presented here there is negligible difference in terms of horizontal thrust predictions where F i and F 0 are the shear forces at the current and previous time step, respectively, and k s is the shear stiffness (in N/m).
For the case when the coefficient of friction μ is modelled as strain hardening or strain softening, it is useful to denote μ = g(β), so that, μ is a function g of an internal state variable β.Then, define the increment of the internal state variable dβ as a scalar product of du p to the power of half: so that where β i and β 0 are the values of the internal state variable at the current and previous time step respectively.Solving of the equations in the stress-update algorithm assumes a dependency of the friction coefficient on the plastic shear displacements instead of the total shear displacements in equation ( 8) for calibration convenience (this very slightly delays the reaching of the asymptote for the exponential decay function in equation ( 8)).By assuming that the plastic strain increment is in the direction of current shear stress the state variable β i can be expressed as and the shear force as Define the yield surface as where F n is the magnitude of the normal force at the contact and μ i is the corresponding friction coefficient at the current time step.
In the stress-update algorithm, the plastic multiplier Λ is sought so that the shear force F i is on the yield surface shear stress being inside the yield surface.As a first approxi-mation, it can be taken as Λ = du/F 0 .If this shear stress is outside the yield surface, Λ is iteratively increased until it lies inside the yield surface.With a pair of bracketing range for Λ, the remaining unknowns in the yield surface (equation ( 16)) are the shear force F i and the current friction coefficient μ i .The first can be calculated from equation ( 15) and the second from μ i = g(β i ), where β i can in turn be calculated from equation ( 14).Bracketing algorithms can be used for this purpose.
Continuously yielding model, equations ( 7) and ( 8) Fig. 1.(a) Stress-deformation of load-unload cycles for uniaxial compression of a three-block column (after figure 14 of Talesnick et al. (2007), with kind permission from Springer Science and Business Media).(b) Contact model derived from Fig. 1(a) after subtracting the deformations of the intact rock material and dividing by three which is the total number of equivalent rock joints, including the topmost and bottommost block interfaces, in the uniaxial test of the three-block column

Fig. 2 .Fig. 4 .Fig. 3 .
Fig. 2. Shear stiffness as a function of shear displacement at different normal stresses.The experimental shear stiffness was taken from the second loading cycle: proposed model (equation (5)) plotted against experimental data deduced from figure 5(c) ofTalesnick (2007) Fig.5.Comparison of thrust build-up of the full-block and between DEM simulations and experiments (blue diamonds) carried out byTalesnick et al. (2007).The forces estimated from strain gauges were obtained from test 180-E(Talesnick, 2006).The magnitude (417 N derived from strain measurements) stated in figure20ofTalesnick et al. (2007) seems to be inconsistent with the strain plots in figure20, where the value is closer to approximately 349 N.Among the three DEM simulations presented here there is negligible difference in terms of horizontal thrust predictions

0Fig. 7 .Fig. 6 .
Fig. 7. Catenary line fit for contact points obtained from DEM analyses at 10g and 40g for the full-block model.The lines of thrust obtained from strain gauge measurements for different gravitational accelerations (24g and 40g) are not distinguishable (figure 20 in Talesnick et al. (2007)).For experimental measurements, the best-fit lines are assumed to pass through 0•0155 m from the centreline at the midspan.Note that, from Fig.5, similar conclusions can be drawn for the other contact models used here

Fig. 8 .Fig. 9 .
Fig. 8. DEM displacement profile for the full-block and k n = linear model with half stiffness (equation (4)) and k s = linear model (equation (2)) in comparison with experiments carried out byTalesnick et al. (2007).Experimental deflections were measured using linear variable differential transformer (LVDT) (test 180-E for the full-block model)

Table 2 .
Joint properties adopted in the DEM simulations and results of comparison with FLAC a Negative indicates underprediction Boon, Houlsby and Utili

Table 3 .
Boon (2013)the models employed in DEM simulations and results of comparison with experimental measurements (refer toBoon (2013)for plot of deflection profiles)