# A variational formulation of constitutive models and updates in nonlinear finite viscoelasticity

Eduardo Fancello, Laurent Stainier, Jean-Philippe Ponthot

► To cite this version:

Eduardo Fancello, Laurent Stainier, Jean-Philippe Ponthot. A variational formulation of constitutive models and updates in nonlinear finite viscoelasticity. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812987

HAL Id: hal-01812987

https://hal.science/hal-01812987

Submitted on 12 Jun 2018

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
A variational formulation of constitutive
models and updates in nonlinear finite
viscoelasticity

Eduardo Fancello\* — Laurent Stainier\*\* — Jean-Philippe Ponthot\*\*

*Departmento de Engenharia Mecânica
Universidade Federal de Santa Catarina
Florianópolis, SC Brazil
fancello@grante.ufsc.br

\*\*Département AéroSpatiale, Mécanique et mAtériaux
Université de Liège, LTAS-MCT, 1 chemin des chevreuils
B-4000 Liège, Belgique
{l.stainier, jp.ponthot}@ulg.ac.be

ABSTRACT. In [ORT 99], Ortiz and Stainier proposed a general variational approach for elasto-
plastic models in finite deformations regime. This work can be extended to finite viscoelastic
models, as shown in [STA 03]. The striking feature of this approach is its variational charac-
teristic which provides an appropriate mathematical structure and allows the use of theoretical
and numerical facilities like, for instance, error estimation studies or mathematical program-
ming algorithms. The present paper follows the same path and focuses on a general varia-
tional approach for finite viscoelastic models. In addition, a specific group of them (generalized
Kelvin-Maxwell models) is analyzed with detail, together with numerical implementation. Fi-
nally, numerical simulations illustrate the performance of the present approach.

RÉSUMÉ. Dans [ORT 99], Ortiz et Stainier ont proposé une approche variationnelle générale de l'élastoplasticité en grandes déformations. Ces travaux peuvent être étendus au cadre de la viscoélasticité en grandes déformations, comme illustré dans [STA 03]. L'atout principal de cette approche est son caractère variationnel, qui fournit un cadre mathématique particulièrement intéressant. Cet article s'intéresse de plus près au développement d'une approche variationnelle générale des modèles de viscoélasticité en grandes déformations. En particulier, un sous-ensemble de ces modèles (matériaux de Kelvin-Maxwell généralisés) sont étudiés en détails, y compris en regard de leur mise en œuvre numérique. Finalement, des simulations numériques illustrent la performance de l'approche proposée.

KEYWORDS: viscoelasticity, variational formulations, nonlinear constitutive models

MOTS-CLÉS : viscoélasticité, formulations variationnelles, modèles constitutifs non-linéaires
# 1. Introduction

In [ORT 99], Ortiz and Stainier proposed a general variational approach for elasto-plastic models in finite deformations regime. This work can be extended to finite viscoelastic models, as shown in [STA 03]. The striking feature of this approach is its variational characteristic which provides a mathematical structure opening the way to many theoretical and numerical tools and methods. The present paper follows the same path and focuses on a general variational approach for finite viscoelastic models.

# 2. Variational form of constitutive equations

Using conventional notation, let us call $\mathbf{F} = \nabla_0 \mathbf{x}$ the gradient of deformations, and $\mathbf{C} = \mathbf{F}^T \mathbf{F}$ the Cauchy strain tensor, respectively. These values may be decomposed in volumetric and isochoric parts. The isochoric tensors are defined as follows:

$$ \hat{\mathbf{F}} = \frac{1}{J^{1/3}} \mathbf{F}, \quad J = \det(\mathbf{F}), \quad \hat{\mathbf{C}} = \hat{\mathbf{F}}^T \hat{\mathbf{F}} = \frac{1}{J^{2/3}} \mathbf{F}^T \mathbf{F}, \quad (1) $$

We will work in the framework of irreversible thermodynamics, with internal variables. Thus, we define a general set $\mathcal{E} = \{\mathbf{F}, \mathbf{F}^i, \mathbf{Q}\}$ of external and internal variables, where $\mathbf{F}^i$ its inelastic part of the (total) deformation, and $\mathbf{Q}$ contains all the remaining internal variables of the model. In addition, a multiplicative decomposition $\mathbf{F} = \mathbf{F}^e \mathbf{F}^i$ of the gradient of deformations is considered. We assume the existence of a free energy potential $W(\mathcal{E})$ and a dissipative potential $\phi(\hat{\mathbf{F}}; \mathcal{E})$, such that the Piola-Kirchhoff stress tensor, comprised of an equilibrium (elastic) and a dissipative (viscous) components, is derived as follows:

$$ \mathbf{P} = \frac{\partial W}{\partial \hat{\mathbf{F}}} (\mathcal{E}) + \frac{\partial \phi}{\partial \hat{\mathbf{F}}} (\hat{\mathbf{F}}; \mathcal{E}). \quad (2) $$

In addition, another dissipative potential $\psi(\hat{\mathbf{F}}^i, \dot{\mathbf{Q}}; \mathcal{E})$ is included to characterize the inelastic behavior related to the inelastic tensor $\mathbf{F}^i$, such that

$$ \mathbf{T} = -\frac{\partial W}{\partial \hat{\mathbf{F}}^i}(\mathcal{E}) = \frac{\partial \psi}{\partial \hat{\mathbf{F}}^i}(\hat{\mathbf{F}}^i, \dot{\mathbf{Q}}; \mathcal{E}), \quad \mathbf{A} = -\frac{\partial W}{\partial \mathbf{Q}}(\mathcal{E}) = \frac{\partial \psi}{\partial \dot{\mathbf{Q}}^i}(\hat{\mathbf{F}}^i, \dot{\mathbf{Q}}; \mathcal{E}). \quad (3) $$

It was shown in [ORT 99] and [STA 03] that an incremental version of the above equations, constituting an incremental update method for the material state, can be obtained from the following incremental potential:

$$ W(\mathbf{F}_{n+1}; \mathcal{E}_n) = \Delta t \phi(\hat{\mathbf{F}}, \mathcal{E}_n) + \min_{\substack{\mathbf{F}^i \\ \mathbf{Q}_{n+1}}} \left\{ W(\mathcal{E}_{n+1}) - W(\mathcal{E}_n) + \Delta t \psi(\hat{\mathbf{F}}^i, \dot{\mathbf{Q}}; \mathcal{E}_n) \right\}, \quad (4) $$

where $\hat{\mathbf{F}}(\mathbf{F}_{n+1}, \mathcal{E}_n)$, $\hat{\mathbf{F}}^i(\mathbf{F}_{n+1}^i, \mathcal{E}_n)$ and $\dot{\mathbf{Q}}(\mathbf{Q}_{n+1}, \mathcal{E}_n)$ are suitable incremental approximations of the rate variables $\hat{\mathbf{F}}, \hat{\mathbf{F}}^i$ and $\dot{\mathbf{Q}}$ respectively.
### 3. A group of visco-hyperelastic models

#### 3.1. General form

A quite general group of viscoelastic materials can be modelled within the present variational framework. Due to the possibility of obtaining analytical or semi-analytical expression for the constitutive updates, only isotropic models will be considered now. However, no theoretical constraints to include more general behaviors are found. The

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Generalized Kelvin-Maxwell model.</figcaption></figure>

rheological mechanism shown in Figure 1(a) is taken as a basis to include different potentials expressions in (4). The model is based on the following assumptions:

- The elastic part of the Kelvin branch is split in isochoric and volumetric energies. The isochoric part is a isotropic function of $\hat{\mathbf{C}} = \hat{\mathbf{F}}^T \hat{\mathbf{F}}$:

$$ \varphi(\hat{\mathbf{C}}) = \varphi(c_1, c_2, c_3), \qquad (5) $$

where $c_j$ are the eigenvalues of $\hat{\mathbf{C}}$. The volumetric part may be defined using the usual expression $U(J) = \frac{K}{2} [\ln J]^2$. The viscous part of the Kelvin branch is an isotropic function of the symmetric part of the rate of deformation:

$$ \phi(\mathbf{D}) = \phi(d_1, d_2, d_3) \quad \text{with} \quad \mathbf{D} = \text{dev} \left( \text{sym} \left( \dot{\mathbf{F}} \mathbf{F}^{-1} \right) \right), \qquad (6) $$

where $d_j$ are the eigenvalues of $\mathbf{D}$.

- The Maxwell branch, connected in parallel, is based on a multiplicative split of strains in an elastic and an isochoric inelastic (viscous) part:

$$ \hat{\mathbf{F}} = \hat{\mathbf{F}}^e \mathbf{F}^v \implies \hat{\mathbf{F}}^e = \hat{\mathbf{F}} \mathbf{F}^{v-1}, \quad \det \mathbf{F}^v = 1. \qquad (7) $$

A flow rule for the internal variable $\mathbf{F}^v$ can be written as:

$$ \dot{\mathbf{F}}^v = \mathbf{D}^v \mathbf{F}^v = (d_j^v \mathbf{M}_j^v) \mathbf{F}^v, \qquad (8) $$

in which the spectral decomposition of $\mathbf{D}^v = \text{sym}(\dot{\mathbf{F}}^v \mathbf{F}^{v-1})$ in eigenvalues $d_j^v$ and eigenprojections $\mathbf{M}_j^v$, $j = 1, 2, 3$, was used. The scalars $d_j$ are chosen to be the internal variables contained in the set $\dot{\mathbf{Q}} = \{d_1, d_2, d_3\}$. In this case, it is important to
note that (8) is a constraint relating the internal variables $\mathbf{F}^v$ and $\mathbf{Q}$. The elastic and viscous potentials associated to this branch are assumed to be isotropic functions of the elastic deformation and viscous stretching, and thus depend on their eigenvalues:

$$ \varphi^e(\hat{\mathbf{C}}^e) = \varphi^e(c_1^e, c_2^e, c_3^e) \quad \text{and} \quad \psi(\mathbf{D}^v) = \psi(d_1^v, d_2^v, d_3^v), \qquad (9) $$

where $c_j^e$ are the eigenvalues of $\hat{\mathbf{C}}^e$.

– Viscous deformations are incrementally updated by exponential mappings:

$$ \Delta \hat{\mathbf{F}} = \hat{\mathbf{F}}_{n+1} \hat{\mathbf{F}}_n^{-1} = \exp[\Delta t \mathbf{D}] \quad \Rightarrow \mathbf{D} = \frac{\Delta q_j}{\Delta t} \mathbf{M}_j = \frac{1}{2\Delta t} \ln(\Delta \hat{\mathbf{C}}). \qquad (10) $$

$$ \Delta \mathbf{F}^v = \mathbf{F}_{n+1}^v \mathbf{F}_n^{v-1} = \exp[\Delta t \mathbf{D}^v] \quad \Rightarrow \mathbf{D}^v = \frac{\Delta q_j^v}{\Delta t} \mathbf{M}_j^v = \frac{1}{2\Delta t} \ln(\Delta \mathbf{C}^v). \qquad (11) $$

Expressions (10) and (11) show that $\mathbf{D}$ and $\mathbf{D}^v$ are approximated by incremental expressions of $\Delta \hat{\mathbf{C}}$ and $\Delta \mathbf{C}^v$ respectively. The exponential mapping has the particular convenient property of providing a isochoric tensor for any traceless argument [WEB 90, MIE 96].

Taking into account (10) and (11), the minimizing variables $\mathbf{Q}_{n+1}, \mathbf{F}_{n+1}^v$ in (4) are replaced by the new incremental variables $\Delta q_j^v, \mathbf{M}_j^v$:

$$ W(\mathbf{F}_{n+1}; \mathcal{E}_n) = W(\mathbf{C}_{n+1}; \mathcal{E}_n) = \Delta\varphi(\hat{\mathbf{C}}_{n+1}) + \Delta t \phi \left(\frac{\Delta q_j^v}{\Delta t}\right) + \Delta U(\theta_{n+1}) \\ + \min_{\mathbf{M}_j^v, \Delta q_j^v} \left\{ \Delta\varphi^e(\hat{\mathbf{C}}_{n+1}^e) + \Delta t \psi \left(\frac{\Delta q_j^v}{\Delta t}\right) \right\}, \qquad (12) $$

such that

$$ \Delta q_j^v \in K_Q = \{p_j \in \mathbb{R} : p_1 + p_2 + p_3 = 0\}, \qquad (13) $$

$$ \mathbf{M}_j^v \in K_M = \{ \mathbf{N}_j \in \text{Sym} : \mathbf{N}_j \cdot \mathbf{N}_j = 1, \mathbf{N}_i \cdot \mathbf{N}_j = 0, i \neq j \}. \qquad (14) $$

The set $K_Q$ enforces the traceless form of $\mathbf{D}^v$, while the set $K_M$ accounts for usual properties of eigenprojections. Moreover, it is easy to verify that both sets are convex on their respective variables. Given isotropic expressions for energy functions, the minimization in (12) can be performed analytically. A simple extension to this model can be obtained by considering a set of P Maxwell branches, as seen in Figure 1(a).

## 3.2. Hencky and Ogden models

Hencky models are based on quadratic forms of logarithmic strain tensors:

$$ \varphi = \mu \sum_{j=1}^{3} (\epsilon_j)^2, \quad \phi = \eta \sum_{j=1}^{3} (d_j)^2, \qquad (15) $$

$$ \varphi^e = \mu^e \sum_{j=1}^{3} (\epsilon_j^e)^2, \quad \psi = \eta^v \sum_{j=1}^{3} (d_j^v)^2. \qquad (16) $$
In this case, it is particularly convenient to obtain simple uncoupled linear expressions for the minimizing argument $\Delta q_j^v$. In spite of the facility offered by Hencky models in terms of analytical treatment, it is well known that these type of hyperelastic potentials do not fit well the behavior of rubber-like materials. For that case, a more adequate choice may be the Ogden model which has also the capability of generalizing other models like neo-Hookean and Mooney-Rivlin. Ogden models are based on the following potentials:

$$ \varphi = \sum_{j=1}^{3} \sum_{p=1}^{N} \frac{\mu_p}{\alpha_p} ([\exp(\epsilon_j)]^{\alpha_p} - 1), \quad \phi = \sum_{j=1}^{3} \sum_{p=1}^{N} \frac{\eta_p}{\alpha_p} ([\exp(d_j)]^{\alpha_p} - 1), \quad (17) $$

$$ \varphi^e = \sum_{j=1}^{3} \sum_{p=1}^{N} \frac{\mu_p^e}{\alpha_p} ([\exp(\epsilon_j^e)]^{\alpha_p} - 1), \quad \psi = \sum_{j=1}^{3} \sum_{p=1}^{N} \frac{\eta_p^v}{\alpha_p} ([\exp(d_j^v)]^{\alpha_p} - 1). \quad (18) $$

## 4. Numerical example

This example presents a pure shear tests of a single 3D element (Figure 1(b)). Material parameters and load characteristics were taken from an equivalent example in [REE 98], in order to perform some useful comparisons. Thus, the rheological model chosen for this example corresponds to that of Figure 1(b). The lateral displacement $u_x$ follows a sinusoidal law $u_x = U \sin wt$, where $w = 0.3 \ s^{-1}$. The material was assumed to be almost incompressible through the choice of a high value for the bulk modulus $K$. Two different models for $\varphi$ were used: Ogden model and Hencky model. In the case of Ogden, we used the following six-parameter fitting: $\mu_1 = 20$, $\mu_2 = -7$, $\mu_3 = 1.5$, $\alpha_1 = 1.8$, $\alpha_2 = -2$, $\alpha_3 = 7$. For the Hencky model, the value $\mu = \sum_i \frac{1}{2}\mu_i\alpha_i = 30.25$ was used, which is the consistent equivalent shear modulus for small deformations. The Maxwell branch uses Hencky model for both potentials with $\mu^e = 77.77$ and viscous coefficient $\eta^v$ such that $\tau = \frac{\eta^v}{\mu^e} = 17.5$.

The time evolution of Cauchy stresses $\sigma_{xy}$ as a function of shear strain $C_{xy}$, for different shear amplitudes, is shown in Figure 2. In the case of small strains both models (Ogden or Hencky main spring) give identical results, and match quite well equivalent results in [REE 98]. As expected, the behavior of the main spring is determinant on the behavior of the whole system for deformations higher than unity. Comparing the results of the Ogden-based model with those of [REE 98] it is possible to see a close correlation of maximum values of stress for all four cases. However, hysteresis loops clearly look “thinner” as the deformation grows along the cycle. This behavior is in agreement with the fact that the Hencky model used in the Maxwell branch provides a contribution in stress much more lower than a corresponding Ogden model for high deformations. In [REE 98], Ogden model is used for both, main and Maxwell springs.

## 5. Conclusion

We have proposed a general variational formulation of nonlinear finite viscoelasticity models. We focused in particular on generalized Kelvin-Maxwell models. We
<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Stress-strain curves in cyclic shear test*</figcaption></figure>

compared the capacity of Hencky- and Ogden-type models to reproduce observed non-linear viscous behaviour in shear tests. It appears that Ogden models perform better, without really involving additional complexity to the numerical implementation.

**6. References**

[MIE 96] MIEHE C., "Exponential map algorithm for stress updates in anisotropic multiplicative elastoplasticity for single crystals", *International Journal for Numerical Methods in Engineering*, vol. 39, 1996, p. 3367-3390.

[ORT 99] ORTIZ M., STAINIER L., "The variational formulation of viscoplastic constitutive updates", *Computer Methods in Applied Mechanics and Engineering*, vol. 171, num. 3-4, 1999, p. 419-444, Elsevier.

[REE 98] REESE S., GOVINDJEE S., "A theory for finite viscoelasticity and numerical aspects", *International Journal of Solids and Structures*, vol. 35, 1998, p. 3455-3482.

[STA 03] STAINIER L., "Une formulation variationelle des algorithmes de calcul des contraintes pour les modèles élastoviscoplastiques et viscoélastiques en grandes transformations", POTIER-FERRY M., BONNET M., BIGNONNET A., Eds., *Géme Colloque National en Calcul des Structures*, vol. 2, Giens, 2003, p. 477-484.

[WEB 90] WEBER G., ANAND L., "Finite deformations constitutive equations and a time integration procedure for isotropic hyperelastic-viscoplastic solids", *Computer Methods in Applied Mechanics and Engineering*, vol. 79, 1990, p. 173-202.