A novel time-dependent damage constitutive model for hard rock materials

The mechanical behavior of deep hard rocks plays a crucial role in the long-term safety and stability of engineering structures. This work focuses on studying the instantaneous and time-dependent fracture propagation of hard rock materials from a theoretical perspective. A new unified constitutive model, based on a modified Mohr–Coulomb (M-C) criterion, is proposed to accurately represent the short and long-term mechanical properties of hard rocks. To capture the strain hardening and strain softening behaviors, damage is utilized as an internal mechanism-driven variable, controlling the expansion and contraction of the plastic yield surface. Additionally, a combination of time-dependent damage law and viscoplastic theory is employed to account for nonlinear creep deformation characteristics. By considering the time-dependent effects, the model can be applied to both instantaneous loading and creep conditions. The general algorithm format is derived in detail, and an explicit integration algorithm is utilized to update the time-dependent damage evolution. Finally, the proposed model is validated by comparing its predictions with the short and long-term mechanical responses of Beishan granite and Rumei dacite. This comprehensive constitutive model improves our understanding of continuum damage mechanics and provides a scientific basis for analyzing and evaluating the long-term safety and stability of deep hard rock engineering projects. A novel integrated constitutive model is proposed to accurately capture the short-term and long-term mechanical responses of hard rocks. Damage is chosen as an internal mechanism-driven variable to govern the expansion and contraction of the plastic yield surface. An explicit integration algorithm is utilized to effectively update the time-dependent damage evolution A novel integrated constitutive model is proposed to accurately capture the short-term and long-term mechanical responses of hard rocks. Damage is chosen as an internal mechanism-driven variable to govern the expansion and contraction of the plastic yield surface. An explicit integration algorithm is utilized to effectively update the time-dependent damage evolution


Introduction
Hard rock materials, such as granite and basalt, are commonly used as geological barriers in nuclear waste repositories (Chen et al. 2017) and as building materials in hydroelectric dam foundation engineering (Shi et al. 2020).In deep high in-situ stress conditions, the time-dependent behavior of these rock materials can lead to delayed failure of structures.As human activities extend deeper into the Earth, understanding the time-dependent behavior of rock materials becomes crucial for predicting the long-term safety and stability of underground rock engineering.The mechanical behaviors of deep-buried hard rocks are influenced by factors such as high geo-stress level, temperature, excavation disturbance, and time effects (Blümling et al. 2007;Chen et al. 2017;Yang et al. 2014;Zhao et al. 2022).The time-dependent behaviors of surrounding rock can significantly impact the service life and normal operation of the practical engineering projects.Therefore, it is essential to develop a comprehensive and unified creep constitutive model to evaluate and analyze the stability of such engineering projects over time.
Rock is a rate-dependent material, meaning that its strength and failure behavior are influenced by the external loading rate (Hashiba and Fukui 2015).The rheology of rock materials, which refers to their time-dependent behavior, has been a subject of interest of many researchers (Li and Xia 2000;Liu et al. 2018).To understand the creep behaviors of rock, abundant in-situ failure investigations (Jiang et al. 2019) and laboratory tests have been carried out on various types of rocks, including argillite (Fabre and Pellet 2006), sandstone (Yang and Hu 2018), marble (Yang et al. 2015), and granite (Liu and Xu 2015) under constant stress conditions.From these tests, it has been observed that there are two types of creep failure in rock materials: brittle creep and ductile creep (Zhao et al. 2019).Brittle creep failure in rock materials refers to a sudden increase in deformation that occurs when the sample fails under an accelerating strain rate.This type of creep undergoes a transition from steady deformation to spontaneous failure, with the process consisting of three stages: transient creep, steady creep, and accelerated creep.The steady creep rate is influenced by the applied stress level and increases as the stress level increases.Additionally, the accelerated creep failure in hard rocks is caused by the evolution of damage resulting from interactions between microcracks.The time-dependent behavior of rock materials typically occurs when the applied stress exceeds a certain threshold.The subcritical growth of microcracks leads to strength and stiffness degradation, causing the creep failure strength to be generally lower than the peak strength under instantaneous loading.The time-dependent mechanical properties of rock materials vary significantly due to differences in mineral composition, grain size, and orientation of bedding planes (Gupta and Mishra 2020).Furthermore, the response of creep failure depends on the stress level, loading path, and environmental conditions (Saksala and Ibrahimbegovic 2014;Taheri et al. 2020;Yang and Hu 2018).These factors play a role in determining the behavior of rock materials during creep failure.
The failure mechanism of rock materials can be investigated using advanced experimental techniques such as acoustic emission (AE) (Browning et al. 2017;Yang et al. 2021), digital image correlation (DIC) (Munoz and Taheri 2017), X-ray computed tomography (CT) (Gupta and Mishra 2020), and scanning electron microscopy (SEM) (Wang et al. 2020b).These techniques allow for the analysis of microcrack initiation, propagation, coalescence, and the evolution of microstructure and microcracks over time.In the case of hard rocks, the nonlinear creep deformation is primarily caused by the initiation, propagation, and coalescence of microcracks.The mentioned test methodologies provide intuitive analysis of the variations in microstructure and the space-time evolution of microcracks.The time-dependent failure mechanism of rock materials can be explained by the subcritical propagation of microcracks within grains and along grain boundaries due to stress corrosion, which progressively deteriorates the material (Brantut et al. 2013).Therefore, the nonlinear time-dependent deformation of materials is governed by two competitive dissipation mechanisms: viscoplastic accumulation and time-dependent damage evolution.Viscoplastic dissipation is dominant in soft rocks, while time-dependent damage dissipation should be given more attention in the case of hard rocks.
Constitutive modeling is a valuable tool for understanding and describing the failure mechanism of rock materials.By utilizing a reasonable model, construction and design schemes can be optimized based on numerical simulation results.Continuum damage mechanics (CDM) (Murakami 2012) provides a scientific framework for modeling, with internal variables representing the nonlinear behavior of materials (e.g., plastic strain and damage variable).Yield or failure criteria are determined based on experimental results, and the constitutive model is derived within a thermodynamic framework.Incorporating ratedependent properties into the yield criteria allows for linking time-dependent behaviors with the constitutive model.For hard rocks, where macroscopic shear failure is the predominant mode during compressive failure (Choens et al. 2021), equivalent plastic or viscoplastic shear strain (or its work) can be used as an internal variable.Developing a unified constitutive model within the thermodynamic framework is crucial to accurately describe progressive damage and degradation mechanisms (Al-Rub and Darabi 2012).However, despite progress, understanding the timedependent mechanical behavior of deep hard rocks and improving time-dependent constitutive models are ongoing challenges for scientists.The theoretical framework of these models is still maturing in some cases, necessitating further improvements for examining the mechanical behaviors of hard rocks over time.
Based on the description of the two dominant mechanisms, rheological constitutive models in rock mechanics can be classified into five categories: empirical models (Cornelius andScott 1993), component models (Song et al. 2022), viscoelastic (Schapery 1969) /viscoplastic models (Borja et al. 2020), and time-dependent damage models (Yuan et al. 2020).Empirical models describe representative phenomena through segmental function fitting based on experimental results.Component models capture nonlinear rheological deformation behaviors by combining elastic elements, plastic elements, damage bodies, and Newtonian dashpots in series or parallel.However, model parameters that change with stress state need to be determined through numerical fitting.
Viscoelastic models generalize the one-dimensional (1D) to three-dimensional (3D) stress space (Lai and Bakker 1996).Viscoplastic models utilize the overstress concept and well-established formulations such as Perzyna (Perzyna 1966) and Duvaut-Lions (Duvaut and Lions 1976) to describe viscoplastic deformation in rock materials (Choo et al. 2021).Timedependent damage models use a progressive damage variable to effectively represent the gradual growth and propagation of microcracks in rocks.These models successfully reproduce the primary, secondary, and tertiary creep deformation stages of hard rocks (Cheng 2022;Cheng et al. 2022;Zhao et al. 2016).However, few studies have considered the coupling mechanism between damage and deformation in rock models.Pure viscoplastic or time-dependent models can only account for one mechanism of creep failure, sometimes leading to overestimation of creep deformation in hard rock materials.Furthermore, thermodynamic consistency is often overlooked during the modeling process.Multiscale time-dependent damage models (Huang and Shao 2013) and joint creep models (Wang and Cai 2021) have been developed to consider both microscopic and macroscopic crack behaviors, but these exceed the scope of this study and won't be elaborated on here.Therefore, the most effective approach so far is to combine timedependent damage and viscoplastic deformation models, considering the physical mechanisms of rock rheological failure.
The objective of this study is to develop a unified constitutive model for hard rocks within the thermodynamic framework.The model incorporates timedependent damage as an internal state variable to capture the progressive deterioration process of hard rocks.Furthermore, it considers the nonlinear creep deformation behavior by coupling the time-dependent damage with viscoplastic deformation.To calculate and update the time-dependent damage factor, the study employs an explicit integration algorithm proposed by Zhao et al. (2016).Numerical simulations of laboratory tests are conducted using a specific set of model parameters.The proposed model's efficiency and reliability are then assessed by comparing its theoretical predictions with the experimental results obtained from tests on Beishan granite and Rumei dacite.
In this study, the convention is to consider compressive stress and strain tensors as positive.Scalars are consistently represented by italic light-face Greek symbols (e.g., k), while vectors and second-order tensors are indicated by bold letters (e.g., ).The fourth-order tensor is represented by blackboard-bold uppercase characters (e.g., ℂ ).The notation for ten- sor products of any second-order tensors and is denoted as follows: 118 Page 4 of 20 Vol:. ( 1234567890)

Mechanical variables description
In this section, we introduce the description of basic mechanical variables, namely deformation and damage.The study adopts the small deformation assumption, quasi-static state, and isothermal condition.Plastic and damage dissipations are assumed to be the inherent mechanisms of rock deformation and failure.Plastic dissipation is associated with the plastic sliding that occurs within the microstructure, while damage dissipation is caused by the propagation of microcracks.Phase III (accelerated creep): During this stage, the creep strain sharply increases with an accelerating creep strain rate over time.The creep deformation behavior exhibited by hard rock materials starts to accelerate with a distinct time-dependent effect.Eventually, the rock materials fail as macrocracks coalesce and form, either single or multiple cracks.The duration of this phase is typically within a few minutes or less for hard rocks, while it may not occur for soft rocks.

Deformation formulation
To explicitly describe the deformation behavior, the strain components of nonlinear deformation are decoupled and separated.The total strain tensor is expressed as the sum of three components using the superposition concept: elastic strain tensor e , instantaneous plastic strain tensor ip , and viscoplastic strain tensor vp .Therefore, the total strain and its incremental formulation can be written as follows: where G represents the plastic potential function.H and m are the viscoplastic parameters.ip and vp are the instantaneous plastic and viscoplastic multipliers, respectively.It is worth noting that H should have the same dimension as f p to maintain dimensional con- sistency.Empirical parameters H and m are utilized to more accurately capture the viscoplastic behaviors.The symbol ⟨⟩ represents the Macaulay brackets, which take the positive format: ⟨x⟩ = (x + �x�)∕2.
During the creep process, the equivalent plastic and viscoplastic shear strain p is composed of two parts: the instantaneous part ip and the timedependent part vp .It can be expressed as follows: (1) The incremental formulation of the instantaneous and time-dependent equivalent plastic shear strain can be expressed independently using Eq.(5a) and Eq.(5b).
where "dev" is interpreted as the deviatoric part of the second-order tensor.

Damage definition
In Fig. 2a, the natural rock is composed of a matrix and microdefects such as microcracks, voids, and porosity.In the framework of continuous damage mechanics, the damage variable is utilized to quantify the initiation and propagation of microcracks, thereby describing the macroscopic damage.serves as an internal state variable that establishes the relationship between the current configuration (A) and the undamaged configuration ( A 1 ).Note that when equals to zero, the material is undamaged, whereas when equals to one, the rock sample is completely damaged and loses its load-carrying capacity.For simplicity, an isotropic damage model is adopted in this study, with defined as the ratio of crack area to the total area.
where A 2 represents the completed damaged zone.As mentioned before, the progressive deterioration and failure of rock materials are primarily caused by the subcritical growth of microcracks.The time-dependent damage variable is a physical mechanism that contributes to the nonlinear timedependent deformation behavior.As shear cracking plays a significant role in the gradual failure of hard rocks (Yu et al. 2021), it is assumed that the damage variable is closely associated with the equivalent plastic shear strain .Throughout the entire rheological process, the total damage can be divided into two components: viscoplastic/plastic deformation-induced damage pd and time-dependent damage cd .It can be expressed as follows: where pd ( p ) can be further decomposed into two parts: the instantaneous plastic strain-induced damage id ip and viscoplastic strain-induced damage vd ( vp ) .cd can be calculated using the time-depend- ent damage factor .These components characterize the progressive degradation of rock materials due to plastic and viscoplastic deformation, respectively.
The time-dependent damage evolution is a nonequilibrium asymptotic motion process that occurs near the equilibrium state.The parameter is used to characterize the progressive microstructure evolution over time.On the other hand, is defined as a stationary variable associated with the equilibrium state.To simplify the analysis, a linear relationship in the form of ( ) − (t) is postulated to measure the kinetics of microstructure rearrangement.This linear format allows for the quantification of the rate at which the microstructure undergoes rearrangement during the deformation process.
where represents a constant parameter that signifies the rate of microstructure damage evolution.By employing the Laplace transform and convolution theorem, can be expressed using an integration format as shown in Eq. ( 9).For a more comprehensive derivation, please refer to the relevant literature (Pietruszczak et al. 2004).We assume that the damaged area loss its bearing capacity and the entire load is borne by the undamaged zone.The effective stress tensor, denoted as ̃ , represents the actual stress acting on the undam- aged material.On the other hand, the nominal stress tensor, denoted as , is defined by the applied mechanical loading on the entire zone of material.The relation between ̃ and can be determined by employing hypotheses of strain equivalence (Lemaitre 1984), stress equivalence (Al-Rub and Darabi 2012), or energy equivalence (Sidoroff 1981).In the case of isotropic damage, the assumption of strain equivalence, which states that the total strain is equal between the current configuration and the undamaged configuration, is widely adopted due to its simplicity and its ability to meet the modeling requirements.To reduce the complexity of model, the hypothesis of strain equivalence is adopted in this study.Therefore, the relation between ̃ and can be expressed as follows:

Thermodynamic framework
For the sake of theoretical completeness, the constitutive model is derived and verified within the framework of thermodynamics in this section.In accordance with the second law of thermodynamics, the Clausius-Duhem inequality must be satisfied to ensure that the dissipation is nonnegative.
where is mass density, is the free energy density, T is temperature, h is heat flux, and s is entropy.
Since the assumption of an isothermal condition is made, the dissipation potential can be further expressed as follows: The free Helmholtz energy density depends on the stress state, damage degree, and plastic deformation accumulation.Considering the coupling between damage and viscoplastic / plastic strain, can split into three components: elastic energy density e , instantaneous plastic energy density ip , and viscoplastic energy density vp , namely: The elastic energy density e is expressed in Eq. (14a) as a quadratic form of the elastic strain tensor.On the other hand, the plastic or viscoplastic dissipation potential function is typically determined based on assumed mathematical formulations that should ideally satisfy the superposition principle.The expressions for ip and vp , presented in Eqs. ( 14b) and (14c) respectively, take into account the coupling between the damage variable and the equivalent plastic strain p .
where is a material parameter that affects the nonlinear dissipation of the system.The term ℂ( ) cor- responds to the forth stiffness tensor of the damaged material, which depends on the damage variable.
where and K are respectively shear modulus and bulk modulus.The fourth tensor and are defined as follows: K ijkl = I ijkl − J ijkl and J ijkl = 1 3 ij kl .Here, denotes the Kronecker Delta and fourth-order unit tensor can be expressed as I ijkl = 1 2 ik jl + il jk .By substituting Eqs.(14a)-(14c) into Eq.( 12), the nonnegative dissipation inequality can be expanded as follows: Therefore, the problem of development the material constitutive model is transformed into a problem of determining the forms of the free energy function and the corresponding dissipation constraints. (13) Page 7 of 20 118 Vol.: (0123456789) The constitutive relation and thermodynamic driving forces can be subsequently obtained from the following equations.
where p represents the thermodynamic force conjugate to ip and Y is thermodynamic force conjugate to damage.The conditions 0 ≤ ≤ 1.0 and a posi- tive value guarantee that both p and Y are always positive.Therefore, by substituting Eqs. ( 2), ( 3), ( 17), ( 18), and (19) into Eq.( 16), the nonnegative dissipation condition is satisfied.

Yield function and time-dependent damage criterion
Rock engineering is typically subjected to compressive stress conditions, and the main mechanism that induces failure is compressive shear.Based on the experimental measurements of hard rocks, such as dacite Wang et al. (2019), it was found that the M-C criterion is suitable for describing the behavior of the rock under the compressive stress condition.The M-C criterion in the shear stress and normal stress plane can be expressed by where represents the friction coefficient, which is tangent function of fraction angle ; c is the cohesion, (MPa); shear stress = ⊗ ⋅ ( − ⊗ ) ; and the normal stress n = ⋅ ⋅ .The Eq. ( 20) in the plane of ( 1 , 3 ) can be expressed as: where k is the slope of strength envelope in the plane.
For hard rocks, brittle failure is characterized by the accumulation of plastic strain and the evolution of damage.Additionally, the rapid increase in the rate of damage is the main reason for the post-peak strain softening observed in hard rocks.Taking inspiration from the CWFS model proposed by Hajiabdolmajid et al. (2002), the damage variable is chosen as an internal mechanism-driven variable to control the change of the yield surface.This allows for capturing both strain hardening and strain softening behaviors.Therefore, the plastic hardening/softening function ( ) is incor- porated into the yield criterion to control the contraction and expansion of the yield surface, which can be expressed in the following format: where c is the criterial damage parameter.Parameters b and n are used to characterize the rate of nonlin- ear strain hardening or softening behaviors.The transition from brittle failure to ductile deformation behaviors can also be accounted for by adjusting these two parameters.Selecting b = 2 and n = 1 uni- formly for hard rock materials can reduce the difficulty of model parameter identification (Wang et al. 2022b).Moreover, the rate of hardening or softening increases as c decreases, while the ductile behavior becomes more pronounced as c increases.The crite- � can be obtained by using the extremum method of ( ) .
The criterial damage serves as a quantitative indicator to determine if unstable failure of hard rocks occurs.
For the sake of simplicity, k and c are assumed to follow the same trend as the change in in order to minimize the number of model parameter identifications.Furthermore, when the associated flow law is adopted, the influence of damage can weaken the deformation of plasticity in the direction of the major principal stress to some extent.This assumption can be expressed as follows: where ( ) max represents the maximum of ( ) .Parameters k 0 and c 0 correspond to the slope and cohesion at the peak strength, respectively.
To ensure algorithmic convenience and eliminate the need for additional model parameters, an Geomech.Geophys.Geo-energ.Geo-resour.associated plastic flow rule ( G = f p ) and a viscoplas- tic flow rule are adopted in this approach.
where ̇ in represents the inelastic strain increment, which includes both instantaneous plastic deformation and viscoplastic deformation.
Inspired on previous studies (Wang and Xu 2020), an exponential damage criterion is adopted to describe the progressive degradation of rock material properties.
where max defines the maximum value of the plastic damage, and Y represents the damage conjugate force.The damage conjugate force is defined as the partial derivative of the free energy with respect to the damage variable within the thermodynamical framework, as shown in Eq. ( 19).
where parameter is used to control the evolution rate of damage Y consists of the elastic part Y e and either the plastic part Y ip , and the viscoplastic part Y vp , depending on whether it involves plastic or visco- plastic deformation.

General format of creep model
When the stress level is within the elastic range, only the time-dependent variables need to be updated.The deformation and stress can be calculated using Hooke's law.However, when the stress state exceeds the plastic yield surface, the incremental forms of the ( 24) yield function and damage criterion can be expressed as follows: The incremental damage variable and damage driven force are respectively written as follows: The total deformation of the rock (in Eq. ( 1)) can be divided into instantaneous deformation and timedependent deformation.Therefore, the unified timedependent constitutive model is derived by considering the instantaneous and creep stages.

Instantaneous stage
Under instantaneous loading conditions, the timedependent parameters are given as d = 0 , λcd = 0 , and λcp = 0 .The instantaneous elastic strain and plas- tic strain are updated based on the applied external load.The elastic strain can be obtained using the generalized Hoek's law shown in Eq. ( 33).The instantaneous plastic strain is calculated using the plastic flow rule in Eq. ( 34) when the stress state is outside of the yield surface.Otherwise, the incremental plastic strain is set to zero.Similar to plasticity theory, the instantaneous rate of damage evolution is obtained using Eq. ( 35).
The instantaneous plastic multiplier can be obtained by substituting Eq. ( 31) and Eq.(34) into the plastic consistency condition of Eq. ( 29).The increment of the instantaneous plastic multiplier is derived as follows: where simplified symbolic forms are given by p = f p and In a similar manner, the instantaneous damage multiplier can be derived by substituting Eq. ( 32) and Eq. ( 35) into the plastic consistency condition of Eq. ( 30).The increment of the instantaneous damage multiplier is derived as follows: where simplified symbolic forms are given by Therefore, the incremental formulation of instantaneous strain can be expressed as follows: 2 Creep deformation stage When a rock sample is under the creep process, the stress remains constant ( ̇ = 0 ).Based on Eq. ( 36) and Eq. ( 37), the increments of the instantaneous plastic multiplier and instantaneous damage multiplier are equal to zero ( λip = 0 and λid = 0 ).The increments of the time-dependent plastic strain, time-dependent equivalent plastic shear strain, and time-dependent deformation-induced damage can be expressed as Eq. ( 39), Eq. ( 40), and Eq. ( 41), respectively.Therefore, the increment of the time-dependent plastic multiplier λvp can be derived by substituting Eqs. ( 39)-(41) into Eq.( 29), resulting in the following expression. (36 Therefore, the time-dependent strain increment can be obtained by substituting Eq. ( 42) into Eq.( 39) during the constant stress level.Considering the coupling of time-dependent deformation and damage evolution, the time-dependent damage multiplier can be derived by Eq. ( 30), resulting in the following expression.
where simplified symbolic forms are given by The creep deformation increment can be written as follows: The creep constitutive model, which describes the relationship between strain and stress, can be obtained by superimposing Eq. ( 38) and Eq. ( 44).Following the approach used in previous literature (Shao et al. 2003), the incremental formulation of total strain can be expressed as two parts in Eq. ( 45).These parts correspond to the increment of stress and the increment of the time-dependent damage factor, respectively.The general format of the stress relaxation model can also be derived using a similar method.However, due to space limitations, the derivation of the relaxation model will not be included in this paper.As a result, the creep deformation of rock is primarily influenced by the stress level and the time-dependent damage.where ( 43) Geomech.Geophys.Geo-energ.Geo-resour.Vol:.( 1234567890)

Numerical implementations
The accurate simulation results rely on the determination of time-dependent damage evolution.According to Liang et al. (2021), creep test results show that the rate of time-dependent damage is influenced by the current stress level and environmental conditions, such as temperature, humidity, and chemical concentration during the creep process.In the case of hard rocks, the creep deformation will eventually reach an asymptote, and the sample will not fail as long as the stress level remains below the threshold.Stress-induced creep deformation and failure only occur when the stress level surpasses the threshold.
To model this behavior, Eq. ( 47) assumes an initial equilibrium state, denoted as , using a sigmoid function.Over time, the time-dependent damage evolution gradually converges towards this state. (46b) where the parameter χ , which is associated with the rock lithology and stress level, plays a crucial role in controlling the evolution rate of time-dependent damage.Its effect on can be observed in Fig. 3.
The results demonstrate that the curves exhibit an S-shaped pattern, with a maximum value approaching 1.0.Furthermore, the initial nonlinear stage decreases as the parameter χ decreases.Since the time-depend- ent damage rate is influenced by the actual stress state, the parameter χ is determined by the exter- nal stress level.Specifically, it can be expressed as χ = χ 0 exp −B Δσ∕Δ 0 , where Δσ represents the absolute difference between the current stress state and the failure strength, while Δ 0 is calculated as the absolute difference between the failure strength and the inviscid stress threshold.
The explicit nonlinear integration algorithm can be used to express the integral formulation of the timedependent damage factor in Eq. 9.
During the entire creep test, the stress is initially applied to a predetermined stress level at a constant stress rate.The time-dependent responses are then monitored under this constant stress condition.As a result, the numerical implementation process consists of two parts: the instantaneous part and the time-dependent part, which aligns with the practical test process.At the end of time t k , the basic mechani- cal variables are assumed to be known.To simulate the creep behaviors of rock materials, the required mechanical variables at time t k+1 are updated at the Gauss points using the stress-control regime.The flowchart depicted in Fig. 4 illustrates the step-bystep procedure for implementing the time-dependent model at the material point level.
The influence of the time incremental step dt k on the numerical simulation results is examined in Fig. 5.The findings suggest that there is only slight fluctuation observed during the accelerated creep stage, and the overall trend of the creep failure characteristic remains largely unaffected by dt k .An inter- esting observation is that the creep strain curve initially shifts to the left and then to the right along the time axis as dt k increases.( 48)

Parameter identification
The model consists of seven instantaneous material parameters.The deformation parameters, such as the elastic modulus E and Poisson's ratio v , can be obtained by analyzing the slope of the stress-strain curve during the elastic deformation stage.The strength parameters, including and c , are deter- mined by performing a linear fitting of the peak strength using the M-C criterion.However, the damage parameters max , , and c cannot be directly determined based on test data due to the limitations of possible loading paths.max can be determined by the triaxial cyclic loading test.On the other hand, and c are identified by finding a set of parameters that well matches the test results.Consequently, a large number of numerical trials are required for each type of rock material in order to ensure the accuracy of the simulation results.
In addition to the seven instantaneous material parameters, five extra parameters are determined for the creep model.To better understand the influence of these parameters on the time-dependent deformation behaviors, a sensitivity analysis is conducted at a constant deviatoric stress level.The sensitivity analysis of parameter is presented in Fig. 6, indicating that creep deformation failure is accelerated with an increase of .This parameter characterizes the rate of microstructure degradation due to stress corrosion.A Fig. 5 The effect of time increment step on the creep deformation during numerical simulation Fig. 6 The sensitivity of parameter Fig. 7 The sensitivity of parameter H Geomech. Geophys.Geo-energ.Geo-resour.
(2023) 9:118 1 3 Page 13 of 20 118 Vol.: (0123456789) higher value of represents a faster damage evolution rate, which in turn reduces the overall lifetime of the structure.Figure 7 shows the sensitivity of parameter H , which has little effect on the axial creep strain during a single creep stage.However, it is found that H may influence the creep deformation at the tran- sient creep stage, as observed in the multistage creep numerical results presented in the next section.On the other hand, the creep deformation of the model is almost insensitive to changes in parameter m , as depicted in Fig. 8.Moreover, these two parameters do not affect the creep time significantly.Considering the limitations of the simple stress path, parameters H and m are included in this work to improve the uni- versality of the model.Identifying the time-dependent model parameters directly from creep test data is generally challenging.Therefore, the aforementioned studies confirm that the proposed model can capture the brittle creep characteristics of hard rocks when a suitable set of model parameters is recalibrated through trial and error.

Beishan granite
Beishan granite is considered a suitable material for nuclear waste storage due to its low porosity and uniformity.To understand its mechanical properties and time-dependent deformation characteristics, researchers conducted a series of instantaneous and multi-stage creep tests on cylindrical samples (Ming 2019).The test results revealed that the mean elastic modulus and Poisson's ratio of Beishan granite are 45.3 GPa and 0.14, respectively.The M-C criterion is selected because it exhibits a linear relationship with the peak strength (as shown in Fig. 9).Table 1 presents the commonly used model parameters for Beishan granite.It is important to note that the critical damage of the rock is not constant but varies with the stress conditions, as indicated by the rock test results.
To better adapt the model to complex stress states, the value of c is modified based on the confining pressure.A small value of c is used to account for brittle failure under low confining pressure, while a relatively larger value is employed to describe ductile deformation under high confining pressure.
The comparison between instantaneous simulation and test results (in Fig. 10) indicates that the   proposed model can be better capture the main mechanical characteristics of granite, especially overlap approximately under the high confining pressure.
The slight discrepancy post-peak may be that the proposed model ignores the localization damage effect.Moreover, stress-induced anisotropy is also found in the hard rocks.More researches will be further studied in this topic.
The comparison between the instantaneous simulation and test results, as shown in Fig. 10, suggests that the proposed model can accurately capture the primary mechanical characteristics of Beishan granite.Particularly, there is a close agreement between the proposed model and the test results, especially under high confining pressure conditions.However, a slight discrepancy is observed after the peak strength, which may be due to the fact that the proposed model does not account for the localized damage effect.Additionally, it has been found that stress-induced anisotropy exists in hard rocks, as discussed by Feng et al. (2021).Further research will be conducted to investigate this aspect in more detail.
Furthermore, the time-dependent model parameters are selected with the following values: 0 = − 20, B = 1.5, H = 1 × 10 8 , m = 0.1, = 7.0 ×10 −6 , and d c = 0.05 .The parameter is influenced by the dif- ference in axial stress during the loading creep condition.The comparison results between the numerical simulation and creep test are presented in Fig. 11.The comparison results demonstrate that the proposed unified creep model can reasonably capture the primary long-term deformation characteristics of the rock.As the time-dependent damage evolves, it gradually converges to a certain value, resulting in stable creep deformation under low deviatoric stress levels.The creep failure of rock is attributed to the accumulation of time-dependent damage, particularly the sharp increase in damage near rock failure during the accelerated creep stage.Therefore, the creep failure of rock requires the stress level to exceed a critical threshold.These findings align with the conclusions drawn from previous acoustic emission tests conducted on granite during the creep process (Wang et al. 2020a).

Rumei Dacite
Dacite is an igneous rock collected from the Rumei dam site.In a previously published literatures (Wang et al. 2022a;Wang et al. 2019), three groups of instantaneous and creep unloading tests were conducted on dacite.The elastic parameters can be referred to from the aforementioned study.Based on the test results, the friction angle and c 0 are determined to be 60.8° and 16.888 MPa, respectively.However, since the initial strength envelope does not strictly meet the M-C criterion, a parameter called k x is introduced to adjust the value of parameter k, where k = k∕k x , in order to reflect the non-linear strength relationship.To capture the difference between peak strength and initial yield strength, the values of k x are selected as 1.1, 1.1, and 1.0 under initial confining pressures of 5 MPa, 10 MPa, and 15 MPa, respectively.It should be noted that the mechanical properties of dacite are influenced by the unloading path and rock heterogeneity.Therefore, the values chosen for the parameter c are 0.25, 0.25, and 0.4 under various initial confining pressures.
The comparison results between the instantaneous numerical simulation and unloading test data depicted in Fig. 12 demonstrate a good agreement, indicating that the proposed model is capable of capturing the non-linear unloading deformation characteristics of dacite.However, it is worth noting that the discrepancies shown in Fig. 12 may be attributed to the neglect of the unloading anisotropic damage effect.The unloading effect induces anisotropic damage evolution, which implies that the constitutive model should also consider the anisotropy observed in hard rocks.It is a topic that could be addressed in future research endeavors.
Besides, in multi-stage unloading creep tests, the parameter is controlled by the lateral stress difference during the unloading creep condition.The time-dependent parameters are determined through repeated numerical trials and are presented in Table 2.It should be noted that different values of are selected to obtain satisfactory simulation results, considering the heterogeneity of the rock samples.However, the regularity of the parameter has not been further discussed due to the limitations of the unloading creep experimental study on hard rocks.In combination with the monitoring of time-dependent microcrack propagation responses using the AE technique in laboratory tests, an accurate answer may be obtained through analysis using a neural network algorithm.The numerical simulation results, as shown in Fig. 13, exhibit good agreement with the unloading creep test data of dacite.The results also indicate that damage evolution is progressive throughout the entire unloading process.The timedependent damage evolution tends to stabilize at the first few unloading stress levels.However, it accumulates unstably over time and eventually leads to sample failure.Therefore, the proposed model effectively captures the main unloading creep failure mechanism of dacite.Moreover, the model can be applicable in describing the time-dependent behavior of other similar hard rock materials.
This study presents a general and flexible framework for the extension of the stress-induced anisotropic damage model based on the hypothesis of energy equivalence (Wang and Xu 2020).The proposed model effectively captures the main mechanical features of the material, although some discrepancies are observed.One of the main reasons for these differences is the use of a linear M-C yield criterion The comparison between instantaneous numerical and unloading test results of dacite and the neglect of anisotropic damage, particularly notable at low confining pressures.Therefore, future efforts will be dedicated to incorporating a nonlinear strength criterion and considering stress-induced anisotropic damage evolution.Furthermore, the proposed model provides a solid foundation for coupling thermal and hydraulic effects.Given that progressive damage degradation primarily occurs in localized zones, future work will focus on developing a unified localization damage model that can replicate the progressive localization damage observed in hard rocks during the entire failure process.

Conclusions
This study introduces a novel unified time-dependent damage constitutive model tailored for hard rocks.
The key novelty of this model lies in the consideration of damage as an internal variable, which allows for the representation of strain hardening or softening behaviors by controlling changes in the yield surface.Furthermore, the incorporation of time-dependent damage within the modified M-C criterion serves as a primary control variable for assessing the creep failure of hard rocks.The creep deformation characteristics are effectively described through the combination of viscoplastic theory and the principle of time-dependent damage.
An explicit integration algorithm is subsequently employed to update the time-dependent damage factor.This approach effectively mitigates the influence of time increments on the numerical results, ensuring accurate and reliable simulations.The proposed model also offers several advantages, including the consideration of engineering design parameters and its comparatively simple format.These features enhance the practical applicability of the model in engineering applications.To validate the efficacy of the proposed model, simulation results are compared with experimental data from instantaneous and creep tests conducted on Beishan granite and Rumei dacite.Through these comparisons, it is confirmed that the model accurately captures the short and long-term failure behavior of hard rocks, which can be attributed to the progressive accumulation of damage over time.
This work establishes a common foundation for the extension of multiple physical fields.While the initial contribution lies in the development of a unified constitutive model, there are still some limitations that need to be addressed in order to make the model applicable under complex stress conditions.To improve the model's accuracy, further refinement is required, particularly in incorporating anisotropic damage and considering the evolution of localization damage, considering the inherent heterogeneity of rock materials.Additionally, the model will need to encompass the effects of structural characteristics, thus expanding its applicability in practical engineering.

Declarations
Competing interests The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 1
Fig. 1 The phase division of creep strain

Fig. 2
Fig. 2 Damage variable definition of rock materials

Fig. 8
Fig.8The sensitivity of parameter m

Fig. 10
Fig. 10 The numerical simulation of granite in the triaxial compression test

Fig. 11
Fig. 11 The numerical and damage evolution of granite in the creep test ( 3 = 5 MPa)

Funding
This study is financially supported by the National Natural Science Foundation of China (Grant No. 52109143), China Postdoctoral Science Foundation (Grant No. 2022M713376), and the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China Institute of Water Resources and Hydropower Research) (Grant No. IWHR-SKL-KF202305).

Table 2
The time-dependent model parameters of dacite under unloading creep test Vol:. (1234567890)