Influence of the scale effect on the parameters of damage plasticity model for mechanical behaviour of rock mass

Abstract

Rock masses are inhomogeneous, fractured, anisotropic, and initially stressed in their natural state. They consist of intact rock or rock monolith and macro-damages, which include fractures, bedding planes, faults and other forms of discontinuity. It has been experimentally confirmed that with the increase of the rock mass scale, the mechanical properties of the rock mass decrease as a consequence of discontinuities. The subject of this paper is the analysis and determination of the scale effect on the damage plasticity model parameters, which describe the mechanical behaviour of the rock mass for different loads and scales. The scale effect is analysed on a constitutive model that can simulate the most complex mechanical behaviours in rock masses (elasto-plastic with damage in stiffness and strength). The analysis is based on the results of experimental tests of rock mass, which were carried out for different scales (scale of the rock mass sample in laboratory and scale of the in situ test – shear test), while the tests were performed at the same microlocation. For the considered experimental tests, appropriate finite element models (FEMs) were formed with boundary conditions that correspond to the ones from the experiments. By simulating the performed experiments with a series of FEM analyses with complex optimization algorithms, the parameters of damage plasticity model were determined for both values of the rock mass scale, proving the scale effect numerically. Based on the obtained results, analysis and discussion of the scale effect for a given rock mass were performed.

Keywords

Scale effect damage plasticity rock mechanics experimental analysis numerical analysis.

Introduction

Rock masses are inhomogeneous, fractured, anisotropic, and initially stressed geological media. They consist of intact rock mass and various types of discontinuities, such as fractures, bedding planes, faults, and joints. The spatial distribution of these discontinuities is typically irregular, resulting in pronounced anisotropy, heterogeneity, and complex mechanical behaviour. As a result, the mechanical response of a rock mass differs significantly from that of intact rock. One of the key phenomena governing this difference is the scale effect, which describes the dependence of the mechanical properties of a rock mass on the size of the investigated domain. Experimental evidence consistently shows that increasing the scale of observation leads to a reduction in mechanical strength due to the increasing influence of discontinuities.

The scale effect has been extensively investigated through laboratory testing on specimens of different sizes and in situ experiments at different observation scales. Early studies established the fundamental trend of strength reduction with increasing scale. Bieniawski (1968) demonstrated a decrease in uniaxial compressive strength of coal with increasing specimen size. Subsequent research extended these observations from intact rock to rock mass behaviour and discontinuities. Muralha and Da Cunha (1990) compared laboratory shear tests with in situ tests on an arch dam foundation in Portugal and reported significantly lower peak shear strength at the field scale. In a similar direction, Barton (2013), among others, analyzed the scale dependence of shear strength of rock discontinuities, showing that both peak and residual strengths decrease with increasing scale. It has also been shown that shear strength along discontinuities is strongly influenced by surface roughness and infill material, while purely frictional resistance is less sensitive to scale effects. Da Cunha and Muralha (1990) further demonstrated a reduction in deformation modulus with increasing scale.

In parallel with experimental investigations, several authors have proposed empirical and semi-empirical relationships to estimate rock mass strength from intact rock properties and discontinuity characteristics. Goldstein et al. (1966) and Pinto da Cunha (1990), among others, contributed early formulations linking rock mass behaviour to measurable small-scale parameters. More recent approaches incorporate scale effects directly into constitutive formulations. Masoumiet al. (2016) directly incorporated the scale effect into a constitutive model proposed by Christensen (2000) through parameters that correlate uniaxial compressive and tensile strength with sample size. Barton and Bandis (1982) introduced a yield function correction for discontinuities, accounting for the scale of discontinuities in situ. Geomechanical classification systems such as RMR, GSI, and Q are widely used in engineering practice and, although not explicitly scale-dependent, implicitly reflect the scale of engineering structures through their empirical correlations with rock mass properties. Geomechanical classifications are often correlated with failure criteria for rock masses, such as the Hoek-Brown constitutive model or the model proposed by Aubertin, Li and Simon (2000), where parameters from the yield function are correlated with parameters from the GSI classification, indirectly introducing the scale effect into the yield condition. Numerical studies by Gaoet al. (2014), Wang et al. (2016), and Chen et al. (2011) further confirm the reduction of rock mass strength with increasing scale and highlight the role of fracture density and distribution.

(1) Regarding the scale effect, there are different approaches in the case of the analysis of the scale effect in concrete, such as the method presented by Van Mier and Man (2008) and by Yu, Jin and Du (2023), in which different experimental and numerical methods for identifying the scale effect in concrete are presented. Like the scale effect observed in rock masses or concrete, the performance of geotechnical composite materials also significantly depends on the interaction between the microstructure and field conditions. Sun et al. (2026) demonstrated that variations in the composition, gradation, and application conditions of modified synchronous grouting materials lead to significant changes in the mechanical characteristics and system behaviour under real field conditions, indicating a pronounced dependence of performance on system scale and heterogeneity. Han et al. (2025) showed that the mechanical response of the grout–soil interface is significantly influenced by specimen dimensions and interface heterogeneity, with both peak and residual shear strengths decreasing as the specimen height increases.

Despite extensive research, the quantification of scale effects remains challenging due to the complexity of discontinuity networks and the limited availability of consistent multi-scale datasets. This makes it difficult to establish unified mechanical parameters valid across different scales.

(2) The present study aims to identify the parameters of a damage–plasticity model through the application of advanced optimization algorithms integrated with a finite element model (FEM) framework. The FEM model was developed to simulate both laboratory (small-scale) and in situ (large-scale) tests conducted on the same rock mass. The scale effect is investigated using a constitutive model capable of capturing the complex mechanical behaviour of rock masses, including elastoplastic response with degradation of stiffness and strength. The analysis relies on experimental data obtained at different scales, encompassing both rock mass specimens and in situ testing at the same micro-location. For each experimental configuration, corresponding FEM models were established with boundary conditions consistent with the testing conditions. The integration of numerical simulations, experimental results, and optimization procedures enabled the identification of damage–plasticity model parameters for both considered scales. The results provide a basis for evaluating and discussing the scale effect in the analyzed rock mass.

Experimental research

Within this section, experimental tests of the rock mass are presented, based on which numerical and scale effect analysis were performed. The results of tests on the rock mass in the exploratory gallery of the future gravity concrete dam “Pošćenje” in Bosnia and Herzegovina are available, namely in situ shear test and uniaxial compression test with unloading on the rock mass sample. Considering the scale effect analysis, these experiments are representative because they were performed at the same location and corespond to different values of the observation scale.

Uniaxial compression test with unloading

Uniaxial compression test with unloading is performed on the rock mass samples according to the recommendations of ISRM.³¹ The values of compressive strength and modulus of elasticity of the rock mass samples are obtained with this test. Often, testing of rock mass samples is carried out in cycles of loading and unloading with deformation measurements so that elastic and plastic deformations can be separated. At the location in question, the experiment was performed with two cycles of loading and unloading, and in the third cycle the sample was brought to failure. In Figure 1, the layout of the test and loading scheme is presented. In Figure 2 the normal stress – normal strain diagram, which was obtained based on the test results of the rock mass sample for all levels of the specified normal stresses is shown.

Figure 1.

(a) Layout of uniaxial compression test, (b) Load scheme of the rock mass sample under uniaxial cyclic compression with loading to failure.

Figure 2.

Uniaxial cyclic compression test: normal stress – normal strain.

In situ shear test

In situ shear test is performed with the aim of determining mechanical properties of the rock mass for which failure occurs. During the test, normal and shear stress are applied. Load on the rock mass is applied via a concrete block. In Figure 3 the layout and scheme of shear test are given.

Figure 3.

(a) Field shear test setup (photograph) and (b) Shear test configuration (schematic).

At the specific location, the test was performed by first applying normal stress in two loading–unloading cycles, followed by the application of constant shear stress. When normal stress is applied (Figure 4), normal displacements ( $δ_{n}$ ) in the direction of normal stress are measured (Figure 5). During the first two loading – unloading cycles, normal displacements were measured and are shown in the diagram in Figure 5. After the second cycle, the normal stress was set to a value of 1.293 MPa, at which point shearing is performed. During the shearing phase, normal displacements were not measured. In the shear phase, shear displacements ( $δ_{s}$ ) in the direction of shear force were measured (Figure 6). The shear failure moment was determined from the shear stress curves as a function of shear displacements.

Figure 4.

Normal stress loading function.

Figure 5.

Normal stress – normal displacement during concrete-rock shear test (results of the experiment for the first two loading–unloading cycles).

Figure 6.

Shear force – shear displacement obtained during concrete-rock shear test (results of the experiment).

Numerical analysis

Constitutive model for mechanical behaviour of rock mass

For the scale effect analysis and modelling of mechanical behaviour of rock mass within this research, a constitutive model with plastic damage (hereinafter referred to as the DP model) (Lubliner et al., 1989; Lee and Fenves, 1998, 2001; Omidi and Lotfi, 2010) was applied. It was numerically implemented in the software package PAK (Živković et al., 2019a, 2019b). By applying this model, it is possible to simulate the complex material behaviours rock masses like elasto-plastic and damage in stress and stiffness.

In the case of an elasto-plastic model with damage, the relationship between total stress tensor and effective stress tensor is given by the damage variable d. The total tensor stress σ is defined asfollows:

σ = (1 - d) \bar{σ}

(1)

where d represents the damage variable (degradation) and $\bar{σ}$ is the effective stress tensor without damage.

The effective stress is calculated using the elastic constitutive matrix without damage, which depends on the initial modulus of elasticity $E_{0}$ and Poisson's ratio $ν_{0}$ .

During deformation, there is a change in material stiffness as a consequence of the development of plastic deformations and stiffness degradation. For the sake of illustration, the stress–strain curve is shown for uniaxial compression test (Figure 7(a)) and uniaxial tension test (Figure 7(b)), where unloading after reaching the peak strength is also analyzed.

Figure 7.

Behaviour of the DP model during uniaxial compression test (a) and uniaxial tension test (b).

For damage in the model, two internal variables are used for compression $κ_{c}$ (Equation 2) and for tension $κ_{t}$ (Equation 3). In case there is no material damage, the variables $κ_{c}$ and $κ_{t}$ have a value of 0, otherwise the variables are greater than 0 and less than 1 (the value 1 means that complete damage has occurred).

κ_{c} = 1 / g_{c} \int_{o}^{ε^{p}} σ_{c} d ε^{p}

(2)

κ_{t} = \frac{1}{g_{t}} \int_{o}^{ε^{p}} σ_{t} d ε^{p}

(3)

where $g_{c}$ and $g_{t}$ are the fracture energies in compression and tension, defined by Equations 4 and 5

g_{c} = \int_{o}^{\infty} σ_{c} d ε^{P}

(4)

g_{t} = \int_{o}^{\infty} σ_{t} d ε^{p}

(5)

The quantities $g_{c}$ and $g_{t}$ represent the areas under the $σ - ε^{P}$ curve for the uniaxial compression test and the uniaxial tension test, respectively.

In the following text, the general notation $χ \in [[c, t]],$ where c refers to compression and t to tension, will be used.

Based on the relation between stress and plastic deformation given by Lubliner et al. (1989), it is possible, by introducing appropriate exponential approximations, to obtain the relation between the total stress and the internal damage variable $κ_{χ}$ (Equation 6):

σ_{χ} = f_{χ} (κ_{χ}) = f_{χ 0} / a_{χ} [(1 + a_{χ}) \sqrt{φ (κ_{χ})} - φ (κ_{χ})]

(6)

where $f_{χ 0}$ is the initial yield stress (without damage), and function $φ_{χ}$ is defined by Equation 7.

φ_{χ} = 1 + a_{χ} (2 + a_{χ}) κ_{χ}

(7)

where $a_{χ}$ is a parameter that affects the shape of the stress–strain curve during uniaxial compression or uniaxial tension after reaching the yield stress.

Considering the relationship between total and effective stress by introducing degradation, the expression for effective stress can be written in the following form (Lubliner et al., 1989; Lee, 1996) (Equations 8 and 9):

σ_{χ}^{-} = f_{χ}^{-} (κ_{χ}) = f_{χ 0} {[(1 / a_{χ}) (1 + a_{χ} - \sqrt{φ (κ_{χ})})]}^{(1 - c_{χ} / b_{χ})} \sqrt{φ (κ_{χ})}

(8)

d_{χ} = 1 - {[(1 / a_{χ}) (1 + a_{χ} - \sqrt{φ (κ_{χ})})]}^{c_{χ} / b_{χ}}

(9)

where $d_{χ}$ is a variable representing material degradation, and $\bar{σ_{χ}}$ is the effective stress.

The parameters of the DP model, which adjust the shape of the stress–strain curve to the experimental curve, are: $D_{c}$ and $D_{t}$ . The parameter $D_{c}$ is defined as the material degradation for the stress corresponding to the uniaxial compressive strength. The parameter $D_{t}$ is defined as the material degradation for the stress corresponding to half of the uniaxial tensile strength. Based on the adopted parameters $a_{χ}$ , $b_{χ}$ , $c_{χ}$ , the value of the parameter $D_{χ}$ can be obtained by Equations 10 11 (Lee, 1996).

D_{c} = 1 - [(1 / a_{c}) (1 + a_{c} - (1 + a_{c}) / 2)]^{c_{c} / b_{c}}

(10)

D_{t} = 1 - {[(1 / 2 a_{t}) (1 + a_{t} - \sqrt{1 + a_{t}^{2}})]}^{c_{t} / b_{t}}

(11)

For multiaxial problems, the degradation variable (Equation 12) is calculated by interpolating between the compressive and tensile degradation variable values (Lee, 1996; Lee and Fenves, 1998)

d = 1 - (1 - d_{c} (κ_{c})) (1 - d_{t} (κ_{t}))

(12)

In the case of cyclic loading, stiffness recovery during the transition from tension to compression is considered (Equation 13).

d = 1 - (1 - d_{c} (κ_{c})) (1 - s d_{t} (κ_{t}))

(13)

where s is an internal variable that depends on the stress state and the stiffness recovery parameter $s_{0}$ .^25,30

The yield function of the DP model is defined by Equation 14 (Lee, 1996).

F (\bar{σ}, κ) = 1 / (1 - α) (α {\bar{I}}_{1} + \sqrt{3 {\bar{J}}_{2}} + β (κ) {\bar{σ}}_{m a x}) - c_{c} (κ)

(14)

where:

$I_{1}^{-}$ is the first invariant of the stress tensor $\bar{σ}$ , $σ_{m a x}^{-}$ is the algebraic maximum of the principal stresses of the tensor $\bar{σ}$ , ${\bar{J}}_{2}$ is the second invariant of the deviatoric part of the stress tensor $\bar{σ}$ , $c_{c} (κ)$ is the cohesion parameter that depends on the internal damage variable, $α$ is the parameter that regulates the slope of the yield function in the principal stress system. $β (κ)$ is the parameter that is a function of the internal damage variable and defines the yield function for different values of the damage variable. This parameter exists if the maximum principal stress ( $σ_{m a x}^{-}$ ) is greater than 0 and depends on the damage.

According to the research (Lee and Fenves, 1998), parameter $β (κ)$ is defined by Equation 15.

β (κ) = {\begin{matrix} {\bar{σ}}_{c} (κ_{c}) / {\bar{σ}}_{t} (κ_{t}) (1 - α) - (1 + α), ({\bar{σ}}_{m a x} > 0) \\ γ, ({\bar{σ}}_{m a x} \leq 0) \end{matrix}

(15)

The parameter $γ$ defines the shape of the yield surface in the case of triaxial compression (Equation 16). This parameter exists if the maximum principal stress $σ_{m a x}^{-}$ is less than 0.

γ = 3 \cdot (1 - K_{c}) / (2 K_{c} - 1)

(16)

The parameter $K_{c}$ represents ratio between the second invariant of the deviatoric part of the stress tensor in tensile and compressive meridians of the principal stresses space. The limitation is defined by Equation 17.

0.5 < K_{c} \leq 1.0

(17)

The change in plastic deformation is defined by Equation 18.

{\dot{ε}}^{p} = \dot{λ} \partial Q / \partial \bar{σ}

(18)

where:

Q (\bar{σ}) = \bar{S} + α_{p} \cdot I_{1}^{-}

(19)

and represents the plastic potential function (Equation 19) which is a linear function of the mean stress and the stress deviator (Lee, 1996). $\bar{S}$ is the norm of the deviatoric stress tensor $\bar{σ}$ , and $α_{p}$ is the dilatancy parameter.

FEM models for experiment simulation

FEM model for uniaxial compression load–unload test

For numerical simulation of a rock mass compression loading-unloading test, a unit-dimension model consisting of a single tetrahedral finite element was created. One-eighth of the problem was modelled by applying symmetrical boundary conditions across three coordinate planes. The model loading was defined as a vertical stress applied to the free top surface of the model. The pressure was specified in accordance with the loading function shown in Figure 1, while the FEM used in the simulation is presented in Figure 8.

Figure 8.

Finite element model (FEM) model for uniaxial loading-unloading compression test.

By using the developed model, numerical simulations were performed, and the calculated stress–strain relationship, was used for the optimization of the constitutive model parameters as shown in the ‘Dp model parameters’ section. During the numerical simulation, the loading of the model was performed in 35 steps. The Figure 9 ilustrates the stress–strain relationship for arbitrary constitutive model parameters.

Figure 9.

Normal stress–strain diagram obtained based on the results of the numerical simulation of the experiment (finite element model (FEM) model).

Based on shown relationship, it can be observed that the model exhibits linear elastic behaviour up to a certain stress level (initial yield stress), after which nonlinearity occurs due to the development of plastic strain accompanied by material hardening. During unloading, the model behaves linearly elastic, with an evident change in the modulus due to the development of material stiffness degradation.

FEM model for in situ shear test

For the numerical simulation of an in situ large-scale shear test, a FEM was developed to replicate real testing conditions. Тhe fractured rock mass is treated as an equivalent continuum, which implicitly accounts for the combined influence of intact rock material and discontinuities, because the analysis based on in situ large-scale shear tests, in which stresses (normal and shear) are transferred through a concrete block into the rock mass, inducing failure within the tested volume. Due to the scale of the experiment, the measured response reflects the integrated behaviour of both intact rock and discontinuities, thus justifying the use of an equivalent continuum model.

The geometry and dimensions of the numerical model are given in Figure 10, while the corresponding FE model is presented in Figure 11. Due to the symmetry of the problem, only one-half of the domain was modelled by applying boundary conditions of symmetry. The model was discretized using 10-node tetrahedral finite elements. The finite element mesh was strategically refined: a higher mesh density was employed in the rock mass zone where the development of plastic strain is anticipated. Conversely, coarser elements were used towards the model boundaries to optimize computational efficiency. A concrete block with dimensions of 80 × 40x 55 cm was modelled on the upper surface, acting as a loading platen through which the load is transferred to the rock mass. These two components are rigidly connected (nodes are merged), thereby simulating material failure without accounting for potential sliding at the concrete-rock interface. The final FE model consists of 5879 nodes and 3802 finite elements. The boundary conditions were defined by restricting horizontal displacements normal to the planes for all nodes on the vertical boundaries, alongside symmetry constraints, whereas the nodes at the base of the model were fully fixed. The DP constitutive model was utilized to simulate the mechanical behaviour of both materials.

Figure 10.

Mesh geometry for concrete-rock shear test with boundary conditions (distances are in cm).

Figure 11.

Finite element model (FEM) model for concrete-rock shear test (isometry).

The loading was applied in two phases: phase 1 – a vertical compressive load was applied by prescribing a normal stress on the top surface of the concrete block; phase 2 – a shear load was applied via prescribed displacements to the concrete block. These displacements were assigned to a single node, which was linked to the remaining nodes on the top surface of the concrete block. Displacement control was specifically implemented to mitigate numerical instabilities often encountered in the constitutive model's response after reaching peak strength. On the FEM model, the normal stress on the surface A'B'C'D’ is applied over 54 increments in accordance with the loading function during the experiment (Figure 4). After establishing a constant normal stress, shear displacement in the y-axis direction is set in equal increments (27 steps) so that the maximum shear displacement corresponds to the maximum value obtained during the experiment, which is 6.045 mm (see Figure 6). As a result of the numerical simulation of shear test, average values of vertical displacements in the corners of the loaded surface A'B'C'D’ are read from the FEM model for each step of normal stresses. As a result of the shearing phase, for each step of shearing displacement, the total shear force in the y-axis direction is obtained, which represents the resultant of all shear forces in the nodes of the surface ABCD.

The results shown in Figures 12 and 13 refer to arbitrary parameters of the constitutive model, whereas the results obtained following parameter optimization are detailed in the ‘Dp model parameters’ section. The shear stress during numerical simulation was calculated following the experimental methodology: as the ratio of the shear force derived from the numerical simulation to the area of the foundation interface.

Figure 12.

Normal stress – vertical displacement for two normal stress loading cycles (finite element model (FEM) model).

Figure 13.

Shear stress–shear displacement from sharing phase of the finite element model (FEM) model.

Based on the shear results obtained from the numerical simulation, it is observed that, for the given shear displacements, the DP model simulates material damage upon reaching the peak shear strength, when degradation occurs, which is followed by a reduction in the shear strength to the residual value.

Algorithm for parameter identification

To define the mechanical behaviour of the constitutive DP model, it is necessary to determine a large number of numerical parameters. Due to their physical meaning, these parameters cannot be directly obtained from experiments, as there are no in situ results that can be directly applied. Therefore, appropriate numerical methods must be used. The proposed algorithm for determining the parameters of the DP model in this paper is based on the use of the FEM model (Section 3.2) and corresponding experimental results (Section 2), one or more combinations of parameters are determined based on which the smallest deviations are obtained between the results from numerical using FEM models and the results obtained from experimental testing.

This process involves performing a large number of numerical simulations and searching for the optimal parameter combination that ensures agreement between experimental test results and FEM model. The problem is formulated as an optimization problem and issolved by using optimization algorithms. The optimization problem is generally defined as follows: It is necessary to determine the minimum of the function $K_{s} (x), s = 1, 2, \dots, S$ so that the following relations are satisfied:

f_{j} (x) \leq 0 j = 1, 2, \dots, J

(20)

l_{k} (x) = 0 k = 1, 2, \dots, K

(21)

x_{i}^{L} \leq x_{i} \leq x_{i}^{U} i = 1, 2, \dots, n

(22)

Solving the optimization problem involves finding an $n$ -dimensional vector $x = \vec{x} = (x_{1}^{*}, x_{2}^{*}, \dots, x_{n}^{*})$ that satisfies conditions referred to by Equations 20, 21 and 22 and optimizes S objective functions (criteria) $K_{1}, K_{2}, \dots, K_{S}$ . The condition is defined by Equation 22 represents the search area for values that are coordinates of a vector $\vec{x}$ (the vector coordinates are the parameters that are the subject of optimization). The value $x_{i}^{L}$ represents the lower and $x_{i}^{U}$ the upper limit of the search area for parameter i. The solution of the problem is determined by both J in equations and K equations. Functions $f_{j} (x)$ and $l_{k} (x)$ are called restriction functions, and they define the admissible search space.

Given that the determination of parameters is a systematic search for the optimal values of a set of parameters, the necessary steps must be defined. In addition, it is determined which steps can be automated and to what extent. The general steps of the optimization algorithm are:

selection of parameter values from the search space,

performing a numerical simulation using the FEM model,

comparison between the results obtained using the FEM model with the experimental results based on optimization criteria,

correction of parameters and, if necessary, re-running the simulation.

The first step is to identify the set of parameters that are significant for optimization from the widest set of parameters. The set of parameters obtained through the optimization process is determined based on the parameter sensitivity analysis (see Chapter 4.3). Based on the results of the sensitivity analysis, the parameters to be included in the optimization algorithm are defined.

After defining the set of parameters, the search space for the parameters in the optimization algorithm is defined. Calculations using the FEM model, comparison of numerical results with experimental data and evaluation of the objective function are repeated iteratively until the minimum value of the objective function is reached. Numerical simulations are carried out each time using a new combination of parameters. These steps in the optimization algorithm can be automated.

Automatic parameter determination represents a mixed optimization problem aimed at minimizing the difference between measured and calculated values. It is based on a large number of possible solutions, where a wide set of equations and inequations must be taken into consideration, as well as logical restrictions defined within the optimization algorithm. Therefore, to solve the mentioned optimization problem, an automatic optimization approach was applied, which supports multi-objective optimization and is based on algorithms that can be parallelized.

Objective functions (criteria) are defined for each experiment separately.

For the uniaxial compression test with unloading, the criterion function is defined by Equation 23.

K_{cp} (p_{1}, p_{2}, \dots, p_{k}) = \sum_{i = 1}^{p} (| {\bar{δ}}_{v, i} - δ_{v, i} |) / p {\bar{δ}}_{v, i}

(23)

where:

$p_{k}$ is the combination of the DP model parameters from the search area ( $k$ represents the number of parameters that are the subject of optimization, which was determined based on the sensitivity analysis for the given experiment) $δ_{v, i}$ is the calculated displacement value at the ith step of the simulation, ${\bar{δ}}_{v, i}$ is the measured displacement value at the ith step of the simulation, p is the total number of simulation steps.

In the shear test, two criteria are defined: one refers to the minimization of results deviation due to normal stresses, and the other refers to the minimization of the deviation between the results obtained in the shearing phase.

The first criterion is defined by Equation 24.

K_{s k} (p_{1}, p_{2}, . ., p_{k}) = \sum_{i = 1}^{p} r_{i} (| {\bar{δ}}_{v, i} - δ_{v, i} |) / p {\bar{δ}}_{v, i}

(24)

where $p_{k}$ is the combination of the DP model parameters from the search area ( $k$ represents the number of parameters that are the subject of optimization, which was determined based on the sensitivity analysis for the given experiment), ${\bar{δ}}_{v, i}$ is the measured normal displacement in the consolidation phase at the ith step of the numerical simulation (obtained by the experiment), $δ_{v, i}$ is the calculated normal displacement in the consolidation phase at the ith step of the numerical simulation (obtained by the FEM model), p is the number of numerical steps of the consolidation phase, $r_{i}$ is the weighting factor that takes into account the importance of individual steps of the simulation (weighting factor value ranges from 0 to 1). The weighting factor is based on an empirical assessment of the influence of individual steps on the specific test and on the intended use of the test results. In this study, the weighting factor value 1.0 was used for all simulation steps.The second criterion is defined by Equation 25.

K_{s s} (p_{1}, p_{2} \dots . p_{k}) = \sum_{i = m}^{n} r_{i} (| {\bar{T}}_{z y, i} - T_{z y, i} |) / ((n - m + 1) {\bar{T}}_{z y, i})

(25)

where $p_{k}$ is the combination of the DP model parameters from the search area ( $k$ represents the number of parameters that are the subject of optimization, which was determined based on the sensitivity analysis for the given experiment), ${\bar{T}}_{z y, i}$ is the measured shear force in the shearing phase at the ith step of the simulation (obtained by the experiment), $T_{z y, i}$ is the calculated shear force in the shear phase at the ith step of the simulation (obtained by the FEM model), n is the total number of simulation steps, m is the ordinal step number at which the shearing phase begins, $r_{i}$ is the weighting factor that takes into account the importance of individual steps of the simulation (ponder value ranges from 0 to 1). The weighting factor is based on an empirical assessment of the influence of individual steps on the specific test and on the intended use of the test results. In this study, the ponder value 1 was used for all simulation steps.

It can be observed that in some tests, a single criterion is defined, while in others, two optimization criteria are used. In the case of a single criterion, the parameter combination from the admissible search space that minimizes the proposed criterion is identified.

In the case of two or more criteria (multi-criteria optimization), there is a set of compromise solutions called the Pareto front (Deb et al., 2002). The set of Pareto optimal solutions contains solutions that are superior to all other solutions from the search space for each criterion, but are inferior to other solutions from this set for other criteria (Srinivas and Deb, 1994). Based on the set of compromise solutions, the decision maker chooses that solution which he considers optimal, taking into account all criteria. In this paper, the solution comprehends a combination of the DP model parameters. To solve the optimization problem of single-criteria and multi-criteria optimization, the parallelized NSGA-II algorithm was used (Deb et al., 2002).

Sensitivity analysis

Sensitivity analysis methodology

Sensitivity analysis shows how much variations of certain parameters affect mechanical behaviour of the rock mass under certain loads. This analysis gives a hierarchical influence of certain parameters on the specific behaviour of the rock mass. This means that certain parameters can have a significantly greater influence on the rock mass behaviour under certain loads compared to other parameters. Based on sensitivity analysis, the number of parameters is reduced, and a set of important parameters is obtained to which more attention is paid in terms of the precise determination of their values. Other parameters are adopted with less regard for their accuracy, because they have no significance on the rock mass behaviour under a certain load. The sensitivity analysis in this research was carried out for the uniaxial cyclic compression test and for the shear test. DAKOTA programme (Adams et al., 2014) and Morris One at a Time (MOAT) (Morris, 1991) methods were used for sensitivity analysis. The MOAT method is intended for computer model research in terms of categorization of input parameters. The idea is to differentiate between the input parameters by the level of influence on the output variable.

According to the MOAT method, each parameter (dimension) of the $k$ -dimensional space of input parameters is evenly divided into p levels, creating a grid of $p^{k}$ points $x \in R^{k}$ , on which the output variable can be evaluated by the simulation. The elementary effect of parameter change is calculated using the finite-difference method (Equation 26).

d_{i} (x) = (y (x + Δ e_{i}) - y (x)) / Δ

(26)

where $e_{i}$ , $i$ th coordinate vector, and for Δ, which represents differential, a higher number is usually taken, considering that the ultimate intention is not to estimate the differential at a given point, but to estimate the sensitivity at the global level.

The idea of the MOAT method is to estimate the distribution of the elementary effects $d_{i}$ that arise as an effect of the input i on the output variable. The distribution obtained by a series of r simulations can be characterized by a modified average (Equation 27) from r samples³⁴ that shows the sensitivity of a certain parameter in relation to other parameters.

μ_{i} = 1 / r * \sum_{j = 1}^{r} | d_{i}^{(j)} |

(27)

The variable $d_{i}$ in this case represents the elementary effect of a certain variable change that is calculated during numerical tests using the FEM model.

The sensitivity analysis of the DP model parameters was carried out on the FEM models of uniaxial cyclic compression and shear tests. As a result of the sensitivity analysis, hierarchically distributed parameters are obtained according to the effect on the output variables in certain experiments.

Results of the sensitivity analysis of uniaxial cyclic compression test

Table 1 provides the search areas for each of the DP model parameters for the sensitivity analysis. Figure 14 shows the sensitivity indicators (mean value) of all analyzed parameters.

Figure 14.

Indicators of parameter sensitivity to total deformation during uniaxial cyclic compression test.

Table 1.

Parameter search area for the sensitivity analysis of uniaxial cyclic compression test.

Parameter	Lower Limit	Upper Limit
$E_{0} (kPa)$	4000,000	8000,000
$ν$	0.2	0.49
$f_{c} (kPa)$	1900	2500
$f_{t} (kPa)$	40	100
$a_{c}$	1	10
$D_{c}$	0.05	0.6
$a_{t}$	0.05	0.6
$D_{t}$	0.05	06
$G_{c} (kN / m)$	15	45
$G_{t} (kN / m)$	3	10
$α_{p}$	0.05	0.25
$α$	0.05	0.25
$γ$	4	8
$s_{0}$	0	0.95

It can be concluded that parameters $E_{0}$ , $f_{c}$ , $a_{c}$ , $D_{c}$ and $G_{c}$ have the greatest influence on the total deformation obtained during the simulation of the uniaxial cyclic compression test. This conclusion is logical because the above parameters, according to the definition of the DP model, define the material behaviour under uniaxial compression. The parameters related to tension ( $f_{t}$ , $a_{t}$ , $D_{t}$ , $G_{t}$ ) are not important because in uniaxial compression test, tension stress does not occur predominantly. Parameters: $α_{p}$ , $α$ and $γ$ obviously have no influence on the total deformation because they have influence in complex stress states, unlike the linear stress state occurring in uniaxial compression test. Poisson's coefficient also has no influence on the calculation results, because the modulus of elasticity has the greatest influence on the total deformation in the z-axis direction.

Results of the sensitivity analysis of shear test in situ

Table 2 gives the search areas for each of the parameters of the DP model for the purposes of the sensitivity analysis of shear test.

Table 2.

Parameter search area for the sensitivity analysis of shear test.

Parameter	Lower limit	Upper limit
$E_{0} (kPa)$	2000,000	7000,000
$ν$	0.20	0.29
$f_{c} (kPa)$	1500	6500
$f_{t} (kPa)$	50	500
$a_{c}$	2	5
$D_{c}$	0.01	0.8
$a_{t}$	0	1
$D_{t}$	0.01	0.8
$G_{c} (kN / m)$	20	1000
$G_{t} (kN / m)$	5	70
$α_{p}$	0.01	0.14
$α$	0.05	0.3
$γ$	3.0	7.0
$s_{0}$	0	1

Indicators of the individual parameter influence on the shear test simulation outcomes are presented graphically in Figure 15, from highest impact parameters, to the lowest ones.

Figure 15.

Influence of parameters on shear force in shear tests.

The influence of individual parameters is analysed through the mean value. It can be seen from Figure 15 that the parameters that determine material shear strength have the greatest influence on the output variable. In this case, these are parameters $f_{c}$ , $γ$ , and $α$ . The initial modulus of elasticity also has great importance for the output variable because before reaching plastic deformations, this parameter directly affects shear stress values. The influence of dilatancy is evident because dilatancy is determined based on shear tests. Other parameters have less influence on the output variable.

Dp model parameters

DP model parameters for uniaxial compression test scale

Based on the sensitivity analysis results (Section 4.3.2), the optimization parameters were adopted (Table 3). Out of 14 parameters of the DP model, 5 parameters were selected that have the greatest influence on the total deformation obtained from uniaxial compression test with unloading. Other parameters do not have a significant impact on the output variable during uniaxial compression test with unloading and are not the subject of optimization. Based on the conducted optimization analysis, using the FEM model for uniaxial compression test with unloading (Section 4.1.1) and experimental results, and the proposed search areas (Table 3) for certain parameters, the parameters that give the lowest value of the optimization criteria were adopted (Table 3).

Table 3.

DP model parameters which are the subject of optimization for uniaxial compression test with unloading.

Parameter	Parameter search area	Adopted parameter value from optimization
$E_{0} (kPa)$	24,000,000–35,000,000	28,000,000
$f_{c} (kPa)$	47,000–50,000	47,000
$a_{c}$	15–50	55
$D_{c}$	0–1	0.35
$G_{c} (kN / m)$	100–200	215

DP model parameters for in situ shear test scale

Based on the results of the sensitivity analysis (Section 4.3.3), the parameters that are the subject of optimization were adopted (Table 4). From 14 parameters of the DP model, 10 parameters were selected that have the greater influence on the mean deformation under normal stress and on the shear force in the shearing phase. Search areas (Table 4) in the optimization algorithm are given for each parameter that is the subject of optimization.

Table 4.

DP model parameters which are the subject of optimization.

Parameters	Parameter search area	Adopted parameter value from optimization
$E_{0} (kPa)$	18,000,000–30,000,000	27,900,000
$ν$	0.2–0.3	0.21
$f_{c} (KPa)$	500–1500	622
$a_{c}$	60–95	70
$D_{c}$	0.01–0.9	0.15
$G_{c} (kN / m)$	2–4	3.2
$G_{t} (kN / m)$	0.05–0.30	0.17
$α_{p}$	0.01–011	0.018
$α$	0.30–0.7	0.45
$γ$	1–5	2.57

The result of the optimization analysis is the values of criteria (objective functions) for different combinations of parameters, which are displayed through the Pareto front of optimal solutions (Figure 16). The horizontal axis shows the values of criterion 1, which refers to normal displacements under the action of normal stress. The vertical axis shows the values of criterion 2, which refers to shear forces during shear test. From the Pareto front, the combination of parameters that gives the lowest value according to both criteria is adopted. The authors adopted the solution with criterion values $K_{1} = 0.216431$ and $K_{2} = 0.152934$ . This point lies in the region of points closest to the coordinate origin.

Figure 16.

Pareto front of optimal solutions.

Results and discussion

For the adopted values of the DP model parameters, using the FEM model of uniaxial compression test with unloading, the results of compressive stress – total strain were obtained. These results are shown alongside the experimental results for the same compressive stress levels (Figure 16). The average deviation of the experimental results from the FEM ones is 16.9%.

For the adopted parameters, results of the shear test simulation using the FEM model were calculated. Figure 17 shows the comparison of normal displacements during the experiment and during the numerical simulation for the same values of normal stresses. In Figure 18, the comparison of shear forces during the experiment and during the numerical simulation for the same values of shear displacements is presented.

Figure 17.

Uniaxial compression test with unloading: normal stress – total strain (experimental and finite element model (FEM) results for the adopted DP model parameters).

Figure 18.

Shear test: normal stress – normal displacement from the experiment and the finite element model (FEM) model (comparative presentation).

Table 5 shows the DP model parameters obtained for the sample scale and for the in situ scale of shear test.

Table 5.

Comparative presentation of the obtained parameters for different scales of the rock mass.

Parameter	Uniaxial load–unload compression test	Shear test
$E_{0} (kPa)$	28,000,000	27,900,000
$f_{c} (kPa)$	47,000	622
$a_{c}$	55	70
$D_{c}$	0.35	0.15
$G_{c} (kN / m)$	215	3.2

Based on the results shown in Table 5, ratio between the parameters obtained from uniaxial cyclic compression test and the shear test are shown (Figures 19 and 20).

Figure 19.

Shear test: shear force – shear displacement from the experiment and the finite element model (FEM) model (comparative presentation).

Figure 20.

Ratio between parameters $f_{c}$ and $G_{c}$ obtained from uniaxial cyclic compression test and shear test.

Based on the results obtained, it can be said that certain differences are observed for the analyzed parameters depending on the scale size (Figure 21):

Figure 21.

Ratio between parameters $E_{0}$ , $a_{c}$ and $D_{c}$ obtained from uniaxial cyclic compression test and shear test.

The uniaxial compressive strength $f_{c}$ obtained from uniaxial cyclic compression on the rock mass sample is 75.5 times higher than the value of this parameter obtained from shear test, which evidently shows the influence of the scale effect on the uniaxial compressive strength.

Parameter $a_{c}$ , which defines the beginning of plastic deformation of the rock mass, has a value obtained from a uniaxial cyclic compression test equal to 0.79 of the value from shear test. A higher value of this parameter allows the beginning of plastic deformation to occur at lower stress levels. In this case, it means that for the rock mass that is of the in situ test scale, the beginning of plastic deformation occurs at lower stresses than in the case of the rock mass sample. It can be said that the influence of the scale effect on this parameter is evident. Furthermore, this can be interpreted by the presence of fractures for the in situ scale in relation to the rock mass sample which has practically no fractures.

The parameter $D_{c}$ , which defines degradation at the maximum uniaxial compressive stress, has a 2.3 times higher value obtained from the uniaxial cyclic compression test than the value obtained from the shear test. A lower value of this parameter shows that the rock mass reaches plastic deformations earlier, that is, for a zero value of this parameter, elastic deformations do not occur. The aforementioned analysis shows that during in situ tests, the development of plastic deformations occurs earlier in comparison to uniaxial compression test, which is of a smaller scale. Thus, the influence of the scale effect on the given parameter can be observed.

A pronounced influence of the scale effect on the $G_{c}$ parameter is observed. It has a 67.2 times higher value in the uniaxial cyclic compression test compared to the value obtained from the shear test. This is a consequence of the uniaxial strength of the rock mass sample being many times higher compared to the in situ scale.

It is observed that practically the same value of the initial modulus of elasticity of the rock mass $E_{0}$ was obtained for both scale sizes. During the load action on the rock mass fractures close, and after that the material behaves according to the linear elastic law. In the case of rock masses that are more fractured, there is an increase in irreversible deformations compared to those rocks that are less fractured, but the elastic characteristics practically do not change. It can be said that the scale effect on the modulus of elasticity of the rock mass was not observed in this particular case.

Figures 22 to 26 show the relationships between the values of the identified parameters obtained from the uniaxial cyclic compression test and the shear test, and the corresponding scale values. The characteristic scale for the laboratory test is 16.61 cm², while for the shear test it is 6400 cm².

Figure 22.

Dependence of value of the $E_{0}$ value on the scale value (laboratory and shear test).

Figure 23.

Dependence of value of the $f_{c}$ value on the scale value (laboratory and shear test).

Figure 24.

Dependence of value of the $a_{c}$ value on the scale value (laboratory and shear test).

Figure 25.

Dependence of value of the $D_{c}$ value on the scale value (laboratory and shear test).

Figure 26.

Dependence of value of the $G_{c}$ value on the scale value (laboratory and shear test).

As for the remaining parameters of the DP model that are not listed here, some of them were the subject of optimization analysis for the shear test. However, for the uniaxial cyclic compression test, they did not have a significant impact on the results and were not the subject of optimization. For this reason, it is not possible to compare and interpret those parameters for the considered scale sizes.

Analysis of the influence of the scale effect on the mentioned parameters could be carried out if there were laboratory test results for a complex stress state, for example, a triaxial test. In that case, the parameters $α, γ$ and $f_{t}$ would be obtained from a laboratory triaxial test, for example, which could be compared with the results obtained from in situ test. If transverse deformations are also measured with uniaxial cyclic compression, Poisson's ratio would have an influence in that case. This could be compared and discussed with the results obtained during the shear test.

Conclusion

The applied DP model has numerically adequately described the mechanical behaviour of rock mass. The uniaxial unloading test and shear test were simulated using the DP model. The obtained results are consistent with the observed behaviour of the rock mass during the experiments. A methodology for identifying the parameters of the DP model, based on experimental results, numerical FEM models and optimization algorithms was proposed and successfully applied to experimental tests conducted in the research gallery at the “Pošćenje” dam, at different scales. The obtained parameters of the DP model show good agreement between the experimental results and the numerical models. By implementing an optimization algorithm for calibrating the DP material parameters within the FEM framework, a clear size-effect was obtained consistent with the results obtained of both laboratory and in situ experimental tests. Based on the results obtained for the two scales of the rock mass, a strong influence of the scale on the parameters $f_{c}$ , $G_{c}$ , $a_{c}$ and $D_{c}$ was determined. It was observed in the considered case that the scale effect has no influence on the values of the initial modulus of elasticity. The application of the research results is multiple, and it refers, among other things, to the use of the DP model for mechanical behaviour of rock mass during stress–strain analyses of building structures. The applied constitutive model is capable to simulate complex forms of stress–strain behaviour of rock mass, except for time-dependent effects.

Footnotes

ORCID iDs

Radovanović Slobodan

Rakić Dragan

Author contributions

Slobodan Radovanović: Literature review of previous studies, methodology, analysis of experimental results, sensitivity analysis, analysis and discussion of results, manuscript writing.

Dejan Divac: Methodology concept, analysis and discussion, supervision.

Dragan Rakić: FEM analysis.

Dragoslav Šumarac: Supervision and review.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data availability statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

References

Adams

Ebeida

Eldred

, et al. (2014) Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis (No. SAND2014-4633). Albuquerque, NM (United States): Sandia National Lab.(SNL-NM).

Aubertin

Simon

(2000) A multiaxial stress criterion for short-and long-term strength of isotropic rock media. International Journal of Rock Mechanics and Mining Sciences 37(8): 1169–1193.

Barton

(2013) Shear strength criteria for rock, rock joints, rockfill and rock masses: Problems and some solutions. Journal of Rock Mechanics and Geotechnical Engineering 5(4): 249–261.

Barton

Bandis

(1982) Effect of block size on the shear behavior of jointed rocks. In: 23rd U.S. Symposium on Rock Mechanics (USRMS), Berkley, California, August 1982. ARMA-82:739-760.

Barton

Lien

Lunde

JJRM

(1974) Engineering classification of rock masses for the design of tunnel support. Rock Mechanics 6: 189–236.

Bieniawski

(1968) In situ strength and deformation characteristics of coal. Engineering Geology 2(5): 325–340.

Bieniawski

(1989) Engineering Rock Mass Classification. New York, USA: John Wiley & Sons.

Bieniawski

Bernede

(1979) Suggested methods for determining the uniaxial compressive strength and deformability of rock materials: Part 1. Suggested method for determining deformability of rock materials in uniaxial compression. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 16(2): 138–140. Pergamon.

Brown

Hoek

(1980) Underground Excavations in Rock. CRC Press.

10.

Chen

Yang

Tan

, et al. (2011) Study on mechanical parameters of fractured rock masses. Science China Technological Sciences 54: 140–146.

11.

Christensen

(2000) Yield functions, damage states, and intrinsic strength. Mathematics and Mechanics of Solids 5: 285–300.

12.

Da Cunha

(1990) Scale effects in rock mechanics. In: Proceedings of the First International Workshop on Scale Effects in Rock Masses (ed AP Da Cunha), Loen, Norway, 7-8 June 1990, pp.3-27. Rotterdam: Published for the ISRM Commission on Scale Effects in Rock Mechanics by A.A. Balkema.

13.

Da Cunha

Muralha

(1990) About LNEC experience on scale effects in the deformability of rock masses. In: Proceedings of the First International Workshop on Scale Effects in Rock Masses (ed AP Da Cunha), Loen, Norway, 7-8 June 1990, pp.219-229. Rotterdam: Published for the ISRM Commission on Scale Effects in Rock Mechanics by A.A. Balkema.

14.

Deb

Pratap

Agarwal

, et al. (2002) A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on Evolutionary Computation 6(2): 182–197.

15.

Gao

Stead

Kang

(2014) Numerical investigation of the scale effect and anisotropy in the strength and deformability of coal. International Journal of Coal Geology 136: 25–37.

16.

Goldstein

Goosev

Pvrogovsky

, et al. (1966) Investigation of mechanical properties of cracked rock. In: 1st ISRM Congress, Lisbon, Portugal, 25 September −1 October, ISRM-1CONGRESS-1966-089:521-524.

17.

Gutierrez

Youn

(2015) 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.

18.

Han

Chen

, et al. (2025) Experimental study of the shear mechanical properties of the grout-soil interface considering grout permeation. Journal of Geotechnical and Geoenvironmental Engineering 152(1): 04025169.

19.

Hoek

Brown

(1997) Practical estimates of rock mass strength. International Journal of Rock Mechanics and Mining Sciences 34(8): 1165–1186.

20.

Hoek

Carranza-Torres

Corkum

(2002) Hoek-Brown failure criterion-2002 edition. Proceedings of NARMS-TAC Conference Vol 1 2002:267-273. www.rocscience.com.

21.

Ivars

Pierce

Darcel

, et al. (2011) The synthetic rock mass approach for jointed rock mass modelling. International Journal of Rock Mechanics and Mining Sciences 48(2): 219–244.

22.

Kujundžić

(1977) Basics of Rock Mechanics I. Belgrade: Society of Civil Engineers and Technicians of Yugoslavia (in Serbian).

23.

Lee

(1996) Theory and implementation of plastic-damage model for concrete structures under cyclic and dynamic loading. Berkeley: University of California.

24.

Lee

Fenves

(1998) A plastic-damage concrete model for earthquake analysis of dams. Earthquake Engineering & Structural Dynamics 27(9): 937–956.

25.

Lee

Fenves

(2001) A return-mapping algorithm for plastic-damage models: 3-D and plane stress formulation. International Journal for Numerical Methods in Engineering 50(2): 487–506.

26.

Lubliner

Oliver

Oller

, et al. (1989) A plastic-damage model for concrete. International Journal of Solids and Structures 25(3): 299–326.

27.

Marinos

Hoek

(2007) The Geological Strength Index (GSI): a characterization tool for assessing engineering properties for rock masses. In Proceedings International Workshop on Rock Mass Classification for Underground Mining (ed. Mark, Pakalnis and Tuchman), Information Circular (Vol. 9498, pp. 87-94).

28.

Masoumi

Saydam

Hagan

(2016) Incorporating scale effect into a multiaxial failure criterion for intact rock. International Journal of Rock Mechanics and Mining Sciences 83: 49–56.

29.

Morris

(1991) Factorial sampling plans for preliminary computational experiments. Technometrics 33(2): 161–174.

30.

Muralha

Da Cunha

(1990) About LNEC experience on scale effects in the mechanical behaviour of joints. In: Proceedings of the First International Workshop on Scale Effects in Rock Masses (ed AP Da Cunha), Loen, Norway, 7-8 June 1990, pp.77-88. Rotterdam: Published for the ISRM Commission on Scale Effects in Rock Mechanics by A.A. Balkema.

31.

Natau

(1990) Scale effects in the determination of the deformability and strength of rock masses. In: Proceedings of the First International Workshop on Scale Effects in Rock Masses (ed AP Da Cunha), Loen, Norway, 7-8 June 1990, pp.77-88. Rotterdam: Published for the ISRM Commission on Scale Effects in Rock Mechanics by A.A. Balkema.

32.

Omidi

Lotfi

(2010) Finite element analysis of concrete structures using plastic-damage model in 3-d implementation. International Journal of Civil Engineering 8(3): 187–203.

33.

Srinivas

Deb

(1994) Muiltiobjective optimization using nondominated sorting in genetic algorithms. Evolutionary Computation 2(3): 221–248.

34.

Sun

Han

Zhao

, et al. (2026) Sustainable utilization of shield waste muck for synchronous grouting: Mixture ratio optimization and assessment of environmental and economic performance. Tunnelling and Underground Space Technology 168: part 2.

35.

Van Mier

JGM

Man

H-K

(2008) Some notes on microcracking, softening, localization, and size effects. International Journal of Damage Mechanics 18(3): 283–309.

36.

Wang

Yang

, et al. (2016) Numerical analysis on scale effect of elasticity, strength and failure patterns of jointed rock masses. Geosciences Journal 20: 539–549.

37.

Jin

(2023) Experiment and meso-scale modelling on combined effects of strain rate and specimen size on uniaxial-compressive failures of concrete. International Journal of Damage Mechanics 32(5): 683–714.

38.

Živković

Kojić

Slavković

, et al. (2019a) PAK-S - Program for FE Structural Analysis of Solids and Structures. Published online. Kragujevac, Serbia.

39.

Živković

Kojić

Slavković

, et al. (2019b) User Manual and Examples for PAK-S – Program for FE Structural Analysis of Solids and Structures. Faculty of Engineering, University of Kragujevac.