# Development of PMC Numerical Model for Soil and Stone Mixture

Authors: Zhili Sui; Zhaoguang Li; Yanhua Gao; Yuan Li; Yiwu Zheng; Yuan Meng
Journal: Advances in Civil Engineering
 (2021)

---

## 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 ofX direction, Y direction, and Z direction was 1.0013, 1.0003, and 1.002, respectively.


---

## Body

## 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:(1)Aσ1+Bσ2+Cσ3=1.In equation (1), A, B and C are: (2)A=1−sinϕc2sinϕc⋅1V0,B=sinϕc−sinϕe2sinϕcsinϕe⋅1V0,C=−1+sinϕe2sinϕe⋅1V0,where φc and φe are internal friction angles in compression and tension, respectively, and V0 is the theoretical value of triaxial tensile strength, as shown in Figure 1.Figure 1
Six-sided pyramidal failure surface of the PMC model.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 V0 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:(3)p=σ1+σ2+σ33,q=12σ1−σ22+σ2−σ32+σ3−σ12,tanθ=2σ2−σ1−σ33σ1−σ3,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:(4)q=6sinϕc3−sinϕcp+6ccsinϕc3−sinϕc.In the triaxial compression test (θ = 60°), the PMC strength criterion in the p−q plane can be expressed as follows:(5)q=6sinϕe3−sinϕep+6cesinϕe3−sinϕe,where cc and ce are compressive cohesion and tensile cohesion, respectively, as follows:(6)cc=V0tanϕc,ce=V0tanϕe.Whenφc = φe and cc = ce, φc in the PMC model is the internal friction angle φ and cc 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 Figure2). 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:(7)b=σ2−σ3σ1−σ3.Figure 2
TSW-50 soil true triaxial testing machine.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 coefficientb = 0, the specimens were first consolidated at p0 = 50, 100, 150, and 200 kPa, then kept σ3 unchanged, and σ1 gradually increased until the failure occurred. When the intermediate principal stress coefficient b = 1, the same consolidation processes were operated, then σ3 remained unchanged, and gradually increased σ1 to the failure strength. The curve of stress deviation vs. the axial principal strain (σ1–σ3∼ε1) is shown in Figure 3.Figure 3
σ1–σ3∼ε1 curve of the soil and stone mixture triaxial test: (a) b = 0; (b) b = 1.
(a)(b)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 inp-q plane, as shown in Figure 4.Figure 4
Comparison of PMC strength curve and test data inp-q plane of soil and stone mixture.The correlation coefficientR2 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)R2Soil-stone mixture38.142.76.40.92Through 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.


## 2.1. PMC Failure Criterion

The PMC failure criterion formula is as follows:(1)Aσ1+Bσ2+Cσ3=1.In equation (1), A, B and C are: (2)A=1−sinϕc2sinϕc⋅1V0,B=sinϕc−sinϕe2sinϕcsinϕe⋅1V0,C=−1+sinϕe2sinϕe⋅1V0,where φc and φe are internal friction angles in compression and tension, respectively, and V0 is the theoretical value of triaxial tensile strength, as shown in Figure 1.Figure 1
Six-sided pyramidal failure surface of the PMC model.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 V0 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:(3)p=σ1+σ2+σ33,q=12σ1−σ22+σ2−σ32+σ3−σ12,tanθ=2σ2−σ1−σ33σ1−σ3,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:(4)q=6sinϕc3−sinϕcp+6ccsinϕc3−sinϕc.In the triaxial compression test (θ = 60°), the PMC strength criterion in the p−q plane can be expressed as follows:(5)q=6sinϕe3−sinϕep+6cesinϕe3−sinϕe,where cc and ce are compressive cohesion and tensile cohesion, respectively, as follows:(6)cc=V0tanϕc,ce=V0tanϕe.Whenφc = φe and cc = ce, φc in the PMC model is the internal friction angle φ and cc 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 Figure2). 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:(7)b=σ2−σ3σ1−σ3.Figure 2
TSW-50 soil true triaxial testing machine.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 coefficientb = 0, the specimens were first consolidated at p0 = 50, 100, 150, and 200 kPa, then kept σ3 unchanged, and σ1 gradually increased until the failure occurred. When the intermediate principal stress coefficient b = 1, the same consolidation processes were operated, then σ3 remained unchanged, and gradually increased σ1 to the failure strength. The curve of stress deviation vs. the axial principal strain (σ1–σ3∼ε1) is shown in Figure 3.Figure 3
σ1–σ3∼ε1 curve of the soil and stone mixture triaxial test: (a) b = 0; (b) b = 1.
(a)(b)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 inp-q plane, as shown in Figure 4.Figure 4
Comparison of PMC strength curve and test data inp-q plane of soil and stone mixture.The correlation coefficientR2 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)R2Soil-stone mixture38.142.76.40.92Through 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:(8)εx=1Eσx−μσy+σz,εy=1Eσy−μσx+σz,εz=1Eσz−μσx+σy.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:(9)Δεi=Δεie+Δεip.(2)
The stress-strain constitutive equation for flow rule is as follows:(10)Δεip=λ∂g∂σi,i=1,2,3,whereg is the plastic potential function.(3)
The stress increment can be determined by the following equation:(11)Δσi=SiΔεne,i=1…n,whereSi is the linear equation which obeys Hooke’s Law.(4)
The new stress state also satisfies the yield equation as follows:(12)fσ¯n+Δσ¯n=0,whereΔσ¯i=SiΔε¯n−SiΔε¯np, and then it can be deducted that(13)f∗SnΔε¯n−λf∗Sn∂g∂σ¯n=f0¯n.The expression of coefficientλ can be obtained as follows:(14)λ=f∗SnΔε¯nf∗Sn∂g/∂σ¯n−f0¯n.


### 3.2. Composite Failure Criterion

The failure criterion of the PMC model is a composite failure criterion, as shown in Figure5. Considering three principal stresses, let(15)σ1≥σ2≥σ3.Figure 5
Schematic diagram of the composed strength model.To neglect hardening of materials, the strain increment can be divided into elastic and plastic strain increment as follows:(16)Δε=Δεe+Δεp,Δεip=λ∂g∂σi,i=1,2,3,where λ is a constant and g (σ1, σ2, σ3) is the plastic potential function.The plastic potential functionsgs (σ1, σ2, σ3) and gt (σ3) 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:(17)σ1V01−sinφc2sinφc+σ2V0sinφc−sinφe2sinφcsinφe−σ3V01+sinφe2sinφe=1,where σ1, σ2, and σ3 are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively; ϕc and ϕe are the internal friction angle in compression and tension cases, respectively; and V0 is the theoretical triaxial tensile strength.Let(18)Nc=1−sinϕc2sinϕc,Ne=1−sinϕe2sinϕe.Then,(19)fs=σ1Nc+σ2Nc−Ne−σ3Ne+1−V0=0.And the plastic potential function is as follows:(20)gs=σ1Nc+σ2Nc−Ne−σ3Ne+1.


#### 3.2.2. Tension Cut-Off

Ifft = 0 is used as the tension strength, then(21)ft=σ3−σt=0.Then, the plastic potential function is as follows:(22)gt=σ3.


### 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:(23)σiN=σiI−λSi∂g∂σn.For shear yield correction,σiI is as follows:(24)Δσ1=α1Δε1e+α2Δε2e+Δε3e,Δσ2=α1Δε2e+α2Δε1e+Δε3e,Δσ3=α1Δε3e+α2Δε2e+Δε3e,where α1 and α2 are the material parameters defined based on shear modulus G and bulk modulus K:(25)α1=K+43G,α2=K−23G.The yield function of PMCfσ1,σ2,σ3 can be derived as follows:(26)λ=fσ1I,σ2I,σ3IfS1∂g/∂σi,S2∂g/∂σi,S3∂g/∂σi,i=I,II,III,where σnIn=1,2,3 is the elastic stress in FLAC3D and f(0n) is the constant term of the linear failure criterion.In equation (27), the expression Si∂g/∂σn is as follows:(27)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε1e+α2Δε2e+Δε3e,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε2e+α2Δε1e+Δε3e,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε3e+α2Δε1e+Δε2e.After the material enters the plastic state, the flow criterion is as follows:(28)∂gs∂σ1=1,∂gs∂σ2=0,∂gs∂σ3=−Nψ,where ψ is the dilatancy angle and(29)Nψ=1+sinψ1−sinψ.Then,(30)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1−α2Nψ,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α21−Nψ,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=−α1Nψ+α2.When the model does not consider the dilatancy, ifψ = 0, then Nψ = 1. The plastic volume strain increment is 0.(31)Δεvp=∑Δεip=λ∂g∂σi=0.


#### 3.3.1. For Tension Failure

For tension failure,σiI is as follows:(32)Δσ1=α1Δε1p+α2Δε2e+Δε3e,Δσ2=α1Δε2p+α2Δε1e+Δε3e,Δσ3=α1Δε3p+α2Δε1e+Δε2e,λ would be(33)λ=σ3I−σ′α1,where σnI (n = 1, 2, 3) is the elastic stress in FLAC3D and f(0n) is the constant term of the linear failure criterion.Si∂g/∂σn is as follows:(34)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α1,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α2,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α3.


### 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 Table2.Table 2
PMC model parameters and keywords corresponding table.
Model keywordsModel parametersBulkBulk modulus,KShearShear modulus,Gc_frictionInternal friction angle in compression,φce_frictionInternal friction angle in tension,φev_paraTheoretical triaxial tensile strength,V0TensionEnsile strength,σtDilationDilation,ψ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.


## 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:(8)εx=1Eσx−μσy+σz,εy=1Eσy−μσx+σz,εz=1Eσz−μσx+σy.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:(9)Δεi=Δεie+Δεip.(2)
The stress-strain constitutive equation for flow rule is as follows:(10)Δεip=λ∂g∂σi,i=1,2,3,whereg is the plastic potential function.(3)
The stress increment can be determined by the following equation:(11)Δσi=SiΔεne,i=1…n,whereSi is the linear equation which obeys Hooke’s Law.(4)
The new stress state also satisfies the yield equation as follows:(12)fσ¯n+Δσ¯n=0,whereΔσ¯i=SiΔε¯n−SiΔε¯np, and then it can be deducted that(13)f∗SnΔε¯n−λf∗Sn∂g∂σ¯n=f0¯n.The expression of coefficientλ can be obtained as follows:(14)λ=f∗SnΔε¯nf∗Sn∂g/∂σ¯n−f0¯n.


## 3.2. Composite Failure Criterion

The failure criterion of the PMC model is a composite failure criterion, as shown in Figure5. Considering three principal stresses, let(15)σ1≥σ2≥σ3.Figure 5
Schematic diagram of the composed strength model.To neglect hardening of materials, the strain increment can be divided into elastic and plastic strain increment as follows:(16)Δε=Δεe+Δεp,Δεip=λ∂g∂σi,i=1,2,3,where λ is a constant and g (σ1, σ2, σ3) is the plastic potential function.The plastic potential functionsgs (σ1, σ2, σ3) and gt (σ3) 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:(17)σ1V01−sinφc2sinφc+σ2V0sinφc−sinφe2sinφcsinφe−σ3V01+sinφe2sinφe=1,where σ1, σ2, and σ3 are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively; ϕc and ϕe are the internal friction angle in compression and tension cases, respectively; and V0 is the theoretical triaxial tensile strength.Let(18)Nc=1−sinϕc2sinϕc,Ne=1−sinϕe2sinϕe.Then,(19)fs=σ1Nc+σ2Nc−Ne−σ3Ne+1−V0=0.And the plastic potential function is as follows:(20)gs=σ1Nc+σ2Nc−Ne−σ3Ne+1.


### 3.2.2. Tension Cut-Off

Ifft = 0 is used as the tension strength, then(21)ft=σ3−σt=0.Then, the plastic potential function is as follows:(22)gt=σ3.


## 3.2.1. Shear Failure

PMC criteria can be expressed as follows:(17)σ1V01−sinφc2sinφc+σ2V0sinφc−sinφe2sinφcsinφe−σ3V01+sinφe2sinφe=1,where σ1, σ2, and σ3 are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively; ϕc and ϕe are the internal friction angle in compression and tension cases, respectively; and V0 is the theoretical triaxial tensile strength.Let(18)Nc=1−sinϕc2sinϕc,Ne=1−sinϕe2sinϕe.Then,(19)fs=σ1Nc+σ2Nc−Ne−σ3Ne+1−V0=0.And the plastic potential function is as follows:(20)gs=σ1Nc+σ2Nc−Ne−σ3Ne+1.


## 3.2.2. Tension Cut-Off

Ifft = 0 is used as the tension strength, then(21)ft=σ3−σt=0.Then, the plastic potential function is as follows:(22)gt=σ3.


## 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:(23)σiN=σiI−λSi∂g∂σn.For shear yield correction,σiI is as follows:(24)Δσ1=α1Δε1e+α2Δε2e+Δε3e,Δσ2=α1Δε2e+α2Δε1e+Δε3e,Δσ3=α1Δε3e+α2Δε2e+Δε3e,where α1 and α2 are the material parameters defined based on shear modulus G and bulk modulus K:(25)α1=K+43G,α2=K−23G.The yield function of PMCfσ1,σ2,σ3 can be derived as follows:(26)λ=fσ1I,σ2I,σ3IfS1∂g/∂σi,S2∂g/∂σi,S3∂g/∂σi,i=I,II,III,where σnIn=1,2,3 is the elastic stress in FLAC3D and f(0n) is the constant term of the linear failure criterion.In equation (27), the expression Si∂g/∂σn is as follows:(27)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε1e+α2Δε2e+Δε3e,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε2e+α2Δε1e+Δε3e,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1Δε3e+α2Δε1e+Δε2e.After the material enters the plastic state, the flow criterion is as follows:(28)∂gs∂σ1=1,∂gs∂σ2=0,∂gs∂σ3=−Nψ,where ψ is the dilatancy angle and(29)Nψ=1+sinψ1−sinψ.Then,(30)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α1−α2Nψ,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=α21−Nψ,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3,=−α1Nψ+α2.When the model does not consider the dilatancy, ifψ = 0, then Nψ = 1. The plastic volume strain increment is 0.(31)Δεvp=∑Δεip=λ∂g∂σi=0.


### 3.3.1. For Tension Failure

For tension failure,σiI is as follows:(32)Δσ1=α1Δε1p+α2Δε2e+Δε3e,Δσ2=α1Δε2p+α2Δε1e+Δε3e,Δσ3=α1Δε3p+α2Δε1e+Δε2e,λ would be(33)λ=σ3I−σ′α1,where σnI (n = 1, 2, 3) is the elastic stress in FLAC3D and f(0n) is the constant term of the linear failure criterion.Si∂g/∂σn is as follows:(34)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α1,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α2,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α3.


## 3.3.1. For Tension Failure

For tension failure,σiI is as follows:(32)Δσ1=α1Δε1p+α2Δε2e+Δε3e,Δσ2=α1Δε2p+α2Δε1e+Δε3e,Δσ3=α1Δε3p+α2Δε1e+Δε2e,λ would be(33)λ=σ3I−σ′α1,where σnI (n = 1, 2, 3) is the elastic stress in FLAC3D and f(0n) is the constant term of the linear failure criterion.Si∂g/∂σn is as follows:(34)S1∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α1,S2∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α2,S3∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=α3.


## 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 Table2.Table 2
PMC model parameters and keywords corresponding table.
Model keywordsModel parametersBulkBulk modulus,KShearShear modulus,Gc_frictionInternal friction angle in compression,φce_frictionInternal friction angle in tension,φev_paraTheoretical triaxial tensile strength,V0TensionEnsile strength,σtDilationDilation,ψ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ϕc=ϕ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 ϕc=ϕe (denoted as PMC-1); and condition 3 by the PMC model with parameters ϕc≠ϕ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 conditionsE (GPa)μφc (°)φe (°)cc (MPa)ce (MPa)V0 (MPa)σci (MPa)M-C10.00.3328.0—8.5—16.028.29PMC-110.00.3328.028.08.58.516.028.29PMC-210.00.3328.033.08.510.416.028.29Note. σci is the uniaxial compressive strength.A cube element (size: 1 × 1 × 1) is created, and the fixed boundary surfacex = 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 6
Model boundary condition diagram.To ensure the quasi-static loading, the loading speed is set at 10−8 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−3. 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

Figure7 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−6. 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−6.Figure 7
Vertical stress curve of the PMC-2 model.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.Table4 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 ϕc=ϕ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 ϕc=ϕ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 conditionsM-CPMC-1PMC-1/M-CPeak strength (MPa)42.1442.171.0007X directionElastic strain (×10−6)1234.031234.201.0001Plastic strain (×10−6)772.77770.990.9977Total strain (×10−6)2006.802005.190.9992Y directionElastic strain (×10−6)555.50556.301.0014Plastic strain (×10−6)0.000.00Total strain (×10−6)555.50556.301.0014Z directionElastic strain (×10−6)−3733.91−3734.151.0001Plastic strain (×10−6)−772.76−770.980.9977Total strain (×10−6)−4506.67−4505.130.9997Average value0.9998Coefficient of variation0.0013Note. The positive sign of strain is expansion, and the negative sign is compression. The model does not consider the dilatancy.Table5 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 parameterConstitutive modelAnalytical solutionNumerical simulate solutionAnalytical 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-1Elastic strain1221.61556.61−3722.001234.20556.30−3734.150.98981.00060.9967Plastic strain783.13−783.13770.99−770.981.01571.0158Total strain2004.74556.61−4505.132005.19556.30−4505.130.99981.00061PMC-2Elastic strain1264.18599.22−3851.001267.07599.22−3854.150.997710.9992Plastic strain654.13−654.13650.98−650.981.00481.0048Total strain1918.31599.22−4505.131918.05599.22−4505.131.000111Average value1.00131.00031.0028Coefficient of variation0.00780.00030.0013It 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.


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

The PMC failure criterion is derived from the M-C strength criterion. Whenϕc=ϕ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 ϕc=ϕe (denoted as PMC-1); and condition 3 by the PMC model with parameters ϕc≠ϕ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 conditionsE (GPa)μφc (°)φe (°)cc (MPa)ce (MPa)V0 (MPa)σci (MPa)M-C10.00.3328.0—8.5—16.028.29PMC-110.00.3328.028.08.58.516.028.29PMC-210.00.3328.033.08.510.416.028.29Note. σci is the uniaxial compressive strength.A cube element (size: 1 × 1 × 1) is created, and the fixed boundary surfacex = 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 6
Model boundary condition diagram.To ensure the quasi-static loading, the loading speed is set at 10−8 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−3. 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

Figure7 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−6. 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−6.Figure 7
Vertical stress curve of the PMC-2 model.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.Table4 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 ϕc=ϕ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 ϕc=ϕ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 conditionsM-CPMC-1PMC-1/M-CPeak strength (MPa)42.1442.171.0007X directionElastic strain (×10−6)1234.031234.201.0001Plastic strain (×10−6)772.77770.990.9977Total strain (×10−6)2006.802005.190.9992Y directionElastic strain (×10−6)555.50556.301.0014Plastic strain (×10−6)0.000.00Total strain (×10−6)555.50556.301.0014Z directionElastic strain (×10−6)−3733.91−3734.151.0001Plastic strain (×10−6)−772.76−770.980.9977Total strain (×10−6)−4506.67−4505.130.9997Average value0.9998Coefficient of variation0.0013Note. The positive sign of strain is expansion, and the negative sign is compression. The model does not consider the dilatancy.Table5 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 parameterConstitutive modelAnalytical solutionNumerical simulate solutionAnalytical 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-1Elastic strain1221.61556.61−3722.001234.20556.30−3734.150.98981.00060.9967Plastic strain783.13−783.13770.99−770.981.01571.0158Total strain2004.74556.61−4505.132005.19556.30−4505.130.99981.00061PMC-2Elastic strain1264.18599.22−3851.001267.07599.22−3854.150.997710.9992Plastic strain654.13−654.13650.98−650.981.00481.0048Total strain1918.31599.22−4505.131918.05599.22−4505.131.000111Average value1.00131.00031.0028Coefficient of variation0.00780.00030.0013It 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.



---

*Source: 9967412-2021-06-28.xml*