Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.49.0048 Acta Polytechnica CTU Proceedings 49:48–54, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague THE CONCEPT OF FRACTURED ROCK-MASS MODELING USING DFN-BASED STATISTICAL VOLUME ELEMENTS Martin Lebeda∗, Petr Kabele Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: martin.lebeda@fsv.cvut.cz Abstract. The overall mechanical response of a fractured rock mass is, to a large extent, affected by naturally occurring fractures that exhibit sizes from millimeters to kilometers. Thus, in analysis of underground structures, such as tunnels, it is required that the fractures’ influence on the stress and deformation state in the vicinity of the structure is taken into account. In the present work, we examine the applicability of the statistical volume element (SVE) approach to determining the apparent stiffness tensor of an equivalent continuum representation of the fractured rock, which is then used in the framework of the finite element method. The equivalent continuum properties are determined by volume-averaging the effect of individual fractures that intersect the SVE, while the fractures are represented using the “parallel plate model”. Stochastically generated discrete fracture networks are used to represent the fractures’ geometry. Presently, we solve the problem linearly for an incremental change of the stress state. An application of the concept is demonstrated on simulation of tunnel excavation. Keywords: Fractured rock mass, discrete fracture network (DFN), finite element method (FEM), statistical volume element (SVE). 1. Introduction Rock masses generally exhibit material heterogeneity at multiple scales. At the level of microscale, the het- erogeneity can be associated with variable properties of mineral grains and microstructural defects. On the other hand, on the macroscale the behavior of rock mass is dominated by the stiffness of the rock matrix, presence of different types of discontinuities (generally called fractures) and the filling material or fluids inside the fractures. Considering the typical dimensions of underground structures, such as tunnels, their analysis and design generally correspond with the macroscale. Hence, when these structures are excavated in frac- tured rock, the effect of fractures must be taken into account [2]. The finite element method (FEM) is often used for the analyses. Then, one efficient way is to include only a few largest fractures as discrete discon- tinuities in the model, while the rest of the fractures is substituted by an equivalent continuum with effective or apparent properties of the rock mass. These overall parameters are evaluated by averaging or homoge- nization procedure [3]. In this paper, we focus on the equivalent continuum modeling and follow up on our previous work [4], where stochastically generated discrete fracture network (DFN) models [5, 6] provide a geometrical description of the fracture system in the rock mass. Note that this approach makes it possible to utilize data from geological survey to statistically describe the fractured rock mass. The equivalent continuum methods are commonly used in conjunction with the concept of the represen- tative volume element (RVE) [7–9]. In principle, the RVE is the smallest volume, which encompasses a suf- ficient sample of the substructures (e.g., fractures), so as to provide effective constitutive properties, which are nearly invariable with respect to the RVE place- ment. On the other hand, on the scale of the whole analyzed problem (e.g., tunnel), RVE must be small enough with respect to the spatial variation of macro- scopic stress to represent a material point. Our previ- ous research [4] shows that, for DFNs with power-law fracture size distribution (with minimum fracture size on the order of decimeters), the RVE size should be on the order of tens of meters. However, in the case of tunnel analysis, this RVE size does not meet the second criterion, as the stress shows high gradients in the close proximity of the excavation (units of meters). In this regard, La Pointe et al. [10] questioned the applicability of the continuum simulations on the scale of 10 to 100 meters, while fracture sizes occurring in the rock mass may range from millimeters to kilo- meters. Furthermore, Ostoja-Starzewski [11] pointed out the spatial randomness of the rock material and incorporated the concept of statistical volume element (SVE) on the level of mesoscale, separating the mi- croscale of heterogeneities from the macroscale that is usually connected with the RVE approach. The procedure of determination of overall properties us- ing volume elements smaller than representative is in agreement with widely accepted study presented by Kanit et al. [12], which suggests the evaluation of the effective properties as mean values of apparent properties of rather small volumes considering suffi- cient number of realizations. Against this background, 48 https://doi.org/10.14311/APP.2024.49.0048 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 49/2024 The concept of fractured rock-mass modeling using . . . Fracture set 1 2 3 4 5 6 7 8 9 Fisher distribution Mean dip [°] 85.365 36.871 3.360 87.507 62.308 88.629 78.236 88.629 83.253 Mean strike [°] 222.949 49.597 313.665 316.780 133.438 346.277 1.040 269.482 85.957 κ [-] 35.056 2.704 28.633 25.351 2.913 18.730 3.815 17.300 3.962 Power law distribution amin [m] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 amax [m] 1 000 1 000 1 000 1 000 1 000 1 000 1 000 1 000 1 000 α [-] 3.328 3.062 3.800 3.789 3.382 3.042 3.000 3.607 3.296 P30 [m−3] 0.269 0.220 0.061 0.246 0.427 0.232 0.146 0.581 0.380 Table 1. Parameters of DFNs [1]. we adopt the SVE approach to simulations of tun- nel excavation, while taking into account the spatial, orientational, and dimensional randomness of rock fractures. In this study, we present ongoing work on develop- ment of the 3D FE stress and displacement analysis that incorporates parameters of equivalent continuum obtained by averaging procedure in conjunction with the 3D stochastic DFN. We use the averaging method based on the approach proposed by Oda et al. [9] that treats individual fractures by the “parallel plate model”, which means that the deformation response of a fracture to normal and shear stress is represented by two dimensionless stiffness parameters. At present stage, we have implemented the FE analysis linearly with constant fractures’ stiffness for an incremental change of stress state. Although the linear solution is admissible only for small stress changes, the applica- bility of the SVE concept can be demonstrated. 2. DFN model of fractures’ geometry The description of the rock-mass inherent fracturing state represents the crucial part of underground struc- tures’ analysis. The stochastically generated DFN approach provides an adequate balance between the computational efficiency and the accuracy [5, 6]. The quantitative input information about simulated frac- ture network is often obtained by geological survey that usually consists of structural-geological mapping of fracture traces on rock outcrops or tunnel walls and of inspection of boreholes (observation windows). As the survey provides only a limited information about fractures deep in the rock mass volume, the data are used in the framework of so-called stochastic DFN modeling. In this approach, individual fractures are idealized as polygons, whose size, orientation, and spatial placement are assumed to follow certain prob- abilistic distributions. The distributions’ parameters are identified so as to match (in a statistical sense) the traces’ sizes, directions, and density within the obser- vation windows. Based on the calibrated probabilistic distributions any desired number of DFN realizations can be generated [13]. In the following parts of this study we use frac- ture network models generated by in-house software DFN_toolbox [14]. Parameters of these DFNs were calibrated based on data collected in the past by the Czech Geological Survey at the underground research facility Bukov, which is located at a depth of 550 me- tres below the Earth’s surface. Detailed description of the used measurement procedures and the acquired data are presented in report [1]. Positioning of fractures’ centers in the model vol- ume is ruled by Poisson random generation process. The abundance of fractures is defined by volumetric density P30 (number of fractures per unit volume). Fractures are modeled as squares, whose sizes a (radii of circumscribed circles) are controlled by the power law distribution with parameters amin, which is the minimum fracture size (location parameter), and α, which is the law’s exponent (shape parameter). For practical reasons, only fractures smaller than chosen value amax are generated. The orientation of fractures is determined by Fisher distribution, with parameters mean dip and mean strike (from which the mean unit normal vector µ can be determined), and κ, which is the concentration parameter. The whole network consists of nine fracture sets, where each of them is defined by individual density, orientation, and size parameters. All of the parame- ters are listed in Table 1. Note that steeply dipping fractures prevail (namely sets 1, 4, 6, 8, 9) while sub- horizontal (set 3) and mildly dipping fractures are less abundant. 3. Averaging procedure The overall effective properties of fractured rock mass are calculated by averaging procedure that was orig- inally proposed by Oda et al. [9]. The main Equa- tion (1) of the procedure represents the relation be- tween macroscopic stress σkl and the corresponding strain εij that are connected by the fourth order com- pliance tensor Cijkl: 49 Martin Lebeda, Petr Kabele Acta Polytechnica CTU Proceedings εij = 1 2G ( 1 2(δikδjl + δilδjk) − ν 1 + ν δijδkl ) σkl + 1 E [( 1 kn − 1 ks ) Fijkl + 1 4ks (δikFjl + δjkFil + δilFjk + δjlFik) ] σkl = Cijklσkl, (1) where G is the shear modulus, E is the Young’s modulus, ν is the Poisson’s ratio of the intact rock, kn is the nondimensional parameter related to the fracture’s normal, ks is the nondimensional parameter related to the tangent stiffness, δij is Kronecker’s delta, Fij and Fijkl are so-called second and fourth rank crack tensors, respectively: Fij = 1 V M∑ k=1 S(k)L(k)n (k) i n (k) j , (2) Fijlm = 1 V M∑ k=1 S(k)L(k)n (k) i n (k) j n (k) l n(k) m , (3) where V is the volume of the domain over which averaging is performed, S(k) is the area of k-th fracture inside the sampling volume, L(k) is a typical size of the fracture (equal to the diameter of a circle with the same area), n (k) i are the components of a unit vector normal to the fracture. Note that, the typical fracture size L is evaluated using the whole area of the fracture and not of the part inside the sampling volume because it is a property of the fracture independent on the size of the sampling volume. The dimensionless stiffness parameters kn and ks are, in this study, set while considering fractures’ me- chanical response in the depth of hundreds of meters below the ground level, which corresponds, for exam- ple, to the potential location of a deep underground radioactive waste repository. It is assumed that the rock mass is geostatically compressed to the point that all fractures are closed, and their normal stiffness is nearly infinite. On the contrary, fractures are assumed to be vulnerable to slip in the case of changing stress state due to, for example, tunnel excavation, meaning that the shear stiffness is almost nil. In that sense, the parameters are set as kn = 10 000, ks = 1.0. Further Figure 1. Illustration example of three SVEs corre- sponding to finite elements. In each SVE, only the intersecting portions of the underlying DFN fractures are shown. investigation of this assumption is needed, however the present study deals with the applicability of the SVE concept, which is feasible regardless of stiffnesses choice. 3.1. Statistical volume element The concept of using the SVE to obtain averaged mechanical properties of fractured and heterogenous materials has already been studied by other authors at the level of microscale, for example, [15, 16]. The most frequently used shapes of SVE are square and circle in the case of 2D models and cube or sphere for 3D models. The size of the statistical sampling element should be chosen regarding coefficient of variation of the monitored variable [17]. In this study, we use either cubical SVEs with a fixed side size (in the concept verification stage, Section 4) or SVEs, whose shape, size, and position correspond with those of the finite elements (for the FE simulations, Section 5) to obtain the apparent stiff- ness tensor of the fractured rock mass. To this end, each SVE’s compliance tensor is first calculated using Equations (1) to (3), S(k), L(k), n (k) i in which are determined from all fractures of the underlaying DFN, which intersect the SVE (Figure 1). Subsequently, the stiffness tensor is obtained by inverting the compli- ance. Note that, in general, this procedure yields fully anisotropic tensors. Finally, the stiffness tensors are assigned to the corresponding finite elements, which generally means that each finite element is composed of a unique material (Figure 2). The whole process has been implemented by means of Python API in the open-source FE code OOFEM [18]. 4. The SVE concept verification First, we need to verify that the SVE concept is appli- cable to analyses, such as, FE simulation of a tunnel excavation in a fractured rock modelled by DFN with parameters given in Table 1. The verification is in- spired by the study presented by Kanit et al. [12]. In particular, we investigate, whether, when averaging 50 vol. 49/2024 The concept of fractured rock-mass modeling using . . . (a). Distribution of the apparent modulus Ex obtained by SVEs. For comparison, the value of Ex calculated by averaging over large AV on the same underlaying DFN is 27 100 MPa. (b). Section view of the underlaying DFN (right). Figure 2. Distribution of the apparent modulus Ex obtained by SVEs and section view of the underlaying DFN. is performed on volumes smaller than RVE, similar apparent moduli are obtained when samples are taken from multiple locations of one realization of the DFN as well as by sampling one location in multiple DFN realizations. To this end, we use 10 different statistical realiza- tions of the DFN model with parameters described in the Section 2. The fractures are generated with centers in a cubical domain with the side size of 50 m. As a reference, we firstly calculate the effective mod- uli for all DFN realizations by volume-averaging over a sub-domain of 30 × 30 × 30 m3. The sub-domain is placed at the center of the larger DFN to reduce the boundary effect. Based on a previous study [19], this sub-domain, denoted as large averaging volume (AV) should be large enough to serve as RVE. Secondly, moduli are calculated from 10 SVEs with dimensions of 2 × 2 × 2 m3, which are placed randomly far apart but all in the central portion of one realization of the DFN. Thirdly, moduli are obtained for one SVE with dimensions of 2 × 2 × 2 m3, which is placed at the same location near the center in all 10 realizations of the DFN. The SVE size was chosen so as to roughly correspond with the size of typical finite elements used in the subsequent simulations. Although the complete stiffness tensors are calcu- lated, we use Young’s moduli Ex, Ey, Ez, which cor- respond to the tensor elements D1111, D2222, D3333, respectively, to illustrate the results. The data shown in Figures 3–5 confirm the expected difference between values obtained by using the SVEs and by using the large averaging volume. When moduli are evaluated over the large AV, the results from 10 realizations show coefficient of variation (COV) of 5.2 % or less. In the case of the SVE calculations, the results are more dispersed with COV ranging between 15.8 and Figure 3. Comparison of stiffness Ex. The boxes indi- cate the quartiles, the whiskers indicate the minimum and maximum and the cross indicates the mean. Figure 4. Comparison of stiffness Ey. The boxes indi- cate the quartiles, the whiskers indicate the minimum and maximum and the cross indicates the mean. 51 Martin Lebeda, Petr Kabele Acta Polytechnica CTU Proceedings Figure 5. Comparison of stiffness Ez. The boxes indi- cate the quartiles, the whiskers indicate the minimum and maximum and the cross indicates the mean. 19.7 % and the mean values of apparent moduli are higher. Furthermore, both SVEs approaches (spatial sampling and multiple realizations) are in a fair agree- ment. Especially, in the case of the Ey and Ez, the mean values differ by less than 2 %. In the case of the difference is approximately 7 %. 5. Numerical simulation of tunnel excavation To demonstrate the use of the SVE modeling approach, 3D simulation of tunnel excavation is performed. The modelled rock domain consists of 30×30×10 m3 prism, where the short dimension is parallel to the tunnel axis. The tunnel passes through the center of the domain and has a diameter of 3 m. The displacement of the rock body is constrained by statically determinate supports. The excavation problem is solved in two calculation steps, see Figure 6. In the first step the initial stress state is induced and in the second step the excavation of the tunnel is simulated by removal of elements. To impose the initial stress state, the FE model is loaded by normal stresses that are applied to the outer surface of the rock body. The initial loading stress tensor is based on field measurements from Rožná mine [20], which is a part of the underground research facility Bukov. The data from vertical boreholes showed that the directions of the principal stresses are vertical and horizontal. The vertical principal stress is calculated as: Sv = −ρ · g · h, (4) where ρ is volumetric weight, g is gravity acceleration, h is depth beneath the ground level. The parameters for calculation of the vertical principal stress (Table 2) are adopted from [21]. As the major (a). State prior to excavation. (b). State after excavation. Figure 6. The finite element mesh. Parameter Value Major horizontal principal stress SH [MPa] −22.8 Major horizontal principal stress Sh [MPa] −14.6 Volumetric weight ρ [kg m−3] 28 000 Depth of the mine beneath the gnd level [m] 550 Vertical principal stress Sv [MPa] −15.4 Table 2. The initial loading stress parameters. and minor horizontal principal stresses, SH and Sh, respectively, we use the average of the values measured on different levels of the mine. The values of the initial stress are listed in Table 2. Orientation of the model is set so, that the tunnel axis is parallel with the minor horizontal principal stress, which means that the direction of Sh corresponds to the global axis y. For the purpose of comparison, simulations are cal- culated, firstly, using the uniform effective stiffness tensor, obtained by averaging over the large volume, in all finite elements. Secondly, the concept of SVE is applied with the same underlaying DFNs. In the latter case, the stiffness tensor is calculated by averag- 52 vol. 49/2024 The concept of fractured rock-mass modeling using . . . ing over SVEs that coincide with the finite elements, which yields nonuniform stiffness over the analyzed domain, as illustrated in Figure 2. Specifically, the values of modulus Ex are plotted over individual ele- ments along the front cross-section perpendicular to the tunnel axis located at y = −3.75 m. The figure also shows the section of the underlaying DFN. It is noted that the case of uniform stiffness, the stress state prior to the tunnel excavation is also uni- form within the modelled domain. On contrary, the calculations using the SVE concept exhibit nonuni- form initial stress state, which is caused by different material stiffness in each finite element (see Figure 2). 5.1. Results of the study The calculations were performed for 10 realizations of the stochastic DFN. The main monitored result values are horizontal and vertical convergences of the tunnel walls measured on the cross sections passing through the tunnel center, see Figure 7 and Figure 8. The results show that the absolute values of tunnel walls convergences are lower in the case of the SVE ap- proach, which is in agreement with Section 4, where higher stiffnesses were obtained while using SVEs. However, the variances of the results are comparable regardless of the size of the averaging volume. Further investigation of this phenomenon is needed since it may be caused by various effects, e.g., low number of statistical realizations or the specific spatial, orien- tational, and dimensional characteristics of the DFN fracture sets. 6. Conclusion The applicability of the SVE concept incorporated in the 3D FEM simulation of tunnel excavation using overall rock-mass parameters obtained by averaging procedure based on DFN has been examined. The results of the study lead to the following conclusions: • The variance of the overall moduli obtained by averaging over SVEs (COV = 15.8–19.7 %) is higher compared to averaging over large volume (COV = 5.2 %), which is an expected outcome. • The agreement of the two different SVE evaluation methods (spatial sampling and multiple realizations) is fairly good with the highest difference of 7 %. • Lower convergences of tunnel walls appear while using the SVE approach compared to large averag- ing volume, which is in agreement with the higher stiffnesses obtained in case of SVEs. • On the other hand, the variance of convergences of tunnel walls are comparable regardless of the size of the averaging volume. This phenomenon requires further investigation since it is not an expected outcome. Figure 7. Horizontal convergence of tunnel walls. The boxes indicate the quartiles, the whiskers indicate the minimum and maximum and the cross indicates the mean for 10 realizations of stochastic DFN. Figure 8. Vertical convergence of tunnel walls. The boxes indicate the quartiles, the whiskers indicate the minimum and maximum and the cross indicates the mean for 10 realizations of stochastic DFN. Acknowledgements This paper was financially supported by Czech Tech- nical University in Prague under SGS project no. SGS23/032/OHK1/1T/11. References [1] L. Gvoždík, P. Kabele, J. Říha, et al. Transport of radionuclides from deep geological repository/testing of conceptual and numeric models. Tech. rep., Progeo s. r. o., 2020. [2] L. Jing. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. International Journal of Rock Mechanics and Mining Sciences 40(3):283–353, 2003. https://doi.org/10.1016/S1365-1609(03)00013-3 [3] 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 53 https://doi.org/10.1016/S1365-1609(03)00013-3 https://doi.org/10.1016/S0167-6636(99)00020-4 Martin Lebeda, Petr Kabele Acta Polytechnica CTU Proceedings [4] M. Lebeda, P. Kabele. Analysis of tunnel excavation based on linear DFN-FEM modelling. Acta Poytechnica CTU Proceedings 40:61–68, 2023. https://doi.org/10.14311/APP.2023.40.0061 [5] J. Andersson, A. M. Shapiro, J. Bear. A stochastic model of a fractured rock conditioned by measured information. Water Resources Research 20(1):79–88, 1984. https://doi.org/10.1029/wr020i001p00079 [6] C. Darel, P. Davy, O. Bour, J.-R. De Dreuzy. Discrete fracture network for the Forsmark site. Tech. rep., SKB – Svensk Kärnbränslehantering AB Swedish Nuclear Fuel and Waste Management Co, 2006. [2023-08-22]. http://www.iaea.org/inis/collection/ NCLCollectionStore/_Public/38/013/38013872.pdf [7] M. Kachanov. Continuum model of medium with cracks. Journal of the Engineering Mechanics Division 106(5):1039–1051, 1980. https://doi.org/10.1061/JMCEA3.0002642 [8] M. Cai, H. Horii. A constitutive model of highly jointed rock masses. Mechanics of Materials 13(3):217–246, 1992. https://doi.org/10.1016/0167-6636(92)90004-W [9] 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 [10] P. R. La Pointe, P. C. Wallmann, S. Follin. Continuum modelling of fractured rock masses: Is it useful? In Prediction and performance in rock mechanics and rock engineering, vol. 2, pp. 343–350. 1996. [11] M. Ostoja-Starzewski. Material spatial randomness: From statistical to representative volume element. Probabilistic Engineering Mechanics 21(2):112–132, 2006. https: //doi.org/10.1016/j.probengmech.2005.07.007 [12] T. Kanit, S. Forest, I. Galliet, et al. Determination of the size of the representative volume element for random composites: statistical and numerical approach. International Journal of Solids and Structures 40(13–14):3647–3679, 2003. https://doi.org/10.1016/S0020-7683(03)00143-4 [13] P. Kabele, O. Švagera, et al. Mathematical modeling of brittle fractures in rock mass by means of the DFN method. Tech. rep., Czech Technical University in Prague, 2017. [14] P. Kabele. DFN_toolbox, 2023. [2023-08-22]. https://gitlab.com/pkabele/DFN_toolbox [15] J. M. Garrard, R. Abedi. Statistical volume element averaging scheme for fracture of quasi-brittle materials. Computers and Geotechnics 117:103229, 2020. https://doi.org/10.1016/j.compgeo.2019.103229 [16] K. Acton, J. Garrard, R. Abedi. Geometric partitioning schemes to reduce modeling bias in statistical volume elements smaller than the scale of isotropic and homogeneous size limits. Computer Methods in Applied Mechanics and Engineering 393:114772, 2022. https://doi.org/10.1016/j.cma.2022.114772 [17] D. Trias, J. Costa, A. Turon, J. E. Hurtado. Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers. Acta Materialia 54(13):3471–3484, 2006. https://doi.org/10.1016/j.actamat.2006.03.042 [18] B. Patzák. OOFEM project home page, 2000. [2023-08-22]. http://www.oofem.org [19] M. Lebeda, P. Kabele. The effect of sampling volume size on the apparent stiffness of jointed rock mass. Acta Poytechnica CTU Proceedings 34:38–42, 2022. https://doi.org/10.14311/APP.2022.34.0038 [20] Z. Bukovská, et al. SÚRAO TZ 464/2020. Získání dat z hlubokých horizontů dolu Rožná [In Czech; Data acquisition from the deep horizons of the Rožná mine]. Tech. rep., SÚRAO – Správa úložišť radioaktivních odpadů, 2020. [21] K. Souček, M. Vavro, L. Staš, et al. SÚRAO ZZ 221/2018. Komplexní geologická charakterizace prostorů PVP Bukov – část II geotechnická charakterizace [In Czech; Complex geological characterization of the Bukov URF – part II geotechnical characterization]. Tech. rep., SÚRAO – Správa úložišť radioaktivních odpadů, 2018. 54 https://doi.org/10.14311/APP.2023.40.0061 https://doi.org/10.1029/wr020i001p00079 http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/38/013/38013872.pdf http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/38/013/38013872.pdf https://doi.org/10.1061/JMCEA3.0002642 https://doi.org/10.1016/0167-6636(92)90004-W https://doi.org/10.1007/BF01023618 https://doi.org/10.1016/j.probengmech.2005.07.007 https://doi.org/10.1016/j.probengmech.2005.07.007 https://doi.org/10.1016/S0020-7683(03)00143-4 https://gitlab.com/pkabele/DFN_toolbox https://doi.org/10.1016/j.compgeo.2019.103229 https://doi.org/10.1016/j.cma.2022.114772 https://doi.org/10.1016/j.actamat.2006.03.042 http://www.oofem.org https://doi.org/10.14311/APP.2022.34.0038 Acta Polytechnica CTU Proceedings 49:48–54, 2024 1 Introduction 2 DFN model of fractures' geometry 3 Averaging procedure 3.1 Statistical volume element 4 The SVE concept verification 5 Numerical simulation of tunnel excavation 5.1 Results of the study 6 Conclusion Acknowledgements References