Integration and Application of a Fabric-Based Modified Cam-Clay Model in FLAC^3D

Abstract

In order to consider the effect of fabric anisotropy in the analysis of geotechnical boundary value problems, this study proposes a modified model based on a fabric-based modified Cam-clay model, which can account for the anisotropic response of soil. The major modification of the original model aims to simplify the equations for numerical implementation by replacing the SMP strength criterion with the Lade’s strength criterion. This model comprehensively considers the inherent anisotropy, induced anisotropy, and three-dimensional strength characteristics of soil. The model is first numerically implemented using the elastic trial–plastic correction method, and then it is encapsulated into the FLAC^3D 6.0 software, and tested through conventional triaxial, embankment loading, and tunnel excavation experiments. Numerical simulation results indicate that considering anisotropy and three-dimensional strength in geotechnical engineering analysis is necessary. By accounting for the interaction between microstructure and macroscopic anisotropy, the model can more accurately represent soil behavior, providing significant advantages for geotechnical analysis.

1. Introduction

To analyze and calculate the mechanical response of geotechnical bodies under various complex conditions, numerous successful constitutive models have been proposed [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. However, the validation and application of most models are largely limited to element-level problems or practical engineering cases with structural symmetry and simple boundary conditions. To analyze the complex geotechnical problems frequently encountered in practice, it is inevitable to integrate these constitutive models into numerical analysis software [15,16,17,18,19]. Currently, the Mohr–Coulomb (M-C) model and the modified Cam-clay (MCC) model in elastoplastic frameworks have been integrated into nearly all numerical analysis software for geotechnical engineering and are widely accepted for analyzing practical engineering problems [20,21]. However, as pointed out by some researchers [22,23,24,25], the selection of constitutive models is crucial in the analysis of geotechnical boundary value problems. An inappropriate choice of constitutive model may result in unrealistic predictions. It is well known that the circular yield surfaces of these two models cannot capture the dependence of strength on Lode’s angle. Consequently, they tend to significantly overestimate the tensile strength of geotechnical materials. Given that underground structures are commonly subjected to cyclic tensile and compressive loads during their service life, such as seismic loads [26,27], the strength/yield criterion should be extended to the three-dimensional stress space when analyzing their dynamic response [28,29,30].

It is also noteworthy that, due to gravitational effects during natural sedimentation, soil particles form spatial arrangements of particles and associated voids with certain statistical characteristics, displaying anisotropic features at the macroscopic level (known as inherent anisotropy) [31,32]. Because of inherent anisotropy, the mechanical response of soil is highly sensitive to the loading direction. Specifically, lower strength is observed when the principal stress direction deviates from the deposition direction [33,34,35]. Additionally, non-coaxial behavior has been observed, where a significant difference exists between the directions of the principal stress and the principal strain increment [36,37,38]. On the other hand, under external loading during service, soil particles undergo rearrangement (leading to induced anisotropy), which further results in the accumulation of pore water pressure or plastic deformation [39,40,41,42]. Extensive research has confirmed that neglecting inherent or induced anisotropy in the analysis of underground structures can lead to designs that deviate from actual conditions and pose safety risks [43,44,45,46].

To analyze the anisotropic behavior of soils, a large number of successful models have been proposed, among which the S-CLAY model series and the SANISAND model series are particularly widely accepted and applied [47,48,49,50,51,52]. Although these two types of models have different applicable scopes, they are both based on the same concept of rotational hardening. These models aim to represent the anisotropic response of soils by using yield surfaces that rotate with plastic deformation. These models have been successful on a macroscopic scale because they can clearly reflect the relationship between soil strength (or deformation) and anisotropy. However, a large body of evidence from a microscopic perspective has shown that this macroscopic anisotropic behavior is closely related to the microstructure of soil [53,54,55,56,57]. Rotational hardening models, which are based on phenomenological approaches, struggle to quantitatively explain this relationship. In recent years, thanks to advancements in DEM technology, many multiscale models incorporating soil microstructural information have been developed [1,3,58,59,60,61]. Notably, interesting attempts have been made by Yao and Kong [62], Yao et al. [63], as well as Tian and Yao [64]. In their work, they proposed a modeling approach called the Anisotropic Transformation Stress (ATS) Method. Under the framework of ATS, the fabric tensor representing microscopic anisotropy information can be easily integrated into the constitutive model through a two-step stress transformation without introducing additional model parameters. Furthermore, by combining a specific strength criterion, the generalization of the yield surface can be achieved without introducing additional parameters. Using this method and combining it with the MCC model, Tian [65] successfully proposed a fabric-based MCC model. The validity of this model has been verified at both macroscopic and microscopic scales.

In this study, building upon the model proposed by Tian [65], the fabric evolution law and the strength criterion are further modified, and an improved version of the fabric-based MCC model is introduced. The purpose of these modifications is to simplify the model equations, thereby facilitating numerical implementation. Specifically, the SMP strength criterion used in the original model is replaced with Lade’s strength criterion, which avoids the need to compute the second invariant of the principal stresses. In this model, the inherent anisotropy of the soil is characterized by a fabric tensor defined based on the initial preferential orientation of particles [66,67,68], while the induced anisotropy is quantified through a fabric evolution rule driven by the plastic strain rate [69,70,71]. Subsequently, the model is numerically implemented using an elastic trial–plastic correction approach and integrated into the FLAC^3D software. In this approach, the plastic multiplier of the model can be easily and explicitly calculated using the quadratic equation root-finding formula. Finally, a series of numerical experiments were conducted using the model in the FLAC^3D software. These experiments included conventional triaxial tests, embankment load tests, and tunnel excavation tests. The results demonstrated the effectiveness and validity of the model. The study found that the model could reasonably capture the three-dimensional strength characteristics and anisotropic behavior of the soil, and these features are essential and cannot be ignored in practical engineering applications.

2. Constitutive Relation of Fabric-Based MCC Model

2.1. MCC Model, a Brief Review

The development of the MCC model was a major breakthrough in the elastoplastic modeling of geomaterials [72]. Although its application in practical engineering remains somewhat debated, the underlying concept of the MCC model is remarkably straightforward and serves as a powerful tool [73]. It is particularly well suited for describing the mechanical behavior of normally consolidated or lightly over-consolidated soft clays. Additionally, many constitutive models in the literature have been built around the core principles of MCC plasticity [17,28,62,74,75]. The basic equations for the MCC model can be given as follows:

$$f=g=q^2+M^{2}p\left(p-p_c\right),$$

(1)

$$p={\sigma }_{ii}/3, s_{ij}={\sigma }_{ij}-p{\delta }_{ij}, q=\sqrt{3s_{ij}{s}_{ij}/2},$$

(2)

$$e_N=e+\lambda \mathrm{l}\mathrm{n}p,$$

(3)

where f and g are the yield and plastic potential functions, respectively. ${\sigma }_{ij}$ is effective stress tensor. $s_{ij}$ is deviatoric stress tensor. p and q are the mean and deviatoric stress, respectively. ${\delta }_{ij}$ is Kronecker tensor. $p_c$ is consolidation pressure. M denotes the critical stress ratio. e is void ratio. $e_N$ denotes the intersection of the normal consolidation line (NCL) with the e axis in the e-lnp space. $\lambda$ denotes the slope of compression line.

2.2. MCC Model Enhanced by Fabirc Anisotropy and Lade’s Criterion

To measure the anisotropic characteristics of soil, various definition methods for fabric tensors have been proposed [76]. According to the work of Yao et al. [63], Tian and Yao [64], and Tian et al. [65], if the preferred orientation of particle long axes is used to define the fabric tensor $F_{ij}$, it can be expressed as follows:

$$F_{ij}=\left(\begin{array}{ccc}\Delta & 0& 0\\ 0& \left(1-\Delta \right)/2& 0\\ 0& 0& \left(1-\Delta \right)/2\end{array}\right),$$

(4)

where $\Delta$ represents the degree of the anisotropy of the preferred orientation of the particles [58,77]. Furthermore, the deviatoric part of $F_{ij}$ can be given as follows [68]:

$$P_{ij}={\displaystyle \frac{5}{2}}\left(3F_{ij}-{\delta }_{ij}\right), A=\sqrt{P_{ij}{P}_{ij}},$$

(5)

where $P_{ij}$ denotes the deviatoric part of $F_{ij}$. A is the norm of $P_{ij}$, used to measure the degree of anisotropy.

In order to quantitatively describe the variation in fabric tensor during the loading process, a fabric evolution law should be defined. Tian and Yao [64] and Tian et al. [65] have proposed a unified fabric evolution law driven by plastic volumetric deformation, which has been successfully applied to both sand and clay. In this paper, to facilitate numerical implementation, the fabric evolution rule is modified as follows:

$$\mathrm{d}{F}_{ij}=c_F\left({\delta }_{ij}/3-\beta {\eta }_{ij}-F_{ij}\right)\mathrm{d}\gamma$$

(6)

where $c_F$ is a model parameter used to control the evolution rate of fabric tensor under plastic loading. β is a model constant used to control the final degree of anisotropy. ${\eta }_{ij}=s_{ij}/p$ is stress ratio tensor. $\mathrm{d}\gamma$ is plastic multiplier.

To extend the yield surface to three-dimensional anisotropic space, Yao et al. [78], Yao and Kong [62], and Yao et al. [63] have proposed the anisotropic transformed stress (ATS) method. This method maps the stress into a virtual transformed stress space by combining the fabric tensor and a specific strength criterion, allowing for the description of the effects of fabric anisotropy and Lode’s angle on yield strength within this space. According to the work of Tian [65], the stress is first modified using the fabric tensor to extend the yield surface of the MCC model to account for anisotropy, as follows:

$${\overline{\sigma }}_{ij}={\displaystyle \frac{3}{2}}\left({\sigma }_{ik}{F}_{kj}+F_{ik}{\sigma }_{kj}\right)-s_{lm}{F}_{ml}{\delta }_{ij},$$

(7)

$$\overline{p}={\overline{\sigma }}_{ii}/3, {\overline{s}}_{ij}={\overline{\sigma }}_{ij}-\overline{p}{\delta }_{ij}, \overline{q}=\sqrt{3{\overline{s}}_{ij}{\overline{s}}_{ij}/2},$$

(8)

where ${\overline{\sigma }}_{ij}$ is referred to as the modified stress tensor [64,65], which contains information related to fabric. ${\overline{s}}_{ij}$ is the deviatoric part of ${\overline{\sigma }}_{ij}$. $\overline{p}$ and $\overline{q}$ are the mean and deviatoric stress of ${\overline{\sigma }}_{ij}$.

Furthermore, to extend the yield surface to a three-dimensional stress space, the modified stress should be combined with a specific strength criterion. In the work of Tian [65], the SMP criterion was used to create a three-dimensional yield surface for the MCC model. However, this method requires obtaining all three invariants of the modified stress, which complicates numerical implementation. For simplicity, this paper employs Lade’s criterion to three-dimensionalize the MCC yield surface. This criterion requires only the first and third invariants of the modified stress and its validity has been confirmed as follows [28,64]:

$${\tilde{\sigma }}_{ij}=\left\{\begin{array}{l}\overline{p}{\delta }_{ij}+{\overline{q}}_c{\overline{s}}_{ij}/\overline{q}, \left(\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n} \overline{q}\ne 0\right)\\ {\overline{\sigma }}_{ij}, \left(\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n} \overline{q}=0\right)\end{array}\right.,$$

(9)

$${\overline{q}}_c={\overline{I}}_1\left\{1+{\displaystyle \frac{\overline{J}}{2}}{\left[\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{1}{3}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{-1}\overline{J}\right)\right]}^{-1}\right\}, \overline{J}=\sqrt{27{\overline{I}}_3/{\overline{I}}_1^3},$$

(10)

where ${\tilde{\sigma }}_{ij}$ is referred to as the transformed stress tensor. ${\overline{I}}_1$ and ${\overline{I}}_3$ represent the first and third invariants of ${\overline{\sigma }}_{ij}$, respectively. $\overline{J}$ is a variable determined by current modified stress. ${\overline{q}}_c$ denotes the strength under triaxial compression condition.

By substituting the transformed stress into yield function defined in the real stress space, the yield function $\tilde{f}$ and plastic potential function $\tilde{g}$ of the modified MCC model can be given as follows [65]:

$$\tilde{f}=\tilde{g}={\tilde{q}}^2+{\tilde{M}}^2\tilde{p}\left(\tilde{p}-{\tilde{p}}_c\right),$$

(11)

$$\tilde{p}={\tilde{\sigma }}_{ii}/3, {\tilde{s}}_{ij}={\tilde{\sigma }}_{ij}-\tilde{p}{\delta }_{ij}, \tilde{q}=\sqrt{3{\tilde{s}}_{ij}{\tilde{s}}_{ij}/2},$$

(12)

$$\tilde{M}=M-\beta M\left(M+3\right),$$

(13)

where ${\tilde{s}}_{ij}$ is the deviatoric part of ${\tilde{\sigma }}_{ij}$. $\tilde{p}$ and $\tilde{q}$ are the mean and deviatoric stress of ${\tilde{\sigma }}_{ij}$. $\tilde{M}$ denotes the critical state ratio in modified/transformed stress space [63,64,65]. It is important to note that Equation (11) implies that the flow rule of this model is associated. This is because Hashiguchi [79] emphasized that if a constitutive model is developed within a rate-independent framework, the associated flow rule must hold. ${\tilde{p}}_c$ denotes the hardening variable, and its evolution rule can be given as follows:

$$\mathrm{d}{\tilde{p}}_c={\displaystyle \frac{1+e}{\lambda -\kappa }}{\tilde{p}}_c\mathrm{d}{\epsilon }_v^p,$$

(14)

where $\kappa$ is the slope of swelling line. $\mathrm{d}{\epsilon }_v^p$ denotes the plastic volumetric increment.

Furthermore, based on the elastoplastic theory, the stress–strain relationship can be expressed as follows:

$$\mathrm{d}p=K\left(\mathrm{d}{\epsilon }_v-\mathrm{d}{\epsilon }_v^p\right)=K\left(\mathrm{d}{\epsilon }_v-\mathrm{d}\gamma {\tilde{N}}_v\right),$$

(15)

$$\mathrm{d}q=3G\left(\mathrm{d}{\epsilon }_q-\mathrm{d}{\epsilon }_q^p\right)=3G\left(\mathrm{d}{\epsilon }_q-\mathrm{d}\gamma {\tilde{N}}_q\right),$$

(16)

$$K={\displaystyle \frac{1+e}{\kappa }}p, G={\displaystyle \frac{3\left(1-2v\right)}{2\left(1+v\right)}}K,$$

(17)

$${\tilde{N}}_v={\displaystyle \frac{\mathsf{\partial }\tilde{g}}{\mathsf{\partial }\tilde{p}}}={\tilde{M}}^2\left(2\tilde{p}-{\tilde{p}}_c\right), {\tilde{N}}_q={\displaystyle \frac{\mathsf{\partial }\tilde{g}}{\mathsf{\partial }\tilde{q}}}=2\tilde{q},$$

(18)

where $\mathrm{d}{\epsilon }_v$ and $\mathrm{d}{\epsilon }_q$ are the volumetric and deviatoric strain increments, respectively. $\mathrm{d}{\epsilon }_q^p$ is the plastic deviatoric strain increment. K and G are the bulk and shear modulus, respectively. v is Poisson’s ratio. ${\tilde{N}}_v$ and ${\tilde{N}}_q$ denotes the plastic volumetric and deviatoric flow directions, respectively.

For convenience of comparison,

Table 1. Model parameters of MCC model and this model.

Model	Basic Parameters [72]	Fabric Parameters [63,64,65]
MCC model	$M, \lambda , \kappa , v, e_N$$(\mathrm{or} e_{ref}$)	/
Fabric-based MCC model	$M, \lambda , \kappa , v, e_N$$(\mathrm{or} e_{ref}$)	$\mathsf{\Delta }, \beta , c_F$

Remarks: M is critical state ratio; λ is the slope of critical state line; κ is the slope of swelling line; v is the Poisson’s ratio; $e_N$ is the intersection of consolidation line with the e axis in the e-lnp space. $e_{ref}$ denotes the reference void ratio at a specific stress level. Δ denotes the initial fabric anisotropy. β denotes the final fabric anisotropy. $C_F$ controls the evolution rate of fabric anisotropy.

lists the model parameters of the MCC model and the fabric-based MCC model defined in this study, where M, λ, κ, v, and $e_N$ (or $e_{ref}$) are referred to as the basic parameters in this study, as they are common to all models in the MCC family. Compared to the MCC model, the model in this study also requires the determination of parameters such as Δ, β, and $c_F$. These parameters, first proposed by Yao et al. [63], Tian and Yao [64], and Tian et al. [65], are related to the fabric anisotropy of the material and are therefore referred to as fabric parameters in this study.

3. Integration of Fabric-Based MCC Model in FLAC^3D

The numerical algorithm flowchart of this model is shown within the blue box in Figure 1, where superscripts k and k + 1 represent the k th and k + 1 th loading steps, respectively. The superscript 0 denotes the initial state, and the superscript tr denotes the elastic trial state. Before starting the calculations, it is necessary to assign initial values to the internal variables. The initial value of the plastic state variable ${\tilde{p}}_c^0$ in the transformed stress space can be determined by substituting the initial values of the other internal variables into Equation (11), as follows:

$${\tilde{p}}_c^0=\mathrm{O}\mathrm{C}\mathrm{R}\left[{\left({\tilde{q}}^0\right)}^2/\left({\tilde{M}}^2{\tilde{p}}^0\right)+{\tilde{p}}^0\right].$$

(19)

Specifically, when ${\tilde{q}}^0$ = 0 (isotropically consolidation), ${\tilde{p}}_c^0$ can be given as follows:

$${\tilde{p}}_c^0=\mathrm{O}\mathrm{C}\mathrm{R}{\tilde{p}}^0,$$

(20)

where OCR is the over-consolidation ratio. Specifically, when Δ = 1/3 and β = 0, the proposed model reduces to the standard MCC model (isotropic state), and, in this case, ${\tilde{p}}_c^0=p_c^0$.

The integration algorithm of the model primarily consists of two parts: elastic prediction and plastic correction. First, assuming purely elastic deformation occurs between steps k and k + 1, the elastic trial stress can be updated using Equations (15) and (16), while the internal variables $F_{ij}$ and ${\tilde{p}}_c$ related to plastic deformation are held constant, as follows:

$$p^{\mathrm{t}\mathrm{r}}=p^k+K^{k}d{\epsilon }_v^{k+1}, q^{\mathrm{t}\mathrm{r}}=q^k+3G^{k}d{\epsilon }_q^{k+1},$$

(21)

$$F_{ij}^{\mathrm{t}\mathrm{r}}=F_{ij}^k, {\tilde{p}}_c^{\mathrm{t}\mathrm{r}}={\tilde{p}}_c^k.$$

(22)

Then, use $p^{\mathrm{t}\mathrm{r}}$, $q^{\mathrm{t}\mathrm{r}}$ and $F_{ij}^{\mathrm{t}\mathrm{r}}$ to calculate ${\tilde{p}}^{\mathrm{t}\mathrm{r}}$, ${\tilde{q}}^{\mathrm{t}\mathrm{r}}$ and ${\tilde{\sigma }}_{ij}^{\mathrm{t}\mathrm{r}}$, and substitute these trial variables into the yield function defined in Equation (11) for evaluation. If the trail yield function ${\tilde{f}}^{\mathrm{t}\mathrm{r}}$ is less than 0, it indicates that the elastic prediction is correct, and the values for state variables at step k + 1 can be directly output as follows:

$$p^{k+1}=p^{\mathrm{t}\mathrm{r}}, q^{k+1}=q^{\mathrm{t}\mathrm{r}}, F_{ij}^{k+1}=F_{ij}^k, {\tilde{p}}_c^{k+1}={\tilde{p}}_c^k.$$

(23)

If ${\tilde{f}}^{\mathrm{t}\mathrm{r}}\ge 0$, this indicates that plastic deformation has occurred, and plastic correction is required. To solve for $\mathrm{d}{\gamma }^{k+1}$, it can be substituted as the sole unknown into the yield function. Various integration algorithms, including explicit and implicit methods [80,81,82], have been proposed to solve for $\mathrm{d}{\gamma }^{k+1}$. Here, noting that the yield function ${\tilde{f}}^{\mathrm{t}\mathrm{r}}$ = 0 is a univariate quadratic equation, $\mathrm{d}{\gamma }^{k+1}$ can be explicitly obtained using the following quadratic formula [83]:

$$a{\left(\mathrm{d}{\gamma }^{k+1}\right)}^2+b\mathrm{d}{\gamma }^{k+1}+c=0,$$

(24)

$$a={\left(\tilde{M}{K}^k{\tilde{N}}_v^{\mathrm{t}\mathrm{r}}\right)}^2+{\left(3G^k{\tilde{N}}_q^{\mathrm{t}\mathrm{r}}\right)}^2, b=-K^k{\left({\tilde{N}}_v^{\mathrm{t}\mathrm{r}}\right)}^2-3G^k{\left({\tilde{N}}_q^{\mathrm{t}\mathrm{r}}\right)}^2, c={\tilde{f}}^{\mathrm{t}\mathrm{r}}.$$

(25)

Once the plastic multiplier $\mathrm{d}{\gamma }^{k+1}$ is obtained, it can be substituted into Equations (6) and (14)–(16) to perform plastic correction of the internal variables as follows:

$$p^{k+1}=p^{\mathrm{t}\mathrm{r}}-K^k\mathrm{d}{\gamma }^{k+1}{\tilde{N}}_v^{\mathrm{t}\mathrm{r}}, q^{k+1}=q^{\mathrm{t}\mathrm{r}}-3G^k\mathrm{d}{\gamma }^{k+1}\mathrm{d}{\tilde{N}}_q^{\mathrm{t}\mathrm{r}},$$

(26)

$$F_{ij}^{k+1}=F_{ij}^k+c_F\left({\delta }_{ij}/3-\beta {\eta }_{ij}^k-F_{ij}^k\right)\mathrm{d}{\gamma }^{k+1},$$

(27)

$${\tilde{p}}_c^{k+1}={\tilde{p}}_c^k\left(1+\mathrm{d}{\gamma }^{k+1}{\tilde{N}}^k{\displaystyle \frac{1+e}{\lambda -\kappa }}\right).$$

(28)

The numerical implementation and coding processes are carried out in a C++ programming environment. After completion, the code generates a file with a .dll extension, which is then imported into the “cmodel” directory of the FLAC^3D software to enable program encapsulation. When using the model, it can be invoked with the command “model configure plugin”.

In practical applications, as depicted in the left half of Figure 1, the initial environment necessary for computation is first set up, including model parameters, geometric information, initial boundary conditions, and initial conditions (e.g., stress field). FLAC3D then calculates the strain increments for each element using the explicit finite difference method (FDM). With these strain increments, the stress increments can be calculated using the constitutive relations shown in the blue region. Finally, FLAC3D checks whether the results satisfy the equilibrium equations. If the equilibrium is met, the results are output; otherwise, iterative calculations are performed.

4. Application of Fabric-Based MCC Model in FLAC^3D

To verify the applicability of the model in FLAC^3D, a series of numerical simulations are conducted in this study, including triaxial drained and undrained simulations, embankment loading simulations, and tunnel excavation simulations. The model parameters used for simulation are listed in

Table 2. Model parameters used for simulation.

Case	Basic Parameters	Fabric Parameters
Triaxial tests	$M=1; v=0.3; \lambda =0.1;$ $\kappa =0.03; e_N=1.391$	$\mathsf{\Delta }=0.25~1/3; \beta =0.01~0.07$ $c_F=2000~\mathrm{12,000}$
Embankment loading	$M=0.888; v=0.3; \lambda =0.161;$ $\kappa =0.062; e_{ref}=1.858 (p=1kPa)$	$\mathsf{\Delta }=1/3; \beta =0.02~0.08$ $c_F=500$
Tunnel excavation	$M=1.85; v=0.28; \lambda =0.072;$ $\kappa =0.006; e_{ref}=1.094 (p=1kPa)$	$\mathsf{\Delta }=0.2~0.32; \beta =0.09~0.11$ $c_F=10~5000$

Remarks: the basic model parameters used in the triaxial tests are selected from the typical range of clay parameters provided by Yu [74], while the range of fabric parameters is determined with reference to the work of Tian and Yao [64] and Tian et al. [65].

.

4.1. Triaxial Tests

To verify the accuracy of this model, a series of conventional triaxial tests (element tests) are first conducted in FLAC^3D. The basic model parameters are shown in Table 2, while the fabric parameters are set as Δ = 1/3 and β = $c_F$ = 0. This configuration is chosen to temporarily ignore the effects of anisotropy and focus solely on verifying whether Lade’s strength criterion is successfully incorporated into the model. The sample is subjected to isotropic consolidation with p = 300 kPa and subsequently undergoes triaxial compression (TC) and triaxial extension (TE) loading under both drained and undrained conditions. For comparison, the MCC model in FLAC^3D is used with the same model parameters for simulation. Figure 2 shows the calculation results of the MCC model and this model, where the symbols represent the calculation results of the MCC model, and the lines represent the simulation results of this model. Under TC conditions, the mechanical response predicted by this model aligns almost perfectly with that of the MCC model, confirming the correctness of this model. However, under TE conditions, the MCC model predicts identical ratios for the critical state stresses under compression and extension (i.e., $M_c$ = $M_e$, q/p = constant). This is because the strength envelope of the MCC model forms a circle in the π-plane. Such predictions of the MCC model are evidently inaccurate, as numerous experiments have demonstrated that the tensile strength of soils is lower than their compressive strength [84,85]. In contrast, the proposed model successfully captures the difference between extensile and compressive strengths, indicating that the yield surface of the MCC model has been effectively extended to the three-dimensional stress space.

Furthermore, to simulate the anisotropic characteristics of the soil, fabric parameters are set to Δ = 0.3, β = 0.03, and $c_F$ = 2000, while keeping all other parameters unchanged. As previously mentioned, due to its inherent anisotropy, the mechanical response of soil is closely related to the direction of the applied load. To simulate this phenomenon, different deposition directions are considered, as shown in Figure 3. When the deposition plane is rotated, applying compressive loads along the z-axis effectively represents natural sedimentary soils subjected to loads from different directions. Based on the work of Tian [86], the relationship between the deposition direction δ and the fabric tensor can be defined as follows:

$$B_{ij}=\left(\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\delta & 0& -\mathrm{s}\mathrm{i}\mathrm{n}\delta \\ 0& 1& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\delta & 0& \mathrm{c}\mathrm{o}\mathrm{s}\delta \end{array}\right)$$

(29)

$$F_{ij}^{\prime }=B_{ik}{B}_{jl}{F}_{kl}=\left(\begin{array}{ccc}\left[\left(\Delta +1\right)+\left(3\Delta -1\right)\mathrm{c}\mathrm{o}\mathrm{s}2\delta \right]/4& 0& \left(3\Delta -1\right)\mathrm{s}\mathrm{i}\mathrm{n}2\delta /4\\ 0& (1-\Delta )/2& 0\\ \left(3\Delta -1\right)\mathrm{s}\mathrm{i}\mathrm{n}2\delta /4& 0& \left[\left(\Delta +1\right)-\left(3\Delta -1\right)\mathrm{c}\mathrm{o}\mathrm{s}2\delta \right]/4\end{array}\right),$$

(30)

where $B_{ij}$ denotes the rotation tensor. $F_{ij}^{\prime }$ denotes the fabric tensor after rotation.

Figure 4 shows the drained and undrained compression responses of the soil under different deposition directions. It can be found that the soil shows a strong directional dependence on its mechanical behavior such as failure strength, stress path, and volumetric deformation. At lower strain levels, as the deposition direction deviates from the z-axis, both the strength and volumetric deformation of the soil gradually decrease. This trend is consistent with the behavior observed in experimental studies [33,34,35,40,56,87]. Interestingly, as the strain level increases, the soil under each deposition direction tends to converge to the same critical state. This phenomenon aligns with the anisotropic critical state theory (ACST) proposed by Li and Dafalias [69], which suggests that anisotropic soils converge to a unique critical state at sufficiently large strain levels. This observation confirms the validity of the model.

The simulation results above demonstrate the validity of the model at the macroscopic level. Next, the accuracy of the model will be verified at the microscopic level under drained compression conditions. Extensive observations have demonstrated that fabric anisotropy evolves during the loading process, and the evolution patterns of soil anisotropy vary depending on the initial conditions. However, upon reaching the critical state, fabric anisotropy converges to a unique value that is independent of the initial conditions [42,68,88]. To verify whether the model implemented in FLAC^3D can accurately capture the characteristics of fabric evolution, a series of parameter sensitivity analyses were conducted in this study. As shown in Figure 5, the relationships between the evolution of fabric anisotropy A and various microscopic parameters are presented. Figure 5a (fabric parameters are set to Δ = 0.3, β = 0.03, and $c_F$ = 4000) shows that for deposition directions of 10° and 30°, the degree of anisotropy gradually decreases with increasing strain levels. In contrast, for deposition directions of 60° and 80°, anisotropy initially decreases and then increases. This behavior is attributed to particle rearrangement at the initial stage, which offsets the initial cross-anisotropy of the soil. As loading continues, the long axes of particles gradually align perpendicular to the loading direction, leading to a renewed increase in anisotropy. Li and Li [89] also observed a similar phenomenon. Despite the differences in evolutionary patterns, all tests predict the same critical state fabric, consistent with the assumptions of the ACST. Figure 5b shows the evolutions of A under different initial values of Δ (fabric parameters are set to β = 0.03 and $c_F$ = 4000). It can be observed that, despite varying initial degrees of anisotropy, all tests reach the same critical state fabric at large strain levels. Figure 5c shows the effect of parameter $c_F$ on the evolution of anisotropy (fabric parameters are set to Δ = 0.3 and β = 0.01). It is evident that, according to Equation (6), the rate of anisotropy evolution increases with the increase in parameter a and eventually tends towards a unique stable state. Figure 5d shows the effect of parameter β on the anisotropic behavior (fabric parameters are set to Δ = 0.3 and $c_F$ = 4000). As β increases, the final degree of anisotropy in the soil increases.

4.2. Embankment Loading Problem

To investigate the effects of anisotropy on typical geotechnical engineering problems, the proposed model is used to simulate the embankment loading case provided by the FLAC^3D [90]. The dimensions and boundary conditions of the model are shown in Figure 6. Consistent with the MCC model, the basic parameters of this model are set as M = 0.888, v = 0.3, λ = 0.162, κ = 0.062, and the reference void ratio $e_{ref}$ = 1.858 when p = 1 kPa. The density of the soil is set to 2 × 10³ kg/m, while the permeability is set to 1 × 10⁻⁸ m/s. The foundation is assumed to be lightly over-consolidated with $p_c^0$ = 160 kPa. The simulation consists of two stages: in the first stage, an embankment load of 50 kPa is applied to the 4 m section on the left side of the top boundary of the model, with no fluid flow during this stage, corresponding to undrained conditions; in the second stage, fluid flow is permitted to simulate approximately three years of consolidation, corresponding to drained conditions. Further details about the simulation can be referenced from the help documentation of FLAC^3D [90]. It should be noted that before starting the calculations, for the anisotropic model, the initial value of the hardening parameter ${\tilde{p}}_c$ needs to be determined using Equation (19).

As the first calculation, the fabric parameters are set to Δ = 1/3 and $c_F$ = 500, while β is varied from 0.02 and 0.08. The purpose of this setup is to make the initial state of this model equivalent to the initial isotropic state of the MCC model, thereby considering only the effects of anisotropic evolution. Figure 7 shows the mechanical responses of the saturated soil foundation predicted by the MCC model and this model. The predicted settlements versus time at the ground surface along the centerline of the embankment are presented in Figure 7a for approximately 3 years of consolidation time. It can be observed that the settlement predicted by the MCC model is approximately 0.193 m, while this model predicts a final settlement of 0.209 m (β = 0.02), and the predicted settlement increases with the parameter β. Interestingly, Karstunen et al. [50] and Yildiz et al. [51] have conducted similar comparative validations and found that considering anisotropy led to higher predicted final settlements compared to the MCC model. This highlights the necessity of accounting for anisotropy in practical engineering problems to avoid unrealistic predictions. Figure 7b,c present the surface settlement troughs during the undrained and drained stages, respectively. It can be observed that in the undrained stage, both the MCC model and the anisotropic model predict noticeable surface heave, with the anisotropic model forecasting a more pronounced degree of both surface heave and settlement. During the drained consolidation stage, the anisotropic model predicts higher settlement values. These trends align well with the conclusions drawn from the study by Karstunen et al. [50]. Figure 7d illustrates the variation in lateral displacement with depth at x = 4.5 m during the drainage phase. It can be observed that the anisotropic model predicts greater lateral displacement compared to the MCC model, a similar trend also identified by Yildiz et al. [51].

Figure 8 shows the total displacement contours predicted by the MCC model and this model (Δ = 1/3, β = 0.08, $c_F$ = 500). It can be observed that during both the undrained and drained stages, the anisotropic model predicts greater horizontal and lateral displacements compared to the MCC model. Given that the elastic behavior of both models is identical, these differences can be attributed to the anisotropic plastic behavior.

4.3. Tunnel Excavation Problem

The final case study focuses on a tunnel excavation simulation, using the Qinghuayuan Tunnel shield tunneling project as the background to analyze the impact of construction on ground disturbance [91]. A half model is established as shown in Figure 9, with a domain width of 70 m and a depth of 82 m, and the tunnel axis is located at a depth of 23.32 m with a radius of 6.32 m. Horizontal constraints are applied to the surrounding boundaries, while both horizontal and vertical constraints are applied to the bottom boundary. The top of the model is set as a free boundary. According to the work of Cui [91], the soil layers in this region mainly consist of silty clay and sandy pebbles, with sandy pebbles being the most representative. For simplicity, the region was simplified as a homogeneous pebble layer. The soil parameters were recorded as effective friction angle φ′ = 45°, Poisson’s ratio v = 0.28, unit weight = 20.2 kN/m³, and plastic index I_p = 11.6. Based on these parameters and the work of Nakase et al. [92], the basic parameters for the MCC model and this model can be approximated as follows:

$$M=6\mathrm{s}\mathrm{i}\mathrm{n}\phi \prime /\left(3-\mathrm{s}\mathrm{i}\mathrm{n}\phi \prime \right)\approx 1.85,$$

(31)

$$\lambda =0.02+0.0045I_p\approx 0.072, \kappa =0.00084\left(I_p-4.6\right)\approx 0.006.$$

(32)

The tunnel excavation process using the stress release method is carried out in two steps. First, the elements in the tunnel excavation area (excavation radius = 6.32 m) are assigned a null model, and the reverse forces are applied to the grid nodes of the tunnel excavation surrounding walls to simulate the stress release during the excavation process. In the second step, the elements in the lining area (the blue region in Figure 9) are assigned an elastic constitutive model, the lining elements are activated, and the reverse forces are removed to achieve equilibrium. The thickness of the lining area is 0.5 m. According to the work of Cui [91], Young’s modulus of the lining elements is set to 35.5 GPa, Poisson’s ratio is set to 0.25, and the unit weight is set to 25 kN/m³. It should be noted that during the geostress equilibrium stage, the fabric evolution parameters $c_F$ and β are set to 0 to ensure that the soil is undisturbed prior to excavation. Then, during the excavation phase, the parameters $c_F$ and β are assigned predetermined values to simulate the impact of the excavation process on the microstructure of the soil.

It is assumed that the reference void ratio $e_{ref}$ of the soil at p = 1 kPa is 1.094, and the OCR = 2.5. Using the empirical formula $K_0$ = (1 − sinφ′) × OCR^sinφ′ [93,94], the static lateral earth pressure coefficient $K_0$ is approximately 0.56. For comparison, the basic material parameters (i.e., M, λ, κ, v and $e_{ref}$) of the MCC model and the model in this study are kept consistent. To study the effect of anisotropy on soil disturbance, the fabric parameters of the model in this study are set as Δ = 0.3, $c_F$ = 150, and β varies between 0.09 and 0.11. Figure 10a shows the surface settlement trough observed in the experiment, predicted by the MCC model, and predicted by the model in this study. It can be seen that the MCC model underestimates the maximum settlement value, while the model in this study provides a better representation of the actual settlement. However, although fabric anisotropy has been introduced into the MCC model, it primarily affects the settlement depth without significantly altering the shape of the settlement basin, which differs from the shape observed in practice. Similar patterns were also noted by Karstunen et al. [50] and Yildiz et al. [51]. According to the work of Xiao et al. [23], this may be because the MCC model does not account for the difference in elastic modulus during loading and unloading. Figure 10b shows the horizontal displacements predicted by both models. It can be observed that the anisotropic model predicts larger horizontal displacements, which is consistent with the observations reported by Cui [91]. However, it is worth mentioning that in the anisotropic model proposed by Cui [91], a significant soil arching was observed around the soil surface at x = 20 m, which does not align with the experimental observations. This discrepancy may be due to the fact that the anisotropic model proposed by Cui [91] is established on the M-C model, which uses the ideal plasticity theory. This theory assumes that the soil does not exhibit subsequent hardening behavior, which may lead to an overestimation of the soil arching. Figure 11 shows the total displacement contour plots near the excavation area predicted by the MCC model and the model in this study. It can be observed that the displacement distribution predicted by both models is generally consistent. However, for the model in this study, as parameter β increases, the displacements near the crown and shoulders of the excavation area significantly increase. This indicates that the effect of anisotropy on the disturbance of the soil layers during excavation should not be overlooked.

In order to investigate the effect of the static lateral pressure coefficient $K_0$ on the soil disturbance caused by tunnel excavation under anisotropic conditions, a series of $K_0$ values are selected for the simulations ($K_0$ = 0.3, 0.5, 0.8, and 1.0). Figure 12 shows the simulation results of this model with constant OCR = 2.5 but different $K_0$, while the fabric parameters are set as Δ = 1/3, $c_F$ = 150, and β = 0.09. From Figure 12a, it can be seen that as $K_0$ increases, the vertical settlement at the soil surface decreases. Figure 12b shows the normalized surface settlement under different $K_0$ conditions. It can be observed that $K_0$ not only affects the depth of the settlement trough but also influences its width. Specifically, lower values of $K_0$ predict narrower settlement troughs. This could be due to the fact that as $K_0$ decreases, the mean effective stress p in the soil decreases, and, as can be seen from Equation (17), the stiffness of the soil also decreases accordingly. This trend is consistent with the observations made by Cui [91] and Dolezalova [95]. Figure 13 shows the total displacement contours under different $K_0$ conditions. It can be observed that as $K_0$ increases, the soil displacement near the excavation area significantly decreases. Therefore, it can be concluded that the resistance of soil to disturbance during tunnel excavation is sensitive to $K_0$, and the value of $K_0$ should be determined as accurately as possible during site investigation.

The effect of OCR on tunnel excavation has always been an interesting topic of geotechnical engineering [96,97,98]. To investigate the effect of OCR on tunnel excavation, this study sets OCR of 1.6, 1.9, 2.2 and 2.5 while keeping $K_0$ = 0.56, and the fabric parameters are set as Δ = 1/3, $c_F$ = 150, and β = 0.09. Figure 14 shows the surface settlement and horizontal displacement under different OCR conditions. It can be observed that as OCR increases, both the settlement and horizontal displacement at the surface decrease. This is consistent with the findings observed by Shi et al. [96] in their centrifuge tests. This may be because, as OCR increases, the size of the initial yield surface increases, making it more difficult for the soil to reach the plastic stage, resulting in higher stiffness. The plastic zone contour map in Figure 15 (blue region) confirms this conclusion, showing that as OCR increases, the predicted plastic zone area gradually decreases. Furthermore, it can also be found that the plastic zones predicted by the model are concentrated at the crown and shoulders of the excavation area, which is similar to the plastic zone distribution predicted by Cui [91].

In order to investigate the effect of initial anisotropy on tunnel excavation, this study set up two sets of numerical tests. Both sets of experiments have the same OCR = 2.5 and $K_0$ = 0.56, and the fabric parameters are set as $c_F$ = 150 and β = 0.09. In one set of tests, the angle δ between the deposition plane and the z-axis is varied between 0° and 90° with Δ = 0.3. In the other set, the value of initial anisotropy Δ is varied between 0.2 and 0.32 while δ = 0°. Figure 16a shows the relationship between settlement at various positions and δ. It can be observed that as δ increases, the settlement at each position (surface, shoulder, crown) first increases, reaching a maximum between approximately δ = 60 and 75°, and then decreases as δ continues to increase. Figure 16b shows the relationship between horizontal displacement at the springline and δ. It can be observed that horizontal displacement increases as δ increases. This may be because, as the deposition direction gradually deviates from the z-axis, the stiffness of the soil decreases, which leads to an initial increase in vertical settlement. As δ continues to increase, the horizontal displacement at the springline continues to grow, and due to compression, the vertical settlement starts to decrease. Figure 17 shows the distribution of horizontal displacement near the excavation area under different δ conditions. It can be observed that for different scenarios, the maximum horizontal displacement is concentrated at the springline, and as δ increases, the maximum horizontal displacement increases.

Figure 18 illustrates the relationship between vertical settlement and horizontal displacement at various locations with respect to the initial anisotropy parameter Δ. It is observed that as Δ decreases, the vertical settlement at the surface, crown, and shoulder positions decreases. However, when Δ reduces from 0.22 to 0.2, the vertical settlement slightly increases. On the other hand, the horizontal displacement at the haunch position decreases monotonically with decreasing Δ. This behavior may be attributed to the fact that as Δ decreases, the inclination of the initial yield surface increases, resulting in different distributions of plastic regions. Figure 19 corroborates this finding, showing that with a reduction in Δ, the distribution of plastic regions in the excavation area shifts from the shoulder position toward the crown position. A similar phenomenon was observed by Liu et al. [99] in their DEM experiments when they found that the disturbed zone caused by tunnel excavation rotates as the initial anisotropy changes.

To investigate the influence of parameter $c_F$ on tunnel excavation, the fabric parameters are set as Δ = 0.25, β = 0.09, and $c_F$ values of 10, 100, 1000, and 5000, with OCR = 1.6 and $K_0$ = 0.56. Figure 20 presents the predicted maximum surface settlement and maximum horizontal displacement for different $c_F$ values. It can be observed that as $c_F$ increases, the predicted maximum settlement and maximum horizontal displacement also increase. However, compared to factors such as OCR and K₀, the influence of $c_F$ on the model predictions is relatively limited. This may be attributed to the fact that the preferential orientation of soil particles changes only minimally during excavation.

5. Conclusions

In order to incorporate the effect of the fabric anisotropy and three-dimensional strength characteristics of the soil into the analysis of a geotechnical engineering problem, this study proposed a modified fabric-based MCC model based on the work of Tian [65]. Compared to the original model proposed by Tian [65], the improvements in the modified model mainly lie in the fabric evolution criterion and strength criterion. These modifications aim to simplify the equations for easier numerical implementation. Then, the elastic trial–plastic correction method was used to numerically implement the model, and it was subsequently integrated into the FLAC^3D software. On this basis, the model was used to simulate the mechanical response of soil under triaxial tests, embankment loading tests as well as tunnel excavation tests, leading to the following conclusions:

(1)

In the analysis of geotechnical engineering problems, incorporating the fabric anisotropy and three-dimensional strength of the soil is essential, as these factors significantly affect the mechanical response of the soil.

(2)

The validity of the model presented in this paper has been validated at both the macroscopic and microscopic levels. The model can reasonably reflect the anisotropic characteristics of the soil. Furthermore, at higher strain levels, the anisotropic soil can converge to a unique critical state, which is consistent with the ACST proposed by Li and Dafalias [69].

(3)

The preliminary simulation results of the model provide a reference for its practical application. Although the cases used in the study are relatively simple, the model effectively reflects the impact of factors such as anisotropy, over-consolidation, and the coefficient of lateral earth pressure on soil disturbance. This can serve as a valuable reference for the design and construction of real-world engineering projects.

However, it must be pointed out that accurately calibrating the microscale parameters related to fabric anisotropy remains a significant challenge when using the model. Although some feasible methods for calibrating microscale parameters through laboratory experiments have been proposed [58,100], the samples used in these experiments are generally remolded soils. Moreover, such methods often require advanced experimental equipment, such as electron microscopes or CT scanners, which are not commonly available in standard soil mechanics laboratories. Currently, the precise determination of in situ soil fabric parameters using simple and practical methods remains a major challenge.