A Constitutive Model of Time-Dependent Deformation Behavior for Sandstone

Considering sandstone’s heterogeneity in the mesoscale and homogeneity in the macroscale, it is very difficult to describe its time-dependent behavior under stress. The mesoscale heterogeneity can affect the initiation and propagation of cracks. Clusters of cracks have a strong influence on the formation of macroscale fractures. In order to investigate the influence of crack evolution on the formation of fractures during creep deformation, a time-dependent damage model is introduced in this paper. First, the instantaneous elastoplastic damage model of sandstone was built based on the elastoplastic theory of rock and the micro-heterogeneous characteristics of sandstone. A viscoelastic plastic creep damage model was established by combining the Nishihara model and the elastoplastic damage constitutive model. The proposed models have been validated by the results of corresponding analytical solutions. To help back up the model, some conventional constant strain rate tests and multi-step creep tests were carried out to analyze the time-dependent behavior of sandstone. The results show that the proposed damage model can not only reflect the time-dependent viscoelastic deformation characteristics of sandstone, but also provide a good fit to the viscoelastic plastic deformation characteristics of sandstone’s creep behavior. The damage model can also reproduce the propagation process of mesoscopic cracks in sandstone upon the damage and failure of micro-units. This research can provide an effective tool for studying the propagation of microscopic cracks in sandstone.


Introduction
Many studies have demonstrated that quartz-rich rock can deform and eventually fail under a constant differential stress (σ 1 -σ 3 ) over extended periods of time, a phenomenon known as brittle creep (or time-dependent deformation, abbreviated to TDD in the following) [1], especially for sandstone (a rock often seen in the roofs and floors of coal seams) [2,3]. Sandstone taken from deep coal mines and deformed in creep experiments generally exhibits three regimes of strain vs. time: (1) primary creep, or decreasing strain rate stage; (2) secondary creep, or constant strain rate stage; (3) tertiary creep, or increasing strain rate stage ending in failure [4,5]. However, some researchers have divided the creep behavior of rock into three sections: elastic instantaneous deformation, viscoelastic TDD, and viscoplastic TDD [6][7][8][9][10][11]. This categorization-identifying three kinds of deformation-is convenient for building a constitutive model that can be used to reproduce the creep deformation of rock and further understand the failure evolution of rock.
There has been considerable research on the viscoelastic plastic deformation of rock. Shao et al. [12,13] proposed a damage model to describe the elastic-plastic short term behavior and the elastoplastic time-dependent behavior. Their damage law is described in the framework of irreversible thermodynamics. The proposed constitutive relationship could be very suitable for semi-brittle rock material with the help of non-associated plastic flow. Zhao et al. [11,14] have separated the total creep strain into instantaneous elastic strain, instantaneous plastic strain, time-dependent viscoelastic strain and time-dependent visco-plastic strain based on experimental data. Additionally, Prof. Zhao proposed a nonlinear elastoviscous plastic creep model on the basis of element models. The results reveal that the strain curves from the creep model agree well with experimental results and give a precise description of the three stages of creep. Additionally, numerous papers have been published on the description of full stages of rock creep behavior [10,[15][16][17][18][19][20][21][22][23][24][25]. However, most of the published papers focus mainly on strain-time curve fitting or seek to reproduce the full stages of creep behavior. Even though all of them have achieved good agreement with experimental results, we still have no idea how crack evolution influences the formation of fractures in the rock. As a result, some researchers have studied crack evolution before rock failure, going beyond matching the numerical and laboratory strain curves. Xu et al. [26,27] adopted a Norton-Bailey equation to construct the creep constitutive model, which was able to characterize the time-dependent creep deformation quite well. The creep model could also illustrate the failure of rock mass slope using the Norton-Bailey creep equation and the time-independent damage evolution law. In addition, Fu et al. [28] have constructed a three-dimensional discrete element grain-based model to simulate the fracture evolution of rock around underground excavations. The simulation results agree well with those found in the laboratory. Wang et al. [25,29] built a grain-based stress corrosion model based on the distinct element method to study time-dependent failure during creep tests for brittle rock. The crack propagation and the damage evolution induced by stress concentration before the creep failure are discussed in these papers. In addition, some other researchers have worked on the creep behavior of engineering structures such as tunnels and roadways [30][31][32].
Although, there are some published articles related to the macroscopic failure behavior or viscoelastic plastic properties of rock induced by microscopic damage evolution during time-independent deformation [33][34][35][36][37][38][39][40][41], few studies closely combined the microscopic damage evolution of rock with macroscopic TDD and viscoelastic plastic mechanical properties of creep behavior. Most of the published models are isotropic, and the homogeneity model cannot effectively describe the heterogeneous characteristics of sandstone in the mesoscale. In this work, an instantaneous elastic-plastic damage constitutive model and a viscoelastic plastic time-dependent damage constitutive model of sandstone are established, based on laboratory test data. These models are validated, respectively, by an analytical solution of an elastic-plastic model and by the Nishihara model. Supplementing this work, a series of conventional constant-strain-rate tests and multi-step creep tests were carried out on sandstone, which is very common in coal mine roofs. Then, the simulation results are carefully analyzed in combination with results from the laboratory. This research builds a damage constitutive model that should be able to describe the full stages of creep behavior for sandstone. The damage model can also express crack propagation or damage evolution before creep failure.

Numerical Model
According to the instantaneous deformation of rock or rock-like material from constantstrain-rate tests [13] and the strain separation of rock creep behavior [14], an intact constitutive model of rock should contain two important parts: an instantaneous (time-independent) elastic-plastic part, and a time-dependent viscoelastic plastic part. The instantaneous elastic-plastic model could only be used for simulated constant-strain-rate tests, while the instantaneous elastic-plastic model and time-dependent viscoelastic plastic model combined could perform the whole behavior of rock during creep tests.

Instantaneous Elastic-Plastic Model
In order to describe crack evolution before failure, a damage model is adopted in this paper. Some elastic-plastic equations are added to the damage model to illustrate the instantaneous behavior of rock.

Heterogeneous Damage Model
Sandstone samples are homogeneous in the macroscale, but also composed of different kinds of mineral grains in the mesoscale (as shown in Figure 1A). All samples were cored, as the ISRM suggested method, from the same block of material to a diameter of 25 mm, and a length of 100 mm, resulting in a length to diameter ratio of 2:1. Therefore, the sandstone specimens are heterogeneous in the mesoscale. Accordingly, it is assumed that the numerical model of rock is composed of mesoscopic heterogeneous units so as to simulate the mesoscopic heterogeneity of the rock specimens [38,42,43]. As shown in Figure 1B, in order to illustrate the heterogeneity of rock elements in the mesoscale, the mechanical properties (such as strength and Young's modulus) are assigned randomly by Weibull statistical distribution [44].
where u is a mechanical property of the rock in the mesoscale (strength, elastic modulus, etc.); u 0 is the average value of the mechanical properties for all elements; x is the uniformity coefficient, reflecting the homogeneity of the whole specimen; and f (u) is the statistical distribution density function of rock elements for the rock's mechanical property u. More information on rock heterogeneity models can be found in these articles [37,38].

Instantaneous Elastic-Plastic Model
In order to describe crack evolution before failure, a damage model is adopted in this paper. Some elastic-plastic equations are added to the damage model to illustrate the instantaneous behavior of rock.

Heterogeneous Damage Model
Sandstone samples are homogeneous in the macroscale, but also composed of different kinds of mineral grains in the mesoscale (as shown in Figure 1A). All samples were cored, as the ISRM suggested method, from the same block of material to a diameter of 25 mm, and a length of 100 mm, resulting in a length to diameter ratio of 2:1. Therefore, the sandstone specimens are heterogeneous in the mesoscale. Accordingly, it is assumed that the numerical model of rock is composed of mesoscopic heterogeneous units so as to simulate the mesoscopic heterogeneity of the rock specimens [38,42,43]. As shown in Figure  1B, in order to illustrate the heterogeneity of rock elements in the mesoscale, the mechanical properties (such as strength and Young's modulus) are assigned randomly by Weibull statistical distribution [44].
where u is a mechanical property of the rock in the mesoscale (strength, elastic modulus, etc.); 0 u is the average value of the mechanical properties for all elements; x is the uniformity coefficient, reflecting the homogeneity of the whole specimen; and   f u is the statistical distribution density function of rock elements for the rock's mechanical property u. More information on rock heterogeneity models can be found in these articles [37,38]. The maximum tensile stress criterion and the strength criterion determined by the plastic yield function [13] are adopted as damage judgment criteria. The maximum tensile stress criterion is: where 1  is the minimum principal stress and 0 t f is the uniaxial tensile strength. Under any stress condition, the tensile damage is judged first by the maximum tensile stress criterion. When the strength does not meet the tensile criterion, it must be decided whether the strength meet the strength criterion. The strength criterion used in the model is based on the plastic yield function, which is described in detail in the following section. The maximum tensile stress criterion and the strength criterion determined by the plastic yield function [13] are adopted as damage judgment criteria. The maximum tensile stress criterion is: where σ 1 is the minimum principal stress and f t0 is the uniaxial tensile strength. Under any stress condition, the tensile damage is judged first by the maximum tensile stress criterion. When the strength does not meet the tensile criterion, it must be decided whether the strength meet the strength criterion. The strength criterion used in the model is based on the plastic yield function, which is described in detail in the following section. When the stress state or strain state of the rock element satisfies the given damage threshold, the element begins to be damaged. The elastic damage of the elastic modulus could be expressed by the following: where E 0 and E are, respectively, the elastic moduli before and after the damage and D is the damage variable. However, the damage of the element and its evolution are all assumed to be isotropic in the model, so E 0 , E and D are all scalars. The damage variable D of the mechanical parameters under compressive and tensile stress states satisfies the piecewise function damage constitutive relation, which can be written: where ε 1 and ε 3 are, respectively, the tensile and compressive principal strains; ε t0 and ε c0 are the corresponding maximum tensile and compressive principal strains when the element experiences tensile or shear damage; n is the element damage evolution coefficient, taken to have the value of 2; σ tr and σ cr are the residual stresses at the stress condition of tension and compression, respectively; E 0 is the initial elastic modulus of the element; F 2 is the strength criterion of the element, which is determined by the plastic yield function F 2 = f f (σ ij , σ f ); and σ ij is the principal stress. At the loading stage, as the plastic strain increases, the element can also start to be damaged when the strength criterion is satisfied. The stress-strain relationship can be described by the elastic-plastic constitutive equation. The matrix expression of the equation is usually written in the following form: where C ep is called the elastic-plastic flexibility matrix, C ep = [C e ] + C p ; [C e ] is the flexibility matrix of elastic strain, while C p is the flexibility matrix of plastic strain. After the strength reaches the limit in the compressive stress state, and considering the strain softening (due to the damage after the strength meets the criterion), the whole stiffness matrix of the elastic-plastic strain can be written:
The instantaneous deformation part can be separated further into two parts: instantaneous elastic deformation and instantaneous plastic deformation: According to the generalized Hooke's law, the three-dimensional constitutive relationship of the elastic body's instantaneous deformation can be written: where s ij is the stress deviatoric tensor; σ m is the spherical stress tensor; δ ij is the Kronecker delta (zero when i = j and 1 when i = j); [C e ] is the elastic flexibility matrix of the element; G 0 and K 0 are, respectively, the shear modulus and the bulk modulus; and µ is Poisson's ratio.
According to the flow law of the traditional elastic-plastic theory, the plastic strain increment can be found [16,46] as follows: where Q p is the plastic potential function of rock; dλ p is the plastic strain increment. The plastic strain increment index can be derived from the plastic consistency condition.
where f p is the plastic subsequent yield function and α p is the plastic hardening function. According to Equations (8)-(10), the expression of total elastic-plastic strain can be found to be: Zhou et al. [13], Shao et al. [12], Jia et al. [47] and Pietruszczak et al. [48] considered that the following yield function can be used to represent the plastic yield surface of the elasticplastic model, which can effectively characterize the instantaneous ductile deformation characteristics of rock: where p, q and θ are the mean stress (the compression force is taken as negative in the model), stress deviator and Lode angle, respectively; α p γ p is the plastic hardening function; R c is the normalization coefficient, which is taken to be the strength of the element in this model; C s is the cohesive coefficient of the rock (with a value of 0.01); and A is a constant related to the internal friction of the rock. The plastic hardening function in Equation (13) can be described by a function related to the effective plastic shear strain [47]: where B is a coefficient controlling the plastic hardening rate and γ p is the effective plastic shear strain, which can be illustrated by the following equation [49]: In order to describe the transition of volumetric strain from compression to dilatancy, a non-associated plastic flow rule [12,13,50,51] is needed. Additionally, a plastic potential function not equal to the yield function should be proposed to achieve the non-associated plastic flow rule. Based on experimental data, the plastic potential function can be found to be: where β p is a coefficient controlling the turning point of volume deformation from compaction to dilatancy; it is taken to be 1.1. When α p γ p < β p , the volume deformation is compaction; when α p γ p > β p , the deformation is dilatancy. Substituting Equations (13) and (18) into Equation (11), the plastic strain increment can be found to be: In general, the elastic-plastic constitutive damage model for the element can be described by Equation (12). The elastic-plastic constitutive model of the compression and tension loading for the element is shown in Figure 2.
Some numerical results are compared with the result from analytical solution to validate the availability of the proposed elastic-plastic model.
(1) Analytical solution model In order to verify the reliability of the elastic-plastic numerical model esta The loading and unloading conditions of the elastic-plastic model are expressed as follows:  In order to verify the reliability of the elastic-plastic numerical model established above, a circular hole was excavated (numerically) in the stress field. The stress field distribution around the circular hole obtained by the proposed numerical model was compared with that obtained by an analytical solution. The comparison was used to achieve the validation of the numerical model.
As shown in Figure 3, the size of the region to be investigated is 200 cm × 200 cm, with a borehole radius of 5 cm, and with the number of units set to 40,000. The horizontal stress is σ h , and the vertical stress is σ v . The dashed line A-A in Figure 3 extends from the right side of the central circular hole to the right boundary of the stress field. In addition, the tangential (vertical) and radial (horizontal) stress data of the elements on the dashed line are extracted from the numerically calculated result and from the analytical solution for further comparison. The analytical solution of the elastic-plastic model with the M-C criter hesion into account, can be described by: 1 sin 1 sin 2 cos 1 sin where  is the internal friction angle of the rock (assigned the value of 2 radial stress,   is the circumferential stress and c is the cohesion force, to be 0.12 MPa.
The analytical solution assumes that the stress distribution around th can be divided into a circumferential plastic zone and a circumferential ela side the plastic zone. The formula for stress distribution in the elastic zone  The analytical solution of the elastic-plastic model with the M-C criterion, taking cohesion into account, can be described by: where φ is the internal friction angle of the rock (assigned the value of 29 • ), σ r is the radial stress, σ θ is the circumferential stress and c is the cohesion force, which is taken to be 0.12 MPa.
The analytical solution assumes that the stress distribution around the circular hole can be divided into a circumferential plastic zone and a circumferential elastic zone outside the plastic zone. The formula for stress distribution in the elastic zone is as follows: where σ e is the radial stress of the elastic-plastic boundary, σ e = (1 − sin φ)σ 0 and σ 0 is the far-field stress.
The stress distribution in the plastic zone can be expressed in the following form: where p is the inner wall pressure of the borehole, with a value of 2 MPa, and r 0 is the radius of the circular hole, taken to be 5 cm.
According to the stress continuity condition of the elastic-plastic boundary, the radius of the boundary between the elastic and plastic regions can be written: where σ 0 is the far-field stress, with a value of 10 MPa.
(2) Numerical result validation Conventional triaxial compressive strength tests were conducted on the sandstone to obtain the stress-strain curves for confining pressures of 0, 20, and 40 MPa. The stress-strain curves are shown in Figure 4A. Based on these curves and the definitions of effective deviatoric stress and mean effective stress, the yield stress surface and the failure stress surface for the sandstone are illustrated in Figure 4B. The plastic yield function of Equation (13) was used to fit the data in Figure 4B, which reveals the undetermined parameters in the plastic yield function. The determined parameters in the plastic yield function for the sandstone are listed in Table 1. They can determine the initial yield function and the corresponding strength criterion for the sandstone. After some calculations by repetition, the physical and mechanical parameters in the numerical model are also obtained as shown in Table 1.  Figure 4A. Based on these curves and the definitions of effective deviatoric stress and mean effective stress, the yield stress surface and the failure stress surface for the sandstone are illustrated in Figure 4B. The plastic yield function of Equation (13) was used to fit the data in Figure 4B, which reveals the undetermined parameters in the plastic yield function. The determined parameters in the plastic yield function for the sandstone are listed in Table 1. They can determine the initial yield function and the corresponding strength criterion for the sandstone. After some calculations by repetition, the physical and mechanical parameters in the numerical model are also obtained as shown in Table 1.  The proposed elastic-plastic damage model, which adopts the boundary condition of a stress field with a circular borehole, can be used to get the stress distribution on the dashed line A-A' in Figure 3. The stress distribution around the circular borehole obtained from the proposed damage model and the analytical solution are shown in Figure 5. The stresses for comparison contain radial stress and circumferential stress. The results obtained from the proposed damage model show two different degrees of homogeneity. When the homogeneity of the rock sample is very high, the strengths of elements on the dashed line tend toward the same value. However, the lower the homogeneity of the rock, the more fluctuant the stress on elements around the result of the analytical solution. The  The proposed elastic-plastic damage model, which adopts the boundary condition of a stress field with a circular borehole, can be used to get the stress distribution on the dashed line A-A in Figure 3. The stress distribution around the circular borehole obtained from the proposed damage model and the analytical solution are shown in Figure 5. The stresses for comparison contain radial stress and circumferential stress. The results obtained from the proposed damage model show two different degrees of homogeneity. When the homogeneity of the rock sample is very high, the strengths of elements on the dashed line tend toward the same value. However, the lower the homogeneity of the rock, the more fluctuant the stress on elements around the result of the analytical solution. The fluctuation of the stress distribution is due to the Weibull distribution of the element strength in the model, which can induce a difference of stress distribution on the elements. solution. Additionally, the plastic deformation is extracted from the whole deformation occurring around the borehole and drawn in Figure 6.
When the far-field stress is 10 MPa and the inner pressure in the borehole is 2 MPa, the results based on the damage model show that the radius of the plastic zone centered on the borehole is about 9 cm. At the same time, the radius of the plastic zone obtained from the analytical solution is 7.91 cm, very closed to that from the damage model. To sum up, the elastic damage model and elastic-plastic damage model with heterogeneity established in the present research are in good agreement with the stress field distribution obtained from the analytical solutions, which supports the reliability of the numerical damage model.  As shown in Figure 5A, the comparative results of radial stress show that when the homogeneity is 200, the radial stress obtained from the elastic-plastic damage model is in good agreement with that from the analytical solution. When the homogeneity is 9, the radial stress obtained from the damage model fluctuates around that of the analytical solution. This indicates that the proposed damage model is promising for getting the distribution of radial stress. As shown in Figure 5B, the comparative results show that the circumferential stress extracted from the result of the proposed elastic-plastic damage model differs little from that of the analytical solution when the homogeneity is very high. The difference is the location of the maximum yield stress around the borehole, a difference that comes about because the yield function used in the proposed model is different from that in the analytical solution. Yet the various features of the circumferential stress curves obtained from the two models are almost the same, showing that the proposed damage model can account for the elastic behavior and the plastic behavior as well as the analytical solution. Additionally, the plastic deformation is extracted from the whole deformation occurring around the borehole and drawn in Figure 6.

Time-Dependent Viscoelastic Plastic Model
The proposed viscoelastic plastic model is adapted from the traditional Nishihara model. The proposed model and the validation with an analytical solution are detailed as the follows. When the far-field stress is 10 MPa and the inner pressure in the borehole is 2 MPa, the results based on the damage model show that the radius of the plastic zone centered on the borehole is about 9 cm. At the same time, the radius of the plastic zone obtained from the analytical solution is 7.91 cm, very closed to that from the damage model. To sum up, the elastic damage model and elastic-plastic damage model with heterogeneity established in the present research are in good agreement with the stress field distribution obtained from the analytical solutions, which supports the reliability of the numerical damage model.

Time-Dependent Viscoelastic Plastic Model
The proposed viscoelastic plastic model is adapted from the traditional Nishihara model. The proposed model and the validation with an analytical solution are detailed as the follows.

Nishihara Model
The traditional Nishihara model is composed of three groups of elements, such as elastic elements, viscoelastic elements and visco-plastic elements (as shown in Figure 7). The one-dimensional creep equation of the Nishihara model is as follows: where σ s is the yield stress at a one-dimensional scale. When σ ≥ σ s , the total strain can be regarded as ε ij = ε e + ε ve + ε vp . Furthermore, according to the three-dimensional form of the generalized Hooke's law, the elastic deformation can be expressed as follows: where G 0 and K 0 are shear modulus and bulk modulus, respectively.

Nishihara Model
The traditional Nishihara model is composed of three gro elastic elements, viscoelastic elements and visco-plastic elemen The one-dimensional creep equation of the Nishihara model is a   Viscoelastic and visco-plastic deformations can be described by the equations: where η 1 and η 2 are, respectively, the shear viscosity coefficient and the plastic viscosity coefficient. During transformation of the viscoelastic formula from one dimension to three dimensions, it is assumed that the rock does not exhibit creep behavior under a spherical stress tensor or hydrostatic pressure. It is further assumed that only stress deviation can cause the rock to exhibit creep behavior [13]. The viscoelastic strain in the TDD can be illustrated as follows: where G 1 is the viscoelastic shear modulus, expressed by G 1 = E 1 2(1+µ) , in which µ is Poisson's ratio.
According to the principles of plastic mechanics, the expression of visco-plastic strain in TDD can be expressed as follows [13,52]: where F is the plastic yield function, F 0 is the initial value of the plastic yield function (which is taken to be 1) and Q is the plastic potential function related to the plastic yield surface. An expression of the Heaviside step function <·> in Equation (30) is shown: in which φ(•) is assumed to have the form of a power function. Usually the exponential index n in the power function is equal to 1 [52]. The power function can be written φ F F 0 = F F 0 . The total deformation equation for the three-dimensional form of the Nishihara model can be written as follows: By substituting the above formulas into Equation (32), the total deformation equation of the traditional Nishihara model in three-dimensional form can be found to be as follows:

Improved Nishihara Model
As stated in the Nishihara model, the TDD can be divided into two partitions: viscoelastic deformation and visco-plastic deformation. The three-dimensional form of TDD can be expressed as A three-dimensional expression for the viscoelastic deformation is provided by Equation (29), and an expression of the visco-plastic deformation is provided by Equation (30). Therefore, the yield function and plastic potential function used in Section 2.1.2 can also be used here to establish the non-associated flow rule under visco-plastic strain in the TDD: The effective shear strain in Equation (16) can be described as follows: Based on Equations (29), (30), (35) and (36), an expression of the TDD can be derived: After synthesizing instantaneous deformation and TDD, the total viscoelastic plastic deformation equation can be found from Equation (7) (as shown in Figure 8). The threedimensional expression of the total deformation equation can then be rewritten: where ε c ij (D, σ c ij ) is expressed as the TDD, which depends on the damage factor and the loading stress.

Analytical Solution Verification
(1) Analytical solution of Nishihara model based on D-P criterio According to the circular tunnel of a perfect elastoplastic solution Prager (D-P) yield criterion [53], the mechanical model shown in Fig  The stress in the stress field around the circular hole with a radius distance from a point in the stress field to the center of the circular ho the work of Hou and Niu [53], the plastic stress distribution around derived using the D-P criterion, is as follows: where  and k are test constants related only to the internal frictio cohesion factor C , and whose values can be taken to be: Figure 8. Schematic diagram of the improved Nishihara model.

Analytical Solution Verification
(1) Analytical solution of Nishihara model based on D-P criterion According to the circular tunnel of a perfect elastoplastic solution with the Drucker-Prager (D-P) yield criterion [53], the mechanical model shown in Figure 9 is established. The stress in the stress field around the circular hole with a radius of R 0 is P 0 , and the distance from a point in the stress field to the center of the circular hole is r. According to the work of Hou and Niu [53], the plastic stress distribution around the circular hole, as derived using the D-P criterion, is as follows: where α and k are test constants related only to the internal friction angle ϕ and the cohesion factor C, and whose values can be taken to be: The radius of the plastic zone can be derived from the following equation. (42) According to published data [53,54], some basic mechanical paramet merical model used for model validation can be determined and are show Based on the analytical solution, the stress distribution and the TDD aroun hole can be obtained.  (2) Verification results Using the analytical solution of the Nishihara model based on the D-P the proposed damage model, the relationship between the displacemen around the borehole and the distance R0 to the boundary of the borehole i Figure 10A. Similarly, the stress distribution around the circular hole plott distance to the boundary of the borehole is shown in Figure 10B.
The comparison of displacements shows that when the homogeneity o mechanical properties is up to 200, there is a little difference between the obtained by the damage model and that by the analytical solution. The sma caused by the difference in yield functions used in the model and the analy In fact, the stress distribution law extracted from the result of the damage mo agreement with that from the analytical solution. Especially when the homo the stress distribution for the damage model is almost the same as that for The stress distribution in the elastic zone can be described as follows: Up to this point, the stress distribution around the circular hole based on the D-P criterion can be achieved [53]. The analytical solution of the TDD around the circular hole based on the Nishihara model was deduced in another paper [54]. The deformation u and the TDD u 0 at r = R 0 are, respectively, the following: According to published data [53,54], some basic mechanical parameters in the numerical model used for model validation can be determined and are shown in Table 2. Based on the analytical solution, the stress distribution and the TDD around the circular hole can be obtained. (2) Verification results Using the analytical solution of the Nishihara model based on the D-P criterion and the proposed damage model, the relationship between the displacement of elements around the borehole and the distance R 0 to the boundary of the borehole is displayed in Figure 10A. Similarly, the stress distribution around the circular hole plotted against the distance to the boundary of the borehole is shown in Figure 10B.

Comparisons with Laboratory Tests
In order to illustrate advantages and disadvantages of the proposed damage model, we have conducted some uniaxial compressive tests on the sandstone. Based on a comparison of stress-strain curves and the failure mode from the numerical model and from the laboratory tests, it will be helpful for readers to understand the proposed model for its possible use in the future.

Consistent Strain Rate Tests
Based on the parameters in Table 1 and the viscoelastic plastic damage model, some uniaxial compressive strength (UCS) tests (or constant-strain-rate tests) were carried out. During the laboratorial tests, specimens were deformed uniaxially at a constant strain rate of 10 −5 s −1 until macroscopic failure. The stress-strain curves obtained from the numerical simulation and the laboratory tests are compared in Figure 11. From the comparison of the stress-strain curves, it can be found that the numerical damage model can exhibit the full stress-strain process of the sandstone, and can also show compression and dilation behavior of volumetric deformation. This feature is basically consistent with that illustrated by the Vosges sandstone in ref. [13]. Additionally, the compression and dilation behavior of the volumetric strain in the numerical simulation can verify the validity of the non-associated flow rule adopted in the damage model. The comparison of displacements shows that when the homogeneity of the element mechanical properties is up to 200, there is a little difference between the displacement obtained by the damage model and that by the analytical solution. The small difference is caused by the difference in yield functions used in the model and the analytical solution. In fact, the stress distribution law extracted from the result of the damage model is in good agreement with that from the analytical solution. Especially when the homogeneity is 200, the stress distribution for the damage model is almost the same as that for the analytical solution. Therefore, the results obtained from the proposed damage model match well with those from the analytical solution.

Comparisons with Laboratory Tests
In order to illustrate advantages and disadvantages of the proposed damage model, we have conducted some uniaxial compressive tests on the sandstone. Based on a comparison of stress-strain curves and the failure mode from the numerical model and from the laboratory tests, it will be helpful for readers to understand the proposed model for its possible use in the future.

Consistent Strain Rate Tests
Based on the parameters in Table 1 and the viscoelastic plastic damage model, some uniaxial compressive strength (UCS) tests (or constant-strain-rate tests) were carried out. During the laboratorial tests, specimens were deformed uniaxially at a constant strain rate of 10 −5 s −1 until macroscopic failure. The stress-strain curves obtained from the numerical simulation and the laboratory tests are compared in Figure 11. From the comparison of the stress-strain curves, it can be found that the numerical damage model can exhibit the full stress-strain process of the sandstone, and can also show compression and dilation behavior of volumetric deformation. This feature is basically consistent with that illustrated by the Vosges sandstone in ref. [13]. Additionally, the compression and dilation behavior of the volumetric strain in the numerical simulation can verify the validity of the non-associated flow rule adopted in the damage model. the stress-strain curves, it can be found that the numerical damage model can exh full stress-strain process of the sandstone, and can also show compression and behavior of volumetric deformation. This feature is basically consistent with th trated by the Vosges sandstone in ref. [13]. Additionally, the compression and behavior of the volumetric strain in the numerical simulation can verify the validi non-associated flow rule adopted in the damage model. Some snapshots of the element strength during failure of the numerical model are extracted and shown in Figure 12. All of these snapshots have the same representative characteristic. As shown in the figure, the distribution of cracks is relatively discrete when the cracks begin to appear in the numerical specimen ( Figure 12B). Then, the cracks begin to aggregate when they are very closed to each other ( Figure 12C). At the same time, outlines of the final fracture have begun to form. However, as shown in Figure 12D, some small cracks also appear around the main fracture surface, after the main fracture surface has come into being. Figure 13 shows the failure picture of sandstone in the laboratory tests after the UCS tests, demonstrating that the failure mode obtained from the damage model is very similar to that from the laboratory tests. All of the failure modes from the damage model and the laboratory tests are oblique shear failure, and their inclined angles are almost the same. The damage model can also confirm the phenomenon that some small cracks appear around the main fracture. All of this evidence demonstrates that the proposed damage model has a very good ability to reproduce the instantaneous failure of sandstone, not only as demonstrated in the stress-strain curves, but also as illustrated by the failure mode and fracture evolution. Some snapshots of the element strength during failure of the numerical mo extracted and shown in Figure 12. All of these snapshots have the same represe characteristic. As shown in the figure, the distribution of cracks is relatively discret the cracks begin to appear in the numerical specimen ( Figure 12B). Then, the crack to aggregate when they are very closed to each other ( Figure 12C). At the same tim lines of the final fracture have begun to form. However, as shown in Figure 12D small cracks also appear around the main fracture surface, after the main fracture has come into being. Figure 13 shows the failure picture of sandstone in the laboratory tests after t tests, demonstrating that the failure mode obtained from the damage model is very to that from the laboratory tests. All of the failure modes from the damage model laboratory tests are oblique shear failure, and their inclined angles are almost th The damage model can also confirm the phenomenon that some small cracks around the main fracture. All of this evidence demonstrates that the proposed d model has a very good ability to reproduce the instantaneous failure of sandsto only as demonstrated in the stress-strain curves, but also as illustrated by the failur and fracture evolution.

Multi-Step Creep Tests
Considering the parameters in Tables 1 and 2, some ducted on sandstone using the proposed damage model. creep test were set at 60%, 65%, 70%, 75%, 80% and 85% o stone. The loading stress at every stage was held for 12 h

Multi-Step Creep Tests
Considering the parameters in Tables 1 and 2, some multi-step creep tests were conducted on sandstone using the proposed damage model. The loading stress levels of the creep test were set at 60%, 65%, 70%, 75%, 80% and 85% of the peak strength for the sandstone. The loading stress at every stage was held for 12 h. Then the test was switched to the next loading stage until the sandstone specimen failed. Figure 14 compares the multi-step creep curves obtained from the damage model and from the laboratory tests. The results show that the creep curves obtained from the numerical model, especially for the axial strain, fit well with those from the laboratory tests, which also supports the feasibility of using the model used in the creep test. The failure evolution drawn from the results of the numerical multi-step creep tests is shown in Figure 15. The snapshots of the failure evolution show that numerous cracks appear in the specimen in the multi-step creep tests. In addition, more cracks accumulate and develop next to the small fracture, more than in the numerical UCS test. Especially for the final failure mode, the number of main fractures occurring during the multi-step creep test was more than that during the UCS test, which can be certified by the result of creep tests in the lab, as shown in Figure 16. A small fracture appears around the main fracture during the experimental creep tests. Numerous fragments of sandstone come out at the fracture surface during these tests. By contrast, the fracture surface during the UCS tests does not have so much rock fragmentation, which means that the damage on the fracture surface is less pronounced for this way of loading. Both the fractures in the failure mode during the numerical simulation and the fragments on the fracture surface during the experimental tests can be interpreted to mean that the TDD has an accelerating effect on the generation and expansion of cracks inside the specimen, and also has some influence on the final fracture form of the specimen.
Based on the creep strain in Figure 14, creep strain rates at different loading stress levels for the experimental test and the numerical test can be obtained and are shown in Figure 17. The comparison shows that the creep strain rates obtained from the numerical simulation results have a consistent trend that matches those from the laboratory tests. The creep strain rates gradually increase with an increase of loading stress. the creep strain rate has a strong linear relationship with the corresponding loading stress in a semi-logarithmic graph.
fragments of sandstone come out at the fracture surface during these tests. By contrast, the fracture surface during the UCS tests does not have so much rock fragmentation, which means that the damage on the fracture surface is less pronounced for this way of loading. Both the fractures in the failure mode during the numerical simulation and the fragments on the fracture surface during the experimental tests can be interpreted to mean that the TDD has an accelerating effect on the generation and expansion of cracks inside the specimen, and also has some influence on the final fracture form of the specimen.   the fracture surface during the UCS tests does not have so much rock fragmentation, which means that the damage on the fracture surface is less pronounced for this way of loading. Both the fractures in the failure mode during the numerical simulation and the fragments on the fracture surface during the experimental tests can be interpreted to mean that the TDD has an accelerating effect on the generation and expansion of cracks inside the specimen, and also has some influence on the final fracture form of the specimen.   Based on the creep strain in Figure 14, creep strain rates at different loading stress levels for the experimental test and the numerical test can be obtained and are shown in Figure 17. The comparison shows that the creep strain rates obtained from the numerical simulation results have a consistent trend that matches those from the laboratory tests. The creep strain rates gradually increase with an increase of loading stress. the creep strain rate has a strong linear relationship with the corresponding loading stress in a semi-logarithmic graph. Based on the creep strain in Figure 14, creep strain rates at different loading stress levels for the experimental test and the numerical test can be obtained and are shown in Figure 17. The comparison shows that the creep strain rates obtained from the numerical simulation results have a consistent trend that matches those from the laboratory tests. The creep strain rates gradually increase with an increase of loading stress. the creep strain rate has a strong linear relationship with the corresponding loading stress in a semi-logarithmic graph. Figure 17. Relationships between creep strain rate and applied stress in the multi-step creep tests in the laboratory and in numerical simulation.

Discussion
All of above comparisons show that the proposed damage model can reveal the basic deformation characteristics of the time-dependent creep behavior and can also be used to investigate damage or evolution of cracks and the final fracture patterns, not only for the instantaneous deformation but also the TDD. However, the deformation process for rock is very complicated due to it being anisotropic in the mesoscale, and due to the variability among different specimens. Even though we have adopted the Weibull distribution of strength parameters in the numerical model to describe the heterogeneity of the sandstone Figure 17. Relationships between creep strain rate and applied stress in the multi-step creep tests in the laboratory and in numerical simulation.

Discussion
All of above comparisons show that the proposed damage model can reveal the basic deformation characteristics of the time-dependent creep behavior and can also be used to investigate damage or evolution of cracks and the final fracture patterns, not only for the instantaneous deformation but also the TDD. However, the deformation process for rock is very complicated due to it being anisotropic in the mesoscale, and due to the variability among different specimens. Even though we have adopted the Weibull distribution of strength parameters in the numerical model to describe the heterogeneity of the sandstone in the mesoscale, it is still very difficult to match the macro behavior of the sandstone sample. The proposed damage model in the present paper can reveal the crack evolution and the failure mode, which are very similar to the results of experimental tests. However, we could not get the numerical and laboratory results for the strain-time curves to match at 95% or more. The development of the TDD constitutive model remains a significant problem. Some researchers have proposed models to simulate the rock creep behavior and have found good fitting curves [14,16,19]. However, the good fitting focuses only on the strain-time curves. If we want to have good matching results for both the strain curves and the failure model, more research work will need to be done.
We have adopted a complicated multinomial in the damage model, and there are some parameters to determine before the model can be used. Considering the differences among rock samples and the varieties of rock, it is very difficult to predict accurately what the creep behavior of another rock sample will be. Additionally, the calculated quantities using the damage model will increase because of the long multinomial, which will reduce the computation speed. The damage model limits the amount of calculations and the built model will not be very complicated. The damage model, if used to solve engineering problems in the future, will be challenging because of its high calculation requirement, but could very well be used to better understand the failure of rock samples. Before it can be used in practical engineering problems, the damage model will need some simplification.

Conclusions
In order to explore the TDD characteristics of sandstone, especially for the evolution of cracks before rock failure, a viscoelastic plastic damage model was established in this work. The elastoplastic behavior and the viscoelastic plastic behavior of the damage model were validated by analytical solutions. After the validation of the damage model, uniaxial compressive strength tests and multi-step creep tests were conducted on the sandstone to investigate the merits and the deficiencies of the model. The research results reveal the following: (1) The instantaneous elastic-plastic damage model can display the process of volumetric strain for sandstone from compression to dilation. Many small fractures are formed within a very short period of time in the specimen before the main fracture emerges. The damage model captures the phenomenon that small damage of the specimen can lead to the generation of micro-cracks, the accumulation of which causes the specimen to fracture. The final failure mode provided by the numerical simulation is consistent with that extracted from laboratory test results. (2) The viscoelastic plastic damage model can reproduce the viscoelastic plastic deformation behavior. However, unlike the result of the conventional constant-strain-rate test, cracks of the rock appear through all of the stress loading stages in the multi-step creep tests. The final macroscopic fracture is the result of gradual accumulation of the cracks, a conclusion that is similar to that of the UCS tests. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.