This paper aims mainly at providing an incremental elastoplastic constitutive model for heterogeneous porous rock-like materials in the frame of micromechanics. The studied material is considered to be made up of randomly distributed spherical pores embedded in a pressure-sensitive solid matrix obeying Drucker–Prager yield function. The effective elastic properties of porous rocks are obtained by the use of Mori and Tanaka homogenization scheme, which are on function of the bulk and shear moduli of the solid matrix and of the value of porosity. For the macroscopic nonlinear phase, a limit analysis-based macroscopic criterion is adopted to derive the basic constitutive rule by considering an associated plastic flow rule. In order to capture the typical hardening effects of rocks, an originally proposed hardening function of the solid matrix is also taken into consideration, which is related on the accumulated equivalent plastic strain. In order to verify its accuracy, the proposed micro-macro constitutive model is implemented by a numerical procedure including elastic predictions and plastic corrections and compared to experimental results of triaxial compression tests of sandstone with different confining pressures. It is observed that the numerical simulation is in accord with the experimental data, indicating that the obtained model is able to predict the main mechanical behaviours of rock-like materials.

1. Introduction

In recent years, rock-like heterogeneous materials have been widely encountered in engineerings, such as nuclear waste storage, petroleum industry, and mining constructions, which contain microstructures (cracks, pores, and inclusions) at different scales [1, 2]. Experiments have shown that inners pores or voids exhibit a great influence on the effective strength and mechanical behaviours on such materials [35], which lead to the complex plastic deformation, brittle-ductile transition, and sensitivity to confining pressure. The safety of rock engineering [69] is the most important subject to the researchers. In order to predict the nonlinear behaviours of porous rock-like materials, a series of constitutive models have been constructed by researchers and engineers. Based on the experimental data, a series of elastoplastic models [1012] are proposed to capture the main features of such materials. However, the effects of the voids on the strength and the failure related to its plastic deformation and voids evolution cannot be explicitly reflected.

To this end, a limit analysis-based kinematical approach of porous media with a von Mises type solid matrix containing spherical voids was firstly introduced by Gurson [13] in the framework of micromechanics, in which the influence of porosity at microscale on the macroscopic strength was directly taken into account. Following Gurson’s step, numerous extensions have been developed accounting for void shape effects [14, 15], size effects [16, 17], materials with pressure-sensitive matrix [1820], and nonassociated constitutive laws [21], which can be applied to rock-like materials to predict the effective yield strength. On the contrary, a theoretical dual stress-based model aiming to provide a lower bound criterion has also been attempted firstly in [22, 23] and recently extended by [24, 25] based on the very new stress variational homogenization (SVM) model by solving the saddle point problem. For the cyclic loading cases, homogenized fatigue criteria of porous media [2628] have also been provided by considering the shakedown limit state based on the homogenization theory.

By using the macroscopic strength criteria and Hill’s incremental approach [29], a micro-macro constitutive model was proposed by Shen et al. [30]. Clayey rocks were considered as three-phase materials (matrix, pores, and inclusions). The evolution of microstructures (pores and inclusions) was explicitly taken into account in the final expression to describe the macroscopic responses of the considered material. Similarly, the cases of different yield functions of the solid phase and of double pores at two scales [31, 32] have also been developed. It is remarked that the choice of the harding rule is very essential to the model in the hardening or softening phase. A nonassociated plastic flow rule is usually adopted to obtain an accurate lateral stain prediction [3335], which will introduce an additional parameter. Consequently, we will provide a modified micromechanics-based constitutive model for porous rock-like materials by considering a new hardening law, such that only an associated plastic flow rule is concerned.

This paper is organized in the following order. In Section 2, a limit analysis-based macroscopic criterion for porous materials with Drucker–Prager type matrix proposed by Maghous et al. [19] is briefly recalled. The influences of frictional parameter of the solid phase and of porosity on the effective strength are also studied in this section. By adopting an associated plastic flow rule, a micro-macro elastoplastic constitutive model is obtained in the rate form in Section 3, which accounts for the evolution of porosity and frictional parameter. A new plastic hardening function depending on the equivalent plastic strain is introduced to improve the accuracy of this model. In Section 4, the proposed model is implemented numerically and compared to available triaxial compression test results of sandstone under different confining pressures. Both the macroscopic axial and lateral strains predicted by the present model are validated. Some concluding remarks are provided in Section 5.

2. Homogenized Macroscopic Strength Criterion of Porous Rocks

In the perspective of studying the mechanical characteristics of rock-like porous media in the framework of micromechanics, the material is considered made up of random voids at microscale embedded by isotropic solid matrix, see Figure 1. The volumes of whole material and of the void are noted as and , respectively. As a result, the volume of the solid phase is obtained as , providing the void volume fraction .

In order to describe the pressure-sensitive property of rock-like materials, the solid matrix is considered obeying the Drucker–Prager yield function. The solid matrix is supposed to be made of elastoplastic material obeying Drucker–Prager yield criterion:where donates the equivalent stress with being the deviatoric part and is the mean stress. The parameters and represent the friction coefficient and the yield stress of the solid matrix under hydrostatic tension.

The local associated plastic flow rule can be obtained from the above yield function:where is the microscopic strain rate tensor and denotes the rate form plastic parameter.

Concerning the effective strength of porous media with the consideration of effects of voids, Maghous et al. [19] have established a closed-form macroscopic criterion by the use of a modified secant method. In the present paper, this homogenized criterion is also adopted to characterize effective plastic behaviour of rock-like materials:where (resp. ) is the macroscopic mean stress defined as its counterpart (resp. ) at a microscale,with being the deviatoric part of macroscopic stress tensor.

It is observed that effective strength of the materials is related to the macroscopic stress states, the value of porosity, the friction coefficient, and the tensive strength of the solid matrix. In the following figures, the chosen criterion is plotted with the variation of different parameters. represents the compression zone. The comparison of the chosen criterion (3) is illustrated in Figure 2 with the variation of porosity. The strength of the homogenized material is obviously reduced with the increase of void volume fraction in the compression zone, while it shrinks slightly with the variation of porosity in the other part.

Considering a constant value of , Figure 3 plots the comparison of the chosen criterion with different values of friction coefficients of the solid matrix. In the compression zone, the loading capacity of material is improved with the increase of the friction parameter , while the strength is clearly reduced in the tension part. Moreover, for the pure deviatoric loading condition , the macroscopic strength of the porous media stays constant, indicating that the maximum shear load bearing capacity depends on the value of .

Concluding from Figures 2 and 3, the macroscopic criterion (3) proposed by Maghous et al. can reflect the tension-compression asymmetry of rock-like materials. It can also present the effects of the value of porosity and of friction coefficient on the strength. For the special case that tends to zero, the criterion will become symmetric about axis , which leads to the one porous media with von Mises type matrix.

3. Micromechanics-Based Elastoplastic Constitutive Model of Porous Rocks

3.1. Macroscopic Elastic Properties of the Material

The solid phase of porous rocks is considered as elastoplastic material, and the total local strain rates can be composed of elastic and plastic parts:

As a result, the macroscopic stress-strain relation of the material in the rate form writeswhere denotes the macroscopic plastic stain rate tensor and is the initial elastic stiffness fourth-order tensor of the porous rocks, depending on the elastic parameters of the solid phase:where and are effective bulk modulus and shear modulus of rocks, respectively, is the fourth-order hydrostatic projector, and represents the fourth-order deviatoric projector.

Following the homogenization studies of rock-like materials [31, 36], the effective elastic properties are described by adopting the Mori and Tanaka [37] scheme. The effective elastic parameters are related to those of the solid matrix and also to the value of porosity. Considering randomly distributed spherical voids embedded by the isotropic solid phase, one has the effective bulk and shear modulus of porous rocks:where and are the bulk and shear moduli of the solid matrix, respectively. For completeness, the general expressions of homogenized elastic parameters of composite materials with different phases for isotropic case at both local and overall scales are given in Appendix A.

3.2. Micromechanics-Based Constitutive Formulation of Porous Rocks

Taking into account the variables in adopted macroscopic strength criterion (3), the macroscopic plastic consistent condition readsand the macroscopic associated plastic flow rulewhere denotes the macroscopic nonnegative plastic multiplier in rate form.

Considering Hill’s lemma of macrohomogeneity [38], one has the following relation of the plastic dissipations at micro- and macroscales by considering a representative volume element with the homogeneous strain rate boundary condition:

Using the local plastic flow rule (2), the local plastic dissipation reads

Combining with the yield condition (1), the following relation can be built:

By the use of relation between the macroscopic plastic dissipation and the microscopic one (11), it is derived that

The evolution law of voids is obtained by kinematical compatibility condition, relating to the plastic macroscopic homogenized volumetric strain rate and the local one at microscale:where the volumetric plastic strain rate .

The microscopic equivalent plastic strain rate and volumetric one can be, respectively, computed from the local yield function:where , is the deviatoric part of . Taking the previous equation into equation (14), one has

Likewise, the macroscopic plastic volumetric strain rate can also be obtained from the macroscopic criterion:

Combining equation (15) with equation (18), the evolution of porosity is determined by

Finally, considering the general form of the macroscopic stain-stress relation,and substituting equations (17), (10), (19), and (20) into the consistent condition (9), the expression of the macroscopic plastic multiplier is computed in the following form:in which the macroscopic plastic hardening modulus is

Introducing the macroscopic plastic flow rule into equation (20), the macroscopic incremental constitutive law of porous rocks readswhere is the effective fourth-order tangent stiffness tensor:

It is worthy noticing that this constitutive model is also partly used in the study of Callovo-Oxfordian argillite by considering homogenization schema at three scales [31].

3.3. Consideration of Hardening Effects with a New Hardening Rule

On the contrary, inspired from the previous studies [30, 33, 39] concerning the plastic hardening or softening phase of rock-like materials beyond the elastic limit, the above constitutive law is not accurate enough to capture the mechanical behaviour of rocks. By comparison to experimental data, the available hardening rule of the solid matrix is always related to the equivalent plastic strain :

Developed from the typical hardening rule for rock-like materials governed by the minimum and maximum values of the frictional parameter used in [30], a three-parameter hardening function depending on the equivalent plastic strain is adopted in the present study:where denotes the initial value of the frictional parameter corresponding to the first plastic yield state. Besides, is a dimensionless parameter related to the current local equivalent plastic strain and to the critical value at peak stress . This proposed hardening function is partly inspired from the research on the plastic damage model [40] of rocks, in which controls the damage evolution rate at the microscale. In the present study, this specific form is also adopted to estimate the evolution of the frictional parameter, which is related to actual equivalent plastic strain. Figure 4 plots the curves of on function of with respect to different values of . It is observed that, for before the peak stress , it does not exhibit obvious differences with the variation of . On the contrary, in the postpeak phase , the value of has an important influence on the evolution rate. It is reasonable to assume that the hardening curve tends to more smooth with a smaller value of .

Consequently, taking into account of the hardening effects of thye solid matrix and the other variables in the effective strength criterion, the plastic consistency condition rewrites

Similar to the derivation in the previous section, the corresponding macroscopic plastic multiplier is computed by taking the macroscopic plastic flow rule (10) under the new plastic consistency condition:where the new macroscopic hardening modulus is given by

As a result, the final effective fourth-order tangent stiffness tensor in the elastoplastic constitutive law (23) is modified as

Equation (23) together with (29) and (30) contribute to the main results of the present paper. The derived micro-macro elastoplastic constitutive model for porous rock-like materials depends on the elastic parameters of the solid phase and the porosity, as well as the local hardening constants. The associated plastic flow is used without considering an additional potential function. Comparing with the plastic hardening law adopted in [30, 31], the one proposed in this work depends is also governed by a dimensionless value related to the equivalent plastic strain at peak stress. The capacity to capture the effective mechanical behaviours, especially in the hardening process, is discussed in the following section.

4. Assessment of the Micromechanics-Based Modified Model

In this section, the established micro-macro model is implemented by a numerical iteration procedure containing two phases at each increment: elastic prediction and plastic correction. The obtained numerical results are compared with the experimental data of triaxial compression tests of sandstone with respect to different confining pressures. Notice that the experimental results used in this paper are firstly provided by Liu et al. [41] in the study of brittle-ductile transition phenomenon of sandstone, which is one of the most abundant sedimentary rocks in nature.

According to [41], the porosity of the chosen sandstone . Besides, the bulk and shear moduli and the hardening parameters of the solid matrix are given in Table 1, which are identified by fitting the experimental results in triaxial compression test with confining pressure (Figure 5). For the elastic parameters, one has and . It must be pointed out that the initial threshold frictional parameter related to the first yield point is very sensitive to the loading state and can be computed by solving the equation of macroscopic criterion (3) with the corresponding stresses, which are obtained from experimental data.

Adopting the same parameters in Table 1, the numerical predictions by the proposed constitutive model of triaxial compression tests with different confining pressures , 50, and 70 MPa are also, respectively, plotted in Figures 68 and compared to the corresponding experimental results. The evolutions of deviatoric stresses with both axial (blue) and lateral (red) deformations are produced. Noticing that most parameters presented in this model can be easily obtained from the experimental tests, except for hardening variables.

In general, the numerical simulations are in good agreement with the experimental results of the triaxial compression tests of sandstone, indicating that the present constitutive law is able to predict its main mechanical behaviours. The effects of the confining pressures on the strength of the porous rocks can be clearly reflected by the proposed model without consideration of the variation of the yield stress in classical phenomenological models. In particular, although an associated plastic flow rule is applied in the present model, the lateral deformations can also be accurately computed numerically in all triaxial compression tests, which are usually overestimated by other constitutive models based on nonassociated plastic flow rules.

On the contrary, there are significant differences between numerical and experimental results for axial strain in Figure 6 and in Figure 7. The softening phase after the peak stress in experiments of the sandstone cannot be correctly described by the proposed model. According to Zhao et al. [40], the hardening curve tends to be more smooth with a smaller value of in (26). The errors in the present study are mainly due to the choice of a certain value of in the hardening function, which only simulates the hardening effects.

5. Conclusion

In this paper, a new incremental constitutive model for rock-like materials with the pressure-sensitive solid matrix has been constructed in the framework of micromechanics to capture its elastoplastic behaviours. Derived from a macroscopic strength criterion of porous media, the obtained homogenized constitutive law has been developed by considering a modified hardening function, which is related to the equivalent local plastic strain and a plastic state dimensionless parameter. In this proposed model, the elastic parameters (bulk and shear moduli) of the solid phase and the porosity are taken into account by the use of homogenization technique. The evolution of the microstructure is quantified by the variation of void volume friction and that of local plastic equivalent strain.

In order to verify its prediction ability of porous rocks’ mechanical behaviours, the numerical simulations of this model has been plotted and compared with experimental data of triaxial compression tests on sandstone. A good agreement between the numerical and experimental results has been observed, indicating that the proposed model is able to describe the main features of sandstone.

Finally, to extend the proposed constitutive model by considering the degradation effects of Young’s modulus of the solid matrix is potentially an interesting topic in the outlook.


A. Homogenized Elastic Parameters of Porous Media by Mori and Tanaka scheme

In Section 3, the effective elastic tensor is given in function of bulk and shear moduli of solid phase and the porosity, which is deduced from Mori and Tanaka scheme [37] by considering ellipsoidal inclusions of the same form and orientation in an infinite matrix.

Concerning a representative elementary volume with uniform boundary condition,where is the location vector and the macroscopic strain tensor :and the fourth-order strain localization tensor is introduced as

For a composite material made up of different phases (the volume fraction of phase is noted as ), the initial macroscopic elastic tensor is expressed as

Considering the isotropic case in both local and overall scales and Eshelby’s solution, the tensor to describe the interactions of different phases writeswhere is the elastic stiffness tensor of solid matrix, is the fourth-order unit tensor, and denotes Hill’s tensor which takes the following form:

Consequently, the effective elastic stiffness tensor can be written as

Combining with the general form of the stiffness tensor (7), the above equation can be reduced to produce the effective elastic parameters at the macroscale. In this study, the solid matrix is considered as Phase 0 and pores as Phase 1 with zero-value bulk and shear moduli. Substituting into equation (A.7), we can finally obtain the macroscopic elastic parameters in (8) for porous rock-like materials. For more detailed derivation of the effective elastic parameters, readers are referred to [30].

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they do not have any financial or nonfinancial conflict of interests.


The authors would like to thank the Huzhou South Taihu Innovation Team for the funding support.