# Un schéma direct d'intégration numérique de modèles visco-hyperélastiques par relaxation radiale

Stéphane Lejeunes, Stéphane Méo, Adnane Boukamel

► To cite this version:

Stéphane Lejeunes, Stéphane Méo, Adnane Boukamel. Un schéma direct d'intégration numérique de modèles visco-hyperélastiques par relaxation radiale. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421809

HAL Id: hal-01421809

https://hal.science/hal-01421809v1

Submitted on 22 Dec 2016

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 CC BY 4.0 - Attribution - International License
# Un schéma direct d'intégration numérique de modèles visco-hyperélastiques par relaxation radiale

S. Lejeunes¹, S. Méo², A. Boukamel¹,³

¹ Laboratoire de Mécanique et d'Acoustique de Marseille (CNRS UPR7051)
31 ch Joseph-Aiguier, 13402 Marseille
lejeunes@lma.cnrs-mrs.fr

² Laboratoire de Mécanique et Rhéologie (Université de Tours)
EPU Dept Productique, BP 407, F37204 Cedex 3, Tours
stephane.meo@univ-tours.fr

³ Ecole Centrale Marseille
Technopôle Château Gombert, 13451 Marseille
adnane.boukamel@ec-marseille.fr

**Résumé** — On propose dans cette communication un schéma numérique d'intégration de lois d'évolutions visco-hyperélastiques s'inspirant à la fois des méthodes par retour en exponentiel [5] et des méthodes par retour radial [8, 7]. L'originalité de ce schéma repose à la fois sur la substitution d'un système différentiel tensoriel, traduisant les lois d'évolutions, par un système différentiel scalaire de deux équations à deux inconnues, et sur un développement limité au premier ordre de ces variables. On obtient ainsi une solution analytique approchée en exponentiel.

**Mots clés** — relaxation radiale, visco-hyperélasticité, éléments-finis.

## 1 Introduction

La modélisation des matériaux polymères à l'aide de modèles rhéologiques, s'appuyant sur la notion d'états intermédiaires, est de nos jours un sujet largement documenté dans la littérature. De nombreux modèles sans cesse plus complexes ont été proposés afin de reproduire au mieux le comportement dissipatif de ces matériaux. Pour autant l'implémentation numérique de ce type de modèles conduit bien souvent à des algorithmes lourds à mettre en oeuvre et coûteux en terme de temps de calcul.

L'idée développée dans ce papier consiste à approximer le système différentiel, tensoriel et non-linéaire obtenu à partir de la loi d'évolution de modèles visco-hyperélastiques par un système beaucoup plus simple puisque linéaire et scalaire. L'obtention d'une telle approximation s'inscrit dans le cadre des algorithmes de retour radial et de retour en exponentiel développés en plasticité et viscoplasticité [8, 7, 4] ainsi qu'en viscoélasticité [5]. Ce schéma est qualifié de directe puis-qu'une solution analytique en exponentielle à l'aide d'une série de Magnus peut être déterminée et évaluée avec un coût numérique plus faible que des méthodes itératives. Un exemple numérique, permettant de discuter des approximations effectuées est présenté dans ce papier.
## 2 Equations constitutives

En utilisant le concept des états intermédiaires [6], on peut décomposer le gradient de la transformation F en parties incompressible et compressible tel que : $\bar{F} = J^{-\frac{1}{2}}F$ avec $J = \det F$, on décompose ensuite le gradient incompressible en parties visqueuse et élastique :

$$ \bar{F} = \bar{F}_e \cdot \bar{F}_v. \qquad (1) $$

Dans le cadre d'un formalisme Eulérien, on peut choisir les tenseurs de dilatation de Cauchy-Green gauche $\mathbf{B}$ et $\bar{\mathbf{B}}_e$ comme variable d'états indépendantes. On suppose également une décomposition additive de l'énergie libre tel que :

$$ \psi(\bar{\mathbf{B}}, \bar{\mathbf{B}}_e, J) = \psi_v(\bar{\mathbf{B}}_e) + \psi_h(\bar{\mathbf{B}}) + \psi_{vol}(J) \qquad (2) $$

En partant de (2), la variation de l'énergie libre spécifique s'écrit donc :

$$ \dot{\psi} = \frac{\partial \psi_h}{\partial \bar{\mathbf{B}}} : \frac{\partial \bar{\mathbf{B}}}{\partial \mathbf{B}} : \dot{\mathbf{B}} + \frac{\partial \psi_v}{\partial \bar{\mathbf{B}}_e} : \dot{\mathbf{B}}_e + \frac{\partial \psi_{vol}}{\partial J} j \qquad (3) $$

$$ \dot{\mathbf{B}} = \mathbf{L} \cdot \mathbf{B} + \mathbf{B} \cdot \mathbf{L}^T \qquad j = J(\mathbf{1} : L) \qquad (4) $$

$$ \dot{\bar{\mathbf{B}}}_e = \mathbf{L} \cdot \bar{\mathbf{B}}_e + \bar{\mathbf{B}}_e \cdot \mathbf{L}^T - 2\bar{\mathbf{V}}_e \cdot \bar{\mathbf{D}}_v^0 \cdot \bar{\mathbf{V}}_e - \frac{2}{3}(\mathbf{1} : \mathbf{D})\bar{\mathbf{B}}_e \qquad (5) $$

ou $\mathbf{L} = \dot{\mathbf{F}} \cdot \mathbf{F}^{-1}$ et $\bar{\mathbf{V}}_e$ est le tenseur de déformation pure (i.e. $\bar{\mathbf{F}}_e = \bar{\mathbf{V}}_e \cdot \bar{\mathbf{R}}$). On peut définir également un taux de déformation anélastique objectif<fnref n="1" /> tel que :

$$ \bar{\mathbf{D}}_v^0 = \bar{\mathbf{R}} \cdot \bar{\mathbf{D}}_v \cdot \bar{\mathbf{R}}^T = (\bar{\mathbf{V}}_e^{-1} \cdot \dot{\mathbf{F}} \cdot \mathbf{F}^{-1} \cdot \bar{\mathbf{V}}_e)_{\text{sym}} - (\dot{\bar{\mathbf{V}}}_e \cdot \bar{\mathbf{V}}_e^{-1})_{\text{sym}} \qquad (6) $$

Dans ce cadre Eulérien et dans un contexte purement mécanique, l'inégalité de Clausius-Duhem<fnref n="2" /> peut s'écrire :

$$ \phi = \sigma : \mathbf{D} - J^{-1} \rho_0 \dot{\psi} \geq 0 \qquad (7) $$

avec $\sigma$, $\mathbf{D}$ le tenseur de Cauchy et le taux de déformation. En utilisant (3), (4), (5) et (6) dans (7) et sous l'hypothèse d'un pseudo-potentiel de dissipation quadratique et dépendant uniquement de $\bar{\mathbf{D}}_v^0$, la loi de comportement et la loi complémentaire peuvent s'écrire :

$$ \left\{ \begin{array}{l} \sigma = \sigma_0 + \sigma_v + p \mathbf{1} \\ 2\rho_0 J^{-1} \bar{\mathbf{V}}_e \cdot \frac{\partial \psi_v}{\partial \bar{\mathbf{B}}_e} \cdot \bar{\mathbf{V}}_e = \frac{\partial \phi_v}{\partial \bar{\mathbf{D}}_v^0} \end{array} \right. \qquad (8) $$

$$ \sigma_0 = 2\rho_0 J^{-1} \left( \bar{\mathbf{B}} \cdot \frac{\partial \psi_0}{\partial \bar{\mathbf{B}}} \right)^D \qquad p = \rho_0 \frac{\partial \psi_{vol}}{\partial J} \qquad \sigma_v = 2\rho_0 J^{-1} \left( \bar{\mathbf{B}}_e \cdot \frac{\partial \psi_v}{\partial \bar{\mathbf{B}}_e} \right)^D \qquad (9) $$

On obtient à travers cette écriture la généralisation du modèle rhéologique de Zener aux grandes transformations (voir fig. 1).

## 3 Formulation Variationnelle

L'application de la contrainte de quasi-incompressibilité peut être traitée à l'aide d'une formulation en lagrangien perturbée. La solution $(\mathbf{u}; p)$ du problème d'équilibre doit annuler la forme

<footnote n="1">¹On a supposé ici que toute la partie rotation obtenue de la décomposition polaire de $\bar{\mathbf{F}}$ est due à la transformation élastique, i.e. que $\bar{\mathbf{R}} = \bar{\mathbf{R}}_e$ et ainsi $\bar{\mathbf{R}}_v = \mathbf{1}$</footnote>

<footnote n="2">²En négligeant les effets thermiques</footnote>
intégrale suivante, pour toutes les fonctions tests $\delta v$ et $\delta p$ choisies dans le même espace que $\mathbf{u}$ et $p$ :

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_0} \Pi : \nabla \delta v d\Omega - \int_{\Omega_0} \delta v \cdot \mathbf{f}_{vol} d\Omega - \int_{\delta\Omega} \delta v \cdot \mathbf{f}_{surf} dS \\[1em] \displaystyle \int_{\Omega_0} \left( -(J(\mathbf{u}) - 1) - \frac{1}{k} p \right) \delta q d\Omega \end{array} \right. \qquad (10) $$

où $\Pi$ est le premier tenseur des contraintes de Piola-Kirchoff. Pour écrire cette forme variationnelle, on a implicitement fait le choix de potentiel suivant : $\psi_{vol}(J) = \frac{k}{2}(J-1)^2$, avec $k$ le module de compressibilité. La résolution numérique de ces équations à travers un algorithme de Newton-Raphson nécessite le calcul de la contrainte et de l'opérateur tangent en effectuant un transport des équations (8) dans un formalisme mixte, que l'on peut écrire :

$$ \Pi = \underbrace{J(\boldsymbol{\sigma}_0 + p) \cdot \mathbf{F}^{-\boldsymbol{T}}}_{\Pi_h} + \underbrace{J\boldsymbol{\sigma}_v \cdot \mathbf{F}^{-\boldsymbol{T}}}_{\Pi_v} \qquad (11) $$

$$ \mathbb{T} = \frac{\partial \Pi}{\partial \mathbf{F}} = \mathbb{T}_{\mathbf{h}} + \mathbb{T}_{\mathbf{v}}, \qquad \mathbb{T}_{\mathbf{h}} = \frac{\partial J(\boldsymbol{\sigma}_0 + p) \cdot \mathbf{F}^{-\boldsymbol{T}}}{\partial \mathbf{F}}, \qquad \mathbb{T}_{\mathbf{v}} = \frac{\partial J \boldsymbol{\sigma}_v \cdot \mathbf{F}^{-\boldsymbol{T}}}{\partial \mathbf{F}} \qquad (12) $$

Le calcul de la contrainte et de l'opérateur tangent hyperélastique $\Pi_h$ et $\mathbb{T}_h$ est classique et ne présente pas de difficultés particulières (cf. [3] par exemple), si bien qu'il reste à intégrer le comportement visqueux i.e. l'évaluation de $\Pi_v$ et de $\mathbb{T}_v$ au cours d'une itération de Newton Raphson.

# 4 Intégration du comportement visco-hyperélastique

## 4.1 Loi d'écoulement du modèle de Zener

<figure><img src="image_2.png" /><figcaption>Figure 1 – Modèle de Zener</figcaption></figure>

Pour la suite et dans l'objectif de simplifier les calculs on considère un modèle rhéologique simple (cf. figure 1) avec le choix de potentiels suivants :

$$ \left\{ \begin{array}{l} \psi_h(\mathbf{B}) = C_1(I_1(\mathbf{B}_e) - 3) + C_2 \ln(\frac{I_2(\mathbf{B}_e)}{3}) \\ \varphi_v(\mathbf{D}_v^0) = \frac{\eta}{2}(\mathbf{D}_v^0 : \mathbf{D}_v^0) \end{array} \right. \qquad \psi_v(\mathbf{B}_e) = G(I_1(\mathbf{B}_e) - 3) $$

ou $C_1, C_2$ et $G$ sont les paramètres matériaux à identifier. En utilisant les potentiels précédent dans les équations (8) et en reportant ces dernières dans (5) on obtient :

$$ \dot{\mathbf{B}}_e = \mathbf{L} \cdot \mathbf{B}_e + \mathbf{B}_e \cdot \mathbf{L} - \tau \mathbf{V}_e \cdot \mathbf{B}_e^D \cdot \mathbf{V}_e - \frac{2}{3}(\mathbf{1} : \mathbf{L}) \mathbf{B}_e \qquad (13) $$

avec $\tau = \frac{4G}{J\eta}$. De plus, par définition de $\mathbf{B}_e$, on a : $\mathbf{V}_e \cdot \mathbf{B}_e^D \cdot \mathbf{V}_e = \mathbf{B}_e^D \cdot \mathbf{B}_e$. On obtient finalement la forme suivante de loi d'écoulement du modèle de Zener :

$$ \dot{\mathbf{B}}_e = \mathbf{L} \cdot \mathbf{B}_e + \mathbf{B}_e \cdot \mathbf{L} - \tau \mathbf{B}_e^D \cdot \mathbf{B}_e - \frac{2}{3}(\mathbf{1} : \mathbf{L}) \mathbf{B}_e \qquad (14) $$
## 4.2 Schéma d'intégration avec prédicteur élastique

Pour la suite on se place dans un intervalle de temps $[t_n, t_{n+1}]$, on considère que les variables d'état sont connues à l'instant $t_n$ et que l'évolution du gradient de la transformation dans l'interval est linéaire (i.e. $\dot{F}$ est constant). On définit une prédiction élastique en fixant $\bar{F}_v$ à sa dernière valeur convergée<fnref n="3" /> :

$$ \bar{F}_e^* = \bar{F} \cdot \bar{F}_{v|n}^{-1} \qquad (15) $$

on obtient également la prédiction élastique du tenseur de Cauchy-Green gauche $\bar{B}_e^*$ :

$$ \bar{B}_e^* = \bar{F} \cdot \bar{F}_{|n} \cdot \bar{B}_{e|n} \cdot \bar{F}_{|n}^{-T} \cdot \bar{F}^T \qquad (16) $$

On peut aussi définir un tenseur normalisé de contrainte déviatorique<fnref n="4" /> :

$$ \mathbf{N}^* = \frac{\boldsymbol{\sigma}_{v}^{*}}{\|\boldsymbol{\sigma}_{v}^{*}\|} = \frac{\bar{\mathbf{B}}_{e}^{*D}}{\|\bar{\mathbf{B}}_{e}^{*D}\|} \qquad (17) $$

Cette prédiction élastique est la solution exacte de l'équation (14) dans le cas ou $\tau = 0$. Dans des méthodes de prédiction correction avec un retour en exponentiel comme dans [5, 4], ce prédicteur élastique est relaxé afin de vérifier la loi d'écoulement. Cependant dans ce type de méthode, la contrainte visqueuse est temporairement figée de manière à obtenir une solution en exponentielle. La vérification de la loi d'écoulement conduit ensuite à un système non-linéaire à résoudre en chaque point de Gauss. Pour la suite, on va s'inspirer des méthodes de retour radial utilisées en plasticité [7]. On fait ainsi l'hypothèse que $\bar{B}_e$ est sous la forme :

$$ \bar{\mathbf{B}}_e = b(t)\mathbf{N}^* + s(t)\mathbf{1} \qquad (18) $$

avec $b(t)$ et $s(t)$ qui sont définis par :

$$ b(t) = \bar{\mathbf{B}}_e : \mathbf{N}^* \quad s(t) = \frac{1}{3} \bar{\mathbf{B}}_e : \mathbf{1} \qquad (19) $$

La dérivée par rapport au temps de $\bar{\mathbf{B}}_e$ s'écrit donc :

$$ \dot{\bar{\mathbf{B}}}_e = \dot{b}\mathbf{N}^* + b\dot{\mathbf{N}}^* + \dot{s}\mathbf{1} \qquad (20) $$

Les équations (18) et (20) nous permettent de réécrire la loi d'écoulement sous la forme d'un système différentiel à deux équations et deux inconnues en utilisant les projections<fnref n="5" /> $\dot{\bar{\mathbf{B}}}_e : \mathbf{1}$ et $\dot{\bar{\mathbf{B}}}_{ev} : \mathbf{N}^*$

$$ \left\{ \begin{array}{l} \dot{s} = \frac{2}{3} (\mathbf{L} : \mathbf{N}^*) b - \frac{\tau}{3} b^2 \\ \dot{b} = \left(2 (\mathbf{L} : \mathbf{N}^{*2}) - \frac{2}{3} (\mathbf{1} : \mathbf{L})\right) b + 2 (\mathbf{L} : \mathbf{N}^*) s - \tau \left( (\mathbf{N}^{*2} : \mathbf{N}^*) b^2 + sb \right) \end{array} \right. \qquad (21) $$

L'intégration de ce système différentiel scalaire non-linéaire et à coefficients non constants peut-être réalisé à l'aide d'un intégrateur numérique itératif tel qu'une méthode de Runge-Kutta ou un θ-schéma. Afin de comparer dans la suite nos résultats avec une méthode directe d'intégration, nous avons implémenté une méthode de Runge-Kutta-Feldberg qui permet d'adapter le pas de temps en fonction d'une évaluation à posteriori de l'erreur d'approximation (voir par exemple [1]).

<footnote n="3">³l'indice $\underline{|}_n$ représente une grandeur évaluée au temps $t_n$.</footnote>

<footnote n="4">⁴Par définition on a : $(\mathbf{N}^* : \mathbf{N}^*) = 1$ et $(\mathbf{N}^* : \mathbf{1}) = 0$ de plus $\tilde{\sigma}_{v}^*$ et $\bar{\mathbf{B}}_{e}^{*D}$ commutent (isotropie du comportement)</footnote>

<footnote n="5">⁵par définition de $\mathbf{N}^*$, nous avons $\mathbf{N}^* : \mathbf{1} = \mathbf{N}^* : \mathbf{N}^* = 0$</footnote>
## 4.3 Approximation en exponentielle

Le schéma numérique proposé dans cette communication se base sur une linéarisation du système différentiel (21) autour de $(b_n = b(t_n), s_n = s(t_n))$, tel que :

$$b(t) = b_n + \delta b \quad s(t) = s_n + \delta s \quad (22)$$

En reportant cette approximation dans (21) et en négligeant les termes d'ordres 2, on obtient un système différentiel linéaire à coefficients non-constants :

$$\left\{ \begin{array}{l} \dot{\delta s} \\ \dot{\delta b} \end{array} \right\} = [A(t)] \left\{ \begin{array}{l} \delta s \\ \delta b \end{array} \right\} + \{C(t)\} \quad (23)$$

Avec les conditions initiales suivantes :

$$\delta s(t_n) = \delta b(t_n) = 0 \quad (24)$$

Le système différentiel admet une solution analytique tel que :

$$\left\{ \begin{array}{l} \delta s(t_{n+1}) \\ \delta b(t_{n+1}) \end{array} \right\} = \int_{t_n}^{t_{n+1}} e^{[\Omega(t_{n+1}-s)]} \{C(s)\} ds \quad (25)$$

Ou $\Omega$ peut s'écrire sous la forme infinie d'une série de Magnus (voir [2]), tel que :

$$\Omega(t) = \int_{t_n}^{t} A(s) ds + \frac{1}{2} \int_{t_n}^{t} [A(s), \int_{t_n}^{t} A(s1) ds1] ds + \frac{1}{4} \int_{t_n}^{t} [A(s), \int_{t_n}^{t} [A(s1), \int_{t_n}^{t} A(s2) ds2] ds1] ds + \dots \quad (26)$$

Avec l'approximation d'une évolution linéaire en temps du gradient de la déformation et de L, on peut obtenir une approximation linéaire en temps de $A(t)$ si bien que la série de magnus s'intègre de manière analytique, de plus pour la suite cette dernière sera tronquée au troisième terme. La résolution du système revient donc au calcul de l'intégrale de (25). Cette dernière est réalisée avec un schéma de Gauss à 3 points. Pour l'opérateur tangent nous n'utilisons que la partie élastique de ce dernier on doit calculer ainsi :

$$T_v = \frac{\partial J \sigma_v^* \cdot F^{-T}}{\partial F}$$

# 5 Application

On s'intéresse dans cette section au cisaillement d'un bloc d'élastomère en déformations planes. Il s'agit d'un domaine de 30x30 mm, encastré sur un bord et soumis à un déplacement imposé transverse sur l'autre bord. Le maillage comprend 10 éléments par côté (quadratiques en déplacement et linéaire en pression). Le cisaillement est appliqué avec un signal sinusoïdal de fréquence 3.3 Hz et avec une amplitude de 100% de déformation en glissement. Les paramètres matériaux du modèle visco-hyperélastique de Zener sont les suivants :

$$C_1 = 0.97 \, \text{Mpa} \quad C_2 = -0.42 \, \text{Mpa} \quad G = 0.77 \, \text{Mpa} \quad \eta = 0.10 \, \text{Mpa.s}^{-1} \quad k = 1500 \, \text{Mpa}$$

La figure 2(a) montre les différentes courbes d'hystérésis, obtenues en faisant varier le pas de temps. Cette figure permet d'illustrer la stabilité du schéma d'intégration, on peut constater que pour des grands pas de temps le schéma approche la solution "fine" avec un léger décalage en terme de dissipation. Cette écart s'atténue assez rapidement en réduisant le pas de temps. Le schéma proposé est donc sensible au choix du pas de temps mais dans une mesure qui semble rester raisonnable, de plus la plupart des méthodes d'intégrations itératives sont également sensibles au pas de temps.

Sur la courbe de l'évolution de la contrainte de cisaillement en fonction du temps (figure 2(b)) on observe un très bon accord avec une méthode d'intégration itérative du système non linéaire (RK45). De plus, une comparaison du coût numérique est donnée dans le tableau 1, le schéma proposé permet de diviser le temps de calcul quasiment par 2 pour ce cas test.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2 – Test de cisaillement</figcaption></figure>

<table><caption>Tableau 1 – Comparaison du coût numérique</caption><thead><tr><td>méthode</td><td>temps CPU total (s)</td><td>Nb d'incrément</td><td>temps/incrément (s)</td></tr></thead><tbody><tr><td>relaxation radiale</td><td>85.9</td><td>320</td><td>0.268</td></tr><tr><td>Kutta Feldberg</td><td>170.8</td><td>320</td><td>0.533</td></tr></tbody></table>

## 6 Conclusion

La méthode d'intégration par relaxation radiale présente des performances tout à fait intéressantes en terme de coût de calcul. Ce gain significatif n'est pas forcément obtenu en contre partie d'une perte de précision, puisqu'au final ce schéma semble converger vers la solution recherchée. Les approximations effectuées pour obtenir ce schéma ont été appliquées au modèle de Zener. Pour autant, le même genre de démarche peut être adapté dans le cas d'un modèle plus compliqué qui présenterait une loi d'écoulement de la forme : $\dot{\bar{B}}_e = F(\bar{B}_e, L, F)$.

## Références

[1] A. Boukamel, S. Méo, S. Lejeunes. FE-Implementation of a statistical hyper-visco-plastic model. Proceedings of the 5th European Conference on Constitutive Models for Rubber. *Constitutive Models for Rubber V*, 255–261, 2008.

[2] A. Iserles., A. Marthinsen, S.P. Nørsett On the implementation of the method of magnus series for linear differential equations. *BIT Numerical Mathematics*, 39, 281–304, 1999.

[3] S. Lejeunes. Modélisation de structures lamifiées élastomère-métal à l'aide d'une méthode de réduction de modèles. Thèse de doctorat. Université d'Aix-Mareille II, 2006.

[4] B. Nedjar. Frameworks for finite strain viscoelastic-plasticity based on multiplicative decompositions. Part II : Computational aspects. *Computer Methods in Applied Mechanics and Engineering*, 191, 1563–1593, 2002.

[5] S. Reese, S. Govindjee. A theory of finite viscoelasticity and numerical aspects. *International Journal of Solid and Structures*, 35, 3455–3482, 1998.

[6] F. Sidoroff. Variables internes en viscoélasticité, 2. Milieux avec configuration intermédiaire. *Journal de Mecanique*, 14, 571–595, 1975.

[7] J.C. Simo. A framework for Finite Strain Elastoplasticity Based on Maximum Plastic Dissipation and the Multiplicative Decomposition : Part II. Computational Aspects. *Computer Methods in Applied Mechanics and Enyineering*, 61, 1–31, 1988.

[8] G. Weber, L. Arnand. Finite deformation constitutive equations and a time integration procedure for isotropic hyperelastic-viscoplastic solids. *Computer Methods in Applied Mechanics and Enyineering*, 79, 173–202, 1990.