Based on the theory of porous media mechanics, the non-linear mechanical response of sandstone is described by using the plastic yield equation considering porosity. At the same time, the damage evolution criterion was established based on the existing damage theory, and independent damage variables were introduced into the plastic yield function to describe the damage evolution of rock in the process of deformation and failure, so as to determine the coupling relationship between plastic deformation and damage development. This model can simulate both pre-peak strengthening and post-peak softening behaviors at the same time. The physical model of porous media is shown in Figure 1.

As is shown in Figure 1, the porous media is considered made up of isotropic solid matrix and random pore at microscale. The volumes of solid matrix and of the void are noted as Ω m and Ω 1 , respectively. The volume of the whole porous media is obtained as Ω = Ω m + Ω 1 . Based on the above statement, unit porosity can be calculated by the following formula: f = Ω 1 Ω = Ω 1 Ω m + Ω 1 < 1 , (1)

Compared with metal materials, the pressure sensitivity and volumetric deformation are two crucial characteristics of rock materials. In order to consider these aspects, here we assumed that the solid matrix is made up of elastoplastic material subjected to Drucker-Prager yield criterion. F ( σ ) = σ d + T ( σ m − h ) ≤ 0 , (2)

Based on the assumption of small deformation, in the field of traditional plastic mechanics, the total strain increment of rock can be decomposed into elastic strain (increment) and plastic strain (increment): ε i j = ε i j e + ε i j p ,   D i j = D i j e + D i j p , (3)

Referring to the previous research theory of porous media [6–8], the effective volume modulus and shear modulus of intact and non-destructive rock materials are expressed as follows: κ 0 h o m = 4 ( 1 − f ) κ s μ s 4 μ s + 3 f κ s , μ 0 h o m = ( 1 − f ) μ s 1 + 6 f κ s + 2 μ s 9 κ s + 8 μ s , (4) where the parameters κ s , μ s , and f represent the elastic bulk modulus, shear modulus of solid matrix, and porosity, respectively, in which κ s and μ s can be derived by the elastic modulus E s and Poisson v s of solid phase.

Concerning the assumption of material isotropy and damage extension isotropy, the scalar ω is used to represent the damage variable. Therefore, the effective bulk modulus and shear modulus of damaged rock material can be expressed as: κ ( ω ) = κ 0 h o m ( 1 − ω ) , μ ( ω ) = μ 0 h o m ( 1 − ω ) , (5)

As a result, the macroscopic elastic stress-strain relation of the rock in the incremental form writes [24]: d Σ = C w : ( D − D P ) + ∂ C w ∂ ω : ( ε − ε p ) d ω , (6) where C w denotes damage stiffness matrix, D and D P represent total strain increment and plastic strain increment, respectively. ε and ε p are the total strain and plastic strain. The total elastic strain tensor is given by ε e = ε − ε p .

In general, the M-C yield criterion and D-P yield criterion are used for the yield function. Based on the above experimental results of triaxial compression, the yield surfaces of sandstone under seepage conditions are close to non-linear characteristic of elliptical surface, thus the typical porous plastic criterion proposed by Ref. [25] is used to describe the mechanical response of sandstone in the study. The plastic yield surface equation is defined as: Φ ( Σ , f , T , ω ) = 1 + 2 f / 3 T 2 Σ d 2 + ( 3 f 2 T 2 − 1 ) Σ m 2 + 2 h ( 1 − f ) Σ m − ( 1 − f ) 2 h 2 , (7) where Σ d , Σ m is the generalized shear stress and mean stress of rock on macro scale respectively. Σ d = Σ ′ : Σ ′ Σ m = Σ 11 + Σ 22 + Σ 33 3 , (8)

T is the plastic hardening function, which reflects the pre-peak strengthening and post-peak softening behavior of sandstone [26]. Based on the thermodynamic framework and the work of [27], the plastic thermodynamic potential of sandstone Γ p can be obtained as: Γ p = ( 1 − ω ) [ T m γ p − ( T m − T 0 ) b 1 ⁡ ln b 1 + γ p b 1 ] , (9)

Thus the expression of hardening function T is derived as: T = T ( γ p , ω ) = ∂ Γ p ∂ γ p = ( 1 − ω ) [ T 0 + ( T m − T 0 ) γ p b 1 + γ p ] , (10) where T 0 and T m are two parameters related to the position of the initial plastic yield surface and the final plastic yield surface for the rock. b ₁ represents the controlling parameter of hardening rate for the rock. The influencing factors of the hardening function include plastic shear strain γ p of solid matrix and damage variables ω , indicating that hardening function T increases with increasing plastic shear strain γ p , but decreases with the increasing damage variable ω , indicating the characteristics of increased plastic deformation and post-peak softening for the rock.

In addition, the non-associated plastic potential function G is used to describe the characteristics of rock from compression to expansion, which is defined as follows [25]: G = G ( Σ , f , t ) = 1 + 2 f / 3 T t Σ d 2 + ( 3 f 2 T t − 1 ) Σ m 2 + 2 h ( 1 − f ) Σ m , (11) where t is plastic hardening function related to damage variable and plastic shear strain of rock and can be given by: t = t ( γ p , ω ) = ( 1 − ω ) ( t 0 + ( t m − t 0 ) γ p b 2 + γ p ) , (12) where t 0 and t m are two parameters related to the position of the initial plastic potential function and the final plastic potential function for the rock. b ₂ represents the controlling parameter of hardening rate for the rock. In addition, the plastic strain rate of sandstone D p is calculated based on the non-associated flow rule as follows: D p = d λ p   ∂ G ∂ Σ , (13) where λ ˙ is the plastic multiplier, and it is used to verify the following loading-unloading condition: { d λ p = 0     i f     Φ < 0   o r   i f   Φ = 0   a n d   Φ ˙ < 0 d λ p ≥ 0   i f   Φ = 0   a n d   Φ ˙ = 0 , (14)

Assuming that the change of pore volume on the microscopic scale only depends on pore plastic compression and expansion, and ignore the influence of new pore nucleation, according to function Eq. 1, we can get that: d f = d ( Ω 1 Ω ) = d Ω 1 Ω − Ω 1 Ω d Ω Ω = ( 1 − f ) ( d Ω Ω − d Ω m Ω m ) , (15) Where d Ω Ω is the mean macroscopic volumetric strain rate t r D p , d Ω m Ω m is the volume strain rate of solid matrix t r e ˙ . Based on the assumption that solid matrix obeying Drucker-Prager yield function and non-associated flow rule, the potential is given by ϕ ( σ ) = σ d + t σ m , thus the mesoscopic strain rate can be written as follows: e ˙ = d Λ ∂ ϕ ∂ σ ,   e ˙ ′ = d Λ σ ′ σ d ,   e ˙ m = 1 3 t d Λ , (16) where e ˙ ′ is the deviatoric strain rate tensor with e ˙ = e ˙ ′ + e ˙ m δ . d Λ is the plastic multiplier of the solid matrix. Therefore, the equivalent plastic strain rate γ ˙ p can be calculated as: γ ˙ p = e ˙ ′ : e ˙ ′ = d Λ , (17)

According to the energy-based equivalence condition provided by Ref. [28], the following relation between plastic strain rate of porous medium material and equivalent plastic strain rate of solid matrix can be derived as presented in [29]: γ ˙ p = Σ : D p ( 1 − f ) ( T h + ( t − T ) Σ m 1 − f ) , (18) With the relations (Eqs 16, 17) and tre=t γ p in hand, the variation of porosity in Eq. 15 can be expressed as: d f = ( 1 − f ) ( t r D p − t γ ˙ p ) , (19)

In accordance with the thermodynamic theory, the effect of the damage driving force related to the free energy release rate in the elastic stage on the damage of rock change its internal structure. The internal cracks of the sandstone specimen are mostly closed based on conventional triaxial compression. Therefore, the damage of the sandstone is primarily caused by plastic shear, while the damage caused by the elastic deformation is very small. According to the previous study [21] and ignoring the effect of elasticity, the damage driving force Y d can be obtained as follows: Y d   = − ∂ Γ p ∂ ω = T m γ p − ( T m − T 0 ) b 1 ⁡ ln b 1 + γ p b 1 , (20)

In addition, the damage evolution criterion of the rock is introduced by the Mazars’ research [30], and the damage evolution function f d can be defined as: f d = ω c t h ( B d Y d ) − ω ≤ 0 , (21) where ω c is the maximum threshold of damage variable and B _d is a parameter related to the rate of damage evolution.

The plastic flow and damage evolution of sandstone are coupled processes under loading conditions [27]. To reflect the effect of pore water pressure on mechanical behavior of rock, the plasticity multiplier λ p and damage multiplier λ d can be determined by coupling plastic flow and damage evolution, which can obtain the consistency conditions of plastic deformation and damage variable for rock material as follows: { Φ ˙ ( Σ , f , T , ω ) = ∂ Φ ∂ Σ : d Σ + ∂ Φ ∂ f d f + ∂ Φ ∂ T ∂ T ∂ γ p d γ p + ∂ Φ ∂ ω d ω = 0 f d ˙ ( Y d , ω ) = ∂ f d ∂ Y d d Y d + ∂ f d ∂ ω d ω = 0 , (22)

In addition, based on the plastic flow rule, the increment of plastic deformation and damage variable are defined as: { d ε p = d λ p ∂ G ∂ Σ d ω = d λ d ∂ f d ∂ Y d , (23)

According to Eqs. 6, 7, 10–13, 18–21, the plastic multiplier and the damage evolution multiplier are obtained: { d λ d = ( T m + ( T m − T 0 ) b 1 / b 1 + γ p ) Σ : ∂ G ∂ Σ / ( 1 − f ) ( T h + ( t − T ) Σ m 1 − f ) d λ p d λ p = ∂ Φ ∂ Σ : C w : D / ∂ Φ ∂ Σ : C w : ∂ G ∂ Σ − ∂ Φ ∂ f ( 1 − f ) ( ∂ G ∂ Σ m − t B ) − ∂ Φ ∂ T ∂ T ∂ γ p B − ( ∂ C w ∂ ω : ε e − ∂ Φ ∂ T ∂ T ∂ w ) A B , (24)

The parameter A and B is written as follows: A = T m + ( T m − T 0 ) b 1 b 1 + γ p ⋅ B = Σ : ∂ G / ∂ Σ ( 1 − f ) ( T h + ( t − T ) Σ m 1 − f ) , (25)