Application of Dimension Extending Technique to Unified Hardening Model

Abstract

This paper provides the process of incremental constitutive integration for the unified hardening model combined with the transformation stress method. The dimension-extending technique takes the hardening function of the hardening/softening model as the same position as the stress components, so that the constitutive integration of the plasticity can be reduced to an initial value problem of differential–complementarity equations, which is solved using the Gauss–Seidel algorithm-based Projection–Correction for the mixed complementarity problem. The Gauss–Seidel based Projection–Correction algorithm does not require the calculation of the Jacobean matrix of the potential function, making it relatively easy to implement in programming. The unified hardening model is proposed based on the modified Cam–Clay model and the sub-loading surface model, and the elastic properties are pressure-dependent. Two processing methods, backward Euler integration and exact elastic property, are used for the variable elasticity properties. The constitutive integration of the increased dimensional unified hardening model is reduced to a special mixed complementarity problem and solved by the proposed algorithm, which does not need to calculate the Jacobean matrix of the potential function, and greatly simplifies the derivation process. Several numerical examples are given to verify the feasibility of the incremental constitutive integration in the unified hardening model, including the single integral point and the boundary value problems. The research results have expanded the scope of use of the Gauss–Seidel based Projection–Correction algorithm.

1. Introduction

The modified Cam–Clay model (MCC) [1] is a widely used elasto-plastic constitutive model [2], which can reflect the friction, compression and shear shrinkage of soil [3]. However, it can only reflect the characteristics of normal consolidated or light over-consolidated clay, and cannot reasonably reflect the dilatancy of heavy over-consolidated clay [4]. Based on critical state theory, various constitutive models have been proposed [5,6,7]. The unified hardening (UH) model is one of them. The UH model is proposed by Yao et al. in the framework of the MCC model and combined with the concept of the sub-loading surface model [8], which can reasonably reflect the hardening, softening, dilatancy, and shrinkage characteristics of soil. The constitutive equation of the elasto-plastic model is generally expressed in the rate form [9]. The numerical algorithm for updating the stress of the elasto-plastic constitutive model for geotechnical materials requires integration of the rate form constitutive equation along the strain increment path within a finite time step. This is also one of the core issues in elasto-plastic numerical calculations, which directly affects the accuracy of the calculation results [10]. The common numerical algorithms of stress updating can be roughly divided into two kinds: the explicit algorithm (the forward Euler method) and the implicit algorithm (the backward Euler method). Using an explicit algorithm, if the entire strain increment includes both elastic and elastic–plastic behavior, it is necessary to find the intersection point of the stress increment vector and the yield surface [11]. The gradients of the yield surface and plastic potential surface are calculated based on the stress state at the starting point of the incremental step. The procedure is easy to implement, but the computational accuracy is low, and the solution is prone to drift from the yield surface. Implicit algorithms generally include two parts: elastic prediction and plastic correction, which require obtaining the gradients of yield and plastic potential surfaces based on unknown stresses. Therefore, an iterative method is needed. Although iterative calculation improves computational accuracy, it also increases the difficulty of derivation. The closest point projection method (CPPM) [12] is currently a widely used implicit algorithm with good computational stability. However, the Hessian matrix is needed in the algorithm. For some complex constitutive models or those that introduce stress Lode angle, the theoretical derivation process is relatively cumbersome, especially for highly non-linear constitutive models, which are prone to the singularity of the Jacobean matrix and non-convergence during calculation. Bicanic et al. [13] improved the initial iteration values of the algorithm using the auxiliary projection surface method, but how to effectively construct the auxiliary projection surface has not yet been well addressed. To avoid calculating the Hessian matrix, Simo et al. [14,15] proposed the cutting-plane algorithm (CPA). Although the simplicity of this algorithm is highly attractive for large-scale calculation, the precise linearization cannot be obtained in a closed form [16]. Regarding the comparison results of these two algorithms from theory to numerical computation, Huang et al. [17] have conducted a more detailed analysis.

The rate form elasto-plastic constitutive model defines an initial value problem of an ordinary differential equations (ODEs) about the stress tensor $\mathit{\sigma }$ changing in the elastic domain $E_{\alpha }$. For the perfect-plastic model, $E_{\alpha }$ is invariant in the stress space. The Kuhn–Tucker complementarity condition and the consistency condition provide constraints for solving the ODEs. Based on the above analysis, the constitutive integration of the perfect-plastic model can be reduced to a mixed complementarity problem (MiCP), which is also a finite-dimensional variational inequality (VI). He [18] developed a projection-contraction algorithm for this finite-dimensional VI, designated by the acronym PCA. However, the convergence condition of the PCA is rather less stringent to the relevant functions, with no need of gradients for any functions. Tapping into the idea of the Gauss–Seidel iteration, Zheng et al. [19] proposed the Gauss–Seidel based Projection–Correction (GSPC) algorithm for solving the MiCP. Applying the monotonicity of the mapping of the MiCP, GSPC is proved convergent theoretically for associative plasticity. For non-associative plasticity, the sufficient condition for GSPC to be convergent is also established if the tension part of the Mohr–Coulomb elastic domain is cut off. GSPC converges for any size of strain increment, and its numerical stability is superior to traditional return mapping algorithms [19]. However, the GSPC algorithm is only applicable to perfect-plastic models, which is a significant deficiency of the algorithm. The proposal of the dimension-extending technique [20] can successfully expand the application of GSPC to a hardening/softening model. In this study, the dimension-extending technique is applied to the UH model, and the transformation stress method based on the SMP (Spatially Mobilized Plane) [21] strength criterion is used to extend a 2D model to 3D space. At the same time, two methods for pressure-related variable elasticity treatment were discussed, Euler integral and precise integral.

2. Introduction to the UH Model and Transformed Stress

2.1. UH Model

Concerning the framework of the MCC model and the concept of the sub-loading surface, Yao et al. proposed the UH model, which can reflect the gradual weakening trend of the over-consolidation ratio of over-consolidated soil throughout the entire stress change process [7]. It includes two yield surfaces, the reference yield surface and the current yield surface. The reference yield surface, as shown in MCC, takes the plastic volumetric strain as the hardening parameter, and the current yield surface has the same form as the reference yield function, as shown in Figure 1.

Reference yield surface equation:

$$\overline{F}=\ln \frac{\overline{p}}{{\overline{p}}_{x0}}+\ln (1+\frac{{\overline{q}}^2}{M^2{\overline{p}}^2})-\frac{1}{c_p}{\epsilon }_v^p=0$$

(1)

Current yield surface equation:

$$F=\ln \frac{p}{p_{x0}}+\ln (1+\frac{q^2}{M^2{p}^2})-\frac{1}{c_p}H=0$$

(2)

Hardening parameter (H):

$$H={\displaystyle \int \frac{M_f^4-{\eta }^4}{M^4-{\eta }^4}}d{\epsilon }_v^p$$

(3)

with:

$$M_f=6\left[\sqrt{\frac{k}{R}(1+\frac{k}{R})}-\frac{k}{R}\right]$$

(4)

$$R=\frac{p}{{\overline{p}}_{x0}}(1+\frac{{\eta }^2}{M^2})exp(-\frac{{\epsilon }_v^p}{c_p})$$

(5)

$$k=\frac{M^2}{12(3-M)}$$

(6)

$$c_p=\frac{\lambda -\kappa }{1+e_0}$$

(7)

M is the critical state stress ratio; ${\epsilon }_v^p$ indicates the plastic volumetric strain; $\eta =q/p$ is the stress ratio. $e_0$ is the initial void ratio; $\overline{p}$ and $\overline{q}$ define the mean stress and generalized shear stress on the reference yield surface. ${\overline{p}}_{x0}$ and $p_{x0}$ are the intersection points of the initial reference yield surface and the initial current yield surface with the p-axis, respectively. λ is the slope of the normal compression line, and κ is the slope of the unloading line.

The UH model refers to both normally-consolidated and over-consolidated soils. For normally-consolidated soil, it can degenerate into the MCC model. It can reasonably describe the characteristics of the over-consolidated soil with shear expansion, shear shrinkage, hardening/softening, and stress path dependence [22].

2.2. Transformed Stress Method

The transformed stress method [23] extends the two-dimensional constitutive model to three-dimensions. It converts the strength criterion in the ordinary stress space into the extended Mises strength criterion in the transformed stress space without changing the form of the yield function, as shown in Figure 2. Yao et al. proposed the transformed stress method based on strength theories, such as SMP criterion [24], and Lade criterion [25].

The transformed stress method based on the SMP strength criterion is used to extend the UH model to a three-dimensional model, and the corresponding transformation stress tensor is given by:

$${\tilde{\sigma }}_{ij}=\left\{\begin{array}{cc}p{\delta }_{ij}+\frac{q_c}{q}({\sigma }_{ij}-p{\delta }_{ij})& q\ne 0\hfill \\ {\sigma }_{ij}& q=0\hfill \end{array}\right.$$

(8)

$$q_c=\frac{2I_1}{3\sqrt{(I_1{I}_2-I_3)/(I_1{I}_2-9I_3)}-1}$$

(9)

${\tilde{\sigma }}_{ij}$ is the stress tensor in the transformed stress space, and the characters with ‘~’ atop thereafter denote the variables in the transformed stress space. ${\delta }_{ij}$ is the Kronecker’s delta; $I_1$, $I_2$ and $I_3$ are the first, second and third invariants, respectively.

$$\left\{\begin{array}{l}{I}_1={\sigma }_1+{\sigma }_2+{\sigma }_3\hfill \\ I_2={\sigma }_1{\sigma }_2+{\sigma }_2{\sigma }_3+{\sigma }_3{\sigma }_1\hfill \\ I_3={\sigma }_1{\sigma }_2{\sigma }_3\hfill \end{array}\right.$$

(10)

In the transformed stress space, the equations of the current yield surface and the reference yield surface are given by:

$$\tilde{F}=c_p\ln \frac{\tilde{p}}{{\tilde{p}}_{x0}}+c_p\ln (1+\frac{{\tilde{q}}^2}{M^2{\tilde{p}}^2})-\tilde{H}=0$$

(11)

$$\tilde{\overline{F}}=c_p\ln \frac{\tilde{\overline{p}}}{{\tilde{\overline{p}}}_{x0}}+c_p\ln (1+\frac{{\tilde{\overline{q}}}^2}{M^2{\tilde{\overline{p}}}^2})-{\epsilon }_v^p=0$$

(12)

with

$$\tilde{H}={\displaystyle \int \frac{{\tilde{M}}_f^4-{\tilde{\eta }}^4}{M^4-{\tilde{\eta }}^4}}d{\epsilon }_v^p$$

(13)

${\tilde{M}}_f$ signifies the potential failure stress ratio, varying with the changes in the current stress point, and $\tilde{\eta }$ is the stress ratio under the transformed stress space.

$${\tilde{M}}_f=6\left[\sqrt{\frac{k}{\tilde{R}}(1+\frac{k}{\tilde{R}})}-\frac{k}{\tilde{R}}\right]$$

(14)

$$\tilde{R}=\frac{\tilde{p}}{{\tilde{\overline{p}}}_{x0}}(1+\frac{{\tilde{\eta }}^2}{M^2})exp(-\frac{{\epsilon }_v^p}{c_p})$$

(15)

In the transformed stress space, the associated flow rule is adopted:

$$d{\mathit{\epsilon }}^p=\Lambda \frac{\partial \tilde{F}}{\partial \tilde{\mathit{\sigma }}}$$

(16)

where $d{\mathit{\epsilon }}^p$ is the plastic strain increment and $\Lambda$ is the plastic multiplier.

From Equation (13) it can be observed that, when the stress approaches the limit state, it will cause the denominator in the equation to be zero and the calculation will fail. For the stability of numerical calculations, the hardening Equation (13) is now rewritten in incremental form and organized as:

$$d\tilde{H}=\frac{{\tilde{M}}_f^4-{\tilde{\eta }}^4}{M^4-{\tilde{\eta }}^4}d{\epsilon }_v^p=\Lambda \frac{{\tilde{M}}_f^4-{\tilde{\eta }}^4}{M^4-{\tilde{\eta }}^4}\frac{\partial \tilde{F}}{\partial \tilde{p}}=\Lambda c_p\frac{1}{\tilde{p}}\frac{{\tilde{M}}_f^4-{\tilde{\eta }}^4}{{(M^2+{\tilde{\eta }}^2)}^2}$$

(17)

In terms of format, the transformed stress method adopts the form of the associated flow rule, but in fact it adopts the non-associated flow rule, which is more in line with the actual situation of the soil.

3. Incremental Constitutive Integration

3.1. Mixed Complementarity Problems

Complementarity: for vectors $(x,y)\in {ℝ}^n\times {ℝ}^n$, if they satisfy the following conditions, then they are called complementarity conditions.

$$0\le \mathit{x}\perp \mathit{y}\ge 0$$

(18)

A problem that includes both complementarity conditions and equality constraints is called an MiCP. Further, let the set $\Omega ={ℝ}_{+}^m\times {ℝ}^n$, two mapping operators ${\mathit{f}}_I$ and ${\mathit{f}}_E$ can be defined on the $\Omega$. ${\mathit{f}}_I$ is the complementarity condition, and ${\mathit{f}}_E$ is the newly introduced function related to the equality condition. Then, the MiCP can be written in the following form to find a set of vectors $(x,y)\in {ℝ}_{+}^m\times {ℝ}^n$ such that the following conditions are satisfied:

$$\left\{\begin{array}{l}{\mathit{f}}_E(x,y)=0;y\hspace{1em}free\hfill \\ 0\le \mathit{x}\perp {\mathit{f}}_I(x,y)\ge 0 \hfill \end{array}\right.$$

(19)

MiCP $({\mathit{f}}_I,{\mathit{f}}_E)$ is a special case of $VI(\Omega ,\mathit{f})$ [26], with $\Omega$ being a convex set. $\mathit{f}$ is defined as:

$$\mathit{f}(x,y)=\left(\begin{array}{c}{\mathit{f}}_E(x,y)\\ {\mathit{f}}_I(x,y)\end{array}\right)$$

(20)

3.2. GSPC Algorithm

He [18] designed the PC algorithm to solve VI $(\Omega ,f)$; if Ω is convex and f is monotonic, then PC converges to the solution of VI $(\Omega ,f)$. MiCP can be simplified to a Nonlinear Complementarity Problem (NCP), if g satisfies the following relationship.

$$g(\mathit{x})=f_I[\mathit{x},y(\mathit{x})]$$

(21)

The equation y(x) is the implicit function obtained by solving ${\mathit{f}}_E(\mathit{x},\mathit{y})=0$. If f_E has a strong non-linearity, we still consider the definition of the variational inequality VI $(\Omega ,f)$ in the form of MiCP $({\mathit{f}}_I,{\mathit{f}}_E)$. If MiCP is simplified to NCP, this can be directly solved using the PC algorithm. However, in the field of geotechnical engineering, f_I and f_E often have completely different dimensions and magnitudes, which necessitates a large number of iterative calculations to find the appropriate adjustment factors in the algorithm. To improve the robustness of the algorithm, Zheng proposed the GSPC algorithm and mathematically demonstrated the sufficient conditions for the convergence of GSPC, which is that the Jacobean matrix J is positive semi-definite. The calculation of the Jacobean matrix J is as follows:

$$\mathit{J}=\frac{\partial \mathit{f}(\mathit{x},\mathit{y})}{\partial (\mathit{x},\mathit{y})}\equiv \left[\begin{array}{cc}\frac{\partial {\mathit{f}}_I}{\partial \mathit{x}}& \frac{\partial {\mathit{f}}_I}{\partial \mathit{y}}\\ \frac{\partial {\mathit{f}}_E}{\partial \mathit{x}}& \frac{\partial {\mathit{f}}_E}{\partial \mathit{y}}\end{array}\right]$$

(22)

GSPC can be invoked in this way:

$$(x,y)=\mathrm{G}\mathrm{S}\mathrm{P}\mathrm{C}(x_0,y_0,{\mathrm{F}}_{\mathrm{I}},{\mathrm{F}}_{\mathrm{E}})$$

(23)

where $(x_0,y_0)$ is an initial guess of (x, y) and ${\mathrm{F}}_{\mathrm{I}},{\mathrm{F}}_{\mathrm{E}}$ are function handles of ${\mathit{f}}_I$ and ${\mathit{f}}_E$, respectively. The pseudocode of GSPC is shown in Appendix A.

Many problems in engineering can be reduced to MiCP, such as seepage flow with free surfaces, contact mechanics of multiple blocks [27], and elasto-plastic [28] et al. In this study, we are only concerned with plasticity problems, in which the constitutive integration of plasticity is reduced to a MiCP, a special case of finite-dimensional VI.

3.3. The Dimension-Extending Technique and Constitutive Integration

For hardening/softening behaviors, the yield surface in the stress space will shrink or expand with the development of plastic deformation. The hardening parameter $\varpi$ in the yield equation $\phi ({\mathsf{\sigma }}_{ij},\varpi )$ characterizes this change. $\varpi$ is a vector form with a dimension greater than one and is suited to complex hardening behaviors, such as the mixed hardening model. The current model adopts a scalar form. The dimension-extending technique regards $\varpi$ as the same position as the components of σ, and then the yield function in two-dimensional space, principal stress space, or general stress space is expanded to the “increasing dimensional stress” space. Taking the MCC model as an example, the yield function of the increased dimensionality MCC model is as follows:

$$\stackrel{⌣}{F}(p,q;p_c)=q^2+M^{2}p(p-p_c)$$

(24)

$p$ and $q$ are the mean stress and generalized shear stress, respectively.

The elastic domain enclosed by the yield function in the $(p,q;p_c)$ space thus turns into:

$$\stackrel{⌣}{E}\triangleq \left\{(p,q;p_c)\in R^3|\stackrel{⌣}{F}(p,q;p_c)\le 0\right\}$$

(25)

In the “increasing dimensional stress” space, the stress points of the MCC model can only be on or in the $\stackrel{⌣}{E}$ without jumping out of it, as shown in Figure 3.

The constitutive equation of the model is generally given in the rate form. If we want to find the stress–strain relationship, we need to integrate the rate equations. The definition of the constitutive integration is to calculate the stress state at time $n+1$ by knowing the stress state ${\mathit{\sigma }}_n$ and strain state ${\mathit{\epsilon }}_n$ at the time $n$, as well as the strain increment $\Delta \mathit{\epsilon }$ at this time step [29].

The constitutive integration of the elasto-plastic model can be reduced to the MiCP mentioned earlier, such as the Mohr–Coulomb (MC) model [19] and the MCC model [20]. The UH model is similar to the MCC model, so the hardening parameters can also be extended to the “increasing dimensional stress” space. The elastic properties of the model are pressure-dependent. $K$ and $G$ are the bulk modulus and shear modulus, respectively.

$$K={\overline{K}}_{0}p\hspace{1em}\hspace{1em}G={\overline{G}}_{0}p$$

(26)

$${\overline{K}}_0=\frac{(1+e_0)}{\kappa }\hspace{1em}\hspace{1em}{\overline{G}}_0=\frac{3(1-2v)}{2(1+v)}{\overline{K}}_0$$

(27)

$v$ is the Poisson ratio, and $e_0$ is the initial void ratio.

There are two methods to deal with this variable elastic property; one is Euler integral approximation, and the other is precise integration [9]. This study mainly uses the two methods to discuss the constitutive integration of the UH model in the “increasing dimensional stress space”.

3.3.1. Description of MiCP with Approximate Elastic Properties (AEP)

The rate form constitutive equations of the UH model in the “increasing dimensional stress” space can be expressed as:

$$\left\{\begin{array}{l}\dot{\mathit{\sigma }}={\dot{\mathit{\sigma }}}_e-{\mathit{\sigma }}_p\hfill \\ \dot{\tilde{H}}=\Lambda c_p\frac{1}{\tilde{p}}\frac{{\tilde{M}}_f^4 - {\tilde{\eta }}^4}{{(M^2 + {\tilde{\eta }}^2)}^2}\hfill \end{array}\right.$$

(28)

$${\dot{\mathit{\sigma }}}_e=\mathit{D}\dot{\mathit{\epsilon }}$$

(29)

$${\mathit{\sigma }}_p\triangleq \mathit{D}{\dot{\mathit{\epsilon }}}_p$$

(30)

$${\dot{\mathit{\epsilon }}}_p\triangleq \Lambda {\nabla }_{\tilde{\mathit{\sigma }}}\tilde{F}$$

(31)

where ${\nabla }_{\tilde{\mathit{\sigma }}}\tilde{\mathrm{F}}$ represents the partial derivative of the yield function corresponding to the stress in the transformed stress space. $\mathit{D}$ is the symmetric positive elastic stiffness matrix:

$$\mathit{D}\triangleq \frac{3K(1-\nu )}{1+\nu }\left[\begin{array}{ccc}1& \frac{\nu }{1 - \nu }& \frac{\nu }{1 - \nu }\\ \frac{\nu }{1 - \nu }& 1& \frac{\nu }{1 - \nu }\\ \frac{\nu }{1 - \nu }& \frac{\nu }{1 - \nu }& 1\end{array}\right]$$

(32)

$K$ is defined as in Equation (26). $v$ is the specific volume.

$$\frac{\partial \tilde{F}}{\partial {\tilde{\sigma }}_i}=\frac{\partial \tilde{F}}{\partial \tilde{p}}\frac{\partial \tilde{p}}{\partial {\tilde{\sigma }}_i}+\frac{\partial \tilde{F}}{\partial \tilde{q}}\frac{\partial \tilde{q}}{\partial {\tilde{\sigma }}_i}\hspace{1em}i=1,2,3$$

(33)

$$\frac{\partial \tilde{F}}{\partial \tilde{p}}=c_p\frac{1}{\tilde{p}}\frac{M^2{\tilde{p}}^2-{\tilde{q}}^2}{M^2{\tilde{p}}^2+{\tilde{q}}^2}$$

(34)

$$\frac{\partial \tilde{F}}{\partial \tilde{q}}=c_p\frac{2\tilde{q}}{M^2{\tilde{p}}^2+{\tilde{q}}^2}$$

(35)

$$\frac{\partial \tilde{p}}{\partial {\tilde{\sigma }}_i}=\frac{1}{3}\hspace{1em}i=1,2,3$$

(36)

$$\frac{\partial \tilde{q}}{\partial {\tilde{\sigma }}_i}=\frac{3({\tilde{\sigma }}_i-\tilde{p})}{2\tilde{q}}\hspace{1em}i=1,2,3$$

(37)

Using the backward Euler integration on interval $(n,n+1)$ for Equation (28), then:

$$\left\{\begin{array}{l}{\mathit{\sigma }}_{n+1}={\mathit{\sigma }}_e-{\mathit{\sigma }}_p\hfill \\ {\tilde{H}}_{n+1}={\tilde{H}}_n+\Delta \tilde{H}\hfill \end{array}\right.$$

(38)

$${\mathit{\sigma }}_e\triangleq {\mathit{\sigma }}_n+\mathit{D}\Delta \mathit{\epsilon }$$

(39)

$$\Delta \tilde{H}=\Lambda c_p\frac{1}{{\tilde{p}}_{n+1}}\frac{{\tilde{M}}_{f,n+1}^4-{\tilde{\eta }}_{n+1}^4}{{(M^2+{\tilde{\eta }}_{n+1}^2)}^2}$$

(40)

${\mathit{\sigma }}_p$ is still calculated by Equation (30).

$f_I:R_{+}\times R^4\to R$

$$f_I(\Lambda ,\tilde{\mathit{\sigma }})\triangleq -\tilde{F}(\tilde{\mathit{\sigma }})$$

(41)

${\mathit{f}}_E:R_{+}\times R^4\to R^4$

$${\mathit{f}}_E(\Lambda ,\widehat{\mathit{\sigma }})\triangleq \left\{\begin{array}{l}{\mathit{\sigma }}_{n+1}+{\mathit{\sigma }}_p-{\mathit{\sigma }}_e\hfill \\ {\tilde{H}}_{n+1}-{\tilde{H}}_n-\Delta \tilde{H}\hfill \end{array}\right.$$

(42)

with

$${\widehat{\mathit{\sigma }}}^T=({\sigma }_1,{\sigma }_2,{\sigma }_3,\tilde{H})$$

(43)

$${\tilde{\mathit{\sigma }}}^T=({\tilde{\sigma }}_1,{\tilde{\sigma }}_2,{\tilde{\sigma }}_3,\tilde{H})$$

(44)

Equations (41) and (42) form $f_I$ and ${\mathit{f}}_E$ in MiCP.

The consistency condition is:

$$\Lambda >0,\hspace{1em}-\tilde{F}=0,\hspace{1em}\Lambda (-\tilde{F})=0$$

(45)

The yield function of the UH model corresponds one-to-one to the current stress. In the “increasing dimensional stress” space, the stress point can only be located on the surface:

$$\tilde{E}\triangleq \left\{({\tilde{\sigma }}_1,{\tilde{\sigma }}_2,{\tilde{\sigma }}_3;\tilde{H})\in R^4|\tilde{F}({\tilde{\sigma }}_1,{\tilde{\sigma }}_2,{\tilde{\sigma }}_3;\tilde{H})=0\right\}$$

(46)

Therefore, the MiCP is adjusted to the following form:

$$\left\{\begin{array}{l}{\mathit{f}}_E(\Lambda ,\widehat{\mathit{\sigma }})=0;\widehat{\mathit{\sigma }}\hspace{1em}free\hfill \\ \begin{array}{l}\Lambda f_I(\Lambda ,\tilde{\mathit{\sigma }})=0,\hspace{1em}if\hspace{1em}{f}_I(\Lambda ,\tilde{\mathit{\sigma }})=0 \\ \Lambda >0\end{array}\hfill \end{array}\right.$$

(47)

Equation (47) only includes consistency conditions, therefore the unloading/loading criterion [16] included in the traditional complementary relationship needs to be adjusted. According to the processing method in reference [30], the over-consolidation parameter R can be used to determine whether the behavior material is in the stage of plastic loading or elastic unloading.

If the behavior is in the elastic section, Hooke’s law is needed. If it enters the plastic section, GSPC needs to be invoked for calculation.

Suppose the absence of plastic strain, then the trial stress is given by:

$${\mathit{\sigma }}_{n+1}^{tr}={\mathit{\sigma }}_e$$

(48)

${\mathit{\sigma }}_e$ can be obtained from Equation (39).

The unloading/loading criterion is then:

$$R_n=\frac{p_n}{{\overline{p}}_{x0}}(1+\frac{{\eta }_n^2}{M^2})exp(-\frac{{\epsilon }_{v,n}^p}{c_p})$$

(49)

$$R_{n+1,tr}=\frac{p_{n+1,tr}}{{\overline{p}}_{x0}}(1+\frac{{\eta }_{n+1,tr}^2}{M^2})exp(-\frac{{\epsilon }_{v,n}^p}{c_p})$$

(50)

$p_{n+1,tr}$ and ${\eta }_{n+1,tr}$ are the mean stress and stress ratio of the trial stress, respectively.

If $R_n>R_{n+1,tr}$, the soil is beneath the unloading state. However, if $R_n\le R_{n+1,tr}$ the soil is beneath the loading or the neutral state yields elasto-plastic strain and GSPC is needed.

If in the transformed stress space, then

$${\tilde{R}}_n=\frac{{\tilde{p}}_n}{{\tilde{\overline{p}}}_{x0}}(1+\frac{{\tilde{\eta }}_n^2}{M^2})exp(-\frac{{\epsilon }_{v,n}^p}{c_p})$$

(51)

$${\tilde{R}}_{n+1,tr}=\frac{{\tilde{p}}_{n+1,tr}}{{\tilde{\overline{p}}}_{x0}}(1+\frac{{\tilde{\eta }}_{n+1,tr}^2}{M^2})exp(-\frac{{\epsilon }_{v,n}^p}{c_p})$$

(52)

If ${\tilde{R}}_n>{\tilde{R}}_{n+1,tr}$, the soil is in elastic state and presents elastic behavior abiding to Hooke’s law.

If ${\tilde{R}}_n\le {\tilde{R}}_{n+1,tr}$, GSPC will be called.

3.3.2. The MiCP with Exact Elastic Properties (EEP)

An alternative to deal with the rate form of the constitutive equation is to precisely integrate the equation. If the plastic strain increment is frozen, considering that the strain increment only causes pure elastic behavior, the elastic trail stress can be obtained as:

$$p^{n+1,tr}=p^n\exp ({\overline{K}}_0\Delta {\epsilon }_{kk})$$

(53)

$$s_{ij}^{n+1,tr}=s_{ij}^n+2G_{ave}\Delta e_{ij}$$

(54)

$$G_{ave}=\frac{{\overline{G}}_0{p}^n(\exp ({\overline{K}}_0\Delta {\epsilon }_{kk})-1)}{{\overline{K}}_0\Delta {\epsilon }_{kk}}$$

(55)

${\overline{K}}_0$ and ${\overline{G}}_0$ are defined in Equation (27). $s_{ij}^n$ and $s_{ij}^{n+1,tr}$ are the deviatoric stress at time n and the trial deviatoric stress at time n + 1, respectively. $\Delta {\epsilon }_{kk}$ and $\Delta e_{ij}$ are the incremental volumetric and deviatoric strains, respectively.

Combining Equations (53) and (54), the stress tensor, which is in the original σ-space, is derived by Equation (56):

$${\sigma }_{ij}^{n+1,tr}=s_{ij}^{n+1,tr}+p^{n+1,tr}{\delta }_{ij}$$

(56)

According to the unloading/loading criterion in Section 3.3.1, determine whether GSPC needs to be invoked.

$$\Delta {\epsilon }_{ij}^p=\Lambda {\tilde{r}}_{ij}^{n+1}$$

(57)

$${\tilde{r}}_{ij}=\frac{\partial f}{\partial \tilde{\sigma }}$$

(58)

By combining the precise integral, we can obtain:

$$p^{n+1}=p^n\exp [{\overline{K}}_0(\Delta {\epsilon }_{kk}-\Lambda {\tilde{r}}_{kk}^{n+1})]$$

(59)

$$s_{ij}^{n+1}=s_{ij}^n+2{\overline{G}}_{ave}[\Delta e_{ij}-\Lambda {\tilde{d}}_{ij}^{n+1}]$$

(60)

$${\overline{G}}_{ave}=\frac{{\overline{G}}_0{p}^n[e^{{\overline{K}}_0(\Delta {\epsilon }_{kk}-\Lambda {\tilde{r}}_{kk}^{n+1})}-1]}{{\overline{K}}_0(\Delta {\epsilon }_{kk}-\Lambda {\tilde{r}}_{kk}^{n+1})}$$

(61)

${\tilde{d}}_{ij}^{n+1}$ is the deviatoric part of ${\tilde{r}}_{ij}^{n+1}$ defined as:

$${\tilde{d}}_{ij}^{n+1}={\tilde{r}}_{ij}^{n+1}-\frac{1}{3}{\tilde{r}}_{kk}^{n+1}{\delta }_{ij}$$

(62)

Then, the stress tensor and hardening parameter are obtained as:

$${\sigma }_{ij}^{n+1}=s_{ij}^{n+1}+p^{n+1}{\delta }_{ij}$$

(63)

$${\tilde{H}}^{n+1}={\tilde{H}}^n+\Delta \tilde{H}$$

(64)

The calculation formula $\Delta \tilde{H}$ refers to Equation (40).

For simplicity of notation, let us define

$$\overline{\rho }={\overline{K}}_0(\Delta {\epsilon }_{kk}-\Lambda {\tilde{r}}_{kk}^{n+1})$$

(65)

$$\overline{\theta }=\Delta e_{ij}-\Lambda {\tilde{d}}_{ij}^{n+1}$$

(66)

Then we have

$${\sigma }_{ij}^{n+1}=s_{ij}^n+2{\overline{G}}_{ave}\overline{\theta }+p^n{e}^{\overline{\rho }}{\delta }_{ij}$$

(67)

with

$${\overline{G}}_{ave}=\frac{{\overline{G}}_0{p}^n[e^{\overline{\rho }}-1]}{\overline{\rho }}$$

(68)

Then $f_I$ and ${\mathit{f}}_E$ can be obtained from Equations (41) and (42).

The above process provides the derivation of incremental constitutive integration in the transformed stress space. Of course, if the three-dimensional situation is not considered, the UH model can be directly applied to the two-dimensional model through the dimension-extending technique.

3.3.3. Jacobean Matrix

Under the premise of operational convergence, the definition of the Jacobean matrix only affects the convergence speed and does not change the accuracy of the calculation results. Therefore, in the UMAT (user subroutine to define a material’s mechanical behavior) for ABAQUS, the elastic–plastic stiffness matrix is used for finite analysis. The elastic–plastic stiffness matrix can be written as follows:

$$D_{ijkl}^{ep}=D_{ijkl}^e-\frac{D_{ijmn}^e\frac{\partial f}{\partial {\tilde{\sigma }}_{mn}}\frac{\partial f}{\partial {\sigma }_{st}}{D}_{stkl}^e}{Y}$$

(69)

$$Y=\frac{\partial f}{\partial {\tilde{\sigma }}_{ii}}+\frac{\partial f}{\partial {\sigma }_{ij}}{D}_{ijkl}^e\frac{\partial f}{\partial {\tilde{\sigma }}_{kl}}$$

(70)

$$\frac{\partial f}{\partial {\sigma }_{ij}}=\frac{\partial f}{\partial \tilde{p}}\frac{\partial \tilde{p}}{\partial {\sigma }_{ij}}+\frac{\partial f}{\partial \tilde{q}}\frac{\partial \tilde{q}}{\partial {\sigma }_{ij}}=\frac{1}{3}\frac{\partial f}{\partial \tilde{p}}{\delta }_{ij}+\frac{\partial f}{\partial \tilde{q}}\frac{\partial q_c}{\partial {\sigma }_{ij}}$$

(71)

The calculation of $\partial q_c/\partial {\sigma }_{ij}$ can be directly obtained from Equations (9) and (10), where the derivation of stress components corresponding to principal stress is needed, which can be obtained from the following:

$$\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}=(2-{\delta }_{ij})l_i^k{l}_j^k$$

(72)

where $i,j,k=1,2,3$, ${\sigma }_k$ is the k-th principal stress; $l^k$ is the unit principal direction corresponding to the k-th principal stress and the repeated indicators do not sum. ${\delta }_{ij}$ is the Kronecker’s delta. Alternatively, the implicit differentiation rule can be adopted. The SMP criterion under triaxial compression can be compiled into Equation (73):

$$\frac{I_2}{I_3}=9\frac{I_1+q_c}{I_1^2+I_1{q}_c-2q_c^2}$$

(73)

Total differential for Equation (73):

$$\frac{\partial q_c}{\partial {\sigma }_{ij}}=\frac{-2\frac{\partial I_2}{\partial {\sigma }_{ij}}{q}_c^2+\frac{\partial A}{\partial {\sigma }_{ij}}(q_c+I_1)+A\frac{\partial I_1}{\partial {\sigma }_{ij}}}{4I_2{q}_c-A}$$

(74)

where $A=I_1{I}_2-9I_3$

$$\frac{\partial A}{\partial {\sigma }_{ij}}=\frac{\partial (I_1{I}_2-9I_3)}{\partial {\sigma }_{ij}}=I_1\frac{\partial I_2}{\partial {\sigma }_{ij}}+\frac{\partial I_1}{\partial {\sigma }_{ij}}{I}_2-9\frac{\partial I_3}{\partial {\sigma }_{ij}}$$

(75)

$\partial I_i/\partial {\sigma }_{ij}$ can be directly solved with the invariant of stress expressed by the stress component, and the detailed process can refer to reference [31].

3.4. The Constitutive Integral of the UH Model under Different Paths

Based on the above analysis, we can achieve the calculation of constitutive integrals for 2D and 3D UH models under different paths. Taking the calculation of unit with equal p and b paths in a three-dimensional principal stress space as an example, the derivation process is detailed below.

In the three-dimensional principal stress space, the unknown quantities are:

$${\mathit{s}}^T\triangleq ({\sigma }_1,{\sigma }_2,{\sigma }_3,{\epsilon }_1,{\epsilon }_2,{\epsilon }_3;\tilde{H})$$

(76)

and plastic multiplier $\Lambda$, eight unknown quantities in total.

The Equation (47) only has five independent equations, and the solution is not closed, so three additional equations need to be added. ${\dot{\epsilon }}_1$ as an undiminished quantity of a process can always be specified, and then the unknown quantity to be solved is:

$${{\mathit{s}}^{\mathbf{\prime }}}^T\triangleq ({\sigma }_1,{\sigma }_2,{\sigma }_3,{\epsilon }_2,{\epsilon }_3;\tilde{H})$$

(77)

For the path conditions of equal p and b, the equation that needs to be supplemented is:

$$\left\{\begin{array}{l}({\dot{\sigma }}_2-{\dot{\sigma }}_3)-b({\dot{\sigma }}_1-{\dot{\sigma }}_3)=0\hfill \\ ({\dot{\sigma }}_1+{\dot{\sigma }}_2+{\dot{\sigma }}_3)=0\hfill \end{array}\right.$$

(78)

Then, the backward Euler integration is applied to Equation (78).

$$\left\{\begin{array}{l}(\Delta {\sigma }_2-\Delta {\sigma }_3)-b(\Delta {\sigma }_1-\Delta {\sigma }_3)=0\\ \Delta {\sigma }_1+\Delta {\sigma }_2+\Delta {\sigma }_3=0\end{array}\right.$$

(79)

${\epsilon }_2,{\epsilon }_3$ and $\tilde{H}$ in ${\mathit{s}}^{\mathbf{\prime }}$ are used as initial variables along with stress and solved simultaneously. A better convergence with dimensional processing is required, and then the initial vector is:

$${{\mathit{S}}^{\mathbf{\prime }}}^T\triangleq ({\sigma }_1,{\sigma }_2,{\sigma }_3,{\xi }_2,{\xi }_3;\tilde{\Omega })$$

(80)

with

$${\xi }_2=K_n\Delta {\epsilon }_2$$

(81)

$${\xi }_3=K_n\Delta {\epsilon }_3$$

(82)

$$\tilde{\Omega }=K_n\tilde{H}$$

(83)

$K_n$ is the bulk modulus at the initial moment of this time step.

Then $\mathrm{MiCP}({\mathit{f}}_E,f_I)$ is rewritten as:

$${\mathit{f}}_E(\Lambda ;{\mathit{S}}^{\mathbf{\prime }})=\left\{\begin{array}{l}\mathit{\sigma }+{\mathit{\sigma }}_p-{\mathit{\sigma }}_e\hfill \\ {\tilde{\Omega }}_{n+1}-{\tilde{\Omega }}_n-\Delta \tilde{\Omega }\hfill \\ (\Delta {\sigma }_2-\Delta {\sigma }_3)-b(\Delta {\sigma }_1-\Delta {\sigma }_3)\hfill \\ \Delta {\sigma }_1+\Delta {\sigma }_2+\Delta {\sigma }_3\hfill \end{array}\right.$$

(84)

$$f_I(\Lambda ,{\tilde{\mathit{S}}}^{\mathbf{\prime }})=-\tilde{F}(\tilde{\mathit{\sigma }})$$

(85)

with ${\tilde{\mathit{S}}}^{\mathbf{\prime }}$:

$${{\tilde{\mathit{S}}}^{\mathbf{\prime }}}^T\triangleq ({\tilde{\sigma }}_1,{\tilde{\sigma }}_2,{\tilde{\sigma }}_3,{\xi }_2,{\xi }_3;\tilde{\Omega })$$

(86)

Equations (84) and (85) are ${\mathit{f}}_E$ and $f_I$ in the MiCP.

For other different paths, corresponding equations can also be supplemented based on the path conditions, such as equal b with undrained path.

If it is an equal b with an undrained path, then the equations are:

$$\left\{\begin{array}{l}({\dot{\sigma }}_2-{\dot{\sigma }}_3)-b({\dot{\sigma }}_1-{\dot{\sigma }}_3)=0\hfill \\ ({\dot{\epsilon }}_1+{\dot{\epsilon }}_2+{\dot{\epsilon }}_3)=0\hfill \end{array}\right.$$

(87)

If it is a confining pressure path with equal b, the equations are:

$$\left\{\begin{array}{l}({\dot{\sigma }}_2-{\dot{\sigma }}_3)-b({\dot{\sigma }}_1-{\dot{\sigma }}_3)=0\hfill \\ {\dot{\sigma }}_3=0\hfill \end{array}\right.$$

(88)

The above provides several common control equations for different paths in three-dimensional principal stress space.

Alternatively, we can also embed the algorithm into the finite element calculations, and then the boundary value problems can be solved directly.

4. Numerical Examples

To verify the applicability of the incremental constitutive integral in the UH model and the reliability of the algorithm combined with the Finite Element Method (FEM), the author has carried out verification work in the following four parts respectively. The first part uses AEP to compare individual integration points under different paths with experimental results, and preliminarily verifies the applicability of incremental constitutive integration in the UH model; in the second part, the algorithm is written as UMAT and embedded into the commercial software ABAQUS 2018 to calculate the unit, which verifies the correctness of the UMAT; the third part, considering the undrained loading condition, calculates the ultimate bearing capacity of the pile; the fourth part calculates the ultimate bearing capacity of a circular foundation based on the UMAT with different OCRs.

4.1. Verification of Different Stress Paths

This section takes the AEP for numerical calculations and compares the numerical calculation results under different paths with the experimental results to preliminarily verify the correctness of the incremental constitutive integral in the UH model.

4.1.1. Undrained Path

Firstly, the undrained stress path of clay under different over-consolidation ratios was verified, and compared with the remolded Kaolin soil test conducted by Loudon in 1967. The parameters of clay are shown in

Table 1. Values of material parameters for Kaolin clay.

λ	κ	ν	e0	M
0.24	0.045	0.2	1.27	0.898

. Figure 4 shows the comparison between the numerical calculation results and the experiment ($p_e$ the equivalent stress), and it can be found that the calculation results are in good agreement with the experiment, proving the correctness of the dimension-extending technique for the UH model.

4.1.2. Drain Path

Using the numerical results of triaxial compression and tension of Fujinomori clay [32] under equal p and drained conditions, the feasibility is verified of the dimension-extending technique based on the transformed stress method in three-dimensional space. The parameters of clay are shown in

Table 2. Values of material parameters for Fujinomori clay.

λ	κ	ν	e0	M
0.1046	0.0231	0.3	0.915	1.36

. When Over-consolidation Ratio (OCR) = 1, 2, 4, p = 196 kPa; when OCR = 8, p = 98 kPa. The comparison results are shown in Figure 5 and Figure 6.

From Figure 5 and Figure 6, it can be seen that the peak strength and shear expansion of the soil increase with the increase of over-consolidation, and the predicted results are in good agreement with the experimental results.

4.1.3. Cyclical Loading

Figure 7, Figure 8 and Figure 9 show the comparison between the numerical calculation results and experimental results of Fujinomori clay under cyclic loading paths with p = 392 kPa, p = 196 kPa, and ${\sigma }_r=const$. The parameters of clay are shown in Table 2. From Figure 7, Figure 8 and Figure 9, it can be observed that, under the cyclic loading path, the strain increment gradually decreases, the hysteresis loop gradually decreases, and the stress–strain curve tends to be dense, which is consistent with the calculation trend in reference [33].

4.1.4. True Triaxial Test

Figure 10 shows the comparison between the numerical calculation results of normally consolidated clay under equal p (p = 196 kPa) and b path and the true triaxial test data. The material parameters are shown in Table 2.

From Figure 10, it can be seen that the combination of the dimension-extending technique and transformation stress method applies to the three-dimensional strength criterion considering intermediate principal stress.

4.1.5. Complex Loading Conditions

To further validate the feasibility of the dimension-extending technique for the UH model, complex loading condition is considered, which is the true triaxial test for Fujinomori clay. The loading path on the π-plane is $\overrightarrow{OABOCDOE}$, as shown in Figure 11. $\overrightarrow{OA}$ is the compression path with ${\sigma }_y$ being the major principal stress; $\overrightarrow{AB}$ is the extension path with ${\sigma }_y/{\sigma }_z={\sigma }_1/{\sigma }_3=3$; $\overrightarrow{BO}$ is the unloading path; $\overrightarrow{OC}$ is the compression path with ${\sigma }_x$ being the major principal stress; $\overrightarrow{CD}$ is the path with ${\sigma }_x/{\sigma }_y={\sigma }_1/{\sigma }_3=3$; $\overrightarrow{DO}$ is the unloading path; $\overrightarrow{OE}$ is the path with ${\sigma }_z$ being the major principal stress.

Figure 12 shows the stress–strain curve, from which it can be seen that the numerical calculation results are in good agreement with the experimental data.

4.2. Element Validation

To apply the algorithm to the finite element, the author compiles the UMAT and then performs unit calculation based on AEP and EEP, respectively. The process of the UMAT is shown in Figure 13.

The parameters of clay are shown in

Table 3. Parameters of the UH model and MCC model in element.

λ	κ	ν	e0	M
0.095504	0.008836	0.3	0.88	1.3636

. A 1 m × 1 m × 1 m cube of soil model is established. Three directions in the bottom and the normal in the sides are constrained in the model. The initial stress is 98 kPa.

Figure 14 shows the element validation result. It can be seen that AEP, EEP, and ABAQUS have a good consistency, which verifies the correctness of the incremental constitutive integral in the UH model.

4.3. Undrained Vertical Compressive Capacity of Pile

In this section, the UH model is used to calculate the vertical load pile in the foundation. The concrete solid pile is located in normally consolidated saturated clay and the water table is level with the foundation. The model dimension is shown in Figure 15, and the model parameters are shown in

Table 4. Parameters of the UH model.

λ	κ	ν	e0	M	γ’ (kN/m3)	k (m/s)
0.2	0.04	0.35	2.0	1.2	8	1 × 10−7

. The pile adopts a linear elastic model, the elastic modulus is E = 20 Gpa, the Poisson’s ratio is $v=0.2$, and the friction coefficient between pile and soil is 0.577. If the soil is normally consolidated clay, the UH model can degenerate into the MCC model.

Figure 16 shows the pile-end soil hole compression field of the AEP, EEP and ABAQUS, which represents the built-in MCC model. Figure 17a is the load-displacement curve of the pile, and Figure 17b shows the change of pore pressure along the radial direction at ten meters depth. When OCR = 1, the UH model degenerates into MCC. The calculation results of AEP and EEP are almost identical to the calculation results of ABAQUS using the MCC model, which proves the correctness of the algorithm.

4.4. Circular Foundation

In this section, a circular foundation load plate test will be taken as an example for three-dimensional modeling, considering the symmetry of the problem, so only 1/4 of the finite boundary needs to be taken for analysis. The diameter of the circular plate is 0.4 m and the center of the shape is located at the symmetry point. The dimension of the model is 1 m wide and 2 m high. The parameters of the clay are shown in

Table 5. Parameters of the UH model and MCC model in the circular foundation.

λ	κ	ν	e0	M
0.015	0.003	0.35	1.008	1.19

. Based on the self-weight of the soil and the overlying soil, construct an initial stress field:

$$\left\{\begin{array}{l}{\sigma }_y={\gamma }^{\prime }h+{\sigma }_{y0}\hfill \\ {\sigma }_x=K_y{\sigma }_y\hfill \end{array}\right.$$

(89)

${\sigma }_y$ and ${\sigma }_x$ are vertical stress and horizontal stress, respectively, where ${\gamma }^{\prime },h,{\sigma }_{y0}$ and $K_y$ are the buoyant density, depth, a constant initial surface pressure, and lateral earth pressure coefficient, respectively. $K_y=0.5$, ${\sigma }_{y0}=10 \mathrm{k}\mathrm{P}\mathrm{a}$, ${\gamma }^{\prime }=8 \mathrm{k}\mathrm{N}/{\mathrm{m}}^3$.

Figure 18a is the displacement cloud diagram, and Figure 18b is the load-displacement curve with different over-consolidation ratios. The results indicate that the ultimate bearing capacity of soil increases, with the OCR increasing. Figure 18b shows that, for soils with the same OCR, the calculated results using ABAQUS are the highest. This is due to the fact that the MCC model built into ABAQUS is based on the generalized Mises criterion extended to general stress space, which tends to overestimate the ultimate bearing capacity of the soil.

5. Conclusions

This paper mainly provides the inferential process of the incremental constitutive integration for the UH model combined with the transformation stress method. It gives a detailed derivation of two processing methods for the variable elastic properties of the UH model in the transformed stress space: AEP and EEP. Based on the complementarity theory and the dimension-extending technique, the constitutive integration problem of the UH model is reduced to a MiCP, and solved by the GSPC algorithm, which does not require calculation of the Hessian matrix, and greatly simplifies the calculation process.

In order to verify the correctness of the algorithm, numerical validation was conducted on a single integration point under different paths, including drainage path, undrained path, cyclic loading, and true triaxial test path. The calculation results indicate that the algorithm is correct and robust. Afterwards, in order to apply the algorithm to boundary value problems, a UMAT subroutine was written and compared with the calculation results of ABAQUS to verify the correctness of the subroutine. The numerical calculation results indicate that the dimension-extending technique and GSPC algorithm are applicable in the UH model, and AEP integration and EEP integration have almost the same accuracy in both processing methods. Therefore, for equations that cannot be solved with EEP integration, AEP can also meet the computational requirements. Eventually, the bearing capacity of a circular foundation with different OCRs was calculated, and the bearing capacity of soil increases, with the OCR increasing. The 3D MCC model based on generalized Mises criterion overestimates the soil strength. The research results have expanded the scope of use of the GSPC algorithm. It can provide a framework of GSPC for a similar hardening/softening model combined with the transformed stress method.