Computationally efficient stress reconstruction from full-field strain measurements

Abstract

Stress reconstruction based on experimentally acquired full-field strain measurements is computationally expensive when using conventional implicit stress integration algorithms. The computational burden associated with repetitive stress reconstruction is particularly relevant when inversely characterizing plastic material behaviour via inverse methods, like the nonlinear Virtual Fields Method (VFM). Spatial and temporal down-sampling of the available full-field strain data is often used to mitigate the computational effort. However, for metals subjected to non-linear strain paths, temporal down-sampling of the strain fields leads to erroneous stress states biasing the identification accuracy of the inverse method. Hence, a significant speedup factor of the stress integration algorithm is required to fully exploit the experimental data acquired by Digital Image Correlation (DIC). To this end, we propose an explicit stress integration algorithm that is independent on the number of images (i.e. strain fields) taken into account in the stress reconstruction. Theoretically, the proposed method eliminates the need for spatial and temporal down-sampling of the experimental full-field data used in the nonlinear VFM. Finally, the proposed algorithm is also beneficial in the emerging field of real-time DIC applications.

Appendix A: The stress update algorithm

The algorithm is based on the NICE scheme presented in [53] for general rate-independent plasticity models defined by a set of Differential–Algebraic Equations (DAE). Classical models consist of a yield condition and a set of evolution equations that determine the change of the state variables with loading given as

$$ {\Phi } = {\Phi }\left( {{\varvec{\sigma}},\sigma_{{\text{Y}}} ,{{\varvec{\upkappa}}}} \right) = 0 $$

$$ {\text{d}}{\varvec{\sigma}} = {\varvec{C}}_{{{\text{el}}}} \user2{ }\left( {{\text{d}}{\varvec{\varepsilon}}^{{{\text{tot}}}} - {\text{d}}{\varvec{\varepsilon}}^{{\text{p}}} } \right), $$

$$ {\text{d}}{\varvec{\varepsilon}}^{{\text{p}}} = \frac{{\partial {\Phi }}}{{\partial {\varvec{\sigma}}}}{\text{d}}\lambda , $$

$$ {\text{d}}\varepsilon_{{{\text{eq}}}}^{{\text{p}}} = \frac{{\varvec{\sigma}}}{{{\upsigma }_{{\text{Y}}} }}{\text{d}}{\varvec{\varepsilon}}^{{\text{p}}} , $$

$$ {\text{d}}\sigma_{Y} = H {\text{d}}\varepsilon_{{{\text{eq}}}}^{p} , $$

$$ {\text{d}}{{\varvec{\upkappa}}} = {\text{d}}{{\varvec{\upkappa}}}\left( {{\varvec{\sigma}},\sigma_{{\text{Y}}} ,{{\varvec{\upkappa}}}} \right), $$

where $\Phi$ is the yield function, ${\varvec{\sigma}}$ is the stress tensor, ${\sigma }_{{\text{Y}}}$ is the yield stress, ${{\varvec{C}}}_{{\text{el}}}$ is the elastic stiffness tensor, ${\varvec{\varepsilon}}$ is the total strain tensor, ${{\varvec{\varepsilon}}}^{{\text{p}}}$ is the plastic strain tensor, ${\varepsilon }_{{\text{eq}}}^{{\text{p}}}$ is equivalent plastic strain, $\lambda$ is the plastic multiplier, $H=\frac{{\partial \sigma }_{{\text{Y}}}}{\partial {\varepsilon }_{{\text{eq}}}^{{\text{p}}}}$ is the strain hardening modulus, and ${\varvec{\upkappa}}$ is a vector of the other state variables. If such a set of equation is transformed in incremental form in view of the NICE formulation, the above equations yield the following explicit set of equation:

$$\Phi + \frac{{\partial \Phi }}{{\partial {\sigma }}}\Delta {\sigma } + \frac{{\partial \Phi }}{{\partial \sigma _{{\text{Y}}} }}\Delta \sigma _{{\text{Y}}} + \frac{{\partial \Phi }}{{\partial {\mathbf{\kappa }}}}\Delta {\mathbf{\kappa }} = 0,$$

$$\Delta \sigma ={\varvec{C}}_{ {\text{el} }} \left( \Delta {\varepsilon }^{ {\text{tot}}} - \frac{ \partial \Phi } { \partial \sigma } \Delta \lambda \right),$$

$$\Delta \sigma_{{\text{Y}}} = H\,\frac{ \sigma } {\sigma _{{\text{Y}}} } \frac{{\partial\Phi }}{{\partial{\sigma }}}\Delta \lambda ,$$

$$\Delta {\mathbf{\kappa }} = \frac{{\partial {\mathbf{\kappa }}}}{{\partial \lambda }}\Delta \lambda .$$

To preserve the generality of the derivation, the following vector/matrix notation of the tensor variables is defined

$$\begin{aligned}&\Delta {\varvec{\upvarepsilon}}={\left\{\begin{array}{cccccc}\Delta {\varepsilon }_{11}&\Delta {\varepsilon }_{22}&\Delta {\varepsilon }_{33}&\Delta {\varepsilon }_{12}&\Delta {\varepsilon }_{23}&\Delta {\varepsilon }_{31}\end{array}\right\}}^{\mathbf{T}},\\ &{\varvec{\Sigma}}=\left\{\begin{array}{c}{\varvec{\upsigma}}\\ {\sigma }_{{\text{Y}}}\\ {\varvec{\upkappa}}\end{array}\right\}={\left\{\begin{array}{cccccccccc}{\sigma }_{11}& {\sigma }_{22}& {\sigma }_{33}& {\sigma }_{12}& {\sigma }_{23}& {\sigma }_{31}& {\sigma }_{{\text{Y}}}& {\kappa }_{1}& \dots & {\kappa }_{p}\end{array}\right\}}^{\mathbf{T}},\end{aligned}$$

$$\begin{aligned}&\mathbf{C}=\left[\begin{array}{c}{\left[{\mathbf{C}}_{{\text{el}}}\right]}_{6{\text{x}}6}\\ {\left[0\right]}_{1{\text{x}}6}\\ {\left[0\right]}_{p{\text{x}}6}\end{array}\right]=\left[\begin{array}{cccccc}{C}_{1111}& {C}_{1122}& {C}_{1133}& {C}_{1112}& {C}_{1123}& {C}_{1131}\\ {C}_{2211}& {C}_{2222}& {C}_{2233}& {C}_{2212}& {C}_{2223}& {C}_{2231}\\ {C}_{3311}& {C}_{3322}& {C}_{3333}& {C}_{3312}& {C}_{3323}& {C}_{3331}\\ {C}_{1211}& {C}_{1222}& {C}_{1233}& {C}_{1212}& {C}_{1223}& {C}_{1231}\\ {C}_{2311}& {C}_{2322}& {C}_{2333}& {C}_{2312}& {C}_{2323}& {C}_{2331}\\ {C}_{3111}& {C}_{3122}& {C}_{3133}& {C}_{3112}& {C}_{3123}& {C}_{3131}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0& 0& 0& 0& 0& 0\end{array}\right];\\ \boldsymbol{ }\boldsymbol{ }\boldsymbol{ }&\mathbf{R}=\left\{\begin{array}{c}{\left[{\mathbf{C}}_{{\text{el}}}\bullet \frac{\partial\Phi }{\partial {\varvec{\upsigma}}}\right]}_{6{\text{x}}1}\\ {\left[-H\frac{{\varvec{\upsigma}}}{{\sigma }_{{\text{Y}}}}\bullet \frac{\partial\Phi }{\partial {\varvec{\upsigma}}}\right]}_{1{\text{x}}1}\\ {\left[-\frac{\partial {\varvec{\upkappa}}}{\partial \lambda }\right]}_{p{\text{x}}1}\end{array}\right\}\end{aligned}$$

(A1)

where $\Delta {\varvec{\upvarepsilon}}$ is a vector of strain increments, a vector ${\varvec{\Sigma}}$ contains stress components, the yield stress, and all other state variables, such as the components of the back-stress tensor in kinematic hardening or the porosity in the Gurson model. The vector $\mathbf{R}$ contains all the terms that are multiplied by $\Delta \lambda$. For multiplications between the stacked matrix and the vector quantities, the rules of the special matrix operator apply here (see [53] for more on notation).

Let us briefly recall the basic idea of the NICE scheme and its formulation from [53]. A consistency (yield) condition $\Phi =0$ is the algebraic equation expanded in Taylor series to approximate $\Phi$ at the end of the increment, while all differential equations are converted to difference form, as is known from explicit methods. This effectively suppresses the permanent drift that typically occurs in explicit methods. The model equations can be written in a generalized stacked matrix form as

$$ \Phi + \frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot \Delta {{\varvec{\Sigma}}} = 0,{ }\;\Delta {{\varvec{\Sigma}}} = {\mathbf{C}}\Delta {{\varvec{\upvarepsilon}}} - {\mathbf{R}} \Delta \lambda $$

(A2)

from where the plastic multiplier $\mathrm{\Delta \lambda }$ can be explicitly extracted

$$ \Delta \lambda = \frac{{\Phi + \frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{C}} \cdot \Delta {{\varvec{\upvarepsilon}}}}}{{\frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{R}}}} $$

(A3)

Considering also the Kuhn-Tucker conditions for the elastic–plastic transition, part of the increment $\beta\Delta {\varvec{\upvarepsilon}}$ requires an elastic update, and $(1-\beta )\Delta {\varvec{\upvarepsilon}}$ remains for the plastic update described above. The coefficient $\beta$ is estimated from actual and trial values of the yield function $\Phi$ and ${\Phi }^{{\text{tr}}}$, respectively

$$ \beta = \frac{{0 - {\Phi }}}{{{\Phi }^{{{\text{tr}}}} - {\Phi }}}. $$

(A4)

The scheme NICE is a conditionally stable scheme, and in [53] it was shown that in each subincrement the size of the next stable subincrement can be computed, so that the integration can continue with the maximum possible stable increment. It was also shown that an Algorithmic Tangent Stiffness (ATS) can be updated in each subincrement, but in a stress reconstruction from a measured strain field, this term is not needed. Briefly, the remainder of the strain increment $(1-\beta )\Delta {\varvec{\upvarepsilon}}$ is further subdivided into smaller subincrements, where the size of the current increment is $\left(1-\beta \right) \alpha\Delta {\varvec{\upvarepsilon}}$ with the ratio $\alpha$ being

$$ \alpha = \frac{{\frac{2}{{{\text{max}}\left( { - {{\varvec{\upomega}}}} \right)}}\frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{R}} - \Phi }}{{\left( {1 - \beta } \right)\frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{C}} \cdot \Delta {{\varvec{\upvarepsilon}}}}}; {{\varvec{\upomega}}} = {\text{eig}}\left( { - \frac{{\partial {\mathbf{R}}}}{{\partial {{\varvec{\Sigma}}}}}} \right). $$

(A5)

Thus, the elastic updating is $\Delta {{\varvec{\Sigma}}} = {\mathbf{C}} \cdot \Delta {{\varvec{\upvarepsilon}}}$, while the plastic updating of each subincrement follows the equation

$$ \begin{aligned}&\Delta {{\varvec{\Sigma}}} = {\mathbf{C}} \cdot \left( {1 - \beta } \right) \alpha \Delta {{\varvec{\upvarepsilon}}} - {\mathbf{R}} \Delta \lambda ;\\ &{ }\Delta \lambda = \frac{{\Phi + \frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{C}} \cdot \left( {1 - \beta } \right) \alpha \;\Delta {{\varvec{\upvarepsilon}}}}}{{\frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}}} \cdot {\mathbf{R}}}}.\end{aligned} $$

(A6)

This was a brief summary of the method, which serves as a theoretical background; a precise notation of the formulation, its robustness, accuracy, and efficiency is described in detail in [53]. Let us now extend the explanation to make it suitable for stress reconstruction from measured strain fields. The DIC measurements capture the surface strains, so the equations need to be reformulated for the plane stress state ${\upsigma }_{33}=0$. For this purpose, the variables in Eq. (A1) are redefined in the following reduced form

$$ \begin{gathered} \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} = \left\{ {\begin{array}{*{20}c} {\Delta \varepsilon_{11} } & {\Delta \varepsilon_{22} } & {\Delta \varepsilon_{12} } \\ \end{array} } \right\}^{{\mathbf{T}}} ,{{\varvec{\Sigma}}}_{{{\text{PS}}}} = \left\{ {\begin{array}{*{20}c} {\sigma_{11} } & {\sigma_{22} } & {\sigma_{12} } & {\sigma_{{\text{Y}}} } \\ \end{array} } \right\}^{{\mathbf{T}}} , \hfill \\ {\mathbf{C}}_{{{\text{red}}}} = \left[ {\begin{array}{*{20}c} {\left[ {{\mathbf{C}}_{{{\text{el}}}} } \right]_{{3{\text{x}}3}} } \\ {\left[ 0 \right]_{{1{\text{x}}3}} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {C_{1111} } & {C_{1122} } & {C_{1112} } \\ {C_{2211} } & {C_{2222} } & {C_{2212} } \\ {C_{1211} } & {C_{1222} } & {C_{1212} } \\ 0 & 0 & 0 \\ \end{array} } \right];\user2{ }\hfill \\ {\mathbf{R}}_{{{\text{red}}}} = \left\{ {\begin{array}{*{20}c} {\left[ {{\mathbf{C}}_{{{\text{el}}}} \cdot \frac{\partial \Phi }{{\partial {{\varvec{\upsigma}}}}}} \right]_{{3{\text{x}}1}} } \\ {\left[ { - H\frac{{{\varvec{\upsigma}}}}{{\sigma_{{\text{Y}}} }} \cdot \frac{\partial \Phi }{{\partial {{\varvec{\upsigma}}}}}} \right]_{{1{\text{x}}1}} } \\ \end{array} } \right\} \hfill \\ \end{gathered} $$

(A7)

$$ {\mathbf{C}}_{{3{\text{Col}}}} = \left\{ {\begin{array}{*{20}c} {C_{1133} } & {C_{2233} } & {C_{1233} } & 0 \\ \end{array} } \right\}^{{\mathbf{T}}} , $$

$$ {\mathbf{C}}_{{3{\text{Row}}}} = \left\{ {\begin{array}{*{20}c} {C_{3311} } & {C_{3322} } & {C_{3312} } \\ \end{array} } \right\}, $$

where, for simplicity, we assume that the model does not contain a variable ${\varvec{\upkappa}}$, and subscrips “PS” and “red” stand for “plane stress” and “reduced”, respectively. However, we define here two additional vectors ${\mathbf{C}}_{3{\text{Col}}}$ and ${\mathbf{C}}_{3{\text{Row}}}$, which represent the omitted third column and third row of the original matrix C, respectively. Equation (A2) can now be rewritten as follows

$$ \Phi + \frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}_{{{\text{PS}}}} }} \cdot \Delta {{\varvec{\Sigma}}}_{{{\text{PS}}}} = 0,\Delta {{\varvec{\Sigma}}}_{{{\text{PS}}}} = {\mathbf{C}}_{{{\text{red}}}} \cdot \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} + {\mathbf{C}}_{{3{\text{Col}}}} { }\Delta \varepsilon_{33} - {\mathbf{R}}_{{{\text{red}}}} \Delta \lambda $$

(A8)

to which the zero-normal-stress condition is added

$$ \Delta \sigma_{33} = {\mathbf{C}}_{{3{\text{Row}}}} \cdot \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} + C_{3333} \Delta \varepsilon_{33} - R_{3} \Delta \lambda = 0, $$

(A9)

where ${R}_{3}$ is the third term of the original vector $\mathbf{R}$. The term $\Delta {\varepsilon }_{33}$ can now be extracted from Eq. (A9) and, after substitution into Eq. (A8), yields

$$ \Delta {{\varvec{\Sigma}}}_{{{\text{PS}}}} = \left( {{\mathbf{C}}_{{{\text{red}}}} - \frac{{{\mathbf{C}}_{{3{\text{Col}}}} { } \otimes {\mathbf{C}}_{{3{\text{Row}}}} { }}}{{C_{3333} }}} \right) \cdot \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} - \left( {{\mathbf{R}}_{{{\text{red}}}} - \frac{{{\mathbf{C}}_{{3{\text{Col}}}} { }R_{3} { }}}{{C_{3333} }}} \right) \Delta \lambda , $$

(A10)

or in a condensed form

$$ \Delta {{\varvec{\Sigma}}}_{{{\text{PS}}}} = {\mathbf{C}}_{{{\text{PS}}}} \cdot \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} - {\mathbf{R}}_{{{\text{PS}}}} \Delta \lambda $$

(A11)

The form of the equation is mathematically identical to Eq. (A2), so the same derivation used for Eq. (A6) leads to the plastic update for the plane stress state

$$ \begin{aligned} &\Delta {{\varvec{\Sigma}}}_{{{\text{PS}}}} = {\mathbf{C}}_{{{\text{PS}}}} \cdot \left( {1 - \beta } \right) \alpha \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} - {\mathbf{R}}_{{{\text{PS}}}} \Delta \lambda ;{ }\\ &\Delta \lambda = \frac{{\Phi + \frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}_{{{\text{PS}}}} }} \cdot {\mathbf{C}}_{{{\text{PS}}}} \cdot \left( {1 - \beta } \right) \alpha \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} }}{{\frac{\partial \Phi }{{\partial {{\varvec{\Sigma}}}_{{{\text{PS}}}} }} \cdot {\mathbf{R}}_{{{\text{PS}}}} }}.\end{aligned} $$

(A12)

Without compromising the generality of the approach, we will now present a simplification of the equations for a particular model. For clarity of presentation of the principle, we have chosen the von Mises plasticity model with hardening combined with isotropic linear elasticity. In the calculation, the following definitions of the non-zero terms are

$$ \begin{aligned} &\Phi = \sigma_{eq} - \sigma_{Y}; \sigma_{eq} = \sqrt {\frac{3}{2}{\mathbf{s}} \cdot {\mathbf{s}}},\\ &C_{1111} = C_{2222} = C_{3333} = \frac{{E\left( {1 - \nu } \right)}}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}, \\ &C_{1122} = C_{1133} = C_{2233} = C_{2211} = C_{3311} = C_{3322} = \frac{E \nu }{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}},\\ &C_{1212} = \mu = \frac{E}{{2\left( {1 + \nu } \right)}}, \end{aligned} $$

(A13)

where $\mathbf{s}$ is a deviatoric stress. The term ${\mathbf{C}}_{{\text{PS}}}$ in Eq. (A11) simplifies to the extended well-known stiffness matrix for plane stress state

$$ {\mathbf{C}}_{{{\text{PS}}}} = \user2{ }{\mathbf{C}}_{{{\text{red}}}} - \frac{{{\mathbf{C}}_{{3{\text{Col}}}} { } \otimes {\mathbf{C}}_{{3{\text{Row}}}} { }}}{{C_{3333} }} = \left[ {\begin{array}{*{20}c} {C_{{{\text{dia}}}} } & {C_{{{\text{offdia}}}} } & 0 \\ {C_{{{\text{offdia}}}} } & {C_{{{\text{dia}}}} } & 0 \\ 0 & 0 & \mu \\ 0 & 0 & 0 \\ \end{array} } \right] = \frac{E}{{1 - \nu^{2} }}\left[ {\begin{array}{*{20}c} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & {\frac{1 - \nu }{2}} \\ 0 & 0 & 0 \\ \end{array} } \right] $$

(A13a)

and the term ${\mathbf{R}}_{{{\text{PS}}}}$ is simplified into

$$ {\mathbf{R}}_{{{\text{PS}}}} = {\mathbf{R}}_{{{\text{red}}}} - \frac{{{\mathbf{C}}_{{3{\text{Col}}}} { }R_{3} { }}}{{C_{3333} }} = \left\{ {\begin{array}{*{20}c} {\frac{3\mu }{{{\upsigma }_{{{\text{eq}}}} }}\left( {s_{11} - s_{33} \frac{\nu }{1 - \nu }} \right)} \\ {\frac{3\mu }{{{\upsigma }_{{{\text{eq}}}} }}\left( {s_{22} - s_{33} \frac{\nu }{1 - \nu }} \right)} \\ {\frac{3\mu }{{{\upsigma }_{{{\text{eq}}}} }}s_{12} } \\ { - H} \\ \end{array} } \right\}. $$

(A13b)

Using Eqs. (A13a) and (A13b), we can further simplify some expressions in Eq. (A12) and define two scalar variables ${V}_{1}$ and ${V}_{2}$

$$ \begin{gathered} V_{1} = \frac{{\partial {\Phi }}}{{\partial {{\varvec{\Sigma}}}_{{{\text{PS}}}} }} \cdot {\mathbf{C}}_{{{\text{PS}}}} \cdot \hfill \\ \Delta {{\varvec{\upvarepsilon}}}_{{{\text{PS}}}} = \frac{3}{{2 {\upsigma }_{eq} }}\left( s_{11} \left( {C_{{{\text{dia}}}} \Delta \varepsilon_{11} \hfill + C_{{{\text{offdia}}}} \Delta \varepsilon_{22} } \right)\right.\\ \left.+ s_{22} \left( {C_{{{\text{offdia}}}} \Delta \varepsilon_{11} + C_{{{\text{dia}}}} \Delta \varepsilon_{22} } \right) + 2s_{12} \mu \Delta \varepsilon_{12} \right), \hfill \\ V_{2} = \frac{{\partial {\Phi }}}{{\partial {{\varvec{\Sigma}}}_{{{\text{PS}}}} }} \cdot {\mathbf{R}}_{{{\text{PS}}}} = 3\mu + H - \frac{{9 \mu s_{33}^{2} \left( {1 - 2\nu } \right)}}{{2 {\upsigma }_{eq}^{2} \left( {1 - \nu } \right)}}. \hfill \\ \end{gathered}$$

(A14)

Consequently, the expression for $\mathrm{\Delta \lambda }$ in Eq. (A12) can be simplified to

$$ \Delta \lambda = \frac{{\Phi + \left( {1 - \beta } \right) \alpha V_{1} }}{{V_{2} }}. $$

(A15)

While the elastic–plastic transition coefficient $\beta$ is defined by Eq. (A4), the expression (A5) can be simplified to

$$\alpha =\frac{\frac{{\sigma }_{{\text{Y}}}}{27\mu }{ V}_{2}-\Phi }{(1-\beta ) {V}_{1}}.$$

(A16)

To speed up the calculation, the stable increment size is calculated once per increment (and not in each subincrement), which means that the plastic integration is performed in each increment in ${N}_{{\text{sub}}}={\text{floor}}\left(1/\alpha +1\right)$ equidistant subincrements.

$${N}_{{\text{sub}}}={\text{floor}}\left(\frac{(1-\beta ) {V}_{1}}{\frac{{\sigma }_{{\text{Y}}}}{27\mu }{ V}_{2}-\Phi }+1\right)$$

(A17)

Thus, the final expression for the stress update is based on Eq. (A12)

$$ \Delta \sigma_{11} = \left( {C_{{{\text{dia}}}} \Delta \varepsilon_{11} + C_{{{\text{offdia}}}} \Delta \varepsilon_{22} } \right)\frac{{\left( {1 - \beta } \right)}}{{N_{{{\text{sub}}}} }} - \frac{3\mu }{{{\upsigma }_{eq} }}\left( {s_{11} - s_{33} \frac{\nu }{1 - \nu }} \right)\Delta \lambda ;{ } \Delta \lambda = \frac{{\Phi + \frac{{\left( {1 - \beta } \right)}}{{N_{{{\text{sub}}}} }} V_{1} }}{{V_{2} }} $$

(A18)

$$\begin{gathered} {\Delta }\sigma_{22} = \left({C_{{{\text{offdia}}}} {\Delta }\varepsilon_{11} +C_{{{\text{dia}}}} \Delta \varepsilon_{22} } \right)\hfill\\ \frac{{\left( {1 - \beta } \right)}}{{N_{{{\text{sub}}}} }} -\frac{3\mu }{{{\upsigma }_{eq} }}\left( {s_{22} - s_{33} \frac{\nu}{1 - \nu }} \right)\Delta \lambda ,{ } \hfill \\ {\Delta}\sigma_{12} = \mu { }\Delta \varepsilon_{12} \frac{{\left( {1 -\beta } \right)}}{{N_{{{\text{sub}}}} }} - \frac{3\mu }{{{\upsigma}_{eq} }}s_{12} { }\Delta \lambda ,{ }\Delta \sigma_{{\text{Y}}} =H\Delta \lambda \hfill \\ \end{gathered}$$

The update algorithm is computationally inexpensive and very easy to use. The current increment is divided into an elastic subincrement and a plastic remainder. After the elastic subincrement is performed, all variables are updated to shift the state to the yield surface. Then, the plastic remainder is divided into stable, equidistant plastic subincrements, with the stress subincrements explicitly derived so that the state can be updated after each subincrement.