Development of PMC Numerical Model for Soil and Stone Mixture

Abstract

Soil and stone mixture is widely distributed in China, and its physical and mechanical properties are complex, which has a significant impact on geotechnical engineering. Usually soil and stone mixture shows anisotropic features along or perpendicular to the direction of settled layers, and the strength will be significantly affected by intermediate principal stress. True triaxial tests were carried on, Paul–Mohr–Coulomb (PMC) failure criterion was used for the strength analysis in soil-stone mixture, and related parameters of PMC model were obtained. A user-defined PMC numerical constitutive model was developed for FLAC3D. Composite failure criteria of shear failure and tension cut-off were applied for numerical analysis, and nonassociated flow rule was proposed based on the Mohr–Coulomb model. Verification modelling was applied as well, and deviation between analytical and numerical solutions in strains of X direction, Y direction, and Z direction was 1.0013, 1.0003, and 1.002, respectively.

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

Soil and stone mixture is widely distributed in China, which has a great impact on geotechnical engineering. Due to the river scouring and accumulation settlement, the gravels or stones in the soil-stone mixture are arranged in a certain regularity, resulting that the physical and mechanical properties are different from the isotropic soil material. Soil-stone mixture has strong anisotropy in the layer strike direction and vertical direction. Meanwhile, its triaxial strength characteristics and failure characteristics are significantly affected by the intermediate principal stress [1–3].

At present, a lot of experiments and theoretical research on soil-stone mixture were carried out, and achievements in the aspects of mechanical properties and strength criteria of soil-stone mixture were obtained. In the experimental research, Wu Ming carried out 4 groups of tests of soil-stone mixture and the results showed that the common feature of specimens composed of coarse and fine materials depended on the content of components [4]. Jing et al. carried out the large-scale triaxial shear test of two mass ratios of soil-rock mixture under four different confining pressures, respectively. The results showed that the change of soil-rock mixture ratio has great influence on the dilatancy of soil-rock mixture [5]. Yu et al. also carried out the shear test of soil-rock mixture, and the test results show that the shear strength of soil-rock mixture increases with the increase in stone content, and the shear strength increase is more obvious when the stone content exceeds 40% [6]. Minshuo et al. [7] carried out in situ shear tests. The results showed typical characteristics of complete stress-strain curves, and the stone content is an important factor affecting the strength and failure form of soil and stone mixture.

Wang et al. carried out a series of numerical uniaxial compression tests to obtain the deformation and failure characteristics of the soil-rock mixture with different compositions and structures. Obvious bedding phenomena usually appeared in the direction of rock block inclination, and the ultimate shear strain increases with the rock block proportion increase [8]. Zhao and Liu carried out a triaxial compression simulation test of the soil-rock mixture by using PFC and analyzed influences of rock content, mesoscopic features, and random distribution on macroscopic shear strength characteristics. The results showed that, with the increase in rock content, the internal friction angle of soil-rock mixture increased continuously, while the cohesion decreased firstly and then increased [9]. Jin and Zeng conducted numerical simulation of large-scale triaxial tests on soil-rock mixture using DEM with three-dimensional flexible membrane boundary [10]. Zhan and You conducted the numerical simulation of the large direct shear test of soil-rock mixture, and the rock content significantly affected the mechanical properties of soil-rock mixture. The results showed that properties of soil-rock mixture which are obviously different from those of the soil are about at 40% rock content [11].

In terms of constitutive model research, Vallejo and Lobo-Guerrero [12, 13] used the Hashin formula to predict the elastic parameters of clay with large particle size. By comparing with the test results, it can be found that the strength of soil-stone mixture can be predicted according to the strength of soil matrix and Guth formula.

Lee and Pyo [14] proposed a multilevel elastic damage model based on the combination of micromechanics formula and multilevel damage model to predict the effective elastic behavior and interface gradual weakening of granular composites.

Ma et al. [15] and others derived the equivalent elastic modulus formula of elliptical pebbles with horizontal arrangement of long axis by using the parallel model. The formula was based on the assumption of macroscopic isotropy and strain uniformity. In the study of intermediate principal stress, Haimson and Rudnicki [16] believed that the intermediate principal stress affects 15% of the strength.

Xiaochun and Dongjun [17] found that the influence of intermediate principal stress on the strength would be enhanced under higher confining pressure and in case of hard rock. The research on triaxial strength characteristics and failure characteristics of soil-stone mixture is still needed to reveal and explore. There is a lack of an elastic-plastic constitutive model which can reflect the effect of intermediate principal stress in practical engineering.

Paul [18, 19] proposed the Paul–Mohr–Coulomb failure criterion based on the Mohr–Coulomb model, and it was implemented by Meyer and Labuz [20]. The PMC failure criterion reflects the influence of intermediate principal stress on materials by presenting the difference of internal friction angle in compression and tension [21].

In this paper, the PMC failure criterion is carried out to fit the result of the triaxial test and to research the strength and failure rules. Material parameters of soil and stone mixture under the PMC failure criterion were determined, a user-defined FLAC3D numerical constitutive model based on the PMC strength criterion was developed, and verification analysis was involved as well.

2. Applicability Analysis of PMC Strength Criterion in Soil-Stone Mixture Material

2.1. PMC Failure Criterion

The PMC failure criterion formula is as follows:$$\begin{array}{}\text{(1)}& A{\sigma }_1+B{\sigma }_2+C{\sigma }_3=1.\end{array}$$

In equation (1), A, B and C are: $$\begin{array}{}\text{(2)}& \begin{array}{c}A=\frac{1-\sin {\varphi }_c}{2\sin {\varphi }_c}\cdot \frac{1}{V_0},\\ B=\frac{\sin {\varphi }_c-\sin {\varphi }_e}{2\sin {\varphi }_c\sin {\varphi }_e}\cdot \frac{1}{V_0},\\ C=-\frac{1+\sin {\varphi }_e}{2\sin {\varphi }_e}\cdot \frac{1}{V_0},\end{array}\end{array}$$where φ_c and φ_e are internal friction angles in compression and tension, respectively, and V₀ is the theoretical value of triaxial tensile strength, as shown in Figure 1.

[figure omitted; refer to PDF]

Zeng et al. [22] proposed a new fitting method based on the stress invariant relationship between the principal stress space and the $p-q$ plane to determine the parameters φ_c, φ_e, and V₀ of geotechnical materials in the PMC failure criterion.

The relationship between three-dimensional principal stress and mean stress, shear stress, and Lode angle is as follows:$$\begin{array}{c}\text{(3)}& p=\frac{{\sigma }_1+{\sigma }_2+{\sigma }_3}{3},q=\frac{1}{\sqrt{2}}\sqrt{{{\sigma }_1-{\sigma }_2}^2+{{\sigma }_2-{\sigma }_3}^2+{{\sigma }_3-{\sigma }_1}^2},\\ \tan \theta =\frac{2{\sigma }_2-{\sigma }_1-{\sigma }_3}{\sqrt{3}{\sigma }_1-{\sigma }_3},\end{array}$$where $p$ is the average stress, q is the generalized shear stress, and θ is the Rhode angle (ranging from 0° to 60°).

Based on the new fitting method, in the triaxial compression test (θ = 0°), the PMC failure criterion in the $p-q$ plane can be expressed as follows:$$\begin{array}{}\text{(4)}& q=\frac{6\sin {\varphi }_c}{3-\sin {\varphi }_c}p+\frac{6c_c\sin {\varphi }_c}{3-\sin {\varphi }_c}.\end{array}$$

In the triaxial compression test (θ = 60°), the PMC strength criterion in the $p-q$ plane can be expressed as follows:$$\begin{array}{}\text{(5)}& q=\frac{6\sin {\varphi }_e}{3-\sin {\varphi }_e}p+\frac{6c_e\sin {\varphi }_e}{3-\sin {\varphi }_e},\end{array}$$where c_c and c_e are compressive cohesion and tensile cohesion, respectively, as follows:$$\begin{array}{}\text{(6)}& \begin{array}{c}{c}_c=V_0\tan {\varphi }_c,\\ c_e=V_0\tan {\varphi }_e.\end{array}\end{array}$$

When φ_c = φ_e and c_c = c_e, φ_c in the PMC model is the internal friction angle φ and c_c is the cohesive force c.

2.2. True Triaxial Test Analysis of Soil-Stone Mixture

True triaxial testing of soil-stone mixture materials was applied on the machine of TSW-50 (see Figure 2). Analysis by the PMC criterion was carried on, and strength parameters were obtained. The loading was adjusted by changing the ratio coefficient b of medium principal stress as follows:$$\begin{array}{}\text{(7)}& b=\frac{{\sigma }_2-{\sigma }_3}{{\sigma }_1-{\sigma }_3}.\end{array}$$

[figure omitted; refer to PDF]

There were 8 groups of experiments involved, 1∼4 groups were soil-stone mixture experiments with medium principal stress ratio coefficient of 0 and 5∼8 groups were experiments with the coefficient of 1. When the intermediate principal stress coefficient b = 0, the specimens were first consolidated at $p_0$ = 50, 100, 150, and 200 kPa, then kept σ₃ unchanged, and σ₁ gradually increased until the failure occurred. When the intermediate principal stress coefficient b = 1, the same consolidation processes were operated, then σ₃ remained unchanged, and gradually increased σ₁ to the failure strength. The curve of stress deviation vs. the axial principal strain (σ₁–σ₃∼ε₁) is shown in Figure 3.

[figures omitted; refer to PDF]

The results show that with the increase in stress level, the triaxial failure strengths are obtained. The failure surface and test failure data are fitted in $p$-q plane, as shown in Figure 4.

[figure omitted; refer to PDF]

The correlation coefficient R² of soil-stone mixture test data is 0.92. The true triaxial test results of soil-stone mixture fit the PMC criterion well, and the trend is consistent. The compression friction angle is smaller than the tensile internal friction angle. The difference between φ_c and φ_e (see Table 1) can evaluate the influence of intermediate principal.

Table 1

PMC material parameters of soil-stone mixture.

Geotechnical materials	φc (°)	φe (°)	V0 (kPa)	R2
Soil-stone mixture	38.1	42.7	6.4	0.92

Through the applicability analysis of the PMC strength criterion in soil-stone mixture material, the fitted PMC failure surface is consistent with the test data, which are in line with the expected trend of soil-stone mixture material. The PMC strength criterion can be used to describe the strength and failure behavior of soil-stone mixture material.

3. Development of PMC Numerical Constitutive Model

3.1. Elasto-Plastic Calculation Equation

In general, Hooke’s Law is used to describe the linear elastic behavior of the PMC model as follows:$$\begin{array}{}\text{(8)}& \begin{array}{c}{\epsilon }_x=\frac{1}{E}{\sigma }_x-\mu {\sigma }_y+{\sigma }_z,\\ {\epsilon }_y=\frac{1}{E}{\sigma }_y-\mu {\sigma }_x+{\sigma }_z,\\ {\epsilon }_z=\frac{1}{E}{\sigma }_z-\mu {\sigma }_x+{\sigma }_y.\end{array}\end{array}$$

In the plastic stage, the material constitutive in FLAC3D software follows the incremental forms. The specific calculation steps are as follows:

(1) Strain increments are divided into elastic strains and plastic strains as follows:$$\begin{array}{}\text{(9)}& \Delta {\epsilon }_i=\Delta {\epsilon }_i^e+\Delta {\epsilon }_i^p.\end{array}$$

(2) The stress-strain constitutive equation for flow rule is as follows:$$\begin{array}{}\text{(10)}& \Delta {\epsilon }_i^p=\lambda \frac{\partial g}{\partial {\sigma }_i},\hspace{1em}i=\mathrm{1,2,3},\end{array}$$

  where $g$ is the plastic potential function.

(3) The stress increment can be determined by the following equation:$$\begin{array}{}\text{(11)}& \Delta {\sigma }_i=S_i\Delta {\epsilon }_n^e,\hspace{1em}i=1\dots n,\end{array}$$

  where S_i is the linear equation which obeys Hooke’s Law.

(4) The new stress state also satisfies the yield equation as follows:$$\begin{array}{}\text{(12)}& f{\underset{¯}{\sigma }}_n+\Delta {\underset{¯}{\sigma }}_n=0,\end{array}$$

  where $\mathrm{\Delta }{\underset{¯}{\sigma }}_i=S_i\mathrm{\Delta }{\underset{¯}{\epsilon }}_n-S_i\mathrm{\Delta }{\underset{¯}{\epsilon }}_n^p$, and then it can be deducted that$$\begin{array}{}\text{(13)}& f^{\ast }{S}_n\Delta {\underset{¯}{\epsilon }}_n-\lambda f^{\ast }{S}_n\frac{\partial g}{\partial {\underset{¯}{\sigma }}_n}=f{\underset{¯}{0}}_n.\end{array}$$

The expression of coefficient λ can be obtained as follows:$$\begin{array}{}\text{(14)}& \lambda =\frac{f^{\ast }{S}_n\Delta {\underset{¯}{\epsilon }}_n}{f^{\ast }{S}_n\partial g/\partial {\underset{¯}{\sigma }}_n-f{\underset{¯}{0}}_n}.\end{array}$$

3.2. Composite Failure Criterion

The failure criterion of the PMC model is a composite failure criterion, as shown in Figure 5. Considering three principal stresses, let$$\begin{array}{}\text{(15)}& {\sigma }_1\ge {\sigma }_2\ge {\sigma }_3.\end{array}$$

[figure omitted; refer to PDF]

To neglect hardening of materials, the strain increment can be divided into elastic and plastic strain increment as follows:$$\begin{array}{c}\text{(16)}& \Delta \epsilon =\Delta {\epsilon }^e+\Delta {\epsilon }^p,\Delta {\epsilon }_i^p=\lambda \frac{\partial g}{\partial {\sigma }_i},\hspace{1em}i=\mathrm{1,2,3},\end{array}$$where λ is a constant and $g$ (σ₁, σ₂, σ₃) is the plastic potential function.

The plastic potential functions $g^s$ (σ₁, σ₂, σ₃) and $g^t$ (σ₃) of the PMC model are chosen for shear failure and tension failure, which can be derived according to the Mohr–Coulomb model of FLAC3D.

3.2.1. Shear Failure

 PMC criteria can be expressed as follows:$$\begin{array}{}\text{(17)}& \frac{{\sigma }_1}{V_0}\frac{1-\sin {\phi }_c}{2\sin {\phi }_c}+\frac{{\sigma }_2}{V_0}\frac{\sin {\phi }_c-\sin {\phi }_e}{2\sin {\phi }_c\sin {\phi }_e}-\frac{{\sigma }_3}{V_0}\frac{1+\sin {\phi }_e}{2\sin {\phi }_e}=1,\end{array}$$where σ₁, σ₂, and σ₃ are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively; ${\varphi }_c$ and ${\varphi }_e$ are the internal friction angle in compression and tension cases, respectively; and V₀ is the theoretical triaxial tensile strength.

 Let$$\begin{array}{}\text{(18)}& \begin{array}{c}{N}_c=\frac{1-\sin {\varphi }_c}{2\sin {\varphi }_c},\\ N_e=\frac{1-\sin {\varphi }_e}{2\sin {\varphi }_e}.\end{array}\end{array}$$

Then,$$\begin{array}{}\text{(19)}& f^s={\sigma }_1{N}_c+{\sigma }_2{N}_c-N_e-{\sigma }_3{N}_e+1-V_0=0.\end{array}$$

And the plastic potential function is as follows:$$\begin{array}{}\text{(20)}& g^s={\sigma }_1{N}_c+{\sigma }_2{N}_c-N_e-{\sigma }_3{N}_e+1.\end{array}$$

3.2.2. Tension Cut-Off

If f^t = 0 is used as the tension strength, then$$\begin{array}{}\text{(21)}& f^t={\sigma }_3-{\sigma }^t=0.\end{array}$$

Then, the plastic potential function is as follows:$$\begin{array}{}\text{(22)}& g^t={\sigma }_3.\end{array}$$

3.3. Elastic-to-Plastic Stress Analysis

The stress of material according to the elastic-to-plastic model defined by its elastic component can be written as follows:$$\begin{array}{}\text{(23)}& {\sigma }_i^N={\sigma }_i^I-\lambda S_i\frac{\partial g}{\partial {\sigma }_n}.\end{array}$$

For shear yield correction, ${\sigma }_i^I$ is as follows:$$\begin{array}{}\text{(24)}& \begin{array}{c}\Delta {\sigma }_1={\alpha }_1\Delta {\epsilon }_1^e+{\alpha }_2\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e,\\ \Delta {\sigma }_2={\alpha }_1\Delta {\epsilon }_2^e+{\alpha }_2\Delta {\epsilon }_1^e+\Delta {\epsilon }_3^e,\\ \Delta {\sigma }_3={\alpha }_1\Delta {\epsilon }_3^e+{\alpha }_2\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e,\end{array}\end{array}$$where α₁ and α₂ are the material parameters defined based on shear modulus G and bulk modulus K:$$\begin{array}{}\text{(25)}& \begin{array}{c}{\alpha }_1=K+\frac{4}{3}G,\\ {\alpha }_2=K-\frac{2}{3}G.\end{array}\end{array}$$

The yield function of PMC $f{\sigma }_1,{\sigma }_2,{\sigma }_3$ can be derived as follows:$$\begin{array}{}\text{(26)}& \lambda =\frac{f{\sigma }_1^I,{\sigma }_2^I,{\sigma }_3^I}{f{S}_1\partial g/\partial {\sigma }_i,S_2\partial g/\partial {\sigma }_i,S_3\partial g/\partial {\sigma }_i},\hspace{1em}i=\text{I},\text{II},\text{III},\end{array}$$where ${\sigma }_n^{I}n=\mathrm{1,2,3}$ is the elastic stress in FLAC3D and f(0_n) is the constant term of the linear failure criterion.

 In equation (27), the expression $S_i\partial g/\partial {\sigma }_n$ is as follows:$$\begin{array}{}\text{(27)}& \begin{array}{c}{S}_1\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},={\alpha }_1\Delta {\epsilon }_1^e+{\alpha }_2\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e,\\ S_2\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},={\alpha }_1\Delta {\epsilon }_2^e+{\alpha }_2\Delta {\epsilon }_1^e+\Delta {\epsilon }_3^e,\\ S_3\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},={\alpha }_1\Delta {\epsilon }_3^e+{\alpha }_2\Delta {\epsilon }_1^e+\Delta {\epsilon }_2^e.\end{array}\end{array}$$

 After the material enters the plastic state, the flow criterion is as follows:$$\begin{array}{}\text{(28)}& \begin{array}{c}\frac{\partial g^s}{\partial {\sigma }_1}=1,\\ \frac{\partial g^s}{\partial {\sigma }_2}=0,\\ \frac{\partial g^s}{\partial {\sigma }_3}=-N_{\psi },\end{array}\end{array}$$where ψ is the dilatancy angle and$$\begin{array}{}\text{(29)}& N_{\psi }=\frac{1+\sin \psi }{1-\sin \psi }.\end{array}$$

 Then,$$\begin{array}{}\text{(30)}& \begin{array}{c}{S}_1\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},={\alpha }_1-{\alpha }_2{N}_{\psi },\\ S_2\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},={\alpha }_{2}1-N_{\psi },\\ S_3\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3},=-{\alpha }_1{N}_{\psi }+{\alpha }_2.\end{array}\end{array}$$

 When the model does not consider the dilatancy, if ψ = 0, then N_ψ = 1. The plastic volume strain increment is 0.$$\begin{array}{c}\text{(31)}& \Delta {\epsilon }_v^p={\displaystyle \sum \Delta {\epsilon }_i^p}=\lambda \frac{\partial g}{\partial {\sigma }_i}\\ =0.\end{array}$$

3.3.1. For Tension Failure

For tension failure, ${\sigma }_i^I$ is as follows:$$\begin{array}{}\text{(32)}& \begin{array}{c}\Delta {\sigma }_1={\alpha }_1\Delta {\epsilon }_1^p+{\alpha }_2\Delta {\epsilon }_2^e+\Delta {\epsilon }_3^e,\\ \Delta {\sigma }_2={\alpha }_1\Delta {\epsilon }_2^p+{\alpha }_2\Delta {\epsilon }_1^e+\Delta {\epsilon }_3^e,\\ \Delta {\sigma }_3={\alpha }_1\Delta {\epsilon }_3^p+{\alpha }_2\Delta {\epsilon }_1^e+\Delta {\epsilon }_2^e,\end{array}\end{array}$$

λ would be$$\begin{array}{}\text{(33)}& \lambda =\frac{{\sigma }_3^I-{\sigma }^{\prime }}{{\alpha }_1},\end{array}$$where ${\sigma }_n^I$ (n = 1, 2, 3) is the elastic stress in FLAC3D and f(0_n) is the constant term of the linear failure criterion.

$S_i\partial g/\partial {\sigma }_n$ is as follows:$$\begin{array}{}\text{(34)}& \begin{array}{c}{S}_1\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}={\alpha }_1,\\ S_2\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}={\alpha }_2,\\ S_3\frac{\partial g^s}{\partial {\sigma }_1},\frac{\partial g^s}{\partial {\sigma }_2},\frac{\partial g^s}{\partial {\sigma }_3}={\alpha }_3.\end{array}\end{array}$$

3.4. Development of Constitutive Elastic-to-Plastic Model for PMC

Based on the software of “Fast Lagrangian Analysis of Continua in 3 Dimensions” (FLAC3D, Itasca Consulting Group Inc.), a user-defined PMC model is developed.

The parameters and corresponding keywords of the PMC model defined in this paper are shown in Table 2.

Table 2

PMC model parameters and keywords corresponding table.

Model keywords	Model parameters
Bulk	Bulk modulus, K
Shear	Shear modulus, G
c_friction	Internal friction angle in compression, φc
e_friction	Internal friction angle in tension, φe
v_para	Theoretical triaxial tensile strength, V0
Tension	Ensile strength, σt
Dilation	Dilation, $\psi$

When calling the custom PMC numerical constitutive model for numerical model calculation, it is necessary to load the PMC constitutive model into the FLAC3D. The validity of the numerical constitutive model of PMC was evaluated by modelling in the true triaxial test process.

4. Verification of PMC Numerical Constitutive Model

4.1. Development of Constitutive Elastic-to-Plastic Model for PMC

The PMC failure criterion is derived from the M-C strength criterion. When ${\varphi }_c={\varphi }_e$, the PMC model is the M-C model. In the numerical simulation test, different strength can be obtained by matching different model parameters. Three conditions are simulated by the true triaxial numerical simulation test. Condition 1 is the calculation by the M-C model; condition 2 by the PMC model with parameters ${\varphi }_c={\varphi }_e$ (denoted as PMC-1); and condition 3 by the PMC model with parameters ${\varphi }_c\ne {\varphi }_e$ (denoted as PMC-2). See Table 3 for parameters of each case. Comparing the simulation results of the PMC model and M-C model is helpful to verify the customized PMC model.

Table 3

Parameter table of the triaxial compression test model.

Working conditions	E (GPa)	μ	φc (°)	φe (°)	cc (MPa)	ce (MPa)	V0 (MPa)	${\sigma }_{\text{ci}}$ (MPa)
M-C	10.0	0.33	28.0	—	8.5	—	16.0	28.29
PMC-1	10.0	0.33	28.0	28.0	8.5	8.5	16.0	28.29
PMC-2	10.0	0.33	28.0	33.0	8.5	10.4	16.0	28.29

Note. σ_ci is the uniaxial compressive strength.

A cube element (size: 1 × 1 × 1) is created, and the fixed boundary surface x = 0, y = 0, and z = 0 are applied. 5 MPa and 10 MPa pressures are applied to the element in X and Y direction, respectively. The displacement control method is used to apply pressure in Z direction until yielding. The vertical stress and strain are monitored in the calculation process (Figure 6).

[figure omitted; refer to PDF]

To ensure the quasi-static loading, the loading speed is set at 10⁻⁸ per step, and a total of 450,000 steps are calculated. The MC model and PMC-1 model entered the plastic stage in 372,000 steps, and the PMC-2 model entered the plastic stage in 385,000 steps. After 450,000 steps of calculation, the strain in Z direction reaches 4.5 × 10⁻³. At this time, the triaxial test blocks under each constitutive model have entered the plastic state. By analyzing the strength peak value, elastic strain, and plastic strain of the test block, the user-defined PMC model can be checked and the failure characteristics of geotechnical materials can be explored.

4.2. Analysis of True Triaxial Numerical Simulation Results

Figure 7 shows the vertical stress curve of the PMC-2 model. With the increase in step length, point A is the peak strength, and the elastic strain (point A) can be obtained as ε_z|_A = 3851.00 × 10⁻⁶. The plastic strain is obtained by subtracting the elastic strain from the total strain, and the plastic deformation at point B is ε_z|_B = 654.13 × 10⁻⁶.

[figure omitted; refer to PDF]

Therefore, it can be verified that the numerical constitutive model of self-defined PMC can be called normally, and the calculation results in special cases are consistent with those of the M-C model in special case.

Table 4 shows the comparison of yield strength and strain results of the true triaxial numerical simulation test of the M-C strength criterion and PMC model when ${\varphi }_c={\varphi }_e$. Case 1 calls the built-in M-C numerical constitutive model of FLAC3D software, and case 2 calls the user-defined PMC numerical constitutive model and sets it as the special parameter case of ${\varphi }_c={\varphi }_e$. At this time, the PMC model regresses to the M-C model.

Table 4

Comparison table of numerical simulation results of the true triaxial test.

Working conditions		M-C	PMC-1	PMC-1/M-C
Peak strength (MPa)		42.14	42.17	1.0007
X direction	Elastic strain (×10−6)	1234.03	1234.20	1.0001
	Plastic strain (×10−6)	772.77	770.99	0.9977
	Total strain (×10−6)	2006.80	2005.19	0.9992
Y direction	Elastic strain (×10−6)	555.50	556.30	1.0014
	Plastic strain (×10−6)	0.00	0.00
	Total strain (×10−6)	555.50	556.30	1.0014
Z direction	Elastic strain (×10−6)	−3733.91	−3734.15	1.0001
	Plastic strain (×10−6)	−772.76	−770.98	0.9977
	Total strain (×10−6)	−4506.67	−4505.13	0.9997
Average value				0.9998
Coefficient of variation				0.0013

Note. The positive sign of strain is expansion, and the negative sign is compression. The model does not consider the dilatancy.

Table 5 shows the comparison between the numerical simulation solution and the theoretical analytical solution obtained from the true triaxial numerical simulation test. The elastic-plastic strain of each direction of the cube element is compared transversely, and the ratio of theoretical analytical solution and numerical simulation solution, the average value, and variation coefficient of each ratio are calculated.

Table 5

Comparison table of numerical simulation results of the true triaxial test.

Computation parameter
Constitutive model		Analytical solution			Numerical simulate solution			Analytical solution/numerical simulate solution
		εx (×10−6)	εy (×10−6)	εz (×10−6)	εx (×10−6)	εy (×10−6)	εz (×10−6)	εx (×10−6)	εy (×10−6)	εz (×10−6)
PMC-1	Elastic strain	1221.61	556.61	−3722.00	1234.20	556.30	−3734.15	0.9898	1.0006	0.9967
	Plastic strain	783.13		−783.13	770.99		−770.98	1.0157		1.0158
	Total strain	2004.74	556.61	−4505.13	2005.19	556.30	−4505.13	0.9998	1.0006	1
PMC-2	Elastic strain	1264.18	599.22	−3851.00	1267.07	599.22	−3854.15	0.9977	1	0.9992
	Plastic strain	654.13		−654.13	650.98		−650.98	1.0048		1.0048
	Total strain	1918.31	599.22	−4505.13	1918.05	599.22	−4505.13	1.0001	1	1
Average value								1.0013	1.0003	1.0028
Coefficient of variation								0.0078	0.0003	0.0013

It can be seen from the table that the ratio of each value of the theoretical analytical solution and the numerical simulation solution is equal to 1 or approaches to 1, and the coefficient of variation approaches to 0.

So far, it can be verified that the user-defined PMC numerical constitutive model developed based on the PMC strength criterion can be normally used in numerical calculation, and the numerical settlement results are consistent with the theoretical calculation results.

The true triaxial numerical simulation test results show that the developed PMC strength numerical constitutive model can reflect the influence of intermediate principal stress on the triaxial strength of soil-stone mixture and accurately simulate the strength peak value, elastic strain, and plastic strain of soil and stone mixture.

5. Conclusions

Through the research on the applicability of the PMC strength criterion in soil-stone mixture, the numerical constitutive model based on the PMC failure criterion is developed, and the true triaxial numerical simulation test is carried out. The main conclusions are as follows:

(1) Through the true triaxial test, it is found that the triaxial strength characteristics of soil-stone mixture are significantly affected by medium principal stress. The material parameters of soil-stone mixture under the PMC model are determined.

(2) Combined with the composite failure criterion of shear failure, a user-defined numerical constitutive model of PMC was developed and embedded in the software platform of FLAC3D. The numerical simulation solution is compared with the theoretical analytical solution, and the self-defined PMC elastic-plastic constitutive model is verified.

(3) The self-defined PMC numerical constitutive model can reflect the influence of medium principal stress on the triaxial strength of soil-stone mixture and accurately simulate the strength peak value, elastic strain, and plastic strain of soil-stone mixture, which provides an effective method for numerical simulation analysis of underground engineering of soil-stone mixture stratum.

Acknowledgments

The study was supported by the National Key Research and Development Project of China (2018YFE0101100) and General Project of Science and Technology Plan of Beijing Education Committee (KM202011418001).

References

[1] Z. Ouyang, S. Li, Z. Dai, "On the influence factors of mechanical properties for soil-rock mixture," Journal of Experimental Mechanics, vol. 25 no. 1, pp. 61-67, 2010.

[2] X. Yang, X. Hou, Y. Pei, "The influence of distribution of big stones on the stability of soil-rock mixture," Science Technology and Engineering, vol. 20 no. 31, pp. 12962-12967, 2020.

[3] Q. Jiang, Y. Xu, H. Wang, "Research on shear deformation characteristics of soil-rock mixtures under different stone contents," Journal of Engineering Geology, vol. 28 no. 5, pp. 951-958, DOI: 10.13544/j.cnki.jeg.2020-346, 2020.

[4] W. Ming, "Experimental study on mechanical properties of soil rock mixed heterogeneous filler," Highways, vol. 42 no. 1, pp. 40-49, 1997.

[5] H. Jing, Y. Zhang, X. Gu, "Large-scale triaxial shear test of soil rock mixture," Journal of Xi’an University of Science and Technology, vol. 39 no. 2, pp. 270-275, 2019.

[6] J. Yu, H. Sun, X. Zhang, "Effect of stone content on shear mechanical properties of soil-rock mixture," Journal of Nantong University (Natural Science Edition), vol. 19 no. 3, pp. 83-89, 2020.

[7] W. Minshuo, L. Xiao, H. Jianming, "In-situ direct-shear test on rock-soil aggregate," Geotechnical Engineering Technique, vol. 21 no. 4, pp. 26-31, 2007.

[8] S. Wang, T. Ji, Q. Xue, "Deformation and failure characteristics of soil-rock mixture considering material composition and random structure," Advances in Materials Science and Engineering, vol. 2019,DOI: 10.1155/2019/3165096, 2019.

[9] Y. Zhao, Z. Liu, "Numerical experiments on triaxial compression strength of soil-rock mixture," Advances in Civil Engineering, vol. 2019,DOI: 10.1155/2019/8763569, 2019.

[10] L. Jin, Y. Zeng, "Numerical simulation of large-scale triaxial tests on soil-rock mixture using DEM with three-dimensional flexible membrane boundary," Chinese Journal of Geotechnical Engineering, vol. 40 no. 12, pp. 2296-2303, 2018.

[11] Y. Zhan, H. Yao, "Influence of rock content on shear strength characteristics of soil-rock mixture," Highways, vol. 63 no. 2, pp. 20-25, 2018.

[12] L. E. Vallejo, S. Lobo-Guerrero, "The elastic moduli of clays with dispersed oversized particles," Engineering Geology, vol. 78 no. 1-2, pp. 163-171, DOI: 10.1016/j.enggeo.2004.12.003, 2005.

[13] L. E. Vallejo, G. S. Lobo, L. F. Seminsky, The Shear Strength of Sand-Gravel Mixtures: Laboratory and Theoretical Analysis, 2014.

[14] H. K. Lee, S. H. Pyo, "Multi-level modeling of effective elastic behavior and progressive weakened interface in particulate composites," Composites Science and Technology, vol. 68 no. 2, pp. 387-397, DOI: 10.1016/j.compscitech.2007.06.026, 2008.

[15] H. Ma, M. Gao, J. Zhang, Q. Yu, "Theoretical model developed for equivalent elastic modulus estimation of cobblestone-soil matrix," Rock and Soil Mechanics, vol. 32 no. 12, pp. 3642-3646, 2011.

[16] B. Haimson, J. W. Rudnicki, "The effect of the intermediate principal stress on fault formation and fault angle in siltstone," Journal of Structural Geology, vol. 32 no. 11, pp. 1701-1711, DOI: 10.1016/j.jsg.2009.08.017, 2010.

[17] L. Xiaochun, X. Dongjun, "Law and degree of effect the intermediate principle stress on strength of rock," Rock and Soil Mechanics, vol. 12 no. 1, 1991.

[18] B. Paul, "Generalized pyramidal fracture and yield criteria," International Journal of Solids and Structures, vol. 4 no. 2, pp. 175-196, DOI: 10.1016/0020-7683(68)90010-3, 1968.

[19] B. Paul, "Macroscopic criteria for plastic flow and brittle fracture," Fracture an Advanced Treatise, vol. 2, pp. 313-496, 1968.

[20] J. P. Meyer, J. F. Labuz, "Linear failure criteria with three principal stresses," International Journal of Rock Mechanics and Mining Sciences, vol. 60,DOI: 10.1016/j.ijrmms.2012.12.040, 2013.

[21] Z. Yiwu, Q. Lan, W. Qian, L. Yuan, "Destruction mode transformation type strength criterion based on PMC model," Journal of China Coal Society, vol. 44 no. 5, pp. 1404-1410, 2019.

[22] F. Zeng, Y. Li, J. F. Labuz, "Paul-Mohr-Coulomb failure criterion for geomaterials," Journal of Geotechnical and Geoenvironmental Engineering, vol. 144 no. 2,DOI: 10.1061/(asce)gt.1943-5606.0001829, 2018.