Governing equations

Hence, in order to calculate the stress response of the integration point, one has to know the elastic deformation gradient at the end of the increment $\tnsr F_\text e(t)$. In contrast to the total deformation gradient, this is, however, not known a priori. It depends on the partitioning of the deformation gradient into elastic and plastic part. \begin{equation} \label{eq:constitutive Fe} \tnsr F_\text e = \tnsr F {\tnsr F_\text p}^{-1} \end{equation} While the material initially deforms purely elastic, the proportion of the plastic deformation will increase once the yield point of the material is reached.

The rate with which the plastic deformation gradient $\tnsr F_\text p$ evolves is determined by the plastic velocity gradient $\tnsr L_\text p$ by means of a nonlinear ordinary differential equation. \begin{equation} \label{eq:constitutive Fp} \dot{\tnsr F}_\text p = \tnsr L_\text p \tnsr F_\text p \end{equation}

Equations \eqref{eq:constitutive S} to \eqref{eq:constitutive deltastate} govern the evolution of stress, strain, and microstructure and have to be solved at the level of each material point. Due to the nonlinear coupling of the unknown variables, this is not easily possible. In order to reduce the complexity of the problem, we decouple the evolution of the microstructure (i.e. state variables) from the evolution of the other variables. This can be reasoned by the fact that the microstructure usually evolves slower than stress and strain (akin to a Born–Oppenheimer approximation). Then, the set of equations can be solved in two levels:

• In an inner level, the stress $\tnsr S$ is obtained by solving equations \eqref{eq:constitutive S} to \eqref{eq:constitutive Lp} at constant state $\vctr\omega$.
• In an outer level, the state $\vctr\omega$ is obtained by solving equations \eqref{eq:constitutive dotstate} and \eqref{eq:constitutive deltastate} for a given stress $\tnsr S$.

Inner level of stress integration

	Figure 1: Modified Newton–Raphson scheme employed for the time integration of the stress.

Linearized Backward Euler

In order to solve the set of equations \eqref{eq:constitutive S} to \eqref{eq:constitutive Lp}, we approximate the time derivative in \eqref{eq:constitutive Fp} by a backward finite difference, and we obtain an implicit equation for the plastic deformation gradient at the end of the increment. \begin{equation*} \frac{\tnsr F_\text p(t) - \tnsr F_\text p(t_0)}{\Delta t} = \tnsr L_\text p(t) \tnsr F_\text p(t) \end{equation*} This equation can be solved for $\tnsr F_\text p(t)$ by simple linear algebra. \begin{equation*} \tnsr F_\text p(t) = {\left(\tnsr I - \Delta t\; \tnsr L_\text p(t)\right)}^{-1} \tnsr F_\text p(t_0) \end{equation*} Inserting this into equation \eqref{eq:constitutive Fe} one obtains an equation for the elastic deformation gradient at the end of the increment. \begin{equation*} \tnsr F_\text e(t) = \tnsr F(t) \; {\tnsr F_\text p}^{-1}(t_0) \left(\tnsr I - \Delta t\; \tnsr L_\text p(t)\right) \end{equation*} By this linearization procedure, one obtains a system of three algebraic equations with three unknowns $\tnsr S(t)$, $\tnsr F_\text e (t)$, and $\tnsr L_\text p (t)$: \begin{align} \tnsr S (t) &= \tnsr S(\tnsr F_\text e (t)) \nonumber \\ \tnsr F_\text e(t) &= \tnsr A \left(\tnsr I - \Delta t\; \tnsr L_\text p(t)\right) \label{eq:backwardSystem} \\ \tnsr L_\text p (t) &= \tnsr L_\text p(\tnsr S (t))\nonumber \end{align} where the known quantity $\tnsr F(t) \; {\tnsr F_\text p}^{-1}(t_0)$ is substituted by $\tnsr A$ for brevity and corresponds to a fully elastic predictor.

For later use, the derivative of the elastic deformation gradient with respect to the plastic velocity gradient follows as \begin{align} \label{eq:backwardFeDerivative} \tnsr F_\text e,_{\scriptscriptstyle\tnsr L_\text p} &= \left[\tnsr A \left(\tnsr I-\Delta t \; \tnsr L_\text p\right)\right],_{\scriptscriptstyle\tnsr L_\text p} = -\Delta t \, \tnsr A \tnsrfour I = -\Delta t \, \tnsr A \otimes\tnsr I \end{align} It is useful to rewrite this in index notation: \begin{align*} \frac{\partial {F_\text e}_{\,ij}}{\partial {L_\text p}_{\,kl}} \vctr g^i \otimes \vctr g_k \otimes \vctr g_l \otimes \vctr g^j &= -\Delta t \, A_{i\cdot}^{\phantom i k} \delta^l_j \,\vctr g^i \otimes \vctr g_k \otimes \vctr g_l \otimes \vctr g^j \end{align*}

Exponential Forward Euler

An alternative way to solve the set of equations \eqref{eq:constitutive S} to \eqref{eq:constitutive Lp} is by considering the correct solution \begin{equation*} \tnsr F_\text p(t) = \exp\left(\Delta t\,\tnsr L_\text p(t)\right) \tnsr F_\text p(t_0) \end{equation*}

Inserting this into equation \eqref{eq:constitutive Fe} one obtains an equation for the elastic deformation gradient at the end of the increment. \begin{align*} \tnsr F_\text e(t) &= \tnsr F(t) \; {\tnsr F_\text p}^{-1}(t_0) \exp\left(- \Delta t\; \tnsr L_\text p(t)\right) \\ &= \tnsr F(t) \; {\tnsr F_\text p}^{-1}(t_0) \left( \tnsr I -\Delta t \;\tnsr L_\text p(t) + \frac{\Delta t^2}{2!}{\tnsr L_\text p}^2 (t) + \mathcal O(\Delta t^3)\right) \end{align*} By truncating the series expansion after second order, one obtains a system of three algebraic equations with three unknowns $\tnsr S(t)$, $\tnsr F_\text e (t)$, and $\tnsr L_\text p (t)$: \begin{align} \tnsr S (t) &= \tnsr S(\tnsr F_\text e (t)) \nonumber \\ \tnsr F_\text e(t) &= \tnsr A \left( \tnsr I -\Delta t \;\tnsr L_\text p(t) + \frac{\Delta t^2}{2!}\tnsr L_\text p(t) \right) \label{eq:exponentialSystem} \\ \tnsr L_\text p (t) &= \tnsr L_\text p(\tnsr S (t))\nonumber \end{align} where the known quantity $\tnsr F(t) \; {\tnsr F_\text p}^{-1}(t_0)$ is substituted by $\tnsr A$ for brevity and corresponds to a fully elastic predictor.

For later use, the derivative of the elastic deformation gradient with respect to the plastic velocity gradient follows as \begin{align} \label{eq:exponentialFeDerivative} \tnsr F_\text e,_{\scriptscriptstyle\tnsr L_\text p} &= \left[\tnsr A \exp\left(- \Delta t\; \tnsr L_\text p \right)\right],_{\scriptscriptstyle\tnsr L_\text p}\nonumber \\ &= \tnsr A \left[ -\Delta t \tnsr I \otimes \tnsr I + \frac{\Delta t^2}{2!}\left( \tnsr L_\text p \otimes \tnsr I + \tnsr I \otimes \tnsr L_\text p \right) \right]\nonumber \\ &= -\Delta t \, \tnsr A \otimes\tnsr I + \frac{\Delta t^2}{2}\left( \tnsr A \tnsr L_\text p \otimes \tnsr I + \tnsr A \otimes \tnsr L_\text p\right)\;. \end{align} It is useful to rewrite this in index notation: \begin{align*} \frac{\partial {F_\text e}_{\,ij}}{\partial {L_\text p}_{\,kl}} \vctr g^i \otimes \vctr g_k \otimes \vctr g_l \otimes \vctr g^j &= -\Delta t \, A_{i\cdot}^{\phantom i k} \delta^l_j \\ &\phantom = +\frac{\Delta t^2}{2} \left( A_{i m} {L_\text p}^{m k} \delta^l_j + A_{i\cdot}^{\phantom i k} {L_\text p}^{l}_j \right) \vctr g^i \otimes \vctr g_k \otimes \vctr g_l \otimes \vctr g^j \end{align*}

Consistent solution

The systems \eqref{eq:backwardSystem} and \eqref{eq:exponentialSystem} can be solved by a Newton–Raphson scheme by either minimizing the residuum in $\tnsr S$, $\tnsr F_\text e$, or $\tnsr L_\text p$. Choosing the norm of the residuum in $\tnsr S$ as an objective function has the advantage that the inverse that is needed for the Newton–Raphson procedure is only $6\times6$ compared to $9\times9$ for $\tnsr L_\text p$ or $\tnsr F_\text e$, since the stress tensor is symmetric, while $\tnsr L_\text p$ and $\tnsr F_\text e$ are not. However, it is much harder to guess and correct for $\tnsr S$ than for $\tnsr L_\text p$, since the plastic velocity is usually very sensitive to a change in the stress. For this reason the Newton–Raphson scheme is chosen around $\tnsr L_\text p$ with the advantage of fast convergence at the expense of a higher cost per iteration step.

Outer level of state integration

Fixed-point iteration scheme

Explicit Euler integrator

Adaptive Euler integrator

Fourth-order explicit Runge-Kutta integrator

Fifth-order adaptive Runge-Kutta integrator