Dolomite fracture modeling using the Johnson-Holmquist concrete material model: Parameter determination and validation

Micha? Kucewicz, Pawe? Baranowski, Jerzy Ma?achowski

Faculty of Mechanical Engineering,Institute of Mechanics and Computational Engineering, Military University of Technology, 2 Gen.S.Kaliskiego Street, Warsaw,00-908, Poland

Keywords:

ABSTRACT In this paper, the Johnson-Holmquist concrete (JHC) constitutive model is adopted for modeling and simulating the fracture of dolomite.A detailed step-by-step procedure for determining all required parameters, based on a series of experiments under quasi-static and dynamic regimes, is proposed.Strain rate coefficients, failure surfaces, equations of state and damage/failure constants are acquired based on the experimental data and finite element analyses.The JHC model with the obtained parameters for dolomite is subsequently validated using quasi-static uniaxial and triaxial compression tests as well as dynamic split Hopkinson pressure bar(SHPB)tests.The influence of mesh size is also analyzed.It shows that the simulated fracture behavior and waveform data are in good agreement with the experimental data for all tests under both quasi-static and dynamic loading conditions.Future studies will implement the validated JHC model in small- and large-scale blasting simulations.

1.Introduction

Knowledge of the mechanics of rock fracture is widely used in mining and civil engineering to improve the effectiveness and safety of underground drilling or rock excavation (Phang et al.,1983; Kulatilake et al., 2013; Kuili and Sastry, 2018; Baranowski et al., 2019).Depending on the processes that specify the loading conditions,various fracture mechanisms of rock are proposed from the micro- and macro-scopic points of view.Most rocks undergo brittle fracture, which is sensitive to the stress state under which deformation occurs(Bieniawski,1967; Sun et al.,2014; Peng et al.,2018).Microcrack propagation decreases the mechanical resistance of the material to external loading, and irreversible deformation occurs as a result of dislocation motion (Jaeger et al., 2007; Zhang and Yu, 2017).In uniaxial compression, the major cracking mechanism is longitudinal splitting or initiation of multiple cracking induced by boundary effects such as friction with machine crossheads (Bahat et al., 2001; Basu et al., 2013).By contrast, in triaxial tests, cracking is induced by shearing and is inclined at a small fracture angle(20°-30°)relative to the axial direction(Bahat et al.,2001; Chakraborty et al., 2019).Increasing the radial pressure induces well-defined shear failure, and further increases in pressure will cause failure along several planes.When a high confining pressure acts laterally on a sample, a brittle-to-ductile transition can be observed(Aadn?y and Looyeh, 2019).

The fracture mechanism is also often related to the microstructure of the sample and the presence of imperfections.The strength and elastic properties of rock are mainly affected by inclusions of different materials characterized by lower strength,trans- or inter-granular cracks, pores, voids, etc (Rajabzadeh et al.,2012; Geng et al., 2018).It is shown that tensile strength is the only important factor in the fracture gradient of rock masses(Aadn?y and Looyeh,2019).While the tensile strength is significant for samples with small dimensions and affects the anisotropy of mechanical properties, as the tested volume increases, the influence of inclusions decreases,and a large rock massif can be treated as almost quasi-isotropic (Borovikov and Vanyagin, 1995; Soga et al., 2014).Furthermore, different mechanisms of rock fracture are observed under static, dynamic and shock loading conditions(Grady et al.,1976; Mishra et al., 2017; Wang et al., 2019a, b).

The growing research interest in rock mechanics has focused on appropriate numerical modeling of rock deformation with special consideration of fracture mechanisms.Among modeling methods,the finite element method(FEM)remains the most commonly used due to its universality and simplicity of implementation(Podgórski,2017; Chró′scielewski et al., 2019; Ke?dzierski et al., 2019; Mayer et al., 2019; Migueis et al., 2019; Kurzawa et al., 2020).FEM describes a solid body in a discrete way without any imperfections while incorporating constitutive relationships.Damage and fracture are characterized by a damage index calculated based on damage constants and mathematical functions.Several constitutive models describing the behavior of brittle materials with different levels of complexity are available in commercial finite element(FE)codes.The most popular ones are the Mohr-Coulomb (M－C)(Hackston and Rutter,2016;Chang and Konietzky,2018;Geng et al.,2018), Johnson-Holmquist ceramics (JH-2) (Ai and Ahrens, 2006;Jaime, 2011; Banadaki and Mohanty, 2012; Wang et al., 2018;Baranowski et al., 2020), Karagozian and Case concrete (KCC)(Mardalizad et al., 2017, 2019; Huang et al., 2020), Cap (Sandier et al.,1974; Schwer and Murray, 2002; Jiang and Zhao, 2015), and Johnson-Holmquist concrete (JHC) (Holmquist et al.,1993; Meyer,2011; Ren et al., 2017) models.

The present paper is the part of a wider project aiming at simulation and optimization of parallel cut-hole blasting and fragmentation.Recently, the authors have tested two different constitutive models capable of representing the rock subjected to high-strain rate loads (Baranowski et al., 2020; Kucewicz et al.,2020).However, other material models were also considered.The present study adopts the JHC model, which is described briefly in Section 3, for modeling dolomite.The main benefits of this constitutive model are the three-range equation of state(EOS)and the inclusion of both volumetric and equivalent plastic strains in the damage assessment.However, many authors have noted a significant deficiency of the JHC model in the reproduction of tensile damage (Holmquist et al.,1993; Polanco-Loria et al., 2008;Lu et al.,2012;Islam et al.,2013;Kong et al.,2016;Li and Shi,2016;Ren et al., 2017).Fortunately, available hydrocodes include two solutions to deal with this limitation: modification of the original failure surface with a user-defined interface (Polanco-Loria et al.,2008; Kong et al., 2016; Li and Shi, 2016) and application of FEs allowing for erosion (finite element deletion) (Wang et al., 2007;Kala and Huˇsek, 2016).The first requires high-level programming skills and access to source codes, while the second unphysically reduces the mass of the material based on previously calibrated erosion parameters,which are dependent on the sizes of the mesh,the modeled problem and the damage parameters.

Although rock characterization and its constitutive modeling have been extensively studied, few works have covered the experimental testing, constitutive description and numerical simulation of dolomite.In the present paper, a detailed characterization of the mechanical properties of dolomite, including microand macro-scopic studies, is presented, and two leading mechanisms of cracking are analyzed.Next, a calibration procedure for JHC constitutive model parameters based on dolomite test results is described.Static, dynamic and shock data are applied to faithfully reproduce the strength of dolomite under various loading conditions.JHC constants for static and dynamic simulations are determined separately, and M－C theory is adopted to calculate the cohesion.The damage parameters representing fracture and cracking during deformation are calibrated numerically, and a broad sensitivity study is performed to minimize the influence of mesh size on the simulation results.The evaluated and correlated JHC parameters are validated in triaxial and uniaxial compression tests under static and dynamic strain rates with the commercial FEM hydrocode LS-DYNA (Hallquist, 2006), and satisfactory quantitative and qualitative agreements are obtained.Finally, an erosion criterion is applied in the FE model to reproduce the loss of continuity in the sample.In future investigations, the JHC model will be implemented in small-scale blasting simulations and compared with the KCC and JH-2 constitutive models.

2.Characterization of dolomite

2.1.Microstructure characterization

Dolomite, also known as calcium magnesium carbonate, is a sedimentary mineral with the chemical composition CaMg(CO3)2.The term dolomite also refers to sedimentary carbonate rock composed mostly of dolomite mineral.Dolomite is not widely distributed worldwide due to the relatively strict requirements for its formation.The most common process is dolomitization, in which the calcium ions in calcite are replaced by magnesium ions.This process is strongly dependent on the Ca/Mg ratio(Machel and Mountjoy,1986; Mariˇci′c et al., 2018).The conditions under which dolomite is formed impact its microstructure, composition and hence mechanical properties (Mariˇci′c et al., 2018).As dolomite is formed via long-term processes,its structure is heterogeneous and includes many imperfections, such as voids, pores, veins, cracks,and inclusions.Some examples of tested samples are presented in Fig.1.The material used in the investigations was excavated in Lover Silesia,Poland.In general,dolomite is a very fragile rock,and extreme caution was taken in the preparation process to obtain cylinders that were intact as possible for the tests.

The chemical composition of the dolomite was analyzed by scanning electron microscopy(SEM),and the percentage content of each element is presented in Fig.2.In the fractured sample, the failure mechanism of the components was highly visible.Two main components with different geometrical dimensions of grains were recognized.The largest grains were formed from anhydrite,which contributed the highest fraction of the microstructure(spectrum 4).The average measured size of these grains was greater than 500 μm.The failure of anhydrite resulted from the formation of a sliding plane(gliding)through the grains,which corresponds to the mode II/III of cracking.The second component was dolomite mineral with an average grain size of 40 μm.Its failure mechanism was deduced as dimpling from the SEM images.In dimpling, the connections between grains are strong enough to resist the load,but the grains themselves split.In addition to these two main components, large inclusions of calcite were observed in the samples.In general,calcite is more brittle and less durable than the other mineral components of dolomite.In the tests, the fracture surfaces, corresponding to the weakest sections of the samples, usually crisscrossed these inclusions.Computed tomography (CT) of the samples showed that for the tested cores, the calcite volume reached 10%, and the inclusion size varied from 1 mm to greater than 10 mm.Fig.3 presents the microscopic view of fractured surfaces.

2.2.Assessment and identification of mechanical properties

Numerical simulations were preceded by material tests to determine material strength properties.These results were then applied directly or indirectly in constitutive model calibration,which is one of the most important steps in the numerical model preparation process.Therefore,in the present paper,the following tests were conducted to characterize the material: static and dynamic uniaxial compression tests, triaxial compression tests at confining pressures of 10 MPa,17.5 MPa and 25 MPa, static (SBT)and dynamic Brazilian (DBT) tests for indirect tensile strength determination, and cyclic uniaxial compression tests.

Fig.1.Macro- and micro-scopic views of the tested dolomite samples.CT- Computed tomography; SEM - scanning electron microscopy.

Fig.2.Chemical composition of the investigated dolomite samples.

All the tests were performed according to the standards suggested by the International Society for Rock Mechanics and Rock Engineering(ISRM)(Hudson and Ulusay,2007).The samples for the uniaxial,triaxial and cyclic compression tests were cylinders with a 2:1 height/diameter ratio (100 mm/50 mm).The discs for the Brazilian test had a 2:1 diameter/depth ratio (50 mm/25 mm).Static tests were performed at a strain rate of 0.0024 s－1, using an INSTRON 8802 universal testing machine.The range of confining pressures in the triaxial test was limited by the pressure chamber specification, which was limited to 30 MPa.For triaxial compression, loading was applied in two steps.Hydrostatic pressure was first applied by pressurized oil acting on all faces of the sample until reaching a value of 10 MPa,17.5 MPa or 25 MPa.In the second step,the load was applied by the machine crosshead with forcecontrolled movement.The test ended when the immediate force drop after crossing the maximum strength stopped the axial displacement of the crosshead; thus, the residual strength could not be measured.For calibration purposes,all experimental data of each representative sample were averaged (samples with deviations over 25% were ignored), as listed in Table 1.

Dynamic tests were performed with a split Hopkinson pressure bar (SHPB) system (Jankowiak et al., 2020).Samples with a 1:1 height/diameter ratio were used to fulfill the force equilibrium condition during the entire deformation process.For assessingtensile strength, the same disc dimensions used in the static tests were adopted.To smooth the signal in the loading phase,a copper pulse shaper with a radius of 8 mm was used.The incident and transmission bar strains were measured using a full-bridge straingage circuit.The deformation was captured using a high-speed camera at 95,000 frames per second (fps).

Fig.3.Microscopic view of dolomite microstructure at the fracture surface.

Table 1Experimental results of dolomite.

3.Description of the JHC constitutive model

The JHC model was initially developed to reproduce the behavior of concrete subjected to large strains, high pressures and high strain rates(Holmquist et al.,1993).Concrete and rock are both quasi-brittle materials with similar fracture mechanisms,but there are two major differences.First, rock is created via natural processes, and samples are obtained from the in situ rock fragment/mass using mechanical processes such as cutting.By contrast,concrete is a material created under laboratory conditions from several components (water, aggregate and cement).Second, rock has bedding planes that affect the anisotropy of its mechanical properties, and these planes are not included in many common constitutive models for brittle materials (Shah et al.,1995).

In general, the JHC model distinguishes unsymmetrical strength properties of material under tensile and compressive loadings.The model relates the internal pressure and differential stress in the deformed material state to determine the actual strength.The normalized equivalent stress(σ*)resulting from differential stresses is referred to as the quasi-static uniaxial compressive strength(UCS)and is calculated as follows (Holmquist et al.,1993):

where σ is the actual equivalent stress.The full specific expression for the yield surface is given by the following equation, which is called the constitutive model (Holmquist et al.,1993):

where A is the normalized cohesion, i.e.a theoretical stress that causes material fracture due to shearing; B is the normalized pressure hardening coefficient; N is the pressure hardening exponent; C is the strain rate coefficient;P*is the normalized pressureand is given by in which P ＝ (σ1＋σ2＋σ3)/3 is the hy-

drostatic part of the stress tensor;is the normalized strain rate of the material and is given byin whichis the reference strain rate (mostly assumed as 1 s－1); and D is the accumulated damage coefficient and takes a value ranging from 0 for unconfined material to 1 for fully damaged material.In the above definitions,normalization is performed by dividing the value of a parameter by.

A graphical representation of Eq.(2)is presented in Fig.4.In the JHC model, the damage of the material is understood as a progressive decrease in its stiffness until the critical fracture surface specifying the residual strength is reached.Fully damaged material can only transfer compressive loads.The model reproduces the state at which the material has the form of an aggregate with no tensile strength.The damage parameter D is an extension of the Johnson-Cook (JC) model, and the plastic strain is supplemented with the volumetric strain increment (to consider hydrostatic compressibility) according to the following equation (Holmquist et al.,1993):

where Δεpand Δμpare the equivalent plastic strain and volumetric plastic strain increments in each integration cycle, respectively;andis the total plastic strain to fracture under constant pressure P and is expressed as follows (Holmquist et al.,1993):

where D1and D2are the damage constants, T*is the maximum tensile hydrostatic strength, and EFMIN is the minimum plastic strain before fracture.

The JHC model treats deviatoric and spherical quantities of stress separately; therefore, an EOS can be implemented.The EOS describes the relationship between pressure and volumetric strain and is described by three characteristic ranges, as shown in Fig.5.The behavior of rocks under triaxial loading conditions is similar to that of concrete, and three corresponding ranges during deformation can be highlighted.These observations confirm that the JHC model is reliable for rock simulation.The first phase(OA)represents a fully reversible linear elastic response of the material.It covers thepressure from the negative cutoff for tensile loads－T*(1－D)to the elastic limit Pcrush(pressure at UCS) under compression and is expressed by (Holmquist et al.,1993):

Fig.4.Visualization of failure surfaces and the damage evolution function in the JHC model.

Fig.5.Relationship between pressure and volumetric strain implemented in EOS.

where K is the elastic bulk modulus;and μ ＝ ρ/ρ0－1,in which ρ and ρ0are the current and initial densities of the material,respectively.A simplified expression for K is given by K＝Pcrush/μcrush, where μcrushis the volumetric strain at the elastic limit.

The second range (AB) is a transitional region.The decrease in stiffness is caused by a gradual collapse of air void pores.Under actual deformation,the movement of the grain boundaries related to microcracking phenomena is locked even though the material strength increases (until locking pressure Plockis reached).This range is expressed by(Holmquist et al.,1993):

where Klockis the bulk modulus of the material during the described irreversible deformation.In fact, the range AB is a linear interpolation between the first range(OA)and the third range(BC)and cannot be controlled directly by the user.Testing has confirmed that the slope of this interpolating line (which corresponds to the bulk modulus)should be less than the slope of OA(Holmquist et al.,1993).

The third range(BC)describes the strength of the material with no pores in the structure.It is characterized by a nonlinear relationship and may be described by the following 3rd-order expression (Holmquist et al.,1993):

4.JHC parameter determination

To calibrate the JHC constitutive model,21 constants were determined.Most were obtained based on the tests,while some are strictly numerical parameters.The main aim of the present paper is to provide a detailed description of the calibration methodology for the JHC model to be used in various computational problems.This methodology may be adopted for simulations of middle/high-strength rocks and rock-like materials regardless of the loading conditions.The procedure was implemented with dolomite as an example.

4.1.Strain rate sensitivity

The tests under static and dynamic loading conditions confirmed that the dolomite was sensitive to strain-rate effects.In the JHC model,strain rate hardening is controlled by the factor C, which is a multiplicator of the strain rate logarithmic function in Eq.(2)and scales both intact and damaged material strengths.Rocks generally exhibit a bilinear dynamic increase factor (DIF), which should be separated for lower and higher strain rates.Moreover, the strain rate sensitivity of granular materials such as rock or concrete is non-symmetrical, and significant differences are observed between hardening under tension and under compression.These differences arise from the change from the intergranular to transgranular cracking mechanism,which leads to faster crack propagation and crack shapes that reflect the decreased area of the planes(Mahanta et al.,2017).The applied linear approach is not adequate to reproduce the kinematic hardening under a wide range of loading rates.As a result,the coefficient C must be calibrated separately for the two ranges of strain rates.Literature data for dolomite strength at various strain rates are limited;thus,the calibration in the present study was based mostly on data from the authors’laboratory.In addition,the dynamic response of limestone,which is the mother rock of dolomite, was assumed to be similar to that of dolomite, and the corresponding strain rate sensitivity regions were assumed (Liu et al.,2018).The limiting strain rate above which the strengthening of dolomite should be considered was assumed to be 1 s－1.For these loading conditions,the first“quasi-static range”is valid.The second “dynamicrange”is reliable up to 125 s－1due to the lack of data for higher strain rates.The authors’ initial simulations confrimed that the JHC model does not scale th e tensile strength with respect to the strain rate, but the value of T*does impact the damage calculation (see Eq.(4)).Therefore, a proper value of this parameter should be adopted depending on the predicted strain rates at which the model will be applied.Again,two different values were used for static and dynamic loadings:＝42.4 MPa.The measured static tensile strength was equal to 4.5 MPa; however, this indirect method underestimates the results,and a correction factor of 1.15 from the literature was used(Pittet and Lema?tre,2000).For dynamic results,underestimation has not been described in the literature,and thus the dynamic tensile strength remained unchanged.

The parameter C for compressive loadings is determined from SHPB tests, at which the dominant observed failure mechanism in brittle material is multiple axial cracking (Kong et al., 2018).The coefficient C is introduced to take into consideration a viscous behavior of material, while other important effects are inertia effects at micro- and macro-scale.All of described effects are timedependent which causes an experiment-dependent behavior on the loading history.However, to simplify the model and due to a high conciseness of dolomite, the maximum strengths calculated from well-known basic SHPB equations at constant strain rate plateau were used for calibration.

The strain rate sensitivity of a material can be directly measured from uniaxial tests performed at different strain rates as the increase in its maximum strength.In the first step, all test results from dynamic compression tests were normalized using fc′＝213 MPa(Fig.6a).These points were connected by straight lines to the normalized tensile strength T* in the P*-σ* system.The average maximum tensile hydrostatic pressure from the tests,which was 1/3 of the measured tensile strength, was Tdynamic ＝ 14.4 MPa and Tstatic＝ 1.73 MPa.In the second step, to unify the strength increase,points at a pressure of P ＝ 71 MPawere determined for each line.This value of P is the pressure atfrom the quasi-static compression test (Fig.6a).The obtained points and the corresponding strain rates are shown in Fig.6b and were interpolated with a logarithmic linear function with slope C.Separate hardening coefficients were determined for lower and higher strain rates than 1 s－1.Eventually, Cstatic＝ 0.00553 and Cdynamic＝0.0307 were adopted for static and dynamic simulations,respectively.The main disadvantage of this method is the difficulty of precisely scaling the strength for a wide range of strain rates,and thus low and high strain rates should be investigated separately.Additionally,the value of C is strongly dependent on the number of experimental points used for calibration.

4.2.Calibration of the failure surface

In the next step, the failure surface was determined based on uniaxial and triaxial static compression tests.Experimental data(Handin et al.,1967;Mogi,1971;Cie′slik,2014)were also adopted to increase the pressure range (up to 600 MPa) at which the model could be reliably applied.The strengths in uniaxial and triaxial compression tests of the materials presented in these previous studies were similar to those of the dolomite described in the present paper.In many researches adopting the JHC model, the material’s normalized cohesion A is calculated by simple fitting of the failure surface equation to experimental points rather than using the actual data and theoretical considerations.Usually, the value of normalized A ranges from 0.3 to 0.8(Holmquist et al.,1993;Islam et al.,2013).In the present paper,the M－C linear model was applied to determine A as an intersection point of circular envelope and shear strength axis.The cohesion indirectly affects the residual strength of the material, because it takes a value of 0 for fully damaged material while the other parameters remain unchanged.Thus, a higher value of cohesion increases the difference between the maximum and residual strengths.

The shear strength of the material in the M－C model is given by(Hackston and Rutter, 2016):

where σnis the stress in the normal direction,φ is the slope of the envelope, and a is the determined cohesion of the material.

When fitting the M－C envelope to quasi-static uniaxial and triaxial compression test data, our test results for confining pressures up to 25 MPa were more reliable than the literature data and had a higher priority in determination of the dolomite’s cohesion.The values of σ1and σ3were determined as the maximum registered strength and radial confining pressure from each uniaxial and triaxial compression test, respectively, for all applied results.The diameter of each circle is a difference between those values.Points representing the critical shear stress were approximated with linear function, and the cross point of this function and ordinate corresponding to shear stress after normalization was assumed as normalized cohesion A (Fig.7a).Finally, we have a ＝63 MPa and φ ＝ 38°; thus, the normalized value of A ＝ 0.296 was used for the static strain rate range.

Different approaches of cohesion determination were used for dynamic range of strain rates.As there is a lack of triaxial compression test results for dynamic conditions,the fitting of M－C envelope to static outcomes may result in a significant error in dolomite stiffness.Thus, the value of A was determined using theleast square method to best fit the constitutive equation to static and dynamic uniaxial test results considering determined values of coefficient C.The value of A ＝0.12 was used for dynamic strain rates.

Fig.6.Procedure for determining the strain rate hardening coefficient for dynamic tests.

Fig.7.JHC failure surfaces calibrated for static and dynamic loading conditions: (a)Cohesion from Mohr-Coulomb theory and(b)JHC constitutive equation for strain rates of 0.0024 s－1 and 94 s－1.UC-uniaxial compression;TXC-triaxial compression;SFMAX- maximum yield strength.

Next,to determine the failure surface, the experimental results were converted into data in P*-σ* space according to Eqs.(3) and(11) (Islam et al.,2013):

where J2is the second invariant of the deviatoric stress tensor.The calculated points were normalized by the UCS(Fig.7).The discrepancy in the results confirms that for organic materials such as rocks, the concurrence of test results is limited, and homogenization of the model must be performed.The failure surface applied in the JHC model is strain rate dependent, thus separate sets of parameters were determined for strain rates inthe range of ε˙static＜1 s－1＜ ε˙dynamic.The strain rate of 1 s－1was set as the reference strain rate ε˙0.Experimental points were approximated by the constitutive equation(Eq.(2))with iterative changes in B and N to minimize the error.To determine the failure surface shape for static loadings,triaxial tests were crucial because all the tests were performed at a strain rate of 0.0024 s－1,and Cstaticwas used.On the other hand,for dynamic loading tests,it was more important to fit the failure surface so that it crossed the points from the dynamic uniaxial compression tests,as shown in Fig.7 for strain rate ˙ε ＝94 s－1.This point is directly crossed by the failure surface scaled by the previously determined coefficient Cdynamic.Triaxial tests have limited validity in dynamic failure surface calibration due to the static conditions under which they are performed.Data for dynamic confined compression tests on dolomite are not available.In Fig.7, additional points are presented and they represent the residual strength from tests at which the post-failure behavior is predicted,as described by Kucewicz et al.(2020).The static failure surface was determined in such a way that these points were crossed when D ＝1 in order to ensure that the model reproduced the residual response of the actual material.Ultimately, B ＝1.751 and N ＝0.865 were obtained for lower strain rates, whereas dynamic B ＝ 2.286 and N ＝ 0.95 were fit to the SHPB data,similar to the normalized cohesion.The value of SFMAX ＝5 resulted from the limited increase in rock strength when the deviatoric stress exceeded σ*＝ 5, as shown in Fig.7b.

The maximum tensile hydrostatic pressure T has an influence on damage calculation (Eq.(4)).It was determined separately from static and dynamic Brazilian tests for ˙εstaticand ˙εdynamicranges of strain rates.This experimental method for measurement of tensile strength is burdened with an error related to uneven distribution of tensile strength in the sample.To reduce this error, the averaged tensile strength was increased by 15%, as suggested by Pittet and Lema?tre (2000).The hydrostatic pressure was 1/3 of calculated tensile strength.Finally, Tstatic＝ 1.73 MPa and Tdynamic＝ 14.4 MPa were used for numerical simulations.

4.3.Equation of state determination

The EOS applied in the JHC model controlling the hydrostatic behavior of the material is divided into three parts and should be determined from hydrostatic or uniaxial strain compression tests(Holmquist et al., 1993).Due to the lack of such data, a different approach that is often used for the JH-2 ceramic model was adopted to determine the EOS parameters.This method is based on uniaxial shock data from flyer impact tests and is described below (Wang et al., 2018; Baranowski et al., 2020).

The first range, which describes a reversible static range, was obtained from the results of uniaxial compression tests.The dolomite bulk modulus K was calculated using Eq.(6).In JHC model,the increase of material stiffness due to lateral confinement is not supported by default.The elastic modulus E,which in this equation is assumed as the slope of the stress-strain curve,was taken as the average from the uniaxial and triaxial tests.This approach compensated for the effect of increasing dolomite stiffness when the additional radial confining pressure was present.In the triaxial tests,the bulk modulus increased significantly until reaching values of σ2＝ σ3＝ 60 MPa.The volumetric strain was not directly measured during uniaxial compression tests, but calculated from Eq.(5) using an abovementioned averaged E.A more detailed description of this procedure was reported by Kucewicz et al.(2020).Eventually, Pcrush＝71.3 MPa, equivalent to 1/3 of the UCS, and μcrush＝ 0.00602 were calculated from Eq.(6).The resulting average bulk modulus was K ＝11,709 MPa.This range of EOS is crucial from the dolomite stiffness point of view.The increase of P over the value corresponding towith maintaining the constant bulk modulus has a limited effect for response of JHC model, until the failure surface is correctly determined (directly crosses the point corresponding to).

In the next step,the third range of the EOS was determined.Its implementation is similar to that for the JH-2 constitutive model,thus a similar approach for its determination using shock Hugoniot data was adopted (Shang et al., 2000; Islam et al., 2013).Experimental points giving pressure as a function of volume changedetermined from flyer impact test results(Heard et al.,1973;Grady et al.,1976;Larson,1980)were calculated and approximated by Eq.(2) with acceptable accuracy.For the investigated dolomite, this range was valid for pressures up to 8 GPa.The modified volumetric strain was used as a function variable for fitting the pressure equation to the experimental data.The value of K1corresponds to the unloading modulus of dolomite.The following constants for the polynomial equation were obtained:K1＝ 31 GPa,K2＝ 134.996 GPa and K3＝28,600 GPa(Fig.8).

Finally, the second range was determined by interpolation between the first and third ranges.The second range is limited by the locking volumetric strain, which was determined as an average value for different dolomites in hydrostatic compression tests performed previously (Handin et al.,1967; Heard et al.,1973).The second range is limited by the value of pressure at which all pores in dolomite are closed and the densification of the material is impeded.The tests showed that for dolomite, which is characterized by low porosity, this region can be neglected.To provide continuity of all three regions of the EOS,the slope of this range was decreased by 5%of the initial dolomite stiffness K,and was equal to K2ndrange＝11.123 GPa.The locking pressure Plockand volumetric strain μlockwere determined by calculation of the cross point of second and third EOS ranges.The linear function crossed the shock stage at point given by Plock＝ 116.1 MPa and μlock＝ 0.00626.It is crucial to remember that μlockis a point at P ＝ 0 calculated with linear projection of μplockto μ axis by the linear function characterized by slope K1.The full determined EOS is presented in Fig.8.

Fig.8.Graphical presentation of the EOS for dolomite.

4.4.Damage and failure calibration

The damage function implemented in the JHC model is controlled by three constants, i.e.D1, D2and EFMIN, which should be determined from cyclic uniaxial compression tests.However,most researches adopt the values of D1＝ 0.04 and D2＝ 1 initially proposed by Johnson and Holmquist(Holmquist et al.,1993).These values determined for concrete are constant regardless of the application of the model and other simulation configurations such as mesh size.The tests showed that dolomite is a very brittle rock.It instantly released the internal elastic energy accumulated during deformation after reaching the maximum strength, and then the sample split into several pieces.In this study, a coupled experimental-numerical approach for obtaining the fracture and post-failure strengths similar to that used for the uniaxial and triaxial compression tests is proposed.In many commonly used constitutive models,the fracture of brittle materials is represented by the damage index 0 ＜ D ＜ 1, in which D ＝ 1 is used for fully fractured material.In the JHC material model, the damage is described by Eq.(3).In addition, the stress-strain curves were compared with the actual test results.For materials in which fracture plays a significant role, the influence of mesh size on the results should not be neglected.Thus, in the present study, the values of the damage parameters were normalized depending on the mesh size.

The adopted method of dolomite modeling and representation using the JHC model and FEM employs a homogenization of the actual material and does not include any imperfections and discontinuities.Thus, the fracture pattern results directly from the loading conditions and strength of the calibrated constitutive model.The results of simulations may reproduce the stiffness of dolomite and its generalized mode of failure; however, calibration of D1and D2could be insufficient to qualitatively reproduce the actual rock fracture patterns, which are mainly predefined by existing pre-cracks.

Initially, numerical simulations were performed on a single cubic element for three different loading conditions:uniaxial straincompression,uniaxial compression and triaxial compression with a confining pressure of 25 MPa.Boundary conditions for each case are presented in Fig.9.The tests were performed under static loading conditions as described in the experimental setup.The trial-and-error approach was used, and the parameters D1and D2were iteratively changed until the best reproduction of the actual post-peak behavior was observed.

Fig.9.Boundary conditions for single-element tests: (a) Uniaxial strain test (hydrostatic compression), (b) Uniaxial compression, and (c) Triaxial compression.DOF means the degree of freedom at mentioned direction.

First, the results of uniaxial strain tests in which hydrostatic loading conditions were reproduced are presented in Fig.10a as a volumetric strain-pressure relationship.The stress state in this case resulted from the volume decrease only, and the strength of the material was calculated directly from the EOS.No shear-induced damage was observed.The resultant curve perfectly agreed with the input data and confirmed that the EOS was correctly determined.In this case,the values of D1and D2did not have any impact on the test results.

Next,the results from uniaxial compression tests were analyzed as a stress-strain relationship.Implementation of the JHC failure surface did not provide residual strength under uniaxial loading conditions, which was in agreement with the actual dolomite response.However, the rate at which material softening occurred depended on D1and D2.A significant impact of EFMIN, which corresponds to the minimum amount of plastic strain before damage occurs, was observed.The default value of 0.01 proposed in Holmquist et al.(1993) was determined by cyclic progressive compression tests on concrete.Different from dolomite, which instantly fractures and splits into pieces, concrete maintains residual continuity after exceeding the maximum strength in uniaxial compression tests.Thus, for dolomite, EFMIN should be smaller than the default value.A study of EFMIN for assumed constants D1＝ 0.06 and D2＝ 1.1 is shown in Fig.10b.Values of EFMIN larger than 0.005 resulted in incorrect damage accumulation.The value of D remained constant before D ＝ 1 was reached, despite the increasing effective plastic strain in the elements.As a result, the model behaved as if the residual strength was present; while for uniaxial loading conditions,this value should be zero.No impact of EFMIN on maximum strength was observed in either uniaxial or triaxial compression.Ultimately, a value of EFMIN ＝0.005 was applied for further tests.

The parameters D1and D2for selected EFMIN ＝0.005 under uniaxial loading conditions are presented in Fig.11.Table 2 lists the values of the parameters for each test.As expected, increasing D1decreased the energy dissipation rate (damage accumulated slowly).The same effect was observed for decreasing D2, which is an exponent in Eq.(4).In Fig.11a-c,evolution of the stress-strain curves depending on the damage constants is presented.In each figure, the legend is replaced with arrows representing the increases in D1and D2to clearly show the results.As the intensity of material softening decreased, the maximum strength increased slightly (up to 3%).However, this change was not significant and was neglected; it was caused by the accumulation of plastic volumetric strain after reaching Pcrush.The same parameter sets were tested under triaxial compression for EFMIN ＝0.005, and the results are shown in Fig.12a-c.The outcomes were similar to those of the uniaxial tests.The only difference was the change in the maximum strength of the material, which decreased up to 10% as D1and D2decreased.The final parameters were set as follows:D1＝ 0.045,D2＝ 1.08 and EFMIN ＝ 0.005.The two complete sets of parameters used for static and dynamic simulations are listed in Table 3 together with all of the required data or tests.

Fig.10.(a) Volumetric response of the material in the uniaxial strain test and (b) EFMIN study under uniaxial loading conditions from a single-element test.

Fig.11.Damage parameters for uniaxial single-element compression tests.

Table 2Values of the parameters for uniaxial and triaxial single-element tests.(a),(b)and(c)refer to corresponding curves in Figs.11 and 12.

5.Validation of the JHC model for dolomite

To confirm the reliability and validity of the proposed calibration procedure, the capability of the JHC model for simulation of dolomite is demonstrated.Quasi-static uniaxial and triaxial compression tests were considered.Moreover, the dynamic uniaxial compression test with use of the SHPB was investigated for different strain rates corresponding to experimental values (see Table 1).The JHC parameters obtained for the dolomite are summarized in Table 3.

5.1.Quasi-static tests

The evaluated damage parameters were applied to the full-scale model in more detail to reproduce the uniaxial and triaxial tests.The aim of this test was to validate the stiffness and fracture pattern from the calibrated JHC model.A mesh size sensitivity study was also performed.Cylindrical samples with a radius of 50 mm and a height of 100 mm were located between two rigid walls.The first wall was constrained, while the second was moved with the prescribed velocity of a rigid wall given by (Hanssen et al., 2002;Kucewicz et al., 2018):

where Tloadis the total duration of loading and dmaxis the final displacement of the rigid wall.

This approach permits the use of an explicit integration scheme, which is more efficient for simulating phenomena in which material damage occurs compared with the implicit method (Kucewicz et al., 2018).Simultaneously, gradually increasing the speed reduces the excitation of stress wave effects, prevents an increase in kinetic energy, and maintains the simplicity of modeling fracture and post-peak material softening.The kinetic energy is less than 0.5% of total energy from simulation.To ensure that low-magnitude stress waves did notaffect early fracture, an additional part-stiffness damping method was adopted.Different from uniaxial compression tests, in the triaxial tests, an additional confining pressure was applied to the side walls of the cylindrical sample, and a twostep analysis was performed.In the first step, dynamic relaxation of the model was implemented to reproduce the hydrostatic compression phase until σ1, σ2, or σ3reached a preset value.In the second step, axial loading was performed with moving rigid wall in the same way as in the uniaxial compression test.Cylindrical samples were reproduced with cubic elements, and the investigated sizes were based on application of the calibrated material model to medium- and large-scale modeling.To ensure reliable transmission of the determined parameters from laboratory-to “field”-scale simulations, five mesh sizes,i.e.0.625 mm,1.25 mm,1.75 mm,2.5 mm and 5 mm,were tested.The results were compared with the experimental outputs qualitatively and quantitatively.

Table 3Summary of JHC constants and data required for their determination.

Fig.13 compares the strength curves from the uniaxial and triaxial tests with the experimental results for all investigated mesh sizes.No significant impact of element counts was observed formaximum strength and post-peak softening.An immediate drop in strength was registered under uniaxial compression, and material softening was less intense compared with the single-element tests due to fracture of multiple elements on planes affected by shearing.This shearing was an effect of friction between the sample and the testing machine.In Fig.13a,the corresponding damage patterns are presented (D ＝1).The densification of the FE model clarified that the failure mechanism was multiple shearing of the sample (Basu et al., 2013), and the constant angle of the shear plane was preserved for each mesh size.For an FE size of 0.625 mm,the number of fractured planes was overestimated due to the idealized geometry,mesh symmetry and immediate release of elastic energy when the fracture was initiated.This relaxation of the sample generated some weak tensile waves that increased the calculated damage.As the mesh size increased, the effect of material model idealization and homogenization became more apparent; despite increase in the fully damaged volume, the material strength remained almost insensitive to mesh size.

Fig.12.Damage parameters for triaxial single-element compression tests.

By comparison, the results of the triaxial compression tests indicated greater sensitivity of the model to mesh size.The element size affected the material softening phase that became smoother with increasing mesh size.In this test, no significant difference in maximum strength was observed.The presence of 25 MPa confining pressure increased the maximum strength by approximately 25% compared with the uniaxial test results.The residual strength was approximately 60%of the maximum value,which was assumed to be correct for confining pressures up to 100 MPa(Heard et al.,1973;Cie′slik,2007).The JHC model does not support a change in failure mechanism from brittle to ductile with increasing radial pressure,thus adopting a value of 60%is a compromise that permits relatively reliable reproduction of the dolomite strength response under a wide range of confining pressures.

The general fracture patterns shown in Fig.14 were in acceptable agreement for small mesh sizes.However, for coarser mesh sizes above the element size of 1.75-2 mm,large differences were observed.The number of main shear surfaces through the whole sample varied, potentially due to the applied explicit integration scheme,but the angle of the shear surfaces remained constant in all simulations.Because the implicit method does not support the JHC material model, the results cannot be verified by this method.Above an element size of 2.5 mm, the shear plane cannot be recognized.The volume of fractured material increased with increasing element size.The failure mechanism is multiple cracking(Basu et al., 2013).As in uniaxial compression tests, the results of the triaxial compression tests showed satisfactory agreement with the experimental outputs;therefore,it can be assumed that the JHC parameters were properly determined and calibrated for static loading conditions.

5.2.SHPB dynamic tests

Other researchers have confirmed that during deformation at high strain rates,a leading failure mechanism is an increase in the volume of crushed material, followed by propagation of a large number of short fractures in the sample(e.g.Donzé et al.,1997).To reproduce the experimental setup, a whole SHPB stand was modeled, as shown in Fig.15.The dolomite sample was placed between two bars with a length of 3190 mm made of C45 highstrength steel.The interaction process between all components was numerically defined using the penalty-based contact algorithm.Due to the use of lubricant in the actual tests, friction was neglected in finite element analysis (FEA).As contact plays a significant role in such simulations,the stiffness factor of contact was increased to minimize penetration.Additional effect presented in dynamic compression that partially results from friction and inertial forces is lateral confinement of sample.To reduce the computational time,both the incident and transmission bars were meshed with cubic elements of 1 mm in size near the gages and contact zones; for the other parts of the bars,10 mm elements were used.The striker was omitted,and the load was applied to the front face of the incident bar to register pressure vs.time relationship in the tests.

Stress uniformity should be preserved in samples during dynamic tests to ensure the validity of the results.Based on one-dimensional wave theory of SHPB tests, the sample is in equilibrium when the criterion given by Eq.(13) is met (Hudson and Ulusay, 2007):

Fig.13.Mesh patterns for the sensitivity study of damage parameters: (a) Uniaxial compression and (b) Triaxial compression at a confining pressure of 25 MPa.

Fig.14.Damage patterns depending on the damage parameter D1.

Fig.15.FEA setup of the SHPB for the dynamic compression test.

where εi, εrand εtare the signals of incident, reflected and transmitted waves,respectively.The force balance was verified for all the tests, and the results for three modeled strain rates, i.e.66 s－1, 94 s－1and 124 s－1, were calculated from Eq.(13) and are shown in Fig.16.Good agreement of the force balance was recorded at strain rates of 66 s－1and 94 s－1,while the stress history at the front of the sample overlapped slightly at 124 s－1.In addition, the force equilibrium at the constant stress stage of deformation was in good agreement, thus the test results were assumed to be correct.The fracture of the cylindrical samples is presented later with the FEA results.

The results from the tests and numerical simulations were calculated in the same way according to Eqs.(14)-(16) (Hudson and Ulusay, 2007; Baranowski and Ma?achowski, 2018):

where σzis the stress in the sample; ˙εzis the strain rate; εzis the axial strain; A0and Asare the cross-sections of the bars and thesample, respectively; E0is the elastic modulus of the transmission bar;c0is the speed of sound in the bar;and Lsis the sample length.

Fig.16.Stress equilibrium in dynamic tests for three different strain rates.

The strain rate in each test using the SHPB was calculated as the average strain rate within the range starting at 75% of the maximum stress and ending at the maximum stress.An identical approach was implemented by Paja?k et al.(2019).In FEA,the strain rate is computed in each element at every time step separately.The strain rate depends on the material strength,its impedance and the shape of the incident wave, and thus local strain rates within the sample may vary up to a dozen percent from the experimental values.The most reliable range at which the calibrated model is applicable depends on the data used for calibration.When the actual strain rate from the simulation was higher than the maximum experimental data used for calibration, the strength of the material was extrapolated and scaled by the JHC strain rate hardening model.This approach limits rapid growth in stiffness by use of the natural logarithm when a small change in the strain rate occurs but may be responsible for the difference between the maximum strength and the fracture pattern described below.

The relationships between stress and strain from the dynamic compression tests are presented in Fig.17 for each investigated strain rate.The maximum strength obtained from FEA was in good agreement with the experimental outcomes.The results are summarized in Table 4.The maximum calculated relative error did not exceed 5%,confirming that the application of two separate material parameter sets for static and dynamic loadings can reproduce the strengthening of dolomite.Slight disagreement in the shape of the compared strength curves was observed at a strain rate of 66 s－1;the material lost its stiffness and continuity distinctly faster than that in the experiment.The fracture, which is represented by the fully damaged material,tended to be noisier in FEA due to the lower slope and maximum value of the incident wave compared with higher strain rates.Local strain rates were measured in this test at three positions of the sample: the two ends and the middle.The results varied from 51 s－1to 115 s－1.As these differences influenced the hardening intensity, the range of maximum strength can vary locally by up to 50 MPa.In other tests, this difference was significantly smaller.

Fig.17 compares the FEA results with the fracture of the samples registered with a fast camera during the experiment.The crack propagation and failure mechanisms were in good agreement.Brittle deformation began near the faces of the sample that were in contact with the steel bars.Subsequently, propagation of longitudinal cracks was observed.The experimental records clearly show that the fracture of the samples was specified by pre-existing microcracks.In addition, as mentioned in Section 2, calcite inclusions in the dolomite structure(visible on the boundary surface in Fig.17 as grey patches) were crushed despite the lack of preexisting cracks in the neighboring area.Thus, calcite inclusions significantly impact the material strength.To reproduce the physical fracture and loss of continuity,erosion based on principal strain was applied.Elements were deleted when 10% of the strain value was reached.As the erosion occurred in fully damaged material,there was no impact of erosion on the results.The discrepancies between the experimental and FEA results are due to the simplification of the discrete model, which, unlike actual samples, contains no imperfections.After taking this homogenization into consideration, satisfactory correlation was reached.

The effect of inertial forces presented in dynamic uniaxial compression tests is strongly associated with inertial lateral confinement of rock.The stress triaxiality through the sample is uneven, thus the highest stress value was measured close to the sample axis,and decreased in radial direction.When the strain rate is low, the stress variation along the radial direction is less prominent and increases with growing of strain rate.In future studies,the amount of strength increase due to described effect will be considered.

6.Conclusions

This paper presents a detailed procedure of parameter determination for the JHC constitutive model to simulate the dolomite.Experimental results and additional literature data from tests performed under quasi-static and dynamic loading conditions at different states of stress were used.Due to the strain rate sensitivity of dolomite, which can be divided into quasi-static and dynamic ranges, two separate sets of parameters were proposed for strain rates below and above 1 s－1.The model was validated under static uniaxial and triaxial compression as well as under uniaxial dynamic compression.Based on the results,the following conclusions can be drawn:

(1) Reliable calibration of the JHC model requires a series of experiments performed under different stress states and strain rates to fully describe the hydrostatic and deviatoric deformation of the rock.

(2) The linear logarithmic strain rate hardening rule applied in the JHC model must be oversimplified to reproduce the response of dolomite rock and its fracture under static and dynamic conditions.Applying the bi-linear model as two separate sets of material constants improved the accuracy of the results from static and dynamic tests.

(3) The calibration of damage parameters is crucial to reasonably reproduce sample fracture.The proposed set of parameters D1, D2and EFMIN provided satisfactory reproduction of damage under quasi-static and dynamic loading conditions.

(4) A strong impact of element mesh size was observed for the qualitative reproduction of dolomite deformation.Whenlarge elements over 2.5 mm were used, the volume of the fractured material began to be overestimated.However,there was no significant impact of mesh size on the quantitative results.

Fig.17.Comparison of fracture (damage) and stress-strain histories from the experiment and FEA studies for three different strain rates.

Table 4Summary of FEA and experimental (EXP) results from uniaxial dynamic compression tests.

Further investigations will include the development of a modeling method that addresses the inability to reliably reproducethe tensile damage in the JHC material model,which has also been reported by other authors(e.g.Holmquist et al.,1993;Polanco-Loria et al., 2008; Lu et al., 2012; Islam et al., 2013; Kong et al., 2016; Li and Shi, 2016; Ren et al., 2017), and application of the model to small- and large-scale blasting simulations.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors would like to thank Prof.Jacek Janiszewski, Dr.Roman Gieleta, Dr.Pawe? Bogusz, Dr.Micha? Stankiewicz, and Dr.Marcin Konarzewski for conducting experimental research.

List of symbols

εiSignal of incident wave

εSignalof reflected wave

r

εtSignal of transmitted wave

εzStrain in sample from SHPB tests

μcrushCrushing volumetric strain

μlockLocking volumetric strain

μplockVolumetric strain under which the material is fully compacted

ρ0Initial density of material

σnStress in normal direction

σsShear stress

σzStress in sample from SHPB tests

σ, σ*Equivalent stress,normalized equivalent stress

ΔεpEquivalent plastic strain increment

ΔμpEquivalent volumetric strain increment

ν Poisson’s ratio

φ Slope of the Mohr-Coulomb envelope

ρ Current density of material

A0Cross-section of SHPB bars

AsCross-section of SHPB sample

A Normalized cohesion

B Pressure hardening coefficient

c0Sound velocity in the SHPB bars

Cstatic／dynamicStrain rate coefficient (for static and dynamic strain rates)

D Accumulated damage coefficient

D1, D2Damage coefficients

E0Elastic modulus of SHPB bars

E Young’s modulus (elastic modulus)

EFMIN Minimum plastic strain before fracture

G Shear modulus

J2Second invariant of deviatoric stress tensor

K Elastic bulk modulus

K1Unloading bulk modulus

K2,K3Pressure constants for EOS

LsLength of SHPB sample

N Pressure hardening exponent

PcrushCrushing pressure

PlockLocking pressure

P, P*Pressure, normalized pressure

smaxMaximum displacement in explicit simulation

TvelDuration of explicit simulation

Tstatic／dynamicMaximum tensile hydrostatic pressure(for static and dynamic strain rates)

t Time