Abstract
A coupled elastic-plasticity-damage constitutive model, AK Model, is applied to predict fracture propagation in rocks. The quasi-brittle material model captures anisotropic effects and the distinct behavior of rocks in tension and compression. Calibration of the constitutive model is realized using experimental data for Carrara marble. Through the Weibull distribution function, heterogeneity effect is captured by spatially varying the elastic properties of the rock. Favorable comparison between model predictions and experiments for single-flawed specimens reveal that the AK Model is reliable and accurate for modelling fracture propagation in rocks.
1. Introduction
An understanding of fracture initiation and propagation in rocks is important in reservoir (hydrocarbon and geothermal) management and energy exploration activities including enhanced geothermal systems, oil and gas extraction, and underground water transport. This has triggered numerous research efforts through experimental studies and numerical modelling of rock fractures. Numerical modelling could be very challenging especially when other field observations like effects of temperature and confining pressure are to be coupled with the fracturing process [1, 2].
While some numerical models assume an elastic rock behavior, other models try to characterize rock behavior into the plastic region. Assumptions such as material isotropy are generally made for simplicity but relaxing such assumptions is necessary to obtain more accurate results. Some models rely on predefining crack locations [1] or extension of preexisting cracks. The coupled damage-plasticity model originally developed by Cicekli et al. [3] and later modified by Abu Al-Rub and Kim [4], referred to as AK Model in this work, captures crack initiation and propagation without the need to predefine crack locations a priori [4]. The model is based on a continuum damage mechanics approach and generally applicable to quasi-brittle materials. Many such models exist with varying degrees in the extent to which important physical phenomena are captured [5–8].
Another crucial consideration in carbonate reservoirs is the heterogeneity of constituent rocks. Studies have shown that the carbonate reservoirs in the Middle-East are very heterogeneous in terms of rock types [9]. Reservoir heterogeneity is also important in carbon capture and sequestration [10]. Moreover, the compressive strength of carbonates has been shown to be a function of the grain size, porosity, and elastic modulus [11]. Hence, for proper reservoir characterization and economic evaluation, fracture propagation, damage, and heterogeneity should all be taken into account. Herein lies the uniqueness of the AK Model. Even with material models abound for describing rock behavior, many do not combine the effects of both damage and heterogeneity and some that do are limited to isotropic and/or elastic properties [12]. This work presents a computational model that considers (1) the plastic deformation in addition to elastic deformation in rocks; (2) damage localization; (3) anisotropy in rock behavior; and (4) heterogeneity of rock sample/unit.
This study focuses on qualitative validation of the model for predicting fracture propagation in rocks and the incorporation of heterogeneity effects. A brief look at previous works on rock fracture propagation is presented next (in Section 2). Subsequently (Section 3), the constitutive model is concisely described without detailing its numerical implementation. Using Carrara marble as a representative rock material, the model is then calibrated (Section 4); and with experiments on single-flawed rock specimens, the capability of the model in predicting fracture patterns is demonstrated in Section 5 (model validation). In Section 6, some methods of incorporating the effects of heterogeneity in modelling geological media are highlighted before the spatial variation of elastic properties is used to study heterogeneity effects through the Weibull distribution function.
2. Fracture Prediction in Rocks
Crack initiation and propagation in rocks have been the subject of various past and current studies. The approach is to study propagation of cracks from precracked (flawed) specimens using experiments (e.g., [13–17]) and numerical simulations (see [18–21]). Bobet and Einstein [21, 22] used both experiments and numerical simulations to study crack initiation, propagation, and coalescence in rock. They obtained good simulation results with the upgraded code FROCK which is based on the Displacement Discontinuity Method (DDM). A later work by Gonçalves Da Silva and Einstein [23] improved the capabilities of the code by introducing a new strain-based criterion as well as a normal stress-dependent criterion for crack development.
AUTODYN is a nonlinear hydrodynamics code compatible with ANSYS; Li and Wong [24] used it to study the influence of flaw inclination angle and loading condition on crack propagation. Tang [12] developed the 2D finite element code, RFPA 2D. The code incorporates the effect of both damage and heterogeneity and its extended 3D version has been used successfully in predicting fracture of specimens in triaxial compression [19]. In another study, Zhang and Wong [25] used the Bonded Particle Method (BPM) to investigate the effect of loading rate on crack behavior of flawed specimens. Crack coalescence mode was observed to change from that dominated by the tensile segment to that dominated by the shear-band.
FROCK is based on DDM which is a boundary element method (BEM) and according to Khair et al. [26], the Finite Element Method (FEM) is superior to BEM in predicting subsurface fractures. While the AUTODYN-based model adopted by [24] captures damage like the RFPA, it relies on Drucker-Prager criterion which has been proven to overestimate intact rock strength [27]. Additionally, RFPA does not consider anisotropy and plastic deformation in rocks. This calls for an elastoplastic-damage model based on FEM. The model by Abu Al-Rub and Kim [4], AK Model, for plain concrete considers damage effects, anisotropy, and plastic deformation. It also adopts the Lubliner yield criterion which is an improvement on the Drucker-Prager criterion. This model is implemented as a user material subroutine in Abaqus, a commercial finite element code.
It is noteworthy that Linear Elastic Fracture Mechanics (LEFM) has also been used to simulate crack propagation in rocks (e.g., see [29] which uses FRANC2D crack propagation simulator). However, LEFM is not able to predict fracture initiation and when a crack is assumed, the stress intensity and fracture toughness on which LEFM is based could be meaningless for the assumed flaw size [30]. Also, the material behavior near the crack tip region (fracture process zone) could be inelastic and nonlinear [30] making LEFM only applicable when the size of process zone (L) is significantly small with respect to the smallest critical dimension of the structure (D) – D/L > 100 [31]. This influenced the decision of [32] in their elastic Abaqus analysis to take stress measurements at a distance from the flaw tip, that is, to avoid stresses in the process zone which make no physical sense. While cohesive zone modelling has been introduced to address these issues in materials like concrete [33], increased number or complexity of fractures might not be easily handled by cohesive zone modelling. The continuum damage mechanics (CDM) approach adopted by the AK Model is generally meritorious because (1) it accounts for localized damage (which actually happens in rock deformation) and microscopic initiation of cracks unlike Fracture Mechanics which considers a clearly defined discrete macroscopic cracks; (2) it is capable of analyzing complex fractures and networks which is not possible in elastic analysis (such as Linear Elastic Fracture Mechanics (LEFM)); (3) there is no need for any special initial assumptions such as initial perturbations; and (4) it has no computational limitations on number of fractures like the Extended Finite Element Method (XFEM) [34].
3. Elastoplastic-Damage Constitutive Model
A coupled elastoplasticity-damage model (AK Model) is adopted here. It was developed by Abu Al-Rub and Kim [4] based on an earlier work by Cicekli et al. [3], hence the name AK. For a full description of the model, please see [3, 4]. The model stands out because it presents a coupled anisotropic damage and plasticity constitutive model that predicts rock’s distinct behavior in compression and tension with the following: (1) a modified continuum damage mechanics framework to include quadratic isotropic and anisotropic variation of the effective (undamaged) stress in terms of the nominal (damaged) stress. The nominal and effective configurations are explained in Section 3.1 and (2) two novel and different damage (power) evolution laws for both tension and compression for a more accurate prediction of rock behavior after damage initiation.
An overview of the model’s main constitutive relations is presented here to give a general idea without detailing its numerical implementation.
3.1. Anisotropic Damage Model
In damage mechanics, the damaged (nominal) configuration of a material is the normal state of the material with imperfections like voids (pores) and cracks. To analyze this, an imaginary state of the material called the effective (undamaged) configuration is assumed (see Figure 1).
The effective (undamaged) configuration of the material considers is to be intact without any voids or cracks. Strain is assumed to be constant in both configurations (i.e., strain equivalence hypothesis) and a damage internal state variable is defined. The damage variable, which is a degradation variable, varies from 0 to 1; a value of zero indicates no damage and one indicates full damage (i.e., fracture). The relationship between the stresses in the damaged and effective configurations is given by where and are the Cauchy stresses in the damaged and effective configurations, respectively. Variables in the effective configuration are denoted by . Thus the damaged elasticity tensor is given in terms of the effective elasticity tensor by
Rock has a distinct behavior in compression and tension. Thus the Cauchy stress tensor is decomposed into two parts using spectral decomposition technique, the positive part (tension) and the negative part (compression). The following relations show the spectral decomposition of the Cauchy stress tensor for both damaged and effective configurations [4, 35]: where represents tensile components and represents compressive components.
Fourth-order projection tensors and damage-effect tensors are introduced as well as the accompanying spectral decomposition into tensile and compressive parts. With this and further simplification, the following relations are obtained [4, 36]:where and are the projection tensors that describe the orientation of the tensile and compressive principal stresses, respectively, where is the identity fourth-order tensor, is the damage-effect tensor, and are the damage variables describing the evolution of cracks due to tensile principal stresses and cracks due to compressive principal stresses, respectively, is the Heaviside step function, are the principal values of , and are the corresponding directions. The model has been validated using several experimental data for concrete; see [4].
3.2. Plasticity Yield Surface
Rock materials show plastic deformation before failure, especially under high confinement pressures. In plastic deformation, there is a need to define three elements: (a) yield criterion to describe the onset of inelastic deformation; (b) flow rule and plastic deformation function for calculating the magnitude and direction of plastic strain rate; and (c) hardening/softening law for evolution of the yield stress. The yield criterion chosen is that developed by Lubliner et al. [37] due to its capability to account for different tensile and compressive behaviors and confinement pressure effects. The nonassociative plasticity flow rule is used to ensure realistic modelling of volumetric expansion of rock under compression. The yield criterion is expressed in the effective configuration (see Figure 1). This is because plasticity evolution is actually driven by the stresses and strains in the intact material. It is expressed as where is the second-invariant of the effective deviatoric stress tensor , where is the Kronecker delta, is the first invariant of the effective Cauchy stress tensor , which incorporates the hydrostatic pressure effect, is the maximum principal effective stress, is the Heaviside step function, and and are dimensionless constants given bywhere is the initial equibiaxial yield strength, is the uniaxial compressive yield strength, and are tensile and compressive isotropic hardening functions, and the equivalent plastic strains are and in tension and compression, respectively. The isotropic hardening functions are given bywhere and are the compressive and tensile initial yield stresses, respectively; , , and are material constants obtained from the effective configuration of the uniaxial stress-strain diagram. The equivalent plastic strains and their rates are expressed aswhere The superimposed dot indicates derivative with respect to time. The maximum and minimum principal values of the plastic strain rate such that are represented, respectively, by and ; indicate a principal value and, as stated earlier, and indicate tensile and compressive variables, respectively; is a weight factor for tension and compression; and represents the Macaulay bracket taken as .
The nonassociative plasticity flow rule is described bywhere is the plastic Lagrange multiplier, is the plastic potential, and is the dilation material constant. The plastic multiplier is obtained using the consistency condition:
3.3. Tensile and Compression Damage Surfaces
The damage growth function adopted in this model incorporates both tensile and compressive damage. It is as follows:where is the tensile or compressive damage isotropic hardening function. When there is no damage, equals the damage threshold . is the damage driving force or the energy release rate expressed as [4, 38]The evolution equation for is as follows:where is the damage multiplier given by
3.4. Tensile and Compressive Damage Evolution Laws
Both exponential and power damage laws could be used for evolution of the damage variables. While the exponential law has less number of material constants, the power law is proven to give more accurate results [4]. Hence, in this study, the power law is adopted for evolution of damage both in tension and in compression. It is expressed for tension and compression, respectively, as where and are material constants. Under uniaxial loading, the tensile damage isotropic hardening function, , and tensile damage threshold, , are, respectively, expressed as where is the tensile yield strength which is almost equal to the ultimate tensile strength for rocks at which tensile damage initiates. The compressive damage isotropic hardening function, , and compressive damage threshold, , are, respectively, where is the uniaxial compressive stress at which damage starts.
4. Model Calibration
Abu Al-Rub and Kim [4] proposed a method of obtaining unique material parameters based on data from cyclic tests. The proposed method was used to obtain the starting values in this work; the parameters had to be adjusted further to obtain a close fit with experimental data. Carrara marble is used as the rock material.
4.1. Carrara Marble as a Representative Rock Material
Various experimental data exist but calibration of the model requires uniaxial cyclic compression and tension data for unique calibration of material constants [4]. Data for cyclic tensile test of rocks is not readily available because it poses a challenge to experimentalists. It is difficult to perform tensile tests on rocks without introducing spurious stresses. Hence indirect tests such as the Brazilian disc test are used [39]. Chen et al. [40] explained the suitability of the Brazilian disc test in determining the tensile strength of both isotropic and anisotropic rocks. However, they argued that elastic isotropic relations cannot and should not be used for analysis of tests on anisotropic rocks. They used analytical methods in addition to experiments to determine the elastic constants and indirect tensile strength of transversely isotropic rocks. To strike a balance between accuracy and available data, it would be strategic to choose data for an approximately isotropic rock for the Brazilian disc test. Even though some studies reveal some level of anisotropy [41], Carrara marble could be reasonably assumed to be isotropic [42].
Wong et al. [43] recently studied the tensile behavior of Carrara marble using the Brazilian disc test. Moreover, Carrara marble has been and is being widely studied by various researchers (e.g., [2, 16, 17, 28, 43–48]). Walton et al. [46] studied the strength, deformability, and dilatancy of carbonate rocks including Carrara marble. Triaxial test data at different confinements were presented. In this work, the tensile material parameters of the AK Model will be calibrated using data from Brazilian disc test by Wong et al. [43]. For the compressive material parameters, triaxial test data presented by Walton et al. [46] will be used. In addition, experimental data for precracked Carrara marble specimens with various flaw (artificially made preexisting crack) geometries exist for validation [2].
4.2. Calibration Based on Data from Uniaxial Compression and Uniaxial Tension Tests
Data for monotonic uniaxial tests are used in this study. The sources of utilized data for Carrara marble have been presented in the previous section. Based on compressive yield strength, the compressive damage threshold was calculated following (19). The tensile stress-strain curve with a tensile strength 6.9 MPa, as in Wong et al. [43], is adopted here. For compression, the stress-strain curve presented by Walton et al. [46] is used. The obtained material properties are unconfined compressive strength = 94.3 MPa; Young’s modulus = 45.3 GPa; Poisson’s ratio = 0.19.
To ensure that the selected values are representative of Carrara marble properties, a brief review of properties reported by other researchers was carried out. According to Evans et al. [49], the compressive yield strength of Carrara marble is approximately 76 MPa with no confinement. The tensile strength of Carrara marble was obtained as 7.5 MPa [45] and with varying Brazilian disc diameter, a range of 6–8 MPa was reported [50]. A compressive strength of 92 MPa was also obtained by other experimenters [51]. According to Siegesmund et al. [52], the properties of Carrara marble are as follows: unconfined compressive strength = 84.6 MPa; Young’s Modulus = 49 GPa, Poisson’s ratio = 0.19; tensile strength = 6.9 MPa. These values confirm that the selected experimental data fall well within range for Carrara marble properties.
4.2.1. Simulation Setup
For both uniaxial tension and compression, a single 1 mm by 1 mm plane stress element is used for calibration. Each element is supported by rollers on the left and bottom edges and a top displacement is imposed either upward (tension) or downward (compression) as shown in Figure 2.
(a) Uniaxial compression
(b) Uniaxial tension
4.2.2. Material Parameters
To obtain the initial material parameters, the calibration procedure outlined by Abu Al-Rub and Kim [4] was adopted using the cyclic compression data in [46]. However, because the confining pressure for this data set was unclear and most likely nonzero, trial and error was further used to make the parameters fit the data for unconfined uniaxial compression. For tension, trial and error was used based on the data by Wong et al. [43]. The material parameters that provide a good fit for the data are detailed in Table 1.
A comparison between the model prediction using these material constants and experimental data is presented in Figure 3(a) for compression and Figure 3(b) for tension. The model gives a good prediction of the compressive behavior of Carrara marble as seen in Figure 3(a) based on data for uniaxial compression. In the case of tension (Figure 3(b)), the data used was extracted from a Brazilian disc test which is indirect tension. Experiments reveal the progressive development of white patches and an initial nonlinear behavior. The current model assumes a linear elastic behavior until yield. Thus specimen behavior on the onset of linear elastic behavior is used to aid calibration by a shift in the strain values. A more detailed approach would have been the simulation of a complete Brazilian test for a more accurate calibration; however, the current fitting yielded satisfactorily results.
(a)
(b)
5. Model Validation
5.1. Simulation Setup
To validate the model, tests on single-flawed Carrara marble were used. Prismatic specimens with a dimension of 152 76 32 mm were tested by Wong and Einstein [16]. The nomenclature associated with the flaw and specimen geometry is presented in Figure 4.
(a)
(b)
A plane stress representation of prismatic Carrara marble specimens with a single flaw is used for validation of the model. As shown in Figure 5, the bottom of the specimen is pinned and displacement is imposed on the top. An example of a finite element mesh is presented in Figure 6 and consists of linear (CPE3) and bilinear (CPE4R) triangular elements (or triangular elements). Due to stress concentration at the flaw tips, the mesh is made finer at those regions. Cases with flaw angles () of 0°, 30°, 45°, 60°, and 75° are considered.
5.2. Predicted Fracture Patterns
Figure 7 shows the various types of cracks observed in all experiments based on the mode of fracture involved [16]. There are tensile cracks (due to tensile forces), shear cracks (due to shear forces), and mixed tensile shear cracks. Also, types 1, 2, and 3 have been identified for both tensile and shear cracks to help in distinguishing various crack geometries.
(a) Type 1 tensile crack (tensile wing crack)
(b) Type 2 tensile crack
(c) Type 3 tensile crack
(d) Mixed tensile shear crack
(e) Type 1 shear crack
(f) Type 2 shear crack
(g) Type 3 shear crack
Newly evolved crack patterns from simulations are compared qualitatively with experiments from Wong [2] as shown in Figure 8. Although there are some cracks that are not captured by the model, overall, there is a very good qualitative match between the simulations and experiments. As shown, CWP means curvilinear white patch and TWC means tensile wing crack [2]. Clearly for , the model is even able to capture the evolution of CWP. Also, TWC, tensile (T), and shear cracks (S) are all captured by the model for different . One can notice that TWC dominates the crack pattern as increases which is in agreement with the model predictions. Therefore, the model considers both white patches and cracks in the form of damage. The force-displacement curves for different single flaw geometries are presented in Figure 9. This figure shows the different softening mechanisms for different flaw geometries.
Double flaws are also considered as shown in Figure 10 for and size of the ligament lengths (L) of 2a and 4a (see Figure 10). The predictions are not as accurate as in the case of the single-flawed geometry. Specifically, Figure 10(a) shows that the simulations do not capture crack coalescence between the two preexisting cracks as compared to experimental observations. This is obvious because the stress distribution in a double-flawed specimen is more complex. However, better agreement is seen in Figure 10(b). The force-displacement curves for the two different double flaw geometries are presented in Figure 9, which shows that as L increases stiffer response is obtained.
(a) Coplanar flaws with an inclination angle of and
(b) Coplanar flaws with an inclination angle of and
It is worthy to mention that experimental results vary from one test to another for the same flaw geometry. This makes it quite difficult to get a definite fracture pattern; most recurrent patterns are taken as representative results. Hence, a qualitative comparison would be sufficient. More accurate results could be obtained if data from the same sample is used for both tension and compression. While the validation experiments use a certain Carrara marble sample, the tensile data was based on a different sample and the data for compression was based on yet another sample. Although all are samples of Carrara marble, there could be slight variations in material behavior.
6. Incorporating the Effects of Material Heterogeneity
Rocks as observed using microscopy have a heterogeneous microstructure; their heterogeneous nature is also visible from outcrops. Stress-induced and preexisting defects in the form of voids, microcracks, and weak interfaces contribute to rock heterogeneity and randomness in microstructure thereby influencing rock behavior and crack propagation. Studies have shown that the compressive strength of carbonates is a function of the grain size, porosity, and elastic properties (see [11, 53]). In marbles and limestones, the uniaxial compressive strength increases linearly with the inverse square root of mean grain size [54, 55]. These studies indicate that heterogeneity, in its various manifestations, is very important in the study and modelling of fracture propagation in rocks [56]. The importance of capturing heterogeneity in hydraulic fracture modelling has also been demonstrated in the work of Sirat et al. [57].
De Marsily et al. [58] presented a history of heterogeneity and how it applies to hydrogeology. They generalized the methods of treating heterogeneity into (1) averaging, that is, defining equivalent homogeneous properties, and (2) describing spatial variability of rock properties from geologic observations or local measurements. These descriptions are implemented through either continuous geostatistical models or discontinuous facies (distinctive rock unit) models.
Liu et al. [56] pointed out various methods used to incorporate heterogeneity in numerical modelling including (1) randomly assigning different properties through a probabilistic distribution or otherwise; (2) use of a mesh with random geometry but equal element properties; (3) projecting a generated microstructure on a regular element network and then assigning element properties depending on the position of each element; (4) combining a random geometry with a generated microstructure; (5) use of the homogenization theory to obtain effective properties through a representative volume element (RVE).
Randomly assigning properties involves statistical modelling which requires a probability distribution function. In this regard, the Weibull distribution which was originally developed for modelling the breaking strength of materials is widely used [56]. On a different note, studies using the random finite element method (RFEM) developed by Fenton and Griffiths (e.g., [59]) for random behavior of soils have also considered the spatial variation of properties using a lognormal distribution [59]. Results from statistical modelling of rock heterogeneity and the homogenization theory as implemented in the interaction code showed a close match with experiments [56]. This builds some level of confidence in adopting the two methods and also verifies the code. or “rock and tool interaction code” is a numerical approach that was developed based on Rock Fracture Process Analysis (RFPA) and the Finite Element Method (FEM).
In this work, the AK Model is extended to study the effects of heterogeneity on rock fracture propagation while adopting the Weibull distribution function for the elastic properties, specifically Young’s modulus.
6.1. Weibull Distribution Function as a Representative Function for Heterogeneity
The effect of heterogeneity is incorporated in the model using the Weibull distribution. This is partly because it was developed to model the strength of materials [60] which is similar to our application and partly because of its flexibility; it was also used successfully for rock fracture propagation in the RFPA code [19] and is generally widely used [19, 56]. The Weibull distribution is given by [19] where represents a mechanical property such as Young’s modulus, is the scale parameter, and the shape of the distribution function is described by the parameter . Tang [12] and Wang et al. [19] considered to be the homogeneity index that increases with an increase in material homogeneity. According to [19] where the Weibull distribution is used to assign element properties, signifies material heterogeneity accounting for cracks and pores in the microscale [19]. Wang et al. [19] considered homogeneity index of 0.6 and 1.1 as values for a heterogeneous rock; values of 1.5 and 2 were considered to be relatively homogeneous. Coarse and medium-grained marble were assigned values of 1.5 and 2, respectively [19].
Studies have shown that Young’s modulus varies spatially in rocks [61, 62]. Hence, it is taken as the property to be varied in this study; that is, . The variation of with the shape parameter is illustrated in Figure 11. As the shape parameter increases, Young’s modulus approaches the base value, = 45300 MPa; that is, homogeneity increases. To implement the Weibull distribution, a built-in function in MATLAB was used.
6.2. Effect of Varying the Homogeneity Index (Shape Factor)
The Weibull distribution is used to spatially vary Young’s modulus. To do this, a random distribution of Young’s modulus is generated using the Weibull distribution and each element in FEM is assigned a value from the distribution. Results are presented in Figure 12 for the 60° flawed specimen setup as used for validation in the previous section. For low values (≤7), Young’s modulus varies significantly thereby giving interesting failure patterns. With an increase in (≥15); however, favorable comparison with experimental data is obtained. Only high values of (very homogeneous) show acceptable results in simulations. This confirms earlier studies asserting that Carrara marble could be reasonably assumed to be a homogeneous material [63, 64]. Only high values of show acceptable results in simulations.
(a) m = 2
(b) m = 4
(c) m = 7
(d) m = 15
(e) No heterogeneity effects
(f) Experiment
The force-displacement diagrams (Figure 13) show an increase in the average elastic modulus (more values are close to the base value) with an increase in homogeneity index, . Also, Figure 13 shows that materials with higher heterogeneity are weaker due to increase potential of cracking in different regions.
7. Conclusion
A coupled elastoplasticity-damage model developed for quasi-brittle materials has been calibrated and validated for rocks using Carrara marble data. The model’s prediction of fracture patterns matches closely with experimental results thereby verifying the model as a capable tool for fracture prediction in rocks. Thus, the model, AK Model, could be used to aid better reservoir management, whether hydrocarbon, geothermal, or underground water reservoirs. Although the model is capable of describing the conferment effect, this was not taken into account in this study and will be the focus of future work. It is important to investigate such effects as confinement significantly affects rock behavior. The Weibull distribution function has also been used to incorporate heterogeneity effects in the coupled elastoplasticity-damage model by varying Young’s modulus. Results hint that Carrara marble could be reasonably assumed to be a homogeneous material.
Competing Interests
The authors declare that they have no competing interests.