# An energy-dissipative momentum-conserving time integration algorithm for viscoplastic hyperelastic models based on the variational updates framework

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot

► To cite this version:

Ludovic Noels, Laurent Stainier, Jean-Philippe Ponthot. An energy-dissipative momentum-conserving time integration algorithm for viscoplastic hyperelastic models based on the variational updates framework. IMPLAST'07, Symposium on Plasticity and Impact Mechanics, 2007, Bochum, Germany. hal-01007857

HAL Id: hal-01007857

https://hal.science/hal-01007857

Submitted on 2 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AN ENERGY-DISSIPATIVE MOMENTUM-CONSERVING TIME INTEGRATION ALGORITHM FOR VISCOPLASTIC HYPERELASTIC MODELS BASED ON THE VARIATIONAL UPDATES FRAMEWORK

L. Noels<fnref n="1" />, L. Stainier¹and J.-Ph. Ponthot<fnref n="1" />

**Keywords:** Finite-element, non-linear dynamics, energy dissipation, time integration, variational updates

**Abstract** Traditional time integration algorithms for finite-element discretization are numerically stable only for linear models. To overcome that drawback, a new class of "Energy Momentum Conserving Algorithms" or EMCA and "Energy-Dissipative Momentum-Conserving" or EDMC algorithms, verifying the conservation/dissipation laws in the non-linear range, was recently developed. They consist in a mid-point scheme with an adequate evaluation of the internal forces, which strongly assume the existence of an energy potential. Following the variational visco-plastic constitutive updates formalism, a general energy-dissipative momentum-conserving algorithm for hyperelastic-based formulation is proposed. The main feature of this variational formulation is that the stress tensor always derives from an incremental potential, even if plastic deformations occur, allowing to use it in the above integration algorithms.

## 1 Introduction

The most widely used time integration algorithms for finite-element systems, based on the Newmark family, are numerically stable only for linear models. To overcome that drawback, a new class of Energy Momentum Conserving Algorithms or EMCA, verifying the conservation laws in the non-linear range, was initiated by Simo and Tarnow [13]. These algorithms consist in a mid-point scheme with an adequate evaluation of the internal forces. This adequate evaluation was given for hyperelastic models by Gonzalez. The EMCA was then extended to dynamic finite deformation plasticity based on a hyperelastic model by Meng and Laursen [5], by Armero [1], and to dynamic finite deformation plasticity based on a hypoelastic model by the present authors [8, 9, 10].

Recently, the authors [7] have proposed a general energy-momentum conserving algorithm for hyperelastic-based formulation, using the variational visco-plastic constitutive updates proposed by Ortiz and Stainier [9]. The main feature of this formulation is that the stress tensor always derives from an incremental potential, even if plastic deformations occur. Therefore, in such a framework we can use the formulation based on the second Piola-Kirchhoff stress tensor as proposed by Gonzalez without any modification.

<footnote n="1">¹University of Liège. Belgium; emails: l.noels@ulg.ac.be, l.stainier@ulg.ac.be, jp.ponthot@ulg.ac.be</footnote>
To avoid the lack of convergence due to the presence of high frequency modes, numerical dissipation was introduced in the conserving algorithms by Arnero and Romero for hyperelastic materials. An extension to numerical dissipation based on the use of the variational visco-plastic constitutive updates is proposed in the present work. This algorithm is shown to be stable, robust and accurate.

## 2 Variational constitutive updates

Let us consider a time step increment $\Delta t = t^{n+1} - t^n$. The purpose of this formalism is to define an incremental potential $\Delta D^{\text{eff}}(\vec{x}_0, \bar{C}(\vec{\phi}_h, \theta^e))$ depending on the initial positions $\vec{x}_0$ and on the modified right-Cauchy deformation tensor $\bar{C} = \sqrt[3]{\frac{\theta^{e2}}{\det(\bar{C})}} \bar{C}$, which is written in terms of the classical Cauchy strain tensor, this last one depending on the deformation mapping $\vec{\phi}_h$, and of the volume variation $\theta^e$ in order to use the quasi incompressible formulation proposed by Simo and Taylor [14]. Similarly, the deviatoric strain tensor is defined by $\hat{C} = \sqrt[3]{\frac{1}{\det(\hat{C})}} \hat{C}$.

The deformation gradient $F = F^{\text{el}} F^{\text{pl}}$ is decomposed into the elastic part $F^{\text{el}}$ and the plastic part

$$F^{\text{pl}n+1} = \exp\left(\left[\varepsilon^{\text{pl}n+1} - \varepsilon^{\text{pl}}\right] N\right) F^{\text{pl}n}, \quad (1)$$

where tensor $N$ is the plastic flow direction, and $\varepsilon^{\text{pl}}$ is the equivalent plastic strain. The free Helmholtz energy $A$ is obtained from the sum of the volume part of the elastic potential $\Phi^{\text{vol}}$, of the deviatoric part $\Phi^{\text{dev}}$, and of the plastic stored energy $\Phi^{\text{pl}}$, leading to

$$A(\bar{F}(\vec{\phi}_h, \theta^e), \varepsilon^{\text{pl}}, N) = \Phi^{\text{vol el}}(\theta^e) + \Phi^{\text{dev el}}(\hat{C}^{\text{el}}) + \Phi^{\text{pl}}(F^{\text{pl}}(\varepsilon^{\text{pl}}, N), \varepsilon^{\text{pl}}), \quad (2)$$

where it has been taken into account that the plastic deformation gradient $F^{\text{pl}}$ depends only on the deviatoric deformation. The incremental potential $\Delta D$ is defined on a time step as

$$\Delta D(\bar{F}(\vec{\phi}_h^{n+1}, \theta^{cn+1}), \varepsilon^{pn+1}, N) \equiv \Phi^{\text{vol el}}(\theta^{cn+1}) - \Phi^{\text{vol el}}(\theta^{cn}) + \Delta D_{\text{dev}}(\hat{F}, \varepsilon^{\text{pl}}, N), \text{ with} \quad (3)$$

$$\begin{aligned} \Delta D_{\text{dev}}(\hat{F}, \varepsilon^{\text{pl}}, N) ={}& \Phi_{\text{dev}}^{\text{el}}(\hat{C}^{\text{el}n+1}(\hat{F}^{n+1}, \varepsilon^{\text{pl}n+1}, N)) - \Phi_{\text{dev}}^{\text{el}}(\hat{C}^{\text{el}n}) + \\ & \Phi^{\text{pl}}(F^{\text{pl}n+1}(\varepsilon^{\text{pl}n+1}, N), \varepsilon^{\text{pl}n+1}) - \Phi^{\text{pl}}(F^{\text{pl}n}, \varepsilon^{\text{pl}n}) + \Delta t \Psi^* \left(\frac{\Delta \varepsilon^{\text{pl}}}{\Delta t}\right), \end{aligned} \quad (4)$$

where $\Psi^*$ is a dissipative potential that ensures convexity of $\Delta D$. It can be shown that the stationary point of $\Delta D$ with respect to $\varepsilon^{\text{pl}n+1}$ and $N$ constitutes the solution of the elasto-plastic problem [9, 7]. Therefore the effective incremental potential $\Delta D^{\text{eff}}$ can be defined as the value of the incremental potential at this stationary point. But in order to introduce numerical dissipation, a dissipative volume potential $\Psi^{\text{vol}}$ and deviatoric
potential $\Psi_{dev}^d$ are added, leading to the new incremental effective potential

$$ \begin{aligned} \Delta D^{\text{eff}}(\vec{x}_0, \vec{F}(\vec{\varphi}_h^{n+1}, \theta^{n+1})) &= \min_{\varepsilon^{\text{p}n+1}, N} \Delta D + \Delta t \Psi_{\text{vol}}^{\text{d}} \left( \frac{\Delta \theta^e}{\Delta t} \right) + \Delta t \Psi_{\text{dev}}^{\text{d}} \left( \frac{\Delta \hat{C}^{\text{el}}}{\Delta t} \right), \end{aligned} \quad (5) $$

Since the minimization process is independent of $\theta^e$, this effective incremental potential can be split into a volume part $\Delta D_{vol}^{\text{eff}}$ and a deviatoric part $\Delta D_{dev}^{\text{eff}}$ defined by

$$ \Delta D_{vol}^{\text{eff}}(\theta^{\text{en}+1}) = \Phi_{vol}^{\text{el}}(\theta^{\text{en}+1}) - \Phi_{vol}^{\text{el}}(\theta^{\text{en}}) + \Delta t \Psi_{vol}^{\text{d}} \left( \frac{\Delta \theta^e}{\Delta t} \right), \quad \text{and} \quad (6) $$

$$ \Delta D_{\text{dev}}^{\text{eff}}(\hat{F}^{n+1}) = \min_{\varepsilon^{\text{p}n+1}, N} \left[ \Delta D_{\text{dev}}(\hat{F}, \varepsilon^{\text{p}}, N) \right] + \Delta t \Psi_{\text{dev}}^{\text{d}} \left( \frac{\Delta \hat{C}^{\text{el}}}{\Delta t} \right). \quad (7) $$

From these definitions, the second Piola-Kirchhoff stress tensor can be deduced by explicit derivation

$$ S^{n+1} = S_{\text{dev}}^{n+1} + p^{n+1} J^{n+1} (C^{n+1})^{-1} \quad \text{with} \quad (8) $$

$$ p^c = \frac{\partial \Delta D_{vol}^{\text{eff}}}{\partial \theta}, \quad \text{and} \quad (9) $$

$$ S_{\text{dev}} = 2 J^{-\frac{2}{3}} \text{DEV} \left( (F^{\text{pl}})^{-1} \left[ \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} + \Delta t \frac{\partial \Psi_{\text{dev}}^{\text{d}}}{\partial \hat{C}^{\text{el}}} \right] (F^{\text{pl}})^{-T} \right). \quad (10) $$

Details of this formulation can be found in [11].

### 3 The EDMC time integration algorithm

The EDMC algorithm described by Armero and Romero [2] can be stated in terms of the nodal position $\vec{x}^\xi$ and appropriate external/internal forces $\vec{f}_{\text{ext}}^\xi = \vec{f}_{\text{int}}^\xi$ as

$$ \vec{x}^{\xi n+1} = \vec{x}^{\xi n} + \Delta t \left[ \frac{\dot{\vec{x}}^{\xi n+1} + \dot{\vec{x}}^{\xi n}}{2} + \vec{v}_{\text{diss}}^{\xi n+1} \right], \quad (11) $$

$$ \dot{\vec{x}}^{\xi n+1} = \dot{\vec{x}}^{\xi n} + \Delta t \frac{\ddot{\vec{x}}^{\xi n+1} + \ddot{\vec{x}}^{\xi n}}{2}, \quad \text{and} \quad (12) $$

$$ M^{\xi \mu} \frac{\ddot{\vec{x}}^{\mu n} + \ddot{\vec{x}}^{\mu n+1}}{2} = \vec{f}_{\text{ext}}^{\xi n+1/2} - \vec{f}_{\text{int}}^{\xi n+1/2}, \quad (13) $$

with the first order dissipation velocities

$$ \vec{v}_{\text{diss}}^{\xi n+1/2} = \chi \frac{\|\dot{\vec{x}}^{\xi n+1}\| - \|\dot{\vec{x}}^{\xi n}\|}{\|\dot{\vec{x}}^{\xi n+1}\| + \|\dot{\vec{x}}^{\xi n}\|} \frac{\dot{\vec{x}}^{\xi n+1} + \dot{\vec{x}}^{\xi n}}{2}, \quad (14) $$

where $\chi$ is the control parameter.
In the context of Energy-Momentum Conserving Algorithm, Gonzalez [4] proposed a new expression of the internal forces. Given the incremental updates variational formulation developed in the previous section, Gonzalez expression can be directly applied, leading to

$$ (\vec{f}_{\text{int}}^{\xi})^{n+\frac{1}{2}} = \int_{V_0} \left\{ \frac{F^{n+1} + F^n}{2} \left[ (S_{\text{dev}}^{*})^{n+\frac{1}{2}} + 2(p^{*})^{n+\frac{1}{2}} dG \right] \vec{D}^{\xi} \right\} dV_0, \quad (15) $$

where $S_{\text{dev}}^{*}$ is the consistent deviatoric stress tensor, $p^{*}$ is the consistent pressure and $dG$ is the modified differentiation of $J$ by $C$. This last term is defined by

$$ dG = \frac{\partial J^{n+\frac{1}{2}}}{\partial C} + \left[ \frac{J^{n+1} - J^n - \frac{\partial J^{n+\frac{1}{2}}}{\partial C} : \Delta C}{\|\Delta C\|^2} \right] \Delta C \quad \text{with} \quad (16) $$

$$ \frac{\partial J^{n+\frac{1}{2}}}{\partial C} = \frac{1}{2} \sqrt{\det \left( \frac{C^{n+1} + C^n}{2} \right) \left[ \frac{C^{n+1} + C^n}{2} \right]^{-1}}, \quad \text{and} \quad (17) $$

$$ \Delta C = C^{n+1} - C^n. \quad (18) $$

Instead of using directly Eq. (9), the modified pressure is computed as

$$ p^{*n+\frac{1}{2}} = \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}} + \frac{\Delta D_{\text{vol}}^{\text{eff}}(\theta^{e n+1}, \theta^e) - \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}} \Delta \theta^e}{\Delta \theta^{e2}} \Delta \theta^e, \quad (19) $$

$$ \left[ \frac{\partial \Delta D_{\text{vol}}^{\text{eff}}}{\partial \theta^e} \right]^{n+\frac{1}{2}} = \frac{\partial \Phi_{\text{vol}}^{\text{el}}}{\partial \theta^e} \left( \frac{\theta^{e n+1} + \theta^e}{2 \Delta t} \right) + \Delta t \frac{\partial \Psi_{\text{vol}}^{\text{d}}}{\partial \Delta \theta^e} \left( \frac{\Delta \theta^e}{2} \right), \quad \text{and} \quad (20) $$

$$ \Delta \theta^e = \theta^{e n+1} - \theta^e. \quad (21) $$

Similarly, Eq. (10) becomes

$$ S_{\text{dev}}^{*} \stackrel{n+\frac{1}{2}}{=} 2 \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial C} \right]^{n+\frac{1}{2}} + 2 \frac{\Delta D_{\text{dev}}^{\text{eff}}(\hat{C}^{n+1}, \hat{C}^n) - \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial C} \right]^{n+\frac{1}{2}} : \Delta C}{\|\Delta C\|^2} \Delta C, \quad (22) $$

$$ \left[ \frac{\partial \Delta D_{\text{dev}}^{\text{eff}}}{\partial C} \right]^{n+\frac{1}{2}} = J^{-\frac{2}{3}} \text{DEV} \left( (F^{\text{pl}})^{-1} \frac{\partial \Phi_{\text{dev}}^{\text{el}}}{\partial \hat{C}^{\text{el}}} \left( \frac{(\hat{C}^{\text{el}})^{n+1} + (\hat{C}^{\text{el}})^n}{2} \right) (F^{\text{pl}})^{-T} \right) + \Delta t J^{-\frac{2}{3}} \text{DEV} \left( (F^{\text{pl}})^{-1} \frac{\partial \Psi_{\text{dev}}^{\text{d}}}{\partial \hat{C}^{\text{el}}} \left( \frac{(\hat{C}^{\text{el}})^{n+1} - (\hat{C}^{\text{el}})^n}{2 \Delta t} \right) (F^{\text{pl}})^{-T} \right), \quad (23) $$

$$ \text{and } \Delta C = C^{n+1} - C^n. \quad (24) $$

Since these last expressions are similar to the ones developed for elastic dynamics, the conserving properties are demonstrated in a straightforward way, see [11] for details.
<figure><img src="image_2.png" /><figcaption>Figure 1: Description of Taylor's impact test: a) geometry and comparison of the results ($t = 80$ µs) for the Taylor's impact using: b-d) EDMC scheme with a spectral radius evolving from $ρ_∞ = 0.9$ to $ρ_∞ = 0.65$; e) second order Chung-Hulbert time integration scheme; f) first order Newmark time integration scheme.</figcaption></figure>

## 4 Numerical example: Taylor's impact test

The purpose of this test is to demonstrate the accuracy of the scheme when large plastic deformations are involved. The bar described in Fig. 1a impacts a rigid wall with an initial velocity of 227 m·s⁻¹. We consider rate-independent plasticity with linear hardening. Material properties are reported in Table 1.
<table><caption>Table 1: Material properties for the Taylor's impact test.</caption><thead><tr><td>Property</td><td>Value</td></tr></thead><tbody><tr><td>Density</td><td>ρ = 8930 kg·m<sup>-3</sup></td></tr><tr><td>Bulk modulus</td><td>K<sub>0</sub> = 130000 N·mm<sup>-2</sup></td></tr><tr><td>Shear modulus</td><td>G<sub>0</sub> = 43333 N·mm<sup>-2</sup></td></tr><tr><td>Yield stress</td><td>σ<sub>0</sub><sup>vm</sup> = 400 N·mm<sup>-2</sup></td></tr><tr><td>Hardening modulus</td><td>h = Y<sub>0</sub>/ε<sup>p</sup><sub>0</sub> = 100 N·mm<sup>-2</sup></td></tr></tbody></table>

This simulation is solved with a time step $Δt = 0.4$ µs, and solutions obtained with the EDMC algorithm for different spectral radius are compared to the solutions obtained with the Newmark [6] and the Chung-Hulbert [3] time integration algorithms. Final deformations are illustrated in Fig. 1b-f. All the simulations give the same solution within 5% error, except the modified EDMC one, which leads to a different solution. This demonstrates the accuracy of the EDMC algorithms.

## 5 Conclusions

In this paper, a new Energy-Dissipative Momentum-Conserving time integration algorithm, based on an incremental variational updates formulation of the elasto-plastic behavior, is proposed. Owing to this formulation, an effective potential can always be defined, even when irreversible deformations occur. The main originality of this paper is to add a dissipative function to this incremental potential, which leads to a new definition of a potential for the system.

The resulting algorithm verifies the conservation of linear and angular momentum and ensures that the numerical dissipation is always positive. An impact test also demonstrates the accuracy and robustness of the algorithm when large plastic deformations occur.

## References

[1] F. Armero. An energy-consistent material-point method for dynamic finite deformation plasticity. *Int. J. Numer. Meth. Eng.*, 65:1608-1638, 2006.

[2] F. Armero and I. Romero. On the formulation of high-frequency dissipative time-stepping algorithms for non-linear dynamics. Part I: Low-order methods for two model problems and nonlinear elastodynamics. *Comput. Meth. Appl. Mech. Eng.*, 190:2603-2649, 2001.

[3] J. Chung and G. J. M. Hulbert. A time integration algorithm for structural dynamics with improved numerical dissipations: the generalized-α method. *J. Appl. Mech.*, 60:371-375, 1993.

[4] O. Gonzalez. Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Comput. Meth. Appl. Mech. Eng.*, 190:1763-1783, 2000.

[5] X. Meng and T. Laursen. Energy consistent algorithms for dynamic finite deformation plasticity. *Comput. Meth. Appl. Mech. Eng.*, 191:1639-1675, 2001.

[6] N. Newmark. A method of computation for structural dynamics. *J. Eng. Mech. Div. ASCE*, 85(EM3):67-94, 1959.
[7] F. Noels, L. Stainier, and J.-P. Ponthot. An energy momentum conserving algorithm using the variational formulation of visco-plastic updates. *Int. J. Numer. Meth. Eng.*, 65:904-942, 2006.

[8] L. Noels, L. Stainier, and J.-P. Ponthot. Energy-momentum conserving algorithm for non-linear hypoelastic constitutive models. *Int. J. Numer. Meth. Eng.*, 59:83-114, 2004.

[9] L. Noels, L. Stainier, and J.-P. Ponthot. On the use of large time steps with an energy-momentum conserving algorithm for non-linear hypoelastic constitutive models. *Int. J. Solids Struct.*, 41:663-693, 2004.

[10] L. Noels, L. Stainier, and J.-P. Ponthot. Simulation of complex impact problems with implicit time algorithm. Application to a blade-loss problem. *Int. J. Impact Eng.*, 32:358-386, 2005.

[11] L. Noels, L. Stainier, and J.-P. Ponthot. A first order energy-dissipative momentum conserving scheme for elasto-plasticity using the variational updates formulation. *Comput. Meth. Appl. Mech. Eng.*, 2007. Submitted.

[12] M. Ortiz and L. Stainier. The variational formulation of viscoplastic updates. *Comput. Meth. Appl. Mech. Eng.*, 171:419-444, 1999.

[13] J. C. Simo and N. Tarnow. The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics. *Z. Angew. Math. Phys.*, 43:757-792, 1992.

[14] J. C. Simo and R. L. Taylor. Quasi-incompressible finite elasticity in principal stretches continuum basis and numerical algorithms. *Comput. Meth. Appl. Mech. Eng.*, 85:273-310, 1991.