Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2022.34.0038 Acta Polytechnica CTU Proceedings 34:38–42, 2022 © 2022 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague THE EFFECT OF SAMPLING VOLUME SIZE ON THE APPARENT STIFFNESS OF JOINTED ROCK MASS Martin Lebeda∗, Petr Kabele Czech Technical University in Prague, Faculty of Civil Engineering, Thakurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: martin.lebeda@fsv.cvut.cz Abstract. Overall mechanical properties of a jointed rock mass are strongly affected by discontinuities – fractures – that naturally occur in rocks. Stochastically-generated discrete fracture network (DFN) modeling, which uses a probabilistic approach to describe the spatial distribution of fractures, such as position, size, or orientation, offers an explicit way to describe geometry of the fracture system. Many in-situ measurements and analyses presented in literature indicate that fractures’ sizes can be adequately represented by the power law probability distribution. The parallel plate model of individual fractures combined with an averaging technique makes it possible to estimate the overall compliance or stiffness of jointed rock mass (Oda et al. [1]). In the present study, a series of numerical simulations of jointed rock mass modeled by DFN and Oda’s approach were conducted to analyze the effect of different sizes of the sampling volume on the overall elastic moduli. The results of the numerical study show that the variance as well as the average of the apparent stiffness decrease as the size of sampling element grows. Keywords: Jointed rock mass, discrete fracture network (DFN), volume averaging, representative volume element (RVE), apparent elastic modulus. 1. Introduction Brittle structures, such as joints and faults, are found in almost every rock mass body. They vary greatly in quantity, dimensions, and arrangement. The presence of these fractures strongly affects the deformational response of the rock mass to stress and needs to be taken into account in geotechnical calculations. The finite element method (FEM) formulated on the ba- sis on the governing equations of solid mechanics is often used for these calculations. In such a case, it is feasible to model discretely (using, e.g., interface ele- ments) only the less-frequent largest fractures, while the smaller ones, which occur in a large quantity, can be represented by the means of so-called equivalent continuum. In the latter approach, a solid containing discontinuities is treated as a homogeneous continuum, whose stiffness or compliance is determined so as to correspond to the relations between the overall stress and overall strain of the fractured body. To this end, homogenization or averaging [2] over a certain sam- pling volume of the fractured body can be used. The smallest sampling volume element of the fractured body “for which the usual spatially constant ‘overall modulus’ macroscopic constitutive representation is a sufficiently accurate model to represent mean constitu- tive response” [3], is called the representative volume element (RVE). If it is possible to identify an RVE that is small enough to be idealized as a material point on the scale of the structural problem analyzed by the FEM, then the jointed rock mass can be modeled as a uniform equivalent continuum. If, however, the size of the RVE is larger than the resolution required for the structural problem, it is necessary to employ the concept of the statistical volume element (SVE). The overall properties of the SVE are, again, determined by homogenization or averaging. However, they need to be evaluated over a sufficient number of sampling volumes of the fractured solid and interpreted statisti- cally. The rock mass in the structural problem solved by the FEM is then modeled as a nonuniform equiva- lent continuum. To this end, random fields respecting the statistics of the SVEs can be used. Stochastically-generated discrete fracture network (DFN) models, which uses a probabilistic approach to describe the spatial distribution of fractures, such as position, size, or orientation, offers an explicit way to describe geometry of the fracture system in rock mass. Many in-situ measurements and analyses pre- sented in literature indicate that fractures’ sizes can be adequately represented by the power law proba- bility distribution [4], while various distributions for spherical data can model their orientations [5]. One appealing aspect of the DFN modeling is that the parameters of the probability distributions can be identified based on data from structural-geological survey, e.g. [6]. It follows that stochastic realizations of the DFN may serve as the geometrical description of the fractures in rock mass to be used for the earlier- discussed sampling and determination of the RVE or the SVE properties [7–9]. In this article, we present results of a numerical case study, in which a body of a jointed rock modeled by stochastic DFN realizations is sampled by volume elements of different sizes and the apparent overall Young moduli are evaluated and 38 https://doi.org/10.14311/APP.2022.34.0038 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 34/2022 The Effect of sampling volume size on the apparent stiffness of jointed rock mass analyzed. The DFNs were generated using the power law probability distribution of fractures’ sizes and uniform probability distribution of fractures’ location and orientation. The parallel plate model of fractures and the averaging scheme based on the crack tensors proposed by Oda [1] were employed to calculate the overall moduli. In contrast to the previous studies [7– 9] we carried out the analysis in 3D, adopted some different assumptions about the probabilistic distri- butions of fractures’ sizes and orientations or used a different averaging approach. 2. DFN model of joints’ geometry Geometry of individual fractures and the whole frac- ture system in the modeled rock mass are described by 3D stochastically generated DFN with the following assumptions. 2.1. Assumptions on fractures’ orientation, spatial distribution, and shape Fractures in natural rock can be usually grouped into sets (populations) based on their geological origin and dominant orientation. However, for the sake of sim- plicity, to avoid the issue of overall anisotropy of the rock mass, we assume the distribution of fracture ori- entations as uniform on the sphere (Fisher et al. [5]). Fractures’ centers are generated by Poisson random process, so the probability distribution of center coor- dinates in space is also uniform. The total number of generated fractures is controlled by a given volumetric density P30 (number of fractures’ centroids over a unit volume). Individual fractures are idealized as squares. 2.2. Power law distribution of the fracture sizes Fractures’ sizes are assumed to follow the power law distribution, whose probability density function is expressed as: p(x) = α− 1 xmin ( x xmin )−α for x > xmin and α > 1, (1) where x is the radius of the circle circumscribed to the fracture, xmin is the minimum fracture size (lo- cation parameter), α is the law’s exponent (shape parameter). The parameters used for generating stochastic DFNs and for subsequent calculations are listed in Ta- ble 1. The volumetric density and power law exponent are adapted from our previous work [10], where they were estimated for a highly jointed metamorphic rock based on data from structural-geological mapping of fracture traces on tunnel walls. Thus, the range of fracture sizes covered by the model is appropriate for the scale of an underground structure (e.g. tunnel). The DFN models used in this paper were generated using DFraM software (Kabele et al. [11]). Parameter Value volumetric density P30 [ 1m3 ] 2.56 power law exponent α [-] 3.4 minimum fracture size xmin [m] 0.3 volume of the rock mass [m3] 100·100·100 Table 1. Parameters of DFNs. 3. The model of the overall stress-strain behavior of jointed rock mass Oda [1] proposed the following formula for overall apparent elastic stress-strain relation of jointed rock mass: εij = 1 E [ (1 + ν)δikδjl −νδijδkl + ( 1 kn − 1 ks ) Fijkl + 1 4ks (δikFjl + δjkFil + δilFjk + δjlFik) ] σkl = Cijklσkl, (2) where E and ν are Young’s modulus and Poisson’s ratio of the intact rock, respectively, kn and ks are nondimensional parameters related to the fracture (joint) normal and tangent stiffness, respectively, δij is Kronecker’s delta, and Fij and Fijkl are so-called second and fourth rank crack tensors, respectively: Fij = 1 V M∑ k=1 S(k)r(k)n (k) i n (k) j , (3) Fijlm = 1 V M∑ k=1 S(k)r(k)n (k) i n (k) j n (k) l n (k) m . (4) Here v is the volume of the sample, S(k) is the area of k-th fracture inside the sampling volume, r(k) is a typical size of the fracture equal to a diameter of a circle with the same area and n(k)i are the components of a unit vector normal to the fracture. Note that the typical size r is assumed as a property of the fracture independent on the size of the sampling volume and it is calculated using the original area of the fracture generated in stochastic DFN. Individual fractures are in Oda’s approach repre- sented by the so-called parallel plate model as a set of parallel flat plates connected by two springs that symbolize normal and shear stiffness. The fractures’ stiffnesses influence the evaluation of the compliance tensor through the non-dimensional stiffness param- eters kn and ks, as seen in Eq. (2). The values of the rock’s and joints’ mechanical parameters used in this paper are listed in Tab. 2. The joint’s stiffnesses were determined assuming a small stress change in a rock at the depth of 550 m and the other parameters 39 Martin Lebeda, Petr Kabele Acta Polytechnica CTU Proceedings Parameter Value normal stiffness kn [-] 1.356 shear stiffness ks [-] 1.188 Poisson’s ratio of intact rock ν [-] 0.164 elastic modulus of intact rock E [GPa] 44.0 Table 2. Parameters of Oda’s parallel plate model. (E and ν) were determined with the basis on in-situ measurements on the same rock as described in [10]. 4. Methodology of the analysis 4.1. General considerations We consider that the overall stiffness of the rock mass evaluated by Oda’s method on the basis of the DFN would be eventually used for a finite element analysis of stress and deformation in the vicinity of an un- derground structure, e.g. a tunnel or deep geological repository of radioactive waste. In this scenario, it is plausible to assume that the equivalent continuum model would be employed to account for joints only up to a certain maximum size Amax, while larger frac- tures would be represented as discrete discontinuities with, e.g., interface elements or extended FEM. In the context of the underlying DFN it means that only the abundant smaller joints need to be present in the stochastically generated DFN, while the less frequent larger fractures would be mapped and implemented in a deterministic way. To set the threshold value Amax, we used a sufficiently low value of the power law’s complementary cumulative distribution function (CCDF), which is defined as the probability of oc- currence of a fracture with size greater than Amax. With the parameters listed in Tab. 1 and CCDF equal to 1·10−5, the corresponding limit fracture size is Amax = 51.4 m·51.4 m = 2642 m2. This value seems acceptable also from the geological survey’s point of view, as brittle geological structures of this size are likely to be mapped deterministically by boreholes and trial tunnels. When the DFN is generated by the DFraM pro- gram, centers of individual fractures are placed by a random process inside a given prismatic bounding box. Thus, there is an inner region along the bounding box boundary, which could be intersected by fractures from outside of the box, but these fractures are not generated as their centers lie outside the box. The “thickness” of this boundary zone corresponds to the radius of the circle circumscribed to the largest frac- ture, that is 51.4 m2 · √ 2 = 36.3 m. To avoid distortion of results by this effect, the sampling volumes must be taken from a subdomain of the DFN bounding box excluding these boundary regions. 4.2. Generation of DFN models and sampling volumes To obtain a statistically relevant set of sampling vol- umes for the evaluation of the apparent overall stiff- ness, 10 different stochastic realizations of DFN model with the same probabilistic parameters and the same parameters of rock mass (Tab. 1 and Tab. 2) but dif- ferent random seed were generated. The bounding box of the DFN fractures’ centroids was a cube with side length of 100 m. Cubic sampling volumes of different sizes were cre- ated by clipping the DFN models, as shown in Fig. 1. The samples were always concentric with the DFN bounding box. Considering the DFN box size and the elimination of the boundary effect, the side length of the largest sampling volume was 27.3 m. The edge of the smallest sampling volume was chosen 0.5 m to be slightly less than the diagonal of the smallest fracture (2·xmin). In total, 10 sampling volumes (with edge sizes 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 27.3 meters) were extracted from each realization of the stochastic DFN. Figure 1. Sample volume. 4.3. Averaging procedure – calculation of overall strain and apparent stiffness Averaging procedure – calculation of overall strain and apparent stiffness For each sampling volume, the overall apparent Young moduli were calculated using Eq. (1). Although uniform probability dis- tributions were used for the location and orienta- tion of the fractures in the underlying DFN model, which would suggest isotropy, the apparent mod- uli were calculated for three orthogonal directions: Ex = 1/C1111,Ey = 1/C2222,Ez = 1/C3333. This way of calculation allowed us to analyze not only the magnitude of the stiffness but also its directional dependence. 5. Results Figures 2–4 present the dependence of apparent de- formation moduli calculated in three orthogonal di- 40 vol. 34/2022 The Effect of sampling volume size on the apparent stiffness of jointed rock mass rections on the size of sampling volume. Figure 2. Apparent deformation moduli – x axis. Figure 3. Apparent deformation moduli – y axis. Figure 4. Apparent deformation moduli – z axis. The results show the same trends for moduli in all three directions. As the sampling volume grows, the mean of the apparent moduli obtained from different realizations of the DFN decreases and converges to a certain value. The maxima have the same tendency, while the minima grow. The dependence of the mean values can be well approximated by an exponential curve, as shown by the dotted line. The graphs also show that the smaller are the sample volumes, the greater are the deviations of the moduli obtained from the different realization of the DFN. 6. Preliminary analysis of anisotropy of apparent stiffness From Figures 2–4 it is obvious, that in spite of the uni- form statistical distribution of fractures’ orientation in the underlying DFN, the apparent overall moduli of the rock mass in three orthogonal directions are not equal. This implies that the overall stiffness of the rock mass predicted by the model is not isotropic. Hereafter, we present the results of a preliminary anal- ysis that will allow us to gain some insight into how the detected anisotropy is related to the sampling element size. Various measures of anisotropy have been proposed in the literature, especially in connection with crystal elasticity – see a review e.g. in [12]. However, as in the present study we examine only the equivalent deformation moduli, we define the index of anisotropy ad-hoc as the coefficient of variation of the appar- ent moduli Ex,Ey,Ez evaluated in three orthogonal directions: α = std dev(Ex,Ey,Ez) mean(Ex,Ey,Ez) . (5) It is obvious that for an isotropic material, where all moduli are equal, is zero, whereas anisotropy, man- ifested by a larger variation among the moduli, results in a larger value of α. Figure 5. Sample size dependence of the index of anisotropy α. Figure 5 shows the dependence of the index on the sampling volume element size. As in the case of the elastic moduli, we can see that with increasing size of the sampling volume, the mean and maximum values as well as the scatter of the index decrease. The minima do not exhibit a clear trend and fluctuate close to zero for all volume sizes. The mean value of α for the largest volume is approaching, but is not equal to, zero. This is understandable, since the index α can be only nonnegative and zero mean value would require that in all model realizations a perfectly isotropic response was achieved. These observations indicate that the model converges to an isotropic response (in the sense of Ex ≈ Ey ≈ Ez). However, we point out that these are only preliminary conclusions and a more systematic analysis involving the complete apparent stiffness or compliance tensor and appropriate index of anisotropy should be carried out. 7. Conclusions The concept of using Oda’s [1] averaging procedure in conjunction with 3D stochastic DFN has been ex- amined with the aim to determine the dependence of apparent rock mass stiffness on sampling volume size. The results of the case study allow us to draw some interesting conclusions: • With growing sample sizes, the mean value of the apparent stiffness is not constant but decreases and converges to a certain value. The minima and max- ima converge to the same value. 41 Martin Lebeda, Petr Kabele Acta Polytechnica CTU Proceedings • In spite of assuming directionally and spatially uniform probability distributions of fractures, the power law variability of fracture size is sufficient to produce anisotropic overall deformation moduli. The measure of anisotropy decreases with increasing sampling volume size. • The great variance of stiffnesses as well as anisotropy calculated on smaller sampling volumes (up to about 15 m large side of the sample cube with the parame- ters used in the present study) means that these vol- ume elements cannot be regarded as representative (RVE), but stochastic (SVE). The corresponding overall properties must be then considered not as effective, but apparent [13]. Should these overall properties be used to define equivalent continuum for finite element analysis of an underground struc- ture, the variance must be taken into account, e.g. by using the concept of random fields. It is noted that this may be relevant for analysis of many struc- tures in underground engineering, such as tunnels, where, due to the spatial variability of the stress field and the size of the structure, model resolu- tion on the order of 100 m or even 10−1 m may be required. • On the other hand, if the required finite element model resolution is on the order of 101 m or larger, using the mean stiffness obtained by the averaging over RVE as uniform effective property of the equiv- alent continuum may be acceptable. This could be the case, for example, when variations of the geostatic stress field need to be estimated for hy- drogeological simulations on a larger scale. Acknowledgements This paper was financially supported by Czech Tech- nical University in Prague under No. SGS project SGS2/037/OHK1/1T/11. References [1] M. Oda, T. Yanabe, Y. Ishizuka, et al. Elastic stress and strain in jointed rock masses by means of crack tensor analysis. Rock Mechanics and Rock Engineering 26(2):89– 112, 1993. https://doi.org/10.1007/BF01023618. [2] M. Hori, S. Nemat-Nasser. On two micromechanics theories for determining micro–macro relations in heterogeneous solids. Mechanics of Materials 31(10):667–682, 1999. https://doi.org/10.1016/S0167-6636(99)00020-4. [3] W. J. Drugan, J. R. Willis. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. Journal of the Mechanics and Physics of Solids 44(4):497–524, 1996. https://doi.org/10.1016/0022-5096(96)00007-5. [4] E. Bonnet, O. Bour, N. E. Odling, et al. Scaling of fracture systems in geological media. Reviews of Geophysics 39(3):347–383, 2001. https://doi.org/10.1029/1999RG000074. [5] N. I. Fisher, T. Lewis, B. J. J. Embleton. Statistical analysis of spherical data. Cambridge University Press (CUP), 1993. https://doi.org/10.1017/cbo9780511564345. [6] P. Kabele, O. Švagera, M. Somr, et al. Mathematical modeling of brittle fractures in rock mass by means of the dfn method, 2017. [7] K.-B. Min, L. Jing. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. International Journal of Rock Mechanics and Mining Sciences 40(6):795–816, 2003. https://doi.org/10.1016/S1365-1609(03)00038-8. [8] M. Gutierrez, D.-J. Youn. Effects of fracture distribution and length scale on the equivalent continuum elastic compliance of fractured rock masses. Journal of Rock Mechanics and Geotechnical Engineering 7(6):626–637, 2015. https://doi.org/10.1016/j.jrmge.2015.07.006. [9] H. Ni, W. Xu, W. Wang, et al. Statistical analysis of scale effect on equivalent elastic modulus of jointed rock masses. Arab J Geosci 9(3):206, 2016. https://doi.org/10.1007/s12517-015-2190-z. [10] M. Lebeda. Určení tuhosti porušeného horninového masivu s využitím dfn modelu, 2020. [11] P. Kabele, O. Švagera, M. Somr, et al. Mathematical modeling of brittle fractures in rock mass by means of the dfn method, 2017. [12] S. I. Ranganathan, M. Ostoja-Starzewski. Universal elastic anisotropy index. Physical Review Letters 101(5):055504, 2008. https://doi.org/10.1103/PhysRevLett.101.055504. [13] T. Kanit, F. N’Guyen, S. Forest, et al. Apparent and effective physical properties of heterogeneous materials: Representativity of samples of two materials from food industry. Computer Methods in Applied Mechanics and Engineering 195(33-36):3960–3982, 2006. https://doi.org/10.1016/j.cma.2005.07.022. 42 https://doi.org/10.1007/BF01023618 https://doi.org/10.1016/S0167-6636(99)00020-4 https://doi.org/10.1016/0022-5096(96)00007-5 https://doi.org/10.1029/1999RG000074 https://doi.org/10.1017/cbo9780511564345 https://doi.org/10.1016/S1365-1609(03)00038-8 https://doi.org/10.1016/j.jrmge.2015.07.006 https://doi.org/10.1007/s12517-015-2190-z https://doi.org/10.1103/PhysRevLett.101.055504 https://doi.org/10.1016/j.cma.2005.07.022 Acta Polytechnica CTU Proceedings 34:38–42, 2022 1 Introduction 2 DFN model of joints’ geometry 2.1 Assumptions on fractures’ orientation, spatial distribution, and shape 2.2 Power law distribution of the fracture sizes 3 The model of the overall stress-strain behavior of jointed rock mass 4 Methodology of the analysis 4.1 General considerations 4.2 Generation of DFN models and sampling volumes 4.3 Averaging procedure – calculation of overall strain and apparent stiffness 5 Results 6 Preliminary analysis of anisotropy of apparent stiffness 7 Conclusions Acknowledgements References