Closest-point return algorithm

Let us assume, that at time $t_n$ the total and plastic strain vectors and internal variables are known

$$\{$$ $\varepsilon$$$_n,$$    $\varepsilon$$$^p_n,$$   $\kappa$$$_n\}\ {\rm given\ at}\ t_n$$

By applying an implicit backward Euler difference scheme to the evolution equations (4.2 and 4.4) and making use of the initial conditions the following discrete non-linear system is obtained

$$\mbox{\boldmath$\varepsilon $} _{n+1} = \mbox{\boldmath$\varepsilon $} _n+\Delta \mbox{\boldmath$\varepsilon $}$$ (4.7)

$$\mbox{\boldmath$\sigma$} _{n+1} = \mbox{\boldmath$D$} ( \mbox{\boldmath$\varepsilon $} _{n+1}- \mbox{\boldmath$\varepsilon $} ^p_{n+1})$$ (4.8)

$$\mbox{\boldmath$\varepsilon $} ^p_{n+1} = \mbox{\boldmath$\varepsilon $} ^p_{n}+\sum\lambda^i\partial_{\mbox{\boldmath$\sigma$}} g_i(\mbox{\boldmath$\sigma$}_{n+1}, \mbox{\boldmath$\kappa$}_{n+1})$$ (4.9)

In addition, the discrete counterpart of the Kuhn-Tucker conditions becomes

$$f_i( \mbox{\boldmath$\sigma$} _{n+1}, \mbox{\boldmath$\kappa$} _{n+1}) =$$ 0 (4.10)

$$\lambda^i_{n+1} \ge$$ 0 (4.11)

$$\lambda^i_{n+1} f_i( \mbox{\boldmath$\sigma$} _{n+1}, \mbox{\boldmath$\kappa$} _{n+1}) =$$ 0 (4.12)

In the standard displacement-based finite element analysis, the strain evolution is determined by the displacement increments computed on the structural level. The basic task on the level of a material point is to evaluate the stress evolution generated by strain history. According to this, the strain driven algorithm is assumed, i.e. that the total strain $\varepsilon$$_{n+1}$ is given. Then, the Kuhn-Tucker conditions determine whether a constraint is active. The set of active constraints is denoted as $J_{act}$ and is defined as

$$J_{act}=\{\beta\in\{1,\cdots,m\}\vert f_\beta=0\ \&\ \dot{f}_\beta=0\}$$ (4.13)

Let's start with the definition of the residual of plastic flow

$$\mbox{\boldmath$R$} _{n+1}=- \mbox{\boldmath$\varepsilon $} ^p_{n+1}+ \mbox{\boldmath$\varepsilon $} ^p_{n}+\sum_{j\in J_{act}}\lambda^j_{n+1}\partial_\sigma g_{n+1}$$ (4.14)

By noting that total strain $\varepsilon$$_{n+1}$ is fixed during the increment we can express the plastic strain increment using (4.2) as

$$\Delta \mbox{\boldmath$\varepsilon $} ^p_{n+1} = - \mbox{\boldmath$D$} \Delta \mbox{\boldmath$\sigma$} _{n+1}$$ (4.15)

The linearization of the plastic flow residual (4.14) yields^4.1

$$\mbox{\boldmath$R$} + \mbox{\boldmath$D$} ^{-1}\Delta \mbox{\boldmath$\sigma$} +\sum\lambda\partial_{\sigma\sigma}g\Delta \mbox{\boldmath$\sigma$} +$$

$$+\sum\lambda\partial_{\sigma\kappa}g\cdot(\partial_\sigma\kappa\Delta \mbox{\boldmath$\sigma$} +\partial_\lambda\kappa\Delta\lambda)+\sum\Delta\lambda\partial_{\sigma}g =0$$ (4.16)

From the previous equation, the stress increment $\Delta$$\sigma$ can be expressed as

$$\Delta \mbox{\boldmath$\sigma$} =- \mbox{\boldmath$H$} ^{-1}\left(\mbox{\boldmath$R$}+\sum\Delta\lambda\partial_{\sigma}g+\sum\lambda\partial_{\sigma\kappa}g\partial_{\lambda}\kappa\Delta\lambda\right)$$ (4.17)

where $H$ is algorithmic moduli defined as

$$\mbox{\boldmath$H$} =\left[\mbox{\boldmath$D$}^{-1}+\sum\lambda\partial_{\sigma\sigma}g+\sum\lambda\partial_{\sigma\kappa}g\partial_{\sigma}\kappa\right]$$ (4.18)

Differentiation of active discrete consistency conditions (4.10) yields

$$\mbox{\boldmath$f$} +\partial_{\sigma} \mbox{\boldmath$f$} \Delta \mbox{\boldmath$\sigma$} +\partial_{\kappa} \mbox{\boldmath$f$} (\partial_\sigma\kappa\Delta \mbox{\boldmath$\sigma$} +\partial_\lambda\kappa\Delta\lambda) =0$$ (4.19)

Finally, by combining equations (4.17) and (4.19), one can obtain expression for incremental vector of consistency parameters $\Delta$$\lambda$

$$\left[\mbox{\boldmath$V$}^T\mbox{\boldmath$H$}^{-1}\mbox{\boldmat... ...{\boldmath$f$}-\mbox{\boldmath$V$}^T\mbox{\boldmath$H$}^{-1}\mbox{\boldmath$R$}$$ (4.20)

where the matrices $U$ and $V$ are defined as

$$\mbox{\boldmath$U$} = \left[\partial_{\sigma}\mbox{\boldmath$g$}+\sum\lambda\partial_{\sigma\kappa}g\partial_{\lambda}\kappa\right]$$ (4.21)

$$\mbox{\boldmath$V$} = \left[\partial_{\sigma}\mbox{\boldmath$f$}+\partial_{\kappa}\mbox{\boldmath$f$}\partial_{\sigma}\kappa\right]$$ (4.22)

Before presenting the final return mapping algorithm, the algorithm for determination of the active constrains should be discussed. A yield surface $f_{i,n+1}$ is active if $\lambda^i_{n+1} > 0$. A systematic enforcement of the discrete Kuhn-Tucker condition (4.10), which relies on the solution of return mapping algorithm, then serves as the basis for determining the active constraints. The starting point in enforcing (4.10) is to define the trial set

$$J^{trial}_{act}=\{j\in\{1,\cdots,m\}\vert f^{trial}_{j,n+1} > 0\}$$ (4.23)

where $J_{act}\subseteq J_{act}^{trial}$. Two different procedures can be adopted to determine the final set $J_{act}$. The conceptual procedure is as follows

• Solve the closest point projection with $J_{act}=J_{act}^{trial}$ to obtain final stresses, along with $\lambda^i_{n+1},\ i\in J_{act}^{trial}$.
• Check the sign of $\lambda^i_{n+1}$. If $\lambda^i_{n+1} <0$, for some $i\in J_{act}^{trial}$, drop the $i-$th constrain from the active set and goto first point. Otherwise exit.

In the procedure 2, the working set $J_{act}^{trial}$ is allowed to change within the iteration process, as follows

• Let $J_{act}^{(k)}$ be the working set at the k-th iteration. Compute increments $\Delta\lambda^{i,(k)}_{n+1},\ i\in J_{act}^{(k)}$.
• Update and check the signs of $\Delta\lambda^{i,(k)}_{n+1}$. If $\Delta\lambda^{i,(k)}_{n+1} < 0$, drop the i-th constrain from the active set $J_{act}^{(k)}$ and restart the iteration. Otherwise continue with next iteration.

If the consistency parameters $\Delta\lambda^{i}$ can be shown to increase monotonically within the return mapping algorithm, the the latter procedure is preferred since it leads to more efficient computer implementation.

The overall algorithm is convergent, first order accurate and unconditionally stable. The general algorithm is summarized in Tab. 4.2.2.

Table 4.1: General multisurface closest point algorithm