Elastic Modulus and Strength of Rock-Like Material with Locked-In Stress

Abstract

Locked-in stress is a special stress objectively existing in the rock which has an important influence on the mechanical properties. However, there were few studies on the locked-in stress in rocks. The locked-in stress was simplified into spherical inclusions of stress. Based on the Eshelby inclusion theory, the influence equation of the locked-in stress on the elastic modulus of the rock is derived. Through experiment, the expression of elastic modulus of rock-like material with locked-in stress based on M-T method is determined. Based on Mohr-Coulomb criterion and effective stress principle, the equations for the influence of locked-in stress on strength of rock-like material under uniaxial and conventional triaxial conditions were deduced. Through experimental verification, an empirical formula for the strength of rock-like material containing locked-in stress is obtained. This article provided a method for further studying of the locked-in stress in rock-like material and has certain reference value.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Rocks with randomly distributed inclusions can be considered as natural composite material. The main component of rock can be regarded as the matrix of composite material; cracks, heterogeneous particles, and pores in the rock can be regarded as inclusions in composite material. The elastic properties of composite materials are important issues and many theoretical results have been achieved so far. For engineering, rock is generally in elastic deformation stage before yielding. In this article, the micromechanics theory was used to study the elastic properties of rock containing inclusions. So far, many achievements have been made in the study of locked-in stresses in metals, for example, Voyiadjis’s study on the determination of locked-in stresses in metals (Voyiadjis [1, 2]), but there is little research on locked-in stresses in rocks. Therefore, this article will further study the situation of inclusion with locked-in stresses.

The strength criterion of rock materials is a macroscopic concept. However, microinclusions can significantly affect the strength of the rock material, which can be derived from experiments. In fact, the influence of the mesostructure on failure criterion not only exists in the study of locked-in stress, but also exists in other respects. For example, Chen [3] studied a new equivalent statistical damage constitutive model on rock block mixed up with fluid inclusions. Chen [4] studied the influence of porosity on the strength of sandstone based on the Mohr-Coulomb strength criterion. Wang [5] studied the strength and constitutive relationship of rock under the joint action of pore water pressure and confining pressure based on statistical theory. Lai [6] and Li [7] studied the elastic properties and constitutive relations of frozen soils at different temperatures. Zhou [8] studied the mechanical properties of rocks under the combined action of temperature, pore water, and confining pressure. Ren [9] researched three-dimensional mesodamage model of rock concrete based on the deformation and expansion of elliptic microcrack. Yang [10] studied the joint action mechanism of cyclic loading and pore water pressure on rocks. This paper intends to establish the criteria for the destruction of rock materials affected by inclusion content and locked-in stress. The Mohr-Coulomb criterion (Jaeger [11]) is a commonly used failure criterion for rock materials. It can describe the failure of rock under uniaxial compression and conventional triaxial compression. The cohesion force of rock-like test blocks used in this paper is provided by cement; the shear friction is provided by granular fine sand. The failure characteristics of samples are in accordance with Mohr-Coulomb strength theory.

In this article, the inclusion mechanics of composite materials and the Mohr-Coulomb strength criterion were used to analyze the mechanism of the variation of elastic modulus and peak strength of rock specimens with inclusions content and value of locked-in stress.

2. Mechanical Model of Rock with Locked-in Stress

Locked-in stress in rocks actually exists in various forms and is inhomogeneous, which is determined by the nature of geological evolution. Carrier materials of locked-in stress (inclusions) are of various shapes and sizes. The material outside the locked-in stress carrier (matrix) also has random diversity. Randomness is the difficulty in quantifying mechanical parameters of rock, as well as in the study of locked-in stress. It is only possible to obtain a simple analytical solution by simplifying the problem into ideal inclusion with limited content in the elastic matrix. With reference to the practice in the mesomechanics of composite materials, one cell element in the rock is analyzed. Let the volume of the cell be 1; the inclusion volume is the volume content $\rho$ of the inclusion in the test block; ${\sigma }^L$ is the locked-in stress when no external load is applied (Figure 1). The problem of rock contain locked-in stress is reduced to the problem that the ellipsoidal heterogeneous inclusions contain initial intrinsic strain in the homogeneous matrix, and, from this, we get the relevant conclusions about the coupling of external loads and locked-in stress in ellipsoidal heterogeneous inclusions.

[figure omitted; refer to PDF]

${\sigma }_{\mathrm{11}}^{\mathrm{0}}$ , ${\sigma }_{\mathrm{22}}^{\mathrm{0}}$, and ${\sigma }_{\mathrm{33}}^{\mathrm{0}}$ are the three principal stresses of the external load; ${\sigma }^L$ is the locked-in stress caused by inconsistent deformation of the inclusion and the matrix (Figure 1). There are several reasons for inconsistent deformation, such as nonrecoverable plasticity, slippage, nonuniform thermal expansion, weathering, and martensite deformation. Figure 2 can simply state the process of the generation of locked-in stress. There was a heterogeneous sphere in the isotropic material; initially, the sphere and the external medium were in close contact with each other, and the locked-in stress is zero. The heterogeneous sphere was taken out alone and expanded by $\Delta r$ for some reason (thermal expansion, etc.). Only one cavity remains in the matrix material, from which the heterogeneous spheres have been removed. Due to the different material properties, the change in the matrix under the same conditions was negligible and so did the radius of the spherical cavity in the base. At this time, a normal force was applied to the outside of the heterogeneous sphere to restore it to its original shape and reimposed it into the cavity. Then the normal applied stress was removed and the heterogeneous sphere expands, causing compression between the sphere and the matrix. Due to the constraints of the external matrix, the amount of springback of the heterogeneous inclusion cannot reach $\Delta r$, but there is a difference $\Delta r^{\prime }$ from $\Delta r$. It can be seen, the cavity radius of the base material is expanded by $\Delta r-\Delta r^{\prime }$; and the compression of the heterogeneous sphere relative to free expansion is $\Delta r^{\prime }$. Therefore, the value of the locked-in stress can be obtained according to the size change and the mechanical parameters of each material (Figure 2).

[figure omitted; refer to PDF]

The problem studied in this article is that the rock contains locked-in stress and is subject to external loads; i.e., the initial state is the self-balanced state in Figure 2. An external load is applied on the basis of Figure 2, as shown in Figure 1.

3. Changes of Locked-in Stress with External Loads

3.1. Tensor Solution Based on Inclusion Theory

The locked-in stress in the inclusion is not fixed. When the rock is subjected to an external load, the deformation of the matrix and the inclusion must be linked to each other. When the external load is tensile, it will reduce the constraints of the matrix on the inclusion. If the adhesion between the inclusion body and the base body is weak, the interface will be separated and the medium will not be continuous. In this study, the applied loads are all compressive stresses. In the elastic stage, the compression of the matrix material will transfer the external load to the inclusion, so that the locked-in stress in the inclusion will increase, and react on the matrix.

Based on inclusion theory (Eshelby [12]) and M-T method (Mori [13]), combined with the simplified mechanical model shown in Figure 1, the locked-in stress in the inclusions after external load was applied on the rock was analyzed.

The problem in Figure 1 is a heterogeneous inclusion with intrinsic strain. The cells are affected by the uniform stress ${\sigma }_{ij}^0$ in the far field. When there is no heterogeneous inclusion in the cell, the resulting strain is ${\epsilon }_{ij}^0$. The relationship between stress and strain is$$\begin{array}{c}\text{(1)}& {\sigma }_{ij}^{\mathrm{0}}=C_{ijkl}^{\mathrm{0}}{\epsilon }_{ij}^{\mathrm{0}}{\epsilon }_{ij}^{\mathrm{0}}={C_{ijkl}^{\mathrm{0}}}^{-\mathrm{1}}{\sigma }_{ij}^{\mathrm{0}}.\end{array}$$

$C_{ijkl}^0$ is the elastic constant tensor of the matrix. Due to the presence of inclusions, the actual cell strain is in fact different from ${\epsilon }_{ij}^0$ and will produce a perturbation strain $\tilde{\epsilon }$. The average stress of the matrix part in this cell is$$\begin{array}{}\text{(2)}& {\sigma }_{ij}^{\left(\mathrm{0}\right)}={\sigma }_{ij}^{\mathrm{0}}+\tilde{\sigma }=L_{\mathrm{0}}\left({\epsilon }_{ij}^{\mathrm{0}}+\tilde{\epsilon }\right).\end{array}$$

Obviously, the perturbation part of the stress in the matrix is$$\begin{array}{}\text{(3)}& \tilde{\sigma }=L_0\tilde{\epsilon }.\end{array}$$

Due to the different elastic properties of the material, the average stress and strain in the inclusions are different from the corresponding values in the matrix, and their differences are ${\sigma }^{\prime }$ and ${\epsilon }^{\prime }$, respectively. Therefore the stress in the inclusion is$$\begin{array}{}\text{(4)}& {\sigma }_{ij}^L=C_{ijkl}^{\mathrm{1}}\left({\epsilon }_{ij}^{\mathrm{0}}+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }\right)=C_{ijkl}^{\mathrm{0}}\left({\epsilon }_{ij}^{\mathrm{0}}+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast }\right).\end{array}$$

If there is a locked-in stress in the inclusion, i.e., the intrinsic strain ${\epsilon }_{ij}^{\ast 1}$, then there are$$\begin{array}{}\text{(5)}& {\sigma }_{ij}^L=C_{ijkl}^1\left({\epsilon }_{ij}^0+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast 1}\right)=C_{ijkl}^0\left({\epsilon }_{ij}^0+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast }\right).\end{array}$$

Subtracting the right sides of (5) and (4) yields$$\begin{array}{}\text{(6)}& C_{ijkl}^{\mathrm{1}}\left({\epsilon }_{ij}^{\mathrm{0}}+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast \mathrm{1}}\right)-C_{ijkl}^{\mathrm{1}}\left({\epsilon }_{ij}^{\mathrm{0}}+{\tilde{\epsilon }}_{ij}+{\epsilon }_{ij}^{\prime }\right)=-C_{ijkl}^{\mathrm{1}}{\epsilon }_{ij}^{\ast \mathrm{1}}.\end{array}$$

It can be seen that, for a given inclusion body’s intrinsic strain (or locked-in stress), the amount of change in external stress to the inclusion body’s stress is irrelevant to the original locked-in stress. Therefore, only the situation of (4) needs to be analyzed.

$C_{ijkl}^1$ is the elastic constant tensor of the inclusion and ${\epsilon }_{ij}^{\ast }$ is the equivalent intrinsic strain of the inclusion. According to the Eshelby inclusion theory there are$$\begin{array}{}\text{(7)}& {\epsilon }_{kl}^{\prime }=S_{ijmn}{\epsilon }_{mn}^{\ast }.\end{array}$$

From (7), (4), and (2) we can find$$\begin{array}{}\text{(8)}& {\sigma }^{\prime }=C_{ijkl}^0\left({\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast }\right)=C_{ijkl}^0\left(S_{ijkl}-I_{ijkl}\right){\epsilon }_{ij}^{\ast }.\end{array}$$

The volume-average stress of the cell should be equal to the average stress of the far-field effect and is$$\begin{array}{}\text{(9)}& {\sigma }_{ij}^0=\left(1-\rho \right){\sigma }_{ij}^{\left(0\right)}+\rho {\sigma }_{ij}^L.\end{array}$$

According to (2) and (4) we get$$\begin{array}{}\text{(10)}& \tilde{\sigma }=-\rho {\sigma }^{\prime },\text{(11)}& \tilde{\epsilon }=-\rho \left({\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^{\ast }\right)=-\rho \left(S_{ijkl}-I_{ijkl}\right){\epsilon }_{ij}^{\ast }.\end{array}$$

Substituting (7) and (11) into (4), the solution is$$\begin{array}{}\text{(12)}& {\epsilon }_{ij}^{\ast }={\left\{C_{ijkl}^{\mathrm{0}}+\left(C_{ijkl}^{\mathrm{1}}-C_{ijkl}^{\mathrm{0}}\right)\left[\rho I_{ijkl}+\left(\mathrm{1}-\rho \right)S_{ijkl}\right]\right\}}^{-\mathrm{1}}\left(C_{ijkl}^{\mathrm{0}}-C_{ijkl}^{\mathrm{1}}\right){\epsilon }^{\mathrm{0}}.\end{array}$$

Substituting (12), (11), (1), and (7) into (4) can solve the relationship between the stress in the inclusion and the far field stress of the cell.$$\begin{array}{}\text{(13)}& {\sigma }_{ij}^L=C_{ijkl}^{\mathrm{1}}\left\{I_{ijkl}+\left[S_{ijkl}-\rho \left(S_{ijkl}-I_{ijkl}\right)A\right]\right\}{C_{ijkl}^{\mathrm{0}}}^{-\mathrm{1}}{\sigma }_{ij}^{\mathrm{0}},\end{array}$$with$$\begin{array}{}\text{(14)}& A={\left\{C_{ijkl}^{\mathrm{0}}+\left(C_{ijkl}^{\mathrm{1}}-C_{ijkl}^{\mathrm{0}}\right)\left[\rho I_{ijkl}+\left(\mathrm{1}-\rho \right)S_{ijkl}\right]\right\}}^{-\mathrm{1}}\left(C_{ijkl}^{\mathrm{0}}-C_{ijkl}^{\mathrm{1}}\right).\end{array}$$

Equation (13) is the general solution of the locked-in stress in the inclusion with the change of the external load of the test block.

3.2. The Solution When the Inclusion Is Spherical

Isotropic fourth-order tensor in composite micromechanics$$\begin{array}{}\text{(15)}& H_{ijkl}=\frac{\mathrm{1}}{\mathrm{3}}\alpha {\delta }_{ij}{\delta }_{kl}+\frac{\mathrm{1}}{\mathrm{2}}\beta \left({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl}-{\frac{\mathrm{2}}{\mathrm{3}}}_{ij}{\delta }_{kl}\right),\text{(16)}& {H^{\prime }}_{ijkl}=\frac{\mathrm{1}}{\mathrm{3}}{\alpha }^{\prime }{\delta }_{ij}{\delta }_{kl}+\frac{\mathrm{1}}{\mathrm{2}}{\beta }^{\prime }\left({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl}-{\frac{\mathrm{2}}{\mathrm{3}}}_{ij}{\delta }_{kl}\right),\end{array}$$is abbreviated as

$H=\left(\alpha ,\beta \right)$ , $H^{\prime }=\left({\alpha }^{\prime },{\beta }^{\prime }\right)$, and in accordance with the algorithm, $H+H^{\prime }=\left(\alpha +{\alpha }^{\prime },\beta +{\beta }^{\prime }\right)$, $H^{-\mathrm{1}}=\left(\mathrm{1}/\alpha ,\mathrm{1}/\beta \right)$, $H{H}^{\prime }=\left(\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }\right)$, fourth-order unit tensor can be expressed as $I_{ijkl}=\left(\mathrm{1,1}\right)$ and the elastic constant tensor can be expressed as $C_{ijkl}=\left(\mathrm{3}K,\mathrm{2}G\right)$.

The equation of the tensor $A$ can be obtained by the following calculation process:$$\begin{array}{}\text{(17)}& A={\left\{C_{ijkl}^{\mathrm{0}}+\left(C_{ijkl}^{\mathrm{1}}-C_{ijkl}^{\mathrm{0}}\right)\left[\rho I_{ijkl}+\left(\mathrm{1}-\rho \right)S_{ijkl}\right]\right\}}^{-\mathrm{1}}\left(\mathrm{3}{K}_{\mathrm{0}}-\mathrm{3}{K}_{\mathrm{1}},\mathrm{2}{G}_{\mathrm{0}}-\mathrm{2}{G}_{\mathrm{1}}\right)=\left(\frac{\mathrm{1}}{K_{\mathrm{0}}/\left(K_{\mathrm{0}}-K_{\mathrm{1}}\right)-\left(\rho +\mathrm{9}{K}_p{K}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)},\frac{\mathrm{1}}{G_{\mathrm{0}}/\left(G_{\mathrm{0}}-G_{\mathrm{1}}\right)-\left(\rho +\mathrm{4}{G}_p{G}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)}\right),\text{(18)}& I_{ijkl}+\left[S_{ijkl}-\rho \left(S_{ijkl}-I_{ijkl}\right)A\right]=\left(\mathrm{9}{K}_p{K}_{\mathrm{0}}+\mathrm{1}+\frac{\mathrm{9}{K}_p{K}_{\mathrm{0}}-\mathrm{1}}{K_{\mathrm{0}}/\left(K_{\mathrm{0}}-K_{\mathrm{1}}\right)-\left(\rho +\mathrm{9}{K}_p{K}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)},\mathrm{4}{G}_p{G}_{\mathrm{0}}+\mathrm{1}+\frac{\mathrm{4}{G}_p{G}_{\mathrm{0}}-\mathrm{1}}{G_{\mathrm{0}}/\left(G_{\mathrm{0}}-G_{\mathrm{1}}\right)-\left(\rho +\mathrm{4}{G}_p{G}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)}\right),\text{(19)}& {\sigma }^L=\left(\mathrm{9}{K}_p{K}_{\mathrm{1}}+\frac{K_{\mathrm{1}}}{K_{\mathrm{0}}}+\frac{\mathrm{9}{K}_p{K}_{\mathrm{1}}-K_{\mathrm{1}}/K_{\mathrm{0}}}{K_{\mathrm{0}}/\left(K_{\mathrm{0}}-K_{\mathrm{1}}\right)-\left(\rho +\mathrm{9}{K}_p{K}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)},\mathrm{4}{G}_p{G}_{\mathrm{1}}+\frac{G_{\mathrm{1}}}{G_{\mathrm{0}}}+\frac{\mathrm{4}{G}_p{G}_{\mathrm{1}}-G_{\mathrm{1}}/G_{\mathrm{0}}}{G_{\mathrm{0}}/\left(G_{\mathrm{0}}-G_{\mathrm{1}}\right)-\left(\rho +\mathrm{4}{G}_p{G}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)}\right){\sigma }^{\mathrm{0}},\end{array}$$with$$\begin{array}{c}\text{(20)}& K_p=\frac{\mathrm{1}}{\mathrm{3}\left(\mathrm{4}{G}_{\mathrm{0}}+\mathrm{3}{K}_{\mathrm{0}}\right)},G_p=\frac{\mathrm{3}\left(\mathrm{2}{G}_{\mathrm{0}}+K_{\mathrm{0}}\right)}{\mathrm{10}{G}_{\mathrm{0}}\left(\mathrm{4}{G}_{\mathrm{0}}+\mathrm{3}{K}_{\mathrm{0}}\right)}\end{array}$$ (Shen [14]),$$\begin{array}{c}\text{(21)}& K_{\mathrm{0}}=\frac{E_{\mathrm{0}}}{\mathrm{3}\left(\mathrm{1}-\mathrm{2}{\mu }_{\mathrm{0}}\right)},G_{\mathrm{0}}=\frac{E_{\mathrm{0}}}{\mathrm{2}\left(\mathrm{1}+{\mu }_{\mathrm{0}}\right)}\end{array}$$is the bulk modulus and shear modulus of the matrix, and$$\begin{array}{c}\text{(22)}& K_{\mathrm{1}}=\frac{E_{\mathrm{1}}}{\mathrm{3}\left(\mathrm{1}-\mathrm{2}{\mu }_{\mathrm{1}}\right)},G_{\mathrm{1}}=\frac{E_{\mathrm{1}}}{\mathrm{2}\left(\mathrm{1}+{\mu }_{\mathrm{1}}\right)}\end{array}$$is the bulk modulus and shear modulus of the inclusion (Xu [15]).

Let$$\begin{array}{}\text{(23)}& \alpha =\mathrm{9}{K}_p{K}_{\mathrm{1}}+\frac{K_{\mathrm{1}}}{K_{\mathrm{0}}}+\frac{\mathrm{9}{K}_p{K}_{\mathrm{1}}-K_{\mathrm{1}}/K_{\mathrm{0}}}{K_{\mathrm{0}}/\left(K_{\mathrm{0}}-K_{\mathrm{1}}\right)-\left(\rho +\mathrm{9}{K}_p{K}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)},\text{(24)}& \beta =\mathrm{4}{G}_p{G}_{\mathrm{1}}+\frac{G_{\mathrm{1}}}{G_{\mathrm{0}}}+\frac{\mathrm{4}{G}_p{G}_{\mathrm{1}}-G_{\mathrm{1}}/G_{\mathrm{0}}}{G_{\mathrm{0}}/\left(G_{\mathrm{0}}-G_{\mathrm{1}}\right)-\left(\rho +\mathrm{4}{G}_p{G}_{\mathrm{0}}\left(\mathrm{1}-\rho \right)\right)},\end{array}$$Equation (25) is obtained by (19):$$\begin{array}{}\text{(25)}& {\sigma }^L=\left[\begin{array}{cc}\begin{array}{ccc}\alpha +\frac{\mathrm{4}}{\mathrm{3}}\beta & \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta & \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta \\ \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta & \alpha +\frac{\mathrm{4}}{\mathrm{3}}\beta & \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta \\ \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta & \alpha -\frac{\mathrm{2}}{\mathrm{3}}\beta & \alpha +\frac{\mathrm{4}}{\mathrm{3}}\beta \end{array}& \begin{array}{ccc}& & \\ & & \\ & & \end{array}\\ \begin{array}{ccc}& & \\ & & \\ & & \end{array}& \begin{array}{ccc}\mathrm{2}\beta & & \\ & \mathrm{2}\beta & \\ & & \mathrm{2}\beta \end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\sigma }_{\mathrm{11}}^{\mathrm{0}}\\ {\sigma }_{\mathrm{22}}^{\mathrm{0}}\\ {\sigma }_{\mathrm{33}}^{\mathrm{0}}\end{array}\\ \begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right]\end{array}$$

3.3. Examples

(1) The matrix material is made to have an elastic modulus of 1MPa and a Poisson’s ratio of 0.3. The inclusion material has a Poisson’s ratio of 0.3, an elastic modulus $E_1$ of 0.2, 0.4, 0.6, 0.8, 1, 2, 3, 4, and 5, and an external load applied to the rock is a uniaxial pressure of 6 MPa. The variation of the component of the locked-in stress parallel to the loading surface with the content of the inclusion and the ratio of the elastic modulus of the inclusion to the matrix can be obtained by (25), as shown in Figure 3(a). As the elastic modulus of the inclusion changes from small to large, the locked-in stress generated in the inclusion is also increasing.

[figures omitted; refer to PDF]

(2) The external load applied to the rock is a hydrostatic pressure of 6MPa, and the other conditions are the same as (1), as shown in Figure 3(b). As the elastic modulus of the inclusion changes from small to large, the locked-in stress generated in the inclusion is also increasing.

(3) The matrix material is made to have an elastic modulus of 1MPa and a Poisson’s ratio of 0.3. The inclusion material has a modulus of elasticity of 1 and a Poisson’s ratio of 0.1, 0.2, 0.3, 0.4, and 0.49, and an external load applied to the rock is a uniaxial pressure of 6MPa. As the Poisson’s ratio of the inclusion changes from small to large, the locked-in stress generated in the inclusion is also increasing.

(4) The external load applied to the rock is a hydrostatic pressure of 6MPa, and the other conditions are the same as (3), as shown in Figure 4(b). As the Poisson’s ratio of the inclusion changes from small to large, the locked-in stress generated in the inclusion is also increasing.

[figures omitted; refer to PDF]

When the elastic parameters of the inclusions are the same as the matrix, This model degenerates into a homogeneous medium, and the locked-in stress generated in the inclusions is independent of the inclusion content. As the inclusion content increases, the influence of elastic parameters becomes smaller and smaller.

Consider the case of uniaxial compression; there is ${\sigma }_{22}^0={\sigma }_{33}^0=\mathrm{0}$ at this time. Test data of Geng [16] is cited here. The matrix material has an elastic modulus of 2560 MPa and a Poisson’s ratio of 0.3; the elastic modulus of inclusion is $E_1=0.113T+8.99\mathrm{M}\mathrm{P}\mathrm{a}$, and the Poisson’s ratio is 0.49. The relation between ${\sigma }^L$ and $\rho$ can be obtained from (27). Figure 5 shows the increase of the locked-in stress when the uniaxial pressure is 2MPa, 4MPa, 6MPa, and 8MPa, respectively, and it is related to the initial locked-in stress and inclusion content.

[figure omitted; refer to PDF]

Secondly, consider the case of hydrostatic compression; i.e. ${\sigma }_{11}^0={\sigma }_{22}^0={\sigma }_{33}^0$. Figure 6 shows the increase of the locked-in stress when the hydrostatic pressure is 2 MPa, 4 MPa, and 6 MPa, respectively, and it is related to the initial locked-in stress and inclusion content.

[figure omitted; refer to PDF]

From Figures 5 and 6 it can be seen that the external load will generate additional stress on the inclusions with locked-in stress, which is related to the material properties and external forces. The effect of the locked-in stress shown in the figure on the results is due to the fact that the locked-in stress of the inclusion in the object studied here is related to the temperature, and the volumetric elasticity of the inclusion is also related to the temperature, which is essentially the effect of temperature. There is no connection with the locked-in stress. When the locked-in stress is small, the increase of it is not negligible. The increase is proportional to the increase in uniaxial compression load and hydrostatic pressure load.

The derivation of (1), (2), (3), (4), (5), (7), (8), (9), (10), (11), (12), (13), (14), (15), (16), (17), (18), (19), (23), (24), and (25) is also applicable to other forms of inclusion material. In fact, this increase in the locked-in stress is different from the locked-in stress but is the additional stress in the loading. After unloading this increment disappears but the locked-in stress does not.

4. Theoretical Study on the Change of Elastic Modulus of Rocks with Locked-in Stress

4.1. Solution without Considering Locked-In Stress

(1 ) Micromechanical Analysis of Equivalent Elastic Modulus. The upper limit of inclusion content is 1 and the lower limit is 0. When the content of the inclusion is 0, it is obvious that the elastic modulus of the rock is the elastic modulus of the intact rock. When the content reaches 1, it becomes a sample composed of the inclusion material, and its elastic modulus is consistent with the inclusion material. Therefore, the relationship between the elastic modulus of inclusion-containing specimens and the variation of inclusion content should satisfy these two prerequisites.

First, the relationship between elastic modulus and content is analyzed. When calculating the elastic modulus, both the matrix and inclusions are considered to be in a linear elastic state. Therefore, the problem of the elastic modulus in this paper is simplified to discuss the problem of spherical inclusions in elastic materials. For the elastic solution of inclusion-containing materials, there are currently Voigt upper limits and Reuss lower limits:$$\begin{array}{c}\text{(26)}& \stackrel{-}{E}=\left(\mathrm{1}-\rho \right)E_{\mathrm{0}}+\rho E_{\mathrm{1}}\left(\mathrm{V}\mathrm{o}\mathrm{i}\mathrm{g}\mathrm{t}\hspace{0.17em}\hspace{0.17em}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}\hspace{0.17em}\hspace{0.17em}\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{t}\right)\stackrel{-}{E}={\left(\left(\mathrm{1}-\rho \right){E_{\mathrm{0}}}^{-\mathrm{1}}+\rho {E_{\mathrm{1}}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\left(\mathrm{R}\mathrm{e}\mathrm{u}\mathrm{s}\mathrm{s}\hspace{0.17em}\hspace{0.17em}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\hspace{0.17em}\hspace{0.17em}\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{t}\right).\end{array}$$

Taylor method is$$\begin{array}{}\text{(27)}& \stackrel{-}{C}=C^{\mathrm{0}}+{\displaystyle \sum }_{r=\mathrm{1}}^{N-\mathrm{1}}\rho {\left[{\left(C^r-C^{\mathrm{0}}\right)}^{-\mathrm{1}}+P^r\right]}^{-\mathrm{1}}.\end{array}$$

M-T method is$$\begin{array}{}\text{(28)}& \stackrel{-}{C}=C^{\mathrm{0}}+{\displaystyle \sum }_{r=\mathrm{1}}^{N-\mathrm{1}}\rho {\left[{\left(C^r-C^{\mathrm{0}}\right)}^{-\mathrm{1}}+\left(\mathrm{1}-\rho \right)P^r\right]}^{-\mathrm{1}},\end{array}$$where

 

$\stackrel{-}{C}$ is the macroscopically equivalent elastic tensor of the rock containing the inclusion,

 

$C^0$ is the elastic tensor of the rock,

 

$C^r$ is the elastic tensor of the r-phase inclusions.

In addition, there are differential method, generalized self-consistent method, and Hashin and Shtrikman’s upper and lower limit methods. For cases with a single inclusion, (27) to (28) can be written in the following scalar form.$$\begin{array}{c}\text{(29)}& \stackrel{-}{K}=K_{\mathrm{0}}+\frac{K_{\mathrm{1}}-K_{\mathrm{0}}}{\mathrm{1}+\mathrm{9}{K}_p\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)}\rho \stackrel{-}{G}=G_{\mathrm{0}}+\frac{G_{\mathrm{1}}-G_{\mathrm{0}}}{\mathrm{1}+\mathrm{4}{G}_p\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)}\rho \\ \text{(30)}& \stackrel{-}{G}=G_{\mathrm{0}}+\frac{G_{\mathrm{1}}-G_{\mathrm{0}}}{\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)}\rho \stackrel{-}{K}=K_{\mathrm{0}}+\frac{K_{\mathrm{1}}-K_{\mathrm{0}}}{\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)}\rho \end{array}$$with$$\begin{array}{c}\text{(31)}& {\stackrel{-}{K}}_p=\frac{\mathrm{1}}{\mathrm{3}\left(\mathrm{4}\stackrel{-}{G}+\mathrm{3}\stackrel{-}{K}\right)},{\stackrel{-}{G}}_p=\frac{\mathrm{3}\left(\mathrm{2}\stackrel{-}{G}+\stackrel{-}{K}\right)}{\mathrm{10}\stackrel{-}{G}\left(\mathrm{4}\stackrel{-}{G}+\mathrm{3}\stackrel{-}{K}\right)}.\end{array}$$The equivalent elastic modulus can be obtained according to$$\begin{array}{}\text{(32)}& \stackrel{-}{E}=\frac{\mathrm{9}\stackrel{-}{K}\hspace{0.17em}\stackrel{-}{G}}{\mathrm{3}\stackrel{-}{K}+\stackrel{-}{G}}\end{array}$$

The matrix material is made to have an elastic modulus of 1MPa and a Poisson’s ratio of 0.3; the inclusion material is made to have a Poisson’s ratio of 0.3 and an elastic modulus $E_1$ of 0.2, 0.6. The variation of the equivalent elastic modulus with the content of inclusions and the ratio of the elastic modulus of the inclusions to the matrix can be obtained by (26), (27), (28), (29), (30), and (32) (Figure 7(a)). Let the matrix material have an elastic modulus of 1MPa and a Poisson’s ratio of 0.3; let Poisson’s ratio of the inclusion be 0.3 and the elastic modulus of the inclusion be 2 and 3. The variation of the equivalent elastic modulus with the inclusion content and the ratio of the elastic modulus of inclusion to the matrix can be obtained by (26), (27), (28), (29), (30), and (32) (Figure 7(b)).

[figures omitted; refer to PDF]

When the elastic modulus of inclusion is larger than that of the matrix, the macroelastic modulus increases with the increase of the inclusion content, and the elastic modulus of inclusion increases, when the elastic modulus of inclusion is smaller than that of the matrix. The macroelastic modulus decreases with the increase of inclusion content and increases with the increase of inclusion elastic modulus. Under the same conditions, the value obtained by the Voigt method is the highest, the value obtained by the Reuss method is the lowest, and the M-T method curve is slightly higher than the Taylor method curve. In the case where the elastic modulus of the inclusion is small enough, the curve of the Taylor method will be distorted (Figure 7).

There is no Poisson’s ratio in the expressions of elastic modulus of Voigt upper limit method and Reuss lower limit method, so the change of Poisson’s ratio will not affect the macroscopic elastic modulus, so it will not be discussed here.

Secondly, we discuss the influence of the Poisson’s ratio change on the elastic modulus when the Taylor method and the M-T method are used. Let the elastic modulus and Poisson ratio of the matrix be 1 and 0.3, the elastic modulus of inclusion be 2 and 0.2, and Poisson’s ratio be 0.2, 0.4. Figure 8 shows that as Poisson’s ratio increases, the slope of the curve becomes larger. The macroelastic modulus at the same inclusion content increases. Under the same conditions, curve of the M-T method is slightly higher than that of Taylor method, but they are very closed.

[figure omitted; refer to PDF]

( 2) Test Verification. Using rock-like materials and rubber particles to make a test block rock-like model containing inclusions, Geng [16] studied the influence of locked-in stress on the mechanical properties of rocks. By heating the test piece, the rubber particles expand, creating a closed stress in the inclusion. Uniaxial compression tests of similar model blocks with different locked-in stresses and different inclusion contents were performed with a water bath and MTS-5T uniaxial compression testing machine. Results show that the peak strength and elastic modulus decrease with the increase of locked-in stress and decrease with the increase of the inclusion content. With the decrease of the inclusion content, the effect of the closed stress on the peak strength and elastic modulus gradually decreases. When the content of the inclusion is zero, the influence of the closure stress also disappears. Through the mathematical fitting analysis, the quantitative formula for calculating the relationship between the peak strength, the elastic modulus of the test piece, the locked-in stress, and the inclusion content was obtained. The expression and value curve of joint impact factor are given. Results of this study reveal the importance of the locked-in stress in the mechanical properties of the rock and provide a reference for the further study of the locked-in stress in rocks.

Substituting the mechanical parameters of the matrix material and the inclusion material in the test of Geng [16] into (26), (27), and (28), the resulting equation is plotted as a curve (Figure 9).

[figure omitted; refer to PDF]

Figure 9 depicted the elastic modulus-content relationship obtained by different methods when the locked-in stress is 0. It can be seen that both the M-T method and the Taylor method can well match the experimental data in the 0-0.04 inclusion stage.

The existing research results suggest that the Taylor method is only applicable to the case where the inclusion content is less than 5% (Du [17]). This is because when the content of inclusions exceeds 5%, the interaction of inclusions will become unnegligible.

There is small but negligible deviation between calculated results and the experimental results. Equations (29) and (30) are similar in structure, where the denominator has one more parameter related to inclusion content. In general, the M-T method is more suitable for this article.

In summary, when there is no locked-in stress in the inclusion, the scalar form of the macroelastic modulus of the sample can be obtained by substituting (30) into (32):$$\begin{array}{}\text{(33)}& \stackrel{-}{E}=\frac{\mathrm{9}\left(K_{\mathrm{0}}+\left(\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)/\left(\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)\right)\right)\rho \right)\left(G_{\mathrm{0}}+\left(\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)/\left(\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)\right)\right)\rho \right)}{\mathrm{3}\left(K_{\mathrm{0}}+\left(\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)/\left(\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p\left(K_{\mathrm{1}}-K_{\mathrm{0}}\right)\right)\right)\rho \right)+\left(G_{\mathrm{0}}+\left(\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)/\left(\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p\left(G_{\mathrm{1}}-G_{\mathrm{0}}\right)\right)\right)\rho \right)}\end{array}$$

4.2. Solution with Considering Locked-In Stress

The relationship between the equivalent elastic modulus and the locked-in stress can not be found by the equation $C_1({\epsilon }_{ij}^0+{\epsilon }_{ij}^{\prime }-{\epsilon }_{ij}^h)=C_0({\epsilon }_{ij}^0+{\epsilon }_{ij}^{\prime }-{\epsilon }^{\ast 1})$ derived from Eshelby and the M-T method. The change in the volumetric content of inclusions due to intrinsic strain is negligible relative to the volume content. In the M-T method, the locked-in stress is an “irrelevant quantity”, because both the Eshelby theory and the M-T method have a preset condition; namely, it’s continuous between media. At present, the research on the inclusion of stress in rocks mainly focuses on the inclusion of fluids, such as pore water or pore gas, but such research is rare. Liu [18] studied the effects of pore pressure changes on rock strength characteristics and found that the increase in pore pressure caused a decrease in elastic modulus. As the pore pressure increases, some cracks open and the contact between the particles decreases, resulting in a change in the matrix skeleton unit. Li [19] studied the relationship between the shear modulus of saturated sand and water pressure through a saturated sand liquefaction test under triaxial conditions.

(1) Mesoscopic Analysis of Equivalent Elastic Modulus. The locked-in stress studied in this article is essentially different from the pore pressure caused by gas or fluid which is isotropic. Due to the fluidity, the geometry of the fluid changes with the pore morphology during rock damage, while the solid inclusions remain relatively intact. Therefore, the locked-in stress field of fluid inclusions under external loads is more complicated. But as the locked-in stress increases, the trend of the decrease of the elastic modulus of the test block does not change.

The locked-in stress caused by the expansion of the inclusion body can not only squeeze the crack parallel to the spherical surface of the inclusion body to make it more compact, but also open the crack perpendicular to the spherical surface and increase the damage of the rock. This initial damage makes the mechanical properties of matrix material change, and the shear modulus and bulk modulus have undergone a certain degree of reduction. Assume that the bulk modulus and shear modulus of the matrix after reduction are $K_0^{\prime }$ and $G_0^{\prime }$, respectively. Therefore, the macroelastic modulus of the test block is$$\begin{array}{}\text{(34)}& \stackrel{-}{K}=K_{\mathrm{0}}^{\prime }+\frac{K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }}{\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p^{\prime }\left(K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }\right)}\rho ,\text{(35)}& \stackrel{-}{G}=G_{\mathrm{0}}^{\prime }+\frac{G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }}{\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p^{\prime }\left(G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }\right)}\rho ,\end{array}$$with$$\begin{array}{c}\text{(36)}& K_p^{\prime }=\frac{\mathrm{1}}{\mathrm{3}\left(\mathrm{4}{G}_{\mathrm{0}}^{\prime }+\mathrm{3}{K}_{\mathrm{0}}^{\prime }\right)},G_p^{\prime }=\frac{\mathrm{3}\left(\mathrm{2}{G}_{\mathrm{0}}^{\prime }+K_{\mathrm{0}}^{\prime }\right)}{\mathrm{10}{G}_{\mathrm{0}}^{\prime }\left(\mathrm{4}{G}_{\mathrm{0}}^{\prime }+\mathrm{3}{K}_{\mathrm{0}}^{\prime }\right)}.\end{array}$$

Therefore, the change of the mechanical properties of the sample is the result of the joint action of the inclusion content and the locked-in stress. When studying the statistical damage constitutive model, Li [20] equated damaged rock materials with pores, and the volume content of the initial damage (i.e., additional pore) caused by the locked-in stress is $P$. In 1961, Greertsma and Smit [17] obtained an empirical formula for the relationship between rock bulk modulus and porosity based on experimental results. This paper assumes that this empirical relationship also applies to the relationship between shear modulus and porosity, yielding $$\begin{array}{}\text{(37)}& K_{\mathrm{0}}^{\prime }=\frac{K_{\mathrm{0}}}{\mathrm{1}+aP},\text{(38)}& G_{\mathrm{0}}^{\prime }=\frac{G_{\mathrm{0}}}{\mathrm{1}+aP},\end{array}$$

and $P$ is a function of the locked-in stress ${\sigma }^L$ and the inclusion content $\rho$. The equation is$$\begin{array}{}\text{(39)}& P=f\left({\sigma }^L,\rho \right).\end{array}$$

Since a and b are undetermined coefficients, (37) and (38) can be abbreviated as$$\begin{array}{}\text{(40)}& K_{\mathrm{0}}^{\prime }=\frac{K_{\mathrm{0}}}{\mathrm{1}+f\left({\sigma }^L,\rho \right)},\text{(41)}& G_{\mathrm{0}}^{\prime }=\frac{G_{\mathrm{0}}}{\mathrm{1}+f\left({\sigma }^L,\rho \right)},\end{array}$$

Equation (39) is the problem of expansion damage of spherical inclusions. This formula is discussed in detail below.

(2) Damage Mechanical Analysis of Spherical Inclusions. The initial damage is caused by the squeezing of the locked-in stress. The expansion of the inclusion, in addition to the strain energy created by the expansion of the inclusion in itself, also performs work on the base material. The volume of compression of the inclusion is $\Delta r^{\prime }$, and the amount of expansion of the matrix body is $\Delta r_0=\Delta r-\Delta r^{\prime }$. $\Delta r$ is the intrinsic radius increment of the inclusion (Figure 2). According to the expression in Geng [16] regarding locked-in stress,$$\begin{array}{}\text{(42)}& {\sigma }^L=\left({\alpha }_{\mathrm{1}}-{\alpha }_{\mathrm{0}}\right)\frac{\mathrm{12}{G}_{\mathrm{0}}{K}_{\mathrm{1}}}{\mathrm{3}{K}_{\mathrm{1}}+\mathrm{4}{G}_{\mathrm{0}}}\left(T-T_{\mathrm{0}}\right)\end{array}$$Equations (43) and (44) are obtained:$$\begin{array}{}\text{(43)}& \Delta r^{\prime }=\frac{\mathrm{1}}{\mathrm{3}{K}_{\mathrm{1}}}\frac{\Delta r}{\left({\mu }_{\mathrm{0}}+\mathrm{1}\right)/\mathrm{2}{E}_{\mathrm{0}}+\left(\mathrm{1}-\mathrm{2}{\mu }_{\mathrm{1}}\right)/E_{\mathrm{1}}}=\frac{\mathrm{4}{G}_{\mathrm{0}}}{\mathrm{3}{K}_{\mathrm{1}}+\mathrm{4}{G}_{\mathrm{0}}}\Delta r\text{(44)}& \Delta r_{\mathrm{0}}=\frac{\mathrm{3}{K}_{\mathrm{1}}}{\mathrm{3}{K}_{\mathrm{1}}+\mathrm{4}{G}_{\mathrm{0}}}\Delta r\end{array}$$

When the locked-in stress is caused by the thermal expansion of inclusions, $\Delta r$ can be expressed by (39):$$\begin{array}{}\text{(45)}& \Delta r=r_{\mathrm{1}}\left({\alpha }_{\mathrm{1}}-{\alpha }_{\mathrm{0}}\right)\Delta T\end{array}$$

$\Delta r$ is the increase in radius due to free thermal expansion when the rubber ball is unconstrained. The thermal expansion of the matrix is neglected here since the coefficient of it is extremely small relative to the rubber ball. According to (42), (43), and (39), it can be obtained that, under the constraint of the matrix, the sphere is compressed by$$\begin{array}{}\text{(46)}& \Delta r^{\prime }=\frac{\left({\alpha }_{\mathrm{1}}-{\alpha }_{\mathrm{0}}\right)}{\mathrm{3}{K}_{\mathrm{1}}/\mathrm{4}{G}_{\mathrm{0}}+\mathrm{1}}\Delta T{r}_{\mathrm{1}}.\end{array}$$

From (44), the radius of the inner wall of the matrix of is expanded by$$\begin{array}{}\text{(47)}& \Delta r_{\mathrm{0}}=r_{\mathrm{1}}\left({\alpha }_{\mathrm{1}}-{\alpha }_{\mathrm{0}}\right)\Delta T\left(\frac{\mathrm{3}{K}_{\mathrm{1}}}{\mathrm{3}{K}_{\mathrm{1}}+\mathrm{4}{G}_{\mathrm{0}}}\right).\end{array}$$

Therefore, the internal energy generated in the sphere during the generation of the locked-in stress is$$\begin{array}{}\text{(48)}& W_{\mathrm{1}}=\frac{\mathrm{4}\pi r^{\mathrm{3}}}{\mathrm{3}}\frac{{\left({\sigma }^L\right)}^{\mathrm{2}}}{\mathrm{2}{K}_{\mathrm{1}}}.\end{array}$$

The work done by the expansion of individual rubber particles on the matrix is$$\begin{array}{}\text{(49)}& W_{\mathrm{0}}=\frac{\Delta r_{\mathrm{0}}}{\Delta r^{\prime }}·W_{\mathrm{1}}=\frac{\mathrm{4}\pi {r_{\mathrm{1}}}^{\mathrm{3}}}{\mathrm{3}}\frac{\Delta r_{\mathrm{0}}}{\Delta r^{\prime }}\frac{{\left({\sigma }^L\right)}^{\mathrm{2}}}{\mathrm{2}{K}_{\mathrm{1}}}.\end{array}$$

The ratio of $\Delta r_0$ and $\Delta r^{\prime }$ can be obtained from (40) and (47):$$\begin{array}{}\text{(50)}& \frac{\Delta r_0}{\Delta r^{\prime }}=\frac{3K_1}{4G_0}.\end{array}$$

Substituting (50) into (49) yields$$\begin{array}{}\text{(51)}& W_{\mathrm{0}}=\frac{\pi {r_{\mathrm{1}}}^{\mathrm{3}}{\left({\sigma }^L\right)}^{\mathrm{2}}}{\mathrm{2}{G}_{\mathrm{0}}}.\end{array}$$

The total rubber content is $\rho$, and the volume of the test block is 1. The total work done by all rubbers in the test block on the matrix is$$\begin{array}{}\text{(52)}& W_{\mathrm{0}W}=\frac{\rho }{\mathrm{4}\pi r_{\mathrm{1}}^{\mathrm{3}}/\mathrm{3}}{W}_{\mathrm{0}}=\frac{\mathrm{3}\rho {\left({\sigma }^L\right)}^{\mathrm{2}}}{\mathrm{8}{G}_{\mathrm{0}}}.\end{array}$$

According to the uniaxial compression test results of the matrix material, the required strain energy for destruction of matrix material is$$\begin{array}{}\text{(53)}& W_c=\frac{{{\sigma }_c}^{\mathrm{2}}}{\mathrm{2}{E}_{\mathrm{0}}}.\end{array}$$where

$W_c$ is the strain energy required for the matrix material destruction, and the damage degree of the matrix material under uniaxial compression failure is 1. The damage degree of the matrix material when no external load is applied is 0. Let the relationship between strain energy and damage degree follow power function distribution:$$\begin{array}{}\text{(54)}& P\left(\frac{W_{\mathrm{0}W}}{W_c}\right)=\frac{h}{g}{\left(\frac{W_{\mathrm{0}W}/W_c}{g}\right)}^{h-\mathrm{1}},\end{array}$$where $h$ and $g$ are distribution parameters that can be combined with the test.

$D$ indicate the degree of damage:$$\begin{array}{}\text{(55)}& D={{\displaystyle \int }}_0^{W_{0W}/W_c}P\left(\frac{W_{0W}}{W_c}\right)={\left(\frac{W_{0W}}{W_{c}g}\right)}^h.\end{array}$$

Combining (52) and (53) yields$$\begin{array}{}\text{(56)}& f\left({\sigma }_L,\rho \right)={\left(\frac{W_{\mathrm{0}W}}{W_{c}g}\right)}^h={\left(\frac{\mathrm{3}\left(\mathrm{1}+\mu \right)\rho {\left({\sigma }^L\right)}^{\mathrm{2}}}{\mathrm{2}g}\right)}^h.\end{array}$$

Substituting (56) into (40) and (41) yields$$\begin{array}{}\text{(57)}& K_{\mathrm{0}}^{\prime }=\frac{K_{\mathrm{0}}}{\mathrm{1}+{\left(\mathrm{3}\left(\mathrm{1}+\mu \right)\rho {\left({\sigma }^L\right)}^{\mathrm{2}}/\mathrm{2}g\right)}^h},\text{(58)}& G_{\mathrm{0}}^{\prime }=\frac{G_{\mathrm{0}}}{\mathrm{1}+{\left(\mathrm{3}\left(\mathrm{1}+\mu \right)\rho {\left({\sigma }^L\right)}^{\mathrm{2}}/\mathrm{2}g\right)}^h}.\end{array}$$

Substituting (57) and (58) into (32) yields$$\begin{array}{}\text{(59)}& \stackrel{-}{E}=\frac{\mathrm{9}\left(K_{\mathrm{0}}^{\prime }+\left(\left(K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }\right)/\left(\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p^{\prime }\left(K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }\right)\right)\right)\rho \right)\left(G_{\mathrm{0}}^{\prime }+\left(\left(G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }\right)/\left(\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p^{\prime }\left(G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }\right)\right)\right)\rho \right)}{\mathrm{3}\left(K_{\mathrm{0}}^{\prime }+\left(\left(K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }\right)/\left(\mathrm{1}+\mathrm{9}\left(\mathrm{1}-\rho \right)K_p^{\prime }\left(K_{\mathrm{1}}-K_{\mathrm{0}}^{\prime }\right)\right)\right)\rho \right)+\left(G_{\mathrm{0}}^{\prime }+\left(\left(G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }\right)/\left(\mathrm{1}+\mathrm{4}\left(\mathrm{1}-\rho \right)G_p^{\prime }\left(G_{\mathrm{1}}-G_{\mathrm{0}}^{\prime }\right)\right)\right)\rho \right)}.\end{array}$$

(3) Test Verification. The test data of Geng [16] is cited here; the locked-in stress is caused by the uncoordinated thermal expansion of the matrix and inclusions. Substituting (57) and (58) into (59) yields an expression of the elastic modulus of the test block when the inclusion contains locked-in stress. Then, the three-dimensional nonlinear fitting of the elastic modulus test data of the Geng [16] by (61) results in $h=0.45$, $g=0.0114$, ($R^2=0.95975$). The fitted result is shown in Figure 10. When ${\stackrel{-}{\sigma }}_L$ in (57) and (58) is 0, (59) degenerates to (33).

[figure omitted; refer to PDF]

5. Strength Change Mechanism of Test Block

5.1. Theoretical Analysis of Peak Strength without Considering Locked-In Stress

The failure mode of the sample is shear failure. Therefore, the Mohr-Coulomb failure criterion is used (Figure 11), in the form of$$\begin{array}{}\text{(60)}& \tau =c+{\sigma }_{n}tg\phi ,\text{(61)}& {\sigma }_{\mathrm{n}}=\frac{{\sigma }_{\mathrm{1}}\left(\mathrm{1}+\cos \mathrm{2}\beta \right)}{\mathrm{2}}.\end{array}$$where $\tau$ is the shear strength, $\sigma$ is the normal stress of the shear plane, $c$ is the cohesion of the sample, and $\phi$ is the internal friction angle of the sample.

[figure omitted; refer to PDF]

For the case of uniaxial compressive stress ${\sigma }_{\mathrm{1}}$, the specimen is broken along the surface at a angle of $\beta =\pi /\mathrm{4}+\phi /\mathrm{2}={\mathrm{60.55}}^{°}$ to the cross-section and $\tau$ can be expressed as$$\begin{array}{}\text{(62)}& \tau =\frac{{\sigma }_{\mathrm{1}}\sin \mathrm{2}\beta }{\mathrm{2}}.\end{array}$$

In the case of softer inclusions, the shear modulus of inclusion is much smaller than that of the matrix, and basically no shear stress is assumed in the sample. For normal stress, the presence of inclusions reduces the area of the matrix on the shear plane.

Due to the small deformation of the specimen, the inclusion can be considered to be always spherical with a radius of r and a content of $f_p$. Take a representative cell with a volume of 1 containing a single inclusion with a volume of $f_p$ as an example at the cut section, and the six faces of the cell are parallel or perpendicular to the shear plane (Figure 12). Therefore, according to the volume formula of the sphere, there is a radius $r={\left(3\rho /4\pi \right)}^{1/3}$ of the inclusion. The part of the sheared broken section passing through the cell is a square region with an area of 1. The section of the sheared broken section passing the inclusion is a circle and its area is$$\begin{array}{}\text{(63)}& {\rho }_s=\pi {\left(\frac{3\rho }{4\pi }\right)}^{2/3}.\end{array}$$

[figures omitted; refer to PDF]

The area of the matrix material part on the broken surface is 1-${\rho }_s$. Assuming that the elastic modulus and shear modulus of the inclusions are very small compared to rock-like materials, it can be considered that there is basically no ability to carry compressive and shear stresses, and the bearing capacity is close to that of the holes compared to the matrix made of rock-like materials. Therefore, the average normal stress and the average shear stress of the matrix part on the shear surface are not equal to $\sigma$ and $\tau$. According to Lemaitre’s strain equivalence assumption (Lemaitre [21]), the true stresses on the shear plane are the effective normal stress ${\sigma }^{\prime }$ and the effective shear stress ${\tau }^{\prime }$, and the corresponding $\sigma$ and $\tau$ are called the nominal normal stress and the nominal shear stress here.

Figure 13 shows the Mohr-Coulomb model based on effective normal stress and effective shear stress. It can be seen that, as the content of inclusions increases, the nominal stress required for the shear plane to reach failure is getting smaller and smaller. According to the assumption of strain equivalence, the effective normal stress and the effective shear stress on the failure surface of the sample containing the inclusion are$$\begin{array}{}\text{(64)}& {\tau }^{\prime }=\frac{\tau }{\mathrm{1}-{\rho }_s},\text{(65)}& {\sigma }_n^{\prime }=\frac{\sigma }{\mathrm{1}-{\rho }_s}.\end{array}$$

[figure omitted; refer to PDF]

According to the theory of elasticity, the relationship between the axial stress ${\sigma }_{\mathrm{1}}$, ${\sigma }_n^{\prime }$, and ${\tau }^{\prime }$ of the specimen before the specimen is destroyed always conforms to (61) and (62). Substituting (61) and (62) into (64) and (65) yields$$\begin{array}{}\text{(66)}& {\sigma }_n^{\prime }=\frac{{\sigma }_1\left(1+\cos 2\beta \right)}{2\left(1-{\rho }_s\right)},\text{(67)}& {\tau }^{\prime }=\frac{{\sigma }_1\sin 2\beta }{2\left(1-{\rho }_s\right)}.\end{array}$$

From (66), (67), and (60), $\sigma$ is obtained:$$\begin{array}{}\text{(68)}& \sigma =\frac{\mathrm{2}c\left(\mathrm{1}-{\rho }_s\right)}{\sin \mathrm{2}\beta -\left(\mathrm{1}+\cos \mathrm{2}\beta \right)tg\phi }\end{array}$$

Taking the assumptions of the distribution of the inclusions into account, the above formula is modified. Parameter $k$ is introduced in (68) and results in$$\begin{array}{}\text{(69)}& {\sigma }_c=\frac{\mathrm{2}ck\left(\mathrm{1}-{\rho }_s\right)}{\sin \mathrm{2}\beta -\left(\mathrm{1}+\cos \mathrm{2}\beta \right)tg\phi }\end{array}$$

Equation (69) needs to satisfy two conditions: (1) when ${\rho }_s=\mathrm{0}$, it degenerates to the strength of the matrix; (2) when ${\rho }_s=\mathrm{1}$, there is no matrix present on the failure surface and the intensity is zero. That is, in (69), $k(1-{\rho }_s)$ needs to satisfy the following: (1) when ${\rho }_s=\mathrm{0}$, $k(1-{\rho }_s)=\mathrm{1}$; (2) when ${\rho }_s=\mathrm{1}$, $k(1-{\rho }_s)=\mathrm{0}$.

Let $k=f(\rho )$; it was found that when $k$ meets the following formula, it can meet the above two conditions:$$\begin{array}{}\text{(70)}& k=\frac{\mathrm{1}}{\mathrm{1}+d\rho }.\end{array}$$

From formula (69) and formula (70), ${\sigma }_c$ is obtained:$$\begin{array}{}\text{(71)}& {\sigma }_c={\sigma }_0\frac{\left(1-{\rho }_s\right)}{1+d\rho }.\end{array}$$where parameter $d$ can be obtained in combination with the uniaxial compression test data.

5.2. Theoretical Analysis of Peak Strength Considering Locked-In Stress

Figure 14 is a schematic diagram of the force of a representative unit on the failure surface of the sample when the inclusion contains a locked-in stress. When the initial locked-in stress of the inclusion is not 0, the influence of the stress concentration effect of the inclusion must be taken into account as well as the influence of the initial locked-in stress. The effect of the initial locked-in stress on the failure shear plane of the specimen is to cancel out a part of the normal stress on the shear plane and reduce the maximum shear stress required for the failure. Let the normal stress and shear stress on the failure surface be ${\sigma }_n^{\prime \prime }$ and ${\tau }^{\prime \prime }$:$$\begin{array}{}\text{(72)}& {\sigma }_n^{\prime \prime }=\frac{{\sigma }_{\mathrm{1}}\left(\mathrm{1}+\cos \mathrm{2}\beta \right)}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}-\frac{{\sigma }^L{\rho }_s}{\mathrm{1}-{\rho }_s},\text{(73)}& {\tau }^{\prime \prime }={\tau }^{\prime }=\frac{{\sigma }_{\mathrm{1}}\sin \mathrm{2}\beta }{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}.\end{array}$$

[figure omitted; refer to PDF]

Substituting (72) and (73) into (60) yields$$\begin{array}{}\text{(74)}& {\sigma }_c=\frac{\mathrm{2}c\left(\mathrm{1}-{\rho }_s-{\sigma }^L{\rho }_{s}tg\phi /c\right)}{\sin \mathrm{2}\beta -\left(\mathrm{1}+\cos \mathrm{2}\beta \right)tg\phi }.\end{array}$$

Considering the stress concentration effect around the inclusion, the locked-in stress term of (74) is modified and the correction coefficient $q$ is introduced, which is considered to be related to the inclusion content:$$\begin{array}{}\text{(75)}& q=f\left(\rho \right).\end{array}$$

Equation (74) is amended to$$\begin{array}{}\text{(76)}& {\sigma }_c=\frac{\mathrm{2}c\left(\mathrm{1}-{\rho }_s-tg\phi {\sigma }^L{\rho }_s/qc\right)}{\sin \mathrm{2}\beta -\left(\mathrm{1}+\cos \mathrm{2}\beta \right)tg\phi }.\end{array}$$

Combined with the correction factor $k\hspace{0.17em}\hspace{0.17em}(k=1/\left(1+d\rho \right))$ related to the distribution of inclusions in (71), we can get$$\begin{array}{}\text{(77)}& {\sigma }_c=\frac{\mathrm{2}c\left(\mathrm{1}-{\rho }_s-tg\phi {\sigma }^L{\rho }_s/qc\right)}{\left(\sin \mathrm{2}\beta -\left(\mathrm{1}+\cos \mathrm{2}\beta \right)tg\phi \right)·\left(\mathrm{1}+d\rho \right)}\end{array}$$

Simplify (77) and obtain$$\begin{array}{}\text{(78)}& {\sigma }_c={\sigma }_{\mathrm{0}}\frac{\mathrm{1}-{\rho }_s-{\sigma }^L{\rho }_s\tan \phi /qc}{\mathrm{1}+d\rho }\end{array}$$where $d$ is the parameter determined by the uniaxial compression test of rock-like blocks without locked-in stress and $q$ is the parameter related to inclusion content determined by uniaxial compression test of rock-like blocks with locked-in stress.

5.3. Theoretical Analysis of Peak Strength Considering Confining Pressure

When considering confining pressure, there are$$\begin{array}{}\text{(79)}& {\sigma }_{\mathrm{n}}^{\prime \prime \prime }=\frac{{\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}+\frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}\cos \mathrm{2}\beta -\frac{{\sigma }^L{\rho }_s}{\mathrm{1}-{\rho }_s}\text{(80)}& {\tau }^{\prime \prime \prime }=\frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}\sin \mathrm{2}\beta \end{array}$$

Substituting (79) and (80) into (60) yields$$\begin{array}{}\text{(81)}& \frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}\sin \mathrm{2}\beta =c+\left(\frac{{\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}+\frac{{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}}{\mathrm{2}\left(\mathrm{1}-{\rho }_s\right)}\cos \mathrm{2}\beta -\frac{{\sigma }^L{\rho }_s}{\mathrm{1}-{\rho }_s}\right)tg\phi \end{array}$$

Equation (81) can be converted into$$\begin{array}{}\text{(82)}& {\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}={\sigma }_{\mathrm{0}}\left(\mathrm{1}-{\rho }_s+\frac{tg\phi {\sigma }_{\mathrm{3}}}{c}-\frac{{\sigma }^L{\rho }_{s}tg\phi }{c}\right)\end{array}$$

By observing (82), it is found that when the other parameters are constant, the confining pressure is proportional to the change of the deviatoric stress intensity, which is because that the equation deduced based on the Mohr-Coulomb equation (60) was actually a linear equation, which does not match the actual situation. The relationship between confining pressure and peak deviatoric stress is often a nonlinear relationship in the parabolic or hyperbolic form. Therefore, it is necessary to adjust the confining pressure item in (82). Parameter $h$ is introduced combined with parameters in (78) to get the final expression of peak deviatoric stress:$$\begin{array}{}\text{(83)}& {\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}={\sigma }_{\mathrm{0}}\frac{\left(\mathrm{1}-{\rho }_s+tg\phi ·h·{\sigma }_{\mathrm{3}}/c-{\sigma }^L{\rho }_{s}tg\phi /qc\right)}{\mathrm{1}+d\rho }.\end{array}$$where $h$, $q$, and $d$ are parameters related to confining pressure and can be obtained by fitting the test data.

When the confining pressure is 0, (83) degenerates to (77).

5.4. Test Verification

( 1) Uniaxial Compression without Considering Locked-In Stress. Using (69) to fit the test data of Geng [16] when there is no locked-in stress, $d=1.72\hspace{0.17em}\hspace{0.17em}(R^2=0.9950)$ can be obtained, as shown in Figure 15.

[figure omitted; refer to PDF]

So when the inclusion does not contain a locked-in stress, the expression of the strength of the test block is$$\begin{array}{}\text{(84)}& {\sigma }_c={\sigma }_{\mathrm{0}}\frac{\left(\mathrm{1}-\pi {\left(\mathrm{3}\rho /\mathrm{4}\pi \right)}^{\mathrm{2}/\mathrm{3}}\right)}{\mathrm{1}+\mathrm{1.72}\rho }.\end{array}$$

( 2) Uniaxial Compression with Considering Locked-In Stress. Figure 16 shows the uniaxial and conventional triaxial compressive stress-strain curves of a standard sample of cement mortar without inclusions obtained by laboratory tests. The confining pressure and its corresponding peak strength in the figure are listed in Table 1.

Table 1

The strength of rock-like materials under different confining pressures.

Confining pressure/MPa	0	2	4	6
Peak strength/MPa	9.91	17.92	24.56	28.61

[figures omitted; refer to PDF]

Using the test data in Table 1, the cohesion c and the friction angle $\phi$ of the cement mortar matrix can be obtained by linear fitting according to Mohr-Coulomb equation (85). The fitting curve is shown in Figure 16(b) in the solid line.$$\begin{array}{}\text{(85)}& {\sigma }_{\mathrm{1}}=\frac{\mathrm{1}+\sin \phi }{\mathrm{1}-\sin \phi }{\sigma }_{\mathrm{3}}+\frac{\mathrm{2}c\cos \phi }{\mathrm{1}-\sin \phi }\end{array}$$

The cohesive force of the matrix material is $c=2.7$1MPa; the friction angle is $\phi ={\mathrm{32.60}}^{°}$ by fitting.

Equation (77) is used to fit the strength test data when there is locked-in stress; it is found that when $q$ meets $q=w+l\rho$, (77) can be consistent with the experimental data, and the fitting result is shown in Figure 17, with $l=5.50$, $w=0.018$ ($R^2=0.9946$).

[figure omitted; refer to PDF]

Therefore, the strength expression of the specimen containing locked-in stress inclusions is$$\begin{array}{}\text{(86)}& \sigma ={\sigma }_{\mathrm{0}}\frac{\left(\mathrm{1}-\pi \left(\mathrm{1}+\left(tg\phi ·{\sigma }^L\right)/\left(\mathrm{0.018}+\mathrm{5.5}\rho \right)c\right){\left(\mathrm{3}\rho /\mathrm{4}\pi \right)}^{\mathrm{2}/\mathrm{3}}\right)}{\mathrm{1}+\mathrm{1.72}\rho }\end{array}$$

When ${\sigma }_L$ in (86) is 0, it degenerates to (69).

It can be seen that the inclusion content and the value of locked-in stress are significant for the change of the peak stress and elastic modulus of rock-like materials, and with the increase of locked-in stress, the effect on elastic modulus and peak stress will also grow bigger.

(3 ) Conventional Triaxial Compression Considering Locked-In Stress. From (83), firstly determine the parameter $h$ based on the experimental values of the pure matrix material under different confining pressures. Then on the basis of the test values without considering the locked-in stress under each confining pressure, the parameter $b$ is determined by fitting. On this basis, the parameters $w$, $l$ are determined based on all the data under the corresponding confining pressure. See Table 2 for various parameters under different confining pressures. The result of the fitting is shown in Figure 18.

Table 2

Theoretical results considering confining pressure.

Confining pressure/MPa	parameter
	$$h$$	$$d$$	$$w$$	$$l$$
0	Any value	1.2	0.018	5.50
2	1.21	2.60	0.013	4.2
4	1.14	6.65	0.01	3.44
6	1.00	10.49	0.005	3.38

[figures omitted; refer to PDF]

From Table 2, it can be seen that the parameters $h$, $d$, $w$, $l$ of (83) change with the confining pressure. Among them, $h$ decreases with the increase of confining pressure; $d$ increases with the increase of confining pressure; $w$ decreases with the increase of confining pressure; $l$ decreases with the increase of confining pressure and then increases. The fitting results of these parameters with the confining pressure are shown in (87). The fitting curve is shown in Figure 19.$$\begin{array}{c}\text{(87)}& h=-\mathrm{0.0525}{\sigma }_{\mathrm{3}}+\mathrm{1.33}d=\mathrm{1.2}{\sigma }_{\mathrm{3}}+\mathrm{1.765}\\ w=-\mathrm{0.0021}{\sigma }_{\mathrm{3}}+\mathrm{0.0178}\\ l=\mathrm{0.0775}{{\sigma }_{\mathrm{3}}}^{\mathrm{2}}-\mathrm{0.821}{\sigma }_{\mathrm{3}}+\mathrm{5.5}\end{array}$$

[figures omitted; refer to PDF]

Substituting the value in (87) into (83) yields the strength expression of the specimen containing the locked-in stress when considering confining pressure.

6. Conclusion

The change mechanism of the elastic modulus and peak strength of the sample along with the variation of the inclusion content and the locked-in stress were analyzed, and the following results were obtained.

(1) According to Eshelby’s equivalent inclusion theory, the variation law of locked-in stress with external load in rock was deduced. Results showed that, with the increase of the external load, the locked-in stress would also increase; the increase range was related to the ratio of the elastic modulus of the inclusion to the matrix, the ratio of Poisson’s ratio of the inclusion and the matrix, and the content of inclusions. When the inclusions and the matrix had the same parameters except the elastic modulus, if the ratio of the elastic modulus of the inclusions to the matrix is greater than 1, the locked-in stress decreases with the increase of the content; while the ratio of the elastic modulus of the inclusion to the matrix is less than 1, the locked-in stress increases with the increase of the content.

(2) Composite micromechanics theory was introduced. The relationship between the elastic modulus of the rock and the content of inclusions without considering the locked-in stress was obtained, based on Voigt upper limit and Reuss lower limit method, Taylor method, and M-T method, respectively. Through verification and comparison with the previous test results, the M-T method was used to describe the equivalent elastic modulus of rock-like material containing inclusions.

(3) The concept of strain energy and damage was introduced, and the relationship between strain energy and damage degree were established through a statistical probability model. Based on the M-T method, the theoretical solution of the equivalent elastic modulus of rock with locked-in stress is deduced. Through the three-dimensional nonlinear fitting, the empirical formula of the change law of the elastic modulus of the sample with the content of the inclusion and the locked-in stress value is obtained, which is in good agreement with the test result.

(4) Mohr-Coulomb strength theory was introduced, and the relationship between the peak stress and the inclusion content of the test block without considering the locked-in stress was obtained by combining the effective stress principle. The empirical formula of the rule that the peak stress of the test block containing locked-in stress inclusions varies with the inclusion content and the locked-in stress value was obtained by the three-dimensional nonlinear fitting method.

(5) The concept of effective stress was extended to the conventional three-axis case. Based on the Mohr-Coulomb model, the empirical formula of the variation of the peak deviatoric stress of the test block containing locked-in stress under different confining pressures with the inclusion body content and the locked-in stress value is obtained. The parameters corresponding to the empirical formulas for different confining pressures were obtained by the test data. Through polynomial fitting, the relationship between these parameters and confining pressure was obtained.

The spherical rubber particles used to simulate stress inclusions can only represent a certain type of stress inclusions. There are many types of stress inclusions in rocks. To conduct systematic research, more advanced instruments are needed to conduct a large number of tests.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank Key Laboratory of Ministry of Education for Geomechanics and Geomechanics and Embankment Engineering, Hohai University (no. 201707) and Natural Science Foundation of Jiangsu Province of China (BK20141067) for their financial support.

References

[1] G. Z. Voyiadjis, P. D. Kiousis, C. S. Hartley, "Analysis of residual stresses in cylindrically anisotropic materials," Experimental Mechanics, vol. 25 no. 2, pp. 145-147, DOI: 10.1007/BF02328804, 1985.

[2] G. Z. Voyiadjis, C. S. Hartley, "Residual-stress determination of concentric layers of cylindrically orthotropic materials," Experimental Mechanics, vol. 27 no. 3, pp. 290-297, DOI: 10.1007/BF02318095, 1987.

[3] X. Chen, H. Xu, H. Geng, L. Dong, J. Zhang, "A New Equivalent Statistical Damage Constitutive Model on Rock Block Mixed Up with Fluid Inclusions," Mathematical Problems in Engineering, vol. 2018,DOI: 10.1155/2018/3080173, 2018.

[4] Y. M. Chen, "Study on damage mechanics of effect of porosity on sandstone strength," Sichuan Hydropower, vol. no. 2, pp. 39-43, 1997.

[5] W. Wang, Z. Y. Tain, Q. Z. Zhu, "Study of statistical damage constitutive model for rock considering pore water pressure," Chinese Journal of Rock Mechanics and Engneering, vol. no. s2, pp. 3676-3682, 2015.

[6] Y. Lai, S. Li, J. Qi, "Strength distributions of warm frozen clay and its stochastic damage constitutive model," Cold Regions Science & Technology, vol. 53 no. 2, pp. 200-215, DOI: 10.1016/j.coldregions.2007.11.001, 2008.

[7] S. Li, Y. Lai, S. Zhang, "An improved statistical damage constitutive model for warm frozen clay based on Mohr–Coulomb criterion," Cold Regions Science & Technology, vol. 57 no. No, pp. 154-159, DOI: 10.1016/j.coldregions.2009.02.010, 2009.

[8] Q. C. Zhou, "Study on the mechanical property of a sandstone under geothermal-mechanical and hydraulic-mechanical coupling," Chinese academy of science (Institute of rock & soil mechanics), vol. 2006, 2006.

[9] Z. J. Ren, A Micromechanical damage model for rocks and concretes based on the deformation and growth of elliptic microcracks, 2008.

[10] H. W. Yang, Study on coupling mechanism of rock and pore water under cyclic loading, 2011.

[11] J. C. Jaeger, N. G. W. Cook, R. W. Zimmerman, Fundamentals of rock mechanics[M]., Maklen, 2007.

[12] J. D. Eshelby, "The determination of the elastic field of an ellipsoidal inclusion, and related problems," Proceedings of the Royal Society London A: Mathematical, Physical and Engineering Sciences, vol. 241, pp. 376-396, DOI: 10.1098/rspa.1957.0133, 1957.

[13] T. Mori, K. Tanaka, "Average stress in matrix and average elastic energy of materials with misfitting inclusions," Acta Metallurgica et Materialia, vol. 21 no. 5, pp. 571-574, DOI: 10.1016/0001-6160(73)90064-3, 1973.

[14] G. L. Shen, G. K. Hu, B. Liu, Mechanics of Composite Materials, 2013.

[15] Z. L. Xu, A brief tutorial on elastic mechanics, 2002.

[16] H. S. Geng, Influence of locked-in stress on compressive properties of rock, 2018.

[17] S. Y. Du, Micromechanics of composite materials, 1998.

[18] X. J. Liu, J. K. Shen, L. X. Liang, al. et, "Effects of pore pressure on rock strength properties," Chinese Journal of Rock Mechanics and Engneering, vol. 30 no. s2, pp. 3457-3463, 2011.

[19] W. Y. Li, H. S. Liu, "Influence of Pore Water Pressure on Shear Modulus and Damping Ratio of Saturated Sands," Chinese Journal of Geotechnical Engineering, vol. 5 no. 4, pp. 56-67, 1983.

[20] X. Li, W.-G. Cao, Y.-H. Su, "A statistical damage constitutive model for softening behavior of rocks," Engineering Geology, vol. 143-144,DOI: 10.1016/j.enggeo.2012.05.005, 2012.

[21] J. A. Lemaitre, "How to use damage mechanics," Nuclear Engineering and Design, vol. 80 no. 2, pp. 233-245, DOI: 10.1016/0029-5493(84)90169-9, 1984.