NICE^h: a higher-order explicit numerical scheme for integration of constitutive models in plasticity

Abstract

The article introduces, as a result of further development of the first-order scheme NICE, a simple and efficient higher-order explicit numerical scheme for the integration of a system of ordinary differential equations which is constrained by an algebraic condition (DAE). The scheme is based on the truncated Taylor expansion of the constraint equation with order h of the scheme being determined by the highest exponent in the truncated Taylor series. The integration scheme thus conceived will be named NICE^h, considering both principal premises of its construction. In conjunction with a direct solution technique used to solve the boundary value problem, the NICE^h scheme is very convenient for integrating constitutive models in plasticity. The plasticity models are defined mostly by a system of algebraic and differential equations in which the yield criterion represents the constraint condition. To study the properties of the new integration scheme, which, like the forward-Euler scheme, is characterised by its implementation simplicity due to the explicitness of its formulations, a damage constitutive model (Gurson–Tvergaard–Needleman model) is considered. The general opinion that the implicit backward-Euler scheme is much more accurate than the thus-far known explicit schemes is challenged by the introduction of the NICE^h scheme. The accuracy of the higher-order explicit scheme in the studied cases is significantly higher than the accuracy of the classical backward-Euler scheme, if we compare them under the condition of a similar CPU time consumption.

Appendices

Appendix A: the GTN constitutive model

Gurson [10] first introduced a yield criterion for a perfectly plastic hollow sphere which does not only depend on the von Mises equivalent stress but also on the hydrostatic stress σ_H and void volume fraction f(or porosity). A model for the nucleation of voids, according to which nucleation depends on the equivalent plastic strain of matrix material, was proposed in [15]. In order to find a better agreement with the finite element analyses of void growth, the Gurson model was modified by [16] into the well-known Gurson–Tvergaard–Needleman (GTN) model, which will serve well for demonstration purposes regarding the proposed NICE^h integration scheme.

1.1 Governing equations

Plastic potential Φ of the GTN model is given by

$$ \Upphi ({\varvec{\sigma}},\sigma_{\text{M}} ,f) = \left( {\frac{{\sigma_{\text{eq}} }}{{\sigma_{\text{M}} }}} \right)^{2} + 2\,f\,q_{1} \,\cosh \left( {\frac{{3\,q_{2} \,\sigma_{\text{H}} }}{{2\sigma_{\text{M}} }}} \right) \quad - (1 + q_{3} \,f^{2} ) $$

(21)

where $\sigma_{\text{eq}} = \sqrt {\tfrac{3}{2}{\user2{s}}:{\user2{s}}}$ is the von Mises equivalent stress with s being the deviatoric part of the stress tensor σ, $\sigma_{\text{M}} = \sigma_{\text{M}} (\bar{\varepsilon }_{m}^{\text{p}} )$ is the yield stress of the matrix material as a function of the equivalent plastic strain $\bar{\varepsilon }_{m}^{p}$, $\sigma_{\text{H}} = \tfrac{1}{3}{\varvec{\delta}}:{\varvec{\sigma}}$ is the hydrostatic stress with δ being the Kronecker tensor, f is the porosity and q ₁, q ₂ and q ₃ are the respective parameters of the damage model.

The law governing the porosity evolution considers two mechanisms, void growth and void nucleation, respectively

$$ {\text{d}}f = {\text{d}}f_{\text{nucleation}} + {\text{d}}f_{\text{growth}} $$

(22)

The first term on the right-hand side can be formulated by considering mass conservation

$$ {\text{d}}f_{\text{growth}} = (1 - f)\,{\text{d}}{\user2{\varepsilon }}^{\text{p}} :{\varvec{\delta}} $$

(23)

and the second one represents the nucleation of voids due to microcracking and decohesion of particle–matrix interface. It is related to the rate of plastic strain of the matrix material

$$ \begin{aligned} {\text{d}}f_{\text{nucleation}} & = A_{n} \,{\text{d}}\bar{\varepsilon }_{m}^{\text{p}} \\ A_{n} & = \frac{{f_{n} }}{{s_{n} \sqrt {2\pi } }}\exp \left[ { - \frac{1}{2}\left( {\frac{{\bar{\varepsilon }_{m}^{\text{p}} - \varepsilon_{n} }}{{s_{n} }}} \right)^{2} } \right] \\ \end{aligned} $$

(24)

where A _n follows a normal distribution about the mean nucleation strain ε _n with a standard deviation s _n . Parameter f _n describes the maximum nucleated void volume fraction. In this study, a decrease of strength of the material due to extensive void coalescence is omitted. For comparison reasons, we follow the formulation given in ABAQUS [3], where voids are nucleated only in tension.

For the elastic part of the constitutive model, Hooke’s law is adopted:

$$ {\varvec{\sigma}} = {\user2{C\varepsilon }}^{\text{e}} = {\user2{C}}({\user2{\varepsilon }} - {\user2{\varepsilon }}^{\text{p}} ) $$

(25)

It is assumed above that the total strain ε may be additively decomposed into elastic and plastic part, ε ^e and ε ^p, respectively:

$$ {\user2{\varepsilon }} = {\user2{\varepsilon }}^{\text{e}} + {\user2{\varepsilon }}^{\text{p}} $$

(26)

When using associated plasticity, the plastic strain rate can be derived as

$$ {\text{d}}{\user2{\varepsilon }}^{\text{p}} = \frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}\;{\text{d}}\lambda $$

(27)

where dλ is called plastic rate multiplier. The evolution of the equivalent plastic strain $\bar{\varepsilon }_{m}^{\text{p}}$ in the matrix material is obtained from the definition of equivalent plastic work:

$$ {\text{d}}\bar{\varepsilon }_{m}^{\text{p}} = \frac{{{\varvec{\sigma}}:\,{\text{d}}{\user2{\varepsilon }}^{\text{p}} }}{{(1 - f)\sigma_{\text{M}} }} $$

(28)

The equations given above, (21)–(28), determine the material behaviour of the GTN model. In partitioning of the GTN model according to thus-far presented mathematical structures of the constitutive models, we can separate the yield condition as the problem constraint equation

$$ \Upphi ({\varvec{\sigma}},\sigma_{\text{M}} ,f) = \left( {\frac{{\sigma_{\text{eq}} }}{{\sigma_{\text{M}} }}} \right)^{2} + 2\,f\,q_{1} \,\cosh \left( {\frac{{3\,q_{2} \,\sigma_{\text{H}} }}{{2\sigma_{\text{M}} }}} \right) \quad - (1 + q_{3} \,f^{2} ) = 0 $$

(29)

and assemble the corresponding evolution equations which are obtained by considering (22)–(28)

$$ {\text{d}}{\varvec{\sigma}} = {\user2{C}}\left( {{\text{d}}{\user2{\varepsilon }} - \frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}{\text{d}}\lambda } \right) $$

$$ {\text{d}}\sigma_{\text{M}} = \frac{{\partial \sigma_{\text{M}} }}{{\partial \bar{\varepsilon }_{m}^{\text{p}} }}{\text{d}}\bar{\varepsilon }_{m}^{\text{p}} $$

(30)

$$ {\text{d}}f = (1 - f)\left( {\frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}:{\varvec{\delta}}} \right){\text{d}}\lambda + A_{n} {\text{d}}\bar{\varepsilon }_{m}^{\text{p}} $$

where

$$ {\text{d}}\bar{\varepsilon }_{m}^{\text{p}} = \frac{{{\varvec{\sigma}}:\frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}{\text{d}}\lambda }}{{\left( {1 - f} \right)\,\sigma_{\text{M}} }} $$

(31)

Appendix B

With reference to the NICE^h integration scheme, we show in this section a procedure of deriving explicit expressions for the plastic rate multiplier Δλ. Only derivation for h = 1 and h = 2 is carried out explicitly, because the procedure is the same for any h. The only difference regards the number of terms in the truncated Taylor series expansion and the related degree of the polynomial equation to be solved.

Derivation for h = 1:

$$ \Upphi + \frac{\partial \Upphi }{\partial \sigma }\Updelta \sigma + \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} }}\Updelta \sigma_{\text{M}} + \frac{\partial \Upphi }{\partial f}\Updelta f = 0 $$

(32)

Considering the evolution equations (12) in (32) and then extracting Δλ we can simply get

$$ \Updelta \lambda = \frac{{\Upphi + \frac{\partial \Upphi }{\partial \sigma }\,E\,\Updelta \varepsilon }}{{\frac{\partial \Upphi }{\partial \sigma }E\frac{\partial \Upphi }{\partial \sigma } - \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} }}\frac{{\partial \sigma_{\text{M}} }}{{\partial \bar{\varepsilon }_{m}^{\text{p}} }}\frac{{\sigma \frac{\partial \Upphi }{\partial \sigma }}}{{(1 - f)\sigma_{\text{M}} }} - \frac{\partial \Upphi }{\partial f}\left( {\left( {1 - f} \right)\frac{\partial \Upphi }{{\partial \sigma_{\text{H}} }} + A_{n} \frac{{\sigma \frac{\partial \Upphi }{\partial \sigma }}}{{(1 - f)\sigma_{\text{M}} }}} \right)}} $$

(33)

The similarity of the last equation with (15) is obvious.

Derivation for h = 2:

$$ \begin{aligned} & \Upphi + \frac{\partial \Upphi }{\partial \sigma }\Updelta \sigma + \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} }}\Updelta \sigma_{\text{M}} + \frac{\partial \Upphi }{\partial f}\Updelta f + \frac{1}{2}\left( {\frac{{\partial^{2} \Upphi }}{{\partial \sigma^{2} }}\Updelta \sigma^{2} + \frac{{\partial^{2} \Upphi }}{{\partial \sigma_{\text{M}}^{2} }}\Updelta \sigma_{\text{M}}^{ 2} + \frac{{\partial^{2} \Upphi }}{{\partial f^{2} }}\Updelta f^{2} } \right) \\ & \quad + \frac{{\partial^{2} \Upphi }}{{\partial \sigma \,\partial \sigma_{\text{M}} }}\Updelta \sigma \,\Updelta \sigma_{\text{M}} + \frac{{\partial^{2} \Upphi }}{\partial \sigma \,\partial f}\Updelta \sigma \,\Updelta f + \frac{{\partial^{2} \Upphi }}{{\,\partial \sigma_{\text{M}} \,\partial f}}\Updelta \sigma_{\text{M}} \,\Updelta f = 0. \\ \end{aligned} $$

(34)

When considering the evolution equations (12) in (34), a second-order polynomial equation in Δλ is obtained, whose roots are

$$ \Updelta \lambda_{12} = \frac{{ - b \pm \sqrt {b^{2} - 4ac} }}{2a} $$

(35)

where

$$ \begin{aligned} a & = \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\partial \sigma_{\text{M}}^{ 2} }}\left( {\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{p}^{\text{M}} }}A} \right)^{2} + \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\partial f^{2} }}B^{2} - \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\partial \sigma^{2} }}\left( {\frac{\partial \Upphi }{\partial \sigma }} \right)^{2} \\ & \quad - \frac{{\partial^{2} \Upphi }}{{\partial \sigma \,\partial \sigma_{\text{M}} }}\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{p}^{\text{M}} }}A\frac{\partial \Upphi }{\partial \sigma } + \frac{{\partial^{2} \Upphi }}{\partial \sigma \,\partial f}B\frac{\partial \Upphi }{\partial \sigma } + \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} \,\partial f}}\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{\text{p}}^{\text{M}} }}A\,B \\ \end{aligned} $$

(36)

$$ \begin{aligned} b & = \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} }}\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{\text{p}}^{\text{M}} }}A - \left( {\frac{\partial \Upphi }{\partial \sigma }} \right)^{2} + \frac{\partial \Upphi }{\partial f}B + \frac{\partial \Upphi }{{\partial \sigma \,\partial \sigma_{\text{M}} }}\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{\text{p}}^{\text{M}} }}A \\ & \quad + \frac{\partial \Upphi }{\partial \sigma \,\partial f}B\,E\,\Updelta \varepsilon + \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} \,\partial f}}\frac{{\partial \sigma_{\text{M}} }}{{\partial \varepsilon_{\text{p}}^{\text{M}} }}A\,B \\ \end{aligned} $$

(37)

$$ c = \Upphi + \frac{\partial \Upphi }{\partial \sigma }E\Updelta \varepsilon + \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\partial \sigma^{2} }}(E\Updelta \varepsilon )^{2} $$

(38)

and

$$ \begin{aligned} A & = \frac{{\sigma \frac{\partial \Upphi }{\partial \sigma }}}{{(1 - f)\sigma_{\text{M}} }} \\ B & = (1 - f)\frac{\partial \Upphi }{{\partial \sigma_{\text{H}} }} + A_{n} \frac{{\sigma \frac{\partial \Upphi }{\partial \sigma }}}{{(1 - f)\,\sigma_{\text{M}} }}. \\ \end{aligned} $$

(39)

Appendix C

Groundings for a derivation of the plastic rate multiplier for the GTN model for the second-order NICE scheme, that is NICE², are presented here. According to the truncated Taylor series expansion (9), the derivation starts with

$$ \Upphi + {\text{d}}\Upphi + {\text{d}}^{2} \Upphi = 0 $$

(40)

which leads to

$$ \begin{gathered} \underbrace {{\Upphi + \frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}:{\user2{\Updelta \sigma }} + \frac{\partial \Upphi }{{\partial \sigma_{\text{M}} }}\Updelta \sigma_{\text{M}} { + }\frac{\partial \Upphi }{\partial f}\Updelta f}}_{{{\text{grounding}}\;{\text{for}}\;{\text{first}} {\hbox{-}} {\text{order}}\;{\text{NICE}}^{ 1} }} \hfill\\ +\underbrace{{\left( {{\varvec{\alpha}} \odot \frac{{\partial^{2} \Upphi }}{{\partial {\varvec{\sigma}}\,\partial {\varvec{\sigma}}}}} \right):({\user2{\Updelta \sigma }} \otimes {\user2{\Updelta \sigma }}) + \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\partial f^{2} }}\,\Updelta f^{2} + \frac{1}{2}\frac{{\partial^{2} \Upphi }}{{\,\partial \sigma_{\text{M}}^{2} }}\Updelta \sigma_{\text{M}}^{2} + \frac{{\partial^{2} \Upphi }}{{\partial {\varvec{\sigma}}\,\partial f}} :{\user2{\Updelta \sigma }}\,\Updelta f + \frac{{\partial^{2} \Upphi }}{{\partial {\varvec{\sigma}}\,\partial \sigma_{\text{M}} }} :{\user2{\Updelta \sigma }}\,\Updelta \sigma_{\text{M}} }}_{{{\text{second}}\;{\text{order differential for second}}{\hbox{-}}{\text{order}}\;{\text{NICE}}^{ 2} }} = 0 \hfill \\ \end{gathered} $$

(41)

where

$$ {\varvec{\alpha}} = \alpha_{ijkl} = \left\{ \begin{array}{ll} \frac{1}{2}{\text{when}} & i = k\,{\text{and}}\,j = l \\ 1&{\text{otherwise}} \\ \end{array} \right. $$

Considering the evolution eqautions (30) and (41), one can extract from the latter the increment of the plastic rate multiplier Δλ as a function of the state variables at the beginning of the increment and the prescribed strain increment Δε.

Appendix D: Stability limit of NICE^h

The sections above focus on the demonstration of accuracy of the NICE^h scheme. It is shown there that explicit and implicit integration of the flow equations lead to different results but are comparable with regard to the accuracy. Besides accuracy, also the stability of numerical integration plays an important role. Thus, in this appendix the NICE^h scheme is evaluated regarding its stability limit.

In [17], the stability of integration algorithms is concisely discussed with emphasis being given on the integration of the plastic flow evolution. In order to generalize the stability analyses, generalized midpoint rule is used

$$ ({\user2{\varepsilon }}^{\text{p}} )^{i + 1} = ({\user2{\varepsilon }}^{\text{p}} )^{i} + \left( {(1 - \alpha )\left( {\frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}} \right)^{i} { + }\alpha \left( {\frac{\partial \Upphi }{{\partial {\varvec{\sigma}}}}} \right)^{{i{ + 1}}} } \right)\Updelta \lambda . $$

(42)

With a proper choice of α (α; 0 ≤ α ≤ 1) in (42) different integration rules of the flow equation can be obtained; with α = 0 and α = 1 the CFE and CBE integration rules are restored, respectively. It is proven in the article that regardless of the yield surface shape the flow rule integration is unconditionally stable if α ≥ 0.5, otherwise integration of the flow rule is only conditionally stable [17]. For conditionally stable algorithms it is necessary, that the applied integration step remains lower than the stability limit step.

As already stated, when boundary value problem is solved with a direct solution technique (explicit dynamics), the whole integration must be performed with an increment size, which enables stable integration of the dynamics equations. The same incrementation is then applied also for the integration of the constitutive equations. Thus, the key question is whether the time increment, as prescribed by the global integration scheme, is sufficiently small for stable integration of the constitutive model with the NICE^h scheme. The size of stable increment for integration of constitutive laws for conditionally stable schemes strongly depends on the chosen scheme and yield surface curvature [17]. Thus, for a precise stability analysis the size of stable increment must be calculated for each stress state independently. In this work only an estimation and approximate comparison of the stable increment of the NICE^h scheme and prescribed increment will be presented.

Let us reconsider the second numerical case, i.e. deep drawing of the cylindrical cup with the punch velocity being 5 m/s, from the stability point of view. Specifically, let us observe a finite element on the upper layer at radius R = 60 mm from the initial blank centre. The element strain history during the drawing process, which comprises its passing along the die radius, can be extracted in a discretised manner from its integration point, when Abaqus-default integration algorithm was used. In order to be able to arbitrarily change the integration step and observe the stability limit of the NICE^h scheme the integration of constitutive equations under such prescribed strain path was then performed in Mathematica software [14]. Figure 8 displays the stress versus strain component in the radial direction, σ_rr versus ε_rr, during the whole integration path. The integration was performed with the NICE¹ scheme with 500, 750 and 10,000 increments, respectively.

Fig. 8

Calculated σ_rr versus ε_rr with NICE¹ in cup drawing simulation

As it can be observed, the integration of stresses with 500 increments becomes unstable at some point, with the error regarding stable result starting to grow rapidly. If the number of increments is increased to 750 increments, the performed integration becomes stable. For this case the stability limit, governed by the integration of flow rule, can be thus estimated through the number of increments, which are demanded for stable integration in the cup drawing simulation.

In elasto-plastic engineering simulation with the boundary value problem considered as a dynamics problem, which is solved with a direct solution technique, the simulations are typically run with 10⁴–10⁶ increments. Considering the above estimated number of increments (750), which are required for stable integration of the constitutive model, one can conclude, that the stable increment of the NICE scheme is large enough for engineering simulations.