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

HAL
open science

# Convergence analysis of XFEM methods for a nonlinear viscoelasticity constitutive law

Cédric Pozzolini, Thouraya Baranger, Yves Renard

► To cite this version:

Cédric Pozzolini, Thouraya Baranger, Yves Renard. Convergence analysis of XFEM methods for a nonlinear viscoelasticity constitutive law. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722110

HAL Id: hal-01722110

https://hal.science/hal-01722110

Submitted on 2 Mar 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.

Public Domain
# Convergence analysis of XFEM methods for a nonlinear viscoelasticity constitutive law

Cédric POZZOLINI 1,2, Thouraya N. BARANGER 3, Yves RENARD 4

**Résumé** — The purpose of this note is to describe and compare some extended finite element approximations and their convergence for crack rubber-like materials. The analysis is based on the finite plane deformation visco-hyperelasticity theory for an incompressible Mooney-Rivlin material. Numerical examples show that for well chosen enrichment strategies, the finite element converges with the expected convergence order. Finally we discuss the numerical results we obtain.

**Mots clés** — extended finite element method, non-linear visco-elasticity, fracture mechanics, incompressible, Cut-off.

## 1 Introduction

There exist several numerical methods to overcome problems dealing with linear fracture mechanics, but very few tools exist to solve large deformations fracture mechanics problems. This work proposes a numerical analysis of viscoelastic strain and displacement fields in the vertex of a crack-tip in large strain fracture mechanics using eXtended Finite Element Method [1]. A major issue in fracture mechanics is the prediction of crack kinking, and this is a open question in this framework. Let us recall that XFEM is a way to introduce the discontinuity across the crack and the asymptotic displacement into the finite element space. It has been initially developed for plane elasticity problems and is now the subject of a large literature. Concerning the specific XFEM enrichment, a "jump function" (or Heaviside function) is used in order to represent the discontinuity due to the crack. Following the ideas already presented in [5], we propose two strategies for the crack tip enrichment. In both of them, an enrichment area of fixed size is defined, centered on the crack tip. In the first strategy, each node contained in the enrichment area has all its degrees of freedom enriched with the crack tip asymptotic bending. For the second strategy, the asymptotic planar displacements are introduced in a global way, using a cut-off function. Along this paper, it will be underlined that our XFEM formulation presents a sense of optimality, both about accuracy and computational cost, since they are nearly equal to those of a classical finite element method on a regular non-cracked problem. Analysis of convergence are performed for nonlinear constitutive law : short and long memory visco-hyperelastic Mooney-Rivlin model.

## 2 Short and long memory visco-hyperelastic laws

We consider a finite plane deformation problem in an incompressible isotropic visco-hyperelastic media, $\Omega$. We assume that $\Omega$ is a two-dimensional-lipschitz-crack-bounded-domain, and we assume the following decomposition of the boundary :

$$ \partial\Omega := \Gamma_N \cup \Gamma_D \cup \Gamma_C^+ \cup \Gamma_C^- $$

where $meas(\Gamma_D) \neq 0$, and $\Gamma_C^+, \Gamma_C^-$ are not empty (Fig. 1).

The crack $\Gamma_C$ is represented by its lips $\Gamma_C^+, \Gamma_C^-$, and these curves are considered free from any soli-citations. The crack tip is given by the intersection $\Gamma_C^+ \cap \Gamma_C^-$. We assume that $\partial\Omega \setminus \Gamma_C$ is partitioned into two sets : $\Gamma_N$ where a Neumann force $g : [0, T] \to \mathbb{L}^2(\Gamma_N)$ is applied, and $\Gamma_D$ where a Dirichlet condition

<footnote n="1">¹ Université Paris Est, IFSTTAR-GRETTIA, Champs sur Marne, France</footnote>

<footnote n="2">² INSA-Lyon, LaMCoS UMR5259, F-69621 Villeurbanne, France, cedric.pozzolini@insa-lyon.fr,</footnote>

<footnote n="3">³ Université de Lyon, CNRS, Université Lyon1, LaMCoS UMR5259, F-69621 Villeurbanne, France, thouraya.baranger@univ-lyon1.fr,</footnote>

<footnote n="4">⁴ INSA-Lyon, CNRS, ICJ UMR5208, LaMCoS UMR5259, F-69621 Villeurbanne, France, yves.renard@insa-lyon.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Example of initial crack domain.</figcaption></figure>

$u = 0$ is prescribed. Let $f: [0, T] \to \mathbb{L}^2(\Omega)$ be the body force applied on $\Omega$. In the following $x \in \Omega$ is the space variable, $t \in [0, T]$ is the time variable, and $d(\star)/dt = (\star)$ refers the time derivative of $(\star)$. The geometry transformation subjected to the incompressibility which maps $\Omega$ on to a domain constraint is given by :

$$ \det(\mathbf{F}) = J = 1. \qquad (1) $$

The classical 3D strain energy of the Mooney-Rivlin hyperelastic form :

$$ \mathbf{W}^{\text{iso}} = C_1(\bar{I}_1 - 3) + C_2(\bar{I}_2 - 3); \qquad (2) $$

is often employed in the description of the static behavior of isotropic rubber-like materials. Using the strain invariants

$$ \bar{I}_1 = \mathrm{tr}(\bar{\mathbf{C}}), \quad \bar{I}_2 = \frac{1}{2}(\mathrm{tr}(\bar{\mathbf{C}}) - \mathrm{tr}((\bar{\mathbf{C}})^2)), $$

with $u$ the displacement, the right Cauchy-Green tensor

$$ \mathbf{C} = \mathbf{F}^t \mathbf{F} = J^{2/3} \bar{\mathbf{C}}, $$

and the material deformation gradient

$$ \mathbf{F} = \mathbf{I} + \nabla u = J^{1/3} \bar{\mathbf{F}}. $$

The asymptotic functions obtained for a notch vertex in the case of hyperelastic Mooney-Rivlin model can be restricted to the crack-tip case, and are given by [2] :

$$ y_1(r, \theta) = r\sin^2(\theta/2), \quad y_2(r, \theta) = \sqrt{r\sin(\theta/2)}. \qquad (3) $$

Where $y_1, y_2$ are the components of the cartesian position in the euler configuration, $r, \theta$ are the radial coordinates of material point. More precisely for $\beta$-short memory Kelving model, we assume that the second Piola-Kirchhoff stress tensor $\mathbf{S} = \mathbf{S}^{\text{vol}} + \mathbf{S}^{\text{iso}}$ is given by Mooney-Rivlin potential $\mathbf{W}^{\text{iso}}$ in following form

$$ \mathbf{S}^{\text{iso}} = 2 \frac{\partial \mathbf{W}^{\text{iso}}(\bar{\mathbf{C}})}{\partial \bar{\mathbf{C}}} + \beta \frac{\partial \dot{\mathbf{W}}^{\text{iso}}(\bar{\mathbf{C}})}{\partial \bar{\mathbf{C}}}. \qquad (4) $$

For long memory generalised Maxwell model, we assume that $\mathbf{S}$ is obtained by Prony series :

$$ \mathbf{S} = \mathbf{S}^{\text{iso}} + \sum_{\alpha=1}^{m} \beta_{\alpha} \mathbf{Q}^{\alpha} \qquad (5) $$

where $\mathbf{S}^{\text{iso}} := 2 \frac{\partial \mathbf{W}^{\text{iso}}(\bar{\mathbf{C}})}{\partial \bar{\mathbf{C}}}$, $\beta_{\alpha} > 0$ is the viscosity parameter, the relaxation time is $\tau_{\alpha} > 0$, and $\mathbf{Q}^{\alpha}$ verifies the evolution equation :

$$ \dot{\mathbf{Q}}^{\alpha} + \frac{\mathbf{Q}^{\alpha}}{\tau_{\alpha}} = \beta_{\alpha} \mathbf{S}^{\text{iso}}. \qquad (6) $$

Note that $\mathbf{S}^{\text{vol}} = 0$ because the incompressibility assumption (1).
The variational formulation, corresponding to the quasi-static equilibrium in the plane domain $\Omega$, reads like this :

$$ egin{cases} 
\text{Find } (u, p) : [0, T] \to (\mathbb{V} \times \mathbb{Q}) \text{ such that :} \\ 
\int_{\Omega} \mathbf{FS} : \nabla v d\Omega - \int_{\Omega} p \operatorname{div} v d\Omega = \int_{\Gamma_N} g v d\Gamma + \int_{\Omega} f v d\Omega, \forall v \in \mathbb{V}, \\ 
\int_{\Omega} q \operatorname{div} u d\Omega = 0, \forall q \in \mathbb{Q}, 
\end{cases} 
\quad (7) $$

where $\mathbb{V}$ is the vector space of the kinematically admissible displacement fields :

$$ \mathbb{V} = \{u \in \mathbb{H}^1(\Omega) \mid u = 0 \text{ on } \Gamma_2\}, \quad \mathbb{Q} = \mathbb{L}^2(\Omega). \quad (8) $$

Discretization of the elasticity problem follows the usual steps. Let us consider an affine mesh of the non cracked domain $\bar{\Omega}$. We approximate $(u, p)$ by $(u^h, p^h) : [0, T] \to (\mathbb{V}^h \times \mathbb{Q}^h)$. The subspaces $\mathbb{V}^h$ and $\mathbb{Q}^h$ are finite dimensional spaces that will be defined later. The existence of a stable finite element approximate solution $(u^h, p^h)$ depends on choosing a pair of spaces $\mathbb{V}^h$ and $\mathbb{Q}^h$ such that the an "inf-sup-type" condition holds. The satisfaction of this condition for a couple $(\mathbb{V}^h \times \mathbb{Q}^h)$ is very difficult to prove in practical situations. To obtain solutions of (7) for long memory viscohyperelastic model, a incremental process of Newton's type is applied, which solves a sequence of linearized problems. Then the total discretized problem associated to (7) with short memory viscohyperelastic model (4) is :

$$ egin{cases} 
\text{for } (u_1^h, p_1^h) \text{ given, } \forall n \ge 1, \text{ find } (u_{n+1}^h, p_{n+1}^h) \in \mathbb{V}^h \times \mathbb{Q}^h, \text{ such that} \\ 
\int_{\Omega} (\mathbf{I} + \nabla u_{n+1}^h) ((\Delta t + \beta)\mathbf{S}^{\text{isoh}}_{n+1} - \beta\mathbf{S}^{\text{isoh}}_n) : \nabla v^h d\Omega - \int_{\Omega} p_{n+1}^h v^h d\Omega = \\ 
\qquad \Delta t \left( \int_{\Gamma_N} g_{n+1}^h v^h d\Gamma + \int_{\Omega} f_{n+1}^h v^h d\Omega \right), \quad \forall v^h \in \mathbb{V}^h, \\ 
\int_{\Omega} q^h \operatorname{div} u_{n+1}^h d\Omega = 0, \forall q^h \in \mathbb{Q}^h, 
\end{cases} 
\quad (9) $$

Concerning the long memory viscohyperelastic (5) model the total discretized problem to solve is :

$$ egin{cases} 
\text{for } (u_1^h, p_1^h) \text{ given, } \forall n \ge 1, \text{ find } (u_{n+1}^h, p_{n+1}^h) \in \mathbb{V}^h \times \mathbb{Q}^h \\ 
\sum_{\alpha=1}^{m} \int_{\Omega} (\mathbf{I} + \nabla u_{n+1}^h) \left( (1 + \beta_{\alpha} \exp(-\frac{\Delta t}{\tau_{\alpha}})) \mathbf{S}^{\text{isoh}}_{n+1} : \nabla v^h \right) d\Omega \\ 
- \sum_{\alpha=1}^{m} \int_{\Omega} \left( (\mathbf{I} + \nabla u_{n+1}^h) \mathbf{H}_{n}^{h,\alpha} : \nabla v^h \right) d\Omega \\ 
- \int_{\Omega} p_{n+1}^h \operatorname{div} v^h d\Omega = \int_{\Gamma_N} g_{n+1}^h v^h d\Gamma + \int_{\Omega} f_{n+1}^h v^h d\Omega, \quad \forall v^h \in \mathbb{V}^h, \\ 
\int_{\Omega} q^h \operatorname{div} u_{n+1}^h d\Omega = 0, \forall q^h \in \mathbb{Q}^h, 
\end{cases} 
\quad (10) $$

with $\mathbf{P}_n = (\mathbf{I} - \frac{1}{3}\bar{\mathbf{C}}_n^{h,-1} \otimes \bar{\mathbf{C}}_n^h)$, $\mathbf{H}_n^{h,\alpha} = \exp(-\frac{\Delta t}{\tau_{\alpha}})\mathbf{P}_n : (\exp(-\frac{\Delta t}{\tau_{\alpha}})\mathbf{Q}_n^{h,\alpha} - \beta_{\alpha}\mathbf{S}^{\text{isoh}}_n)$ and

$$ \mathbf{Q}_{n+1}^{h,\alpha} = \beta_{\alpha} \exp(-\frac{\Delta t}{\tau_{\alpha}}) \mathbf{S}^{\text{isoh}}_{n+1} + \mathbf{H}_n^{h,\alpha}. $$

### 3 XFEM enrichments

Classically the displacement discontinuity across the crack is represented using a Heaviside-like function, which is multiplied by the finite element shape functions. An enrichment area of fixed size is defined and some nonsmooth functions are added inside this area.

We perform a convergence analysis for (9) and (10) of two variant of eXtended finite element methods. The first one is a standard X-FEM with a fixed enriched zone [1], the second one is an X-FEM
<figure><img src="image_3.png" /><figcaption>Fig. 2 – Example of cut-off exponential function.</figcaption></figure>

with a Cut-off function (Fig. 2) to ensure a smooth transition between enriched zone and non enriched zone [4]. Unless the X-FEM with a fixed enriched zone increase the computational cost, the X-FEM with cut-off function method gives optimal convergence results without increasing significantly the computational cost and without degrading the condition number of the tangent matrix. The finite element method is defined independently from the crack on a mesh of the non-cracked domain $\Omega$. At the nodes whose corresponding shape function support is cut by the crack, an enrichment function of Heaviside type is considered :

$$H(x) = \begin{cases} +1, & \text{if } (x-x^*) \cdot n \le 0, \\ -1, & \text{elsewhere,} \end{cases} \qquad (11)$$

where $x^*$ denotes the crack tip and $n$ is a given normal to the crack. Moreover, the nodes of the triangle containing the crack tip are enriched with the following singular functions given in polar coordinates :

$$\{F_j\}_{1 \le j \le 2 \times m} = \{\sqrt{r} \sin(k\theta/2), \sqrt{r} \cos(k\theta/2)\}_{1 \le k \le m} \qquad (12)$$

In the proposed variant of XFEM, we intend to enrich a whole area around the crack tip by using a cut-off function denoted $\chi$. Let us recall that the discretized displacement using $\chi$ can be decomposed as :

$$u^h = \sum_i a_i \phi_i + \sum_i b_i H \phi_i + \sum_j c_j F_j \chi \qquad (13)$$

Where $\phi_i$ is the classic FE basis, $H$ Heaviside step up function (11), $F_j$ functions taking account singularities (12). Since the mesh is independent of the crack geometry, the Heaviside function represents the discontinuity of the displacement field along the crack and the singular functions allow to take into account the asymptotic displacement at the crack tip. Such a method enables to discretize the domain without explicitly meshing the crack surfaces, and hence the crack propagation simulations can be done without remeshing. For the model problem, we consider a Lagrange finite element method of second order defined on a regular triangulation of the non-cracked domain $\Omega$. Then $V^h$ is defined by piecewise P2 basis functions, denoted $\phi_i$ (P2 is the set of second degree polynomials).

We use a mixed finite element method to prevent the locking of the finite element approximation in the incompressible limit (1). To interpolate the pressure with $Q^h$, we use a Lagrange finite element method of first order. The main focus in this contribution is to compare various approaches for the numerical evaluation of stress field, and displacement. These approaches have been proposed in the context of both classical and enriched finite elements for linear material [6]. However, their convergence and the robustness have not yet been extended to hyper-viscoelastic material, which are the goals of this study. The convergence results are obtained with tests in the following cracked specimen (Fig. 3).

## 3.1 Singular functions relatively to the crack tip

On supposing that the global notch problem of the Mooney Rivlin hyperelastostatic admits a solution, the asymptotic analysis of displacement field near the crack tip give a $r$-singularity of order 1/2. Then for model (4) and (5) the displacement and the pressure field near the crack tip are of order 1/2.
## 3.2 Enrichments strategy

We perform a convergence analysis for two variant of extended finite element methods. The first one is a standard X-FEM with a fixed enriched zone [1], [3], adding degrees of freedom on each node, for plane elasticity problems. The second one is an X-FEM with a Cut-off function to ensure a smooth transition between enriched zone and non enriched zone [4]. This second strategy is inspired by the so called "XFEM dof gathering with pointwise matching", introduced in [5]. The idea of such method is to reduce the number of degrees of freedom that are introduced in the finite element basis. For this, the support of the nonsmooth added functions is the whole enrichment area and they are not multiplied by the finite element basis functions. We introduce the boundary of the enrichment area, say $\Gamma_e$. It cuts into two sub-domains : the enrichment area, say $\Omega_1$, and the rest of $\Omega$, say $\Omega_2$. The unknowns defined on each domain $\Omega_i$ are denoted by $u_i^h$ and their expressions read

$$ \left\{ \begin{aligned} u_1^h &= \sum_{i \in N_1} a_i \phi_i + \sum_{i \in J_1} b_i H \phi_i + \sum_j c_j F_j \chi \\ u_2^h &= \sum_{i \in N_2} a_i \phi_i + \sum_{i \in J_2} b_i H \phi_i \end{aligned} \right. \qquad (14) $$

where $N_i$ are the sets of degrees of freedom that are located in $\Omega_i$, ($N_1 \cap N_2$ is not empty : it corresponds to the set of nodes that are on the boundary $\Gamma_e$). In a same way, $J_i$ are the sets of degrees of freedom of $\Gamma_e$ that are located in $\Omega_i$ ($J_1 \cap J_2$ is not empty for the same reason). Naturally, a matching condition is needed at the interface between the enrichment area and the rest of the domain, in order to insure the continuity of the function and its derivatives.

## 4 Numerical experiments : Test problems

The numerical experiments presented in this section were performed with the open-source finite element library Getfem++ developed by our team. We test the previous methods on two benchmark problems of a rubber like material with a straight through crack. The elastic structure is occupying a bounded rectangular domain $\Omega = [-7, +7] \times [-25, +25]$, and the material is assumed to be isotropic. In the experiments, the Mooney Rivlin constants are taken as $C_{10} = 0.57559$, and $C_{01} = 0.07585$. The enrichment area diameter is 0.30. The specimen is subjected to a sudden strain that is kept constant over the duration of the test, and the stress decrease is recorded over prolonged period. The initial stress is due to the elastic response of the material. Then, the stress relaxes over time due to the viscous effects in the material (Fig. 4).

<table><thead><tr><td>viscosity</td><td>relaxation time</td></tr></thead><tbody><tr><td>β<sub>1</sub> = 0.112</td><td>τ<sub>1</sub> = 0.2</td></tr><tr><td>β<sub>2</sub> = 0.692</td><td>τ<sub>2</sub> = 2</td></tr><tr><td>β<sub>3</sub> = 0.123</td><td>τ<sub>3</sub> = 20</td></tr><tr><td>β<sub>4</sub> = 0.0832</td><td>τ<sub>4</sub> = 200</td></tr><tr><td>β<sub>5</sub> = 0.0416</td><td>τ<sub>5</sub> = 1779</td></tr></tbody></table>

Parameters of the long memory viscohyperelastic model
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Displacement field for the opening mode using a cut-off function.</figcaption></figure>

Références

[1] N. Moes, J. Dolbow, T. Belytschko, *A finite element method for crack growth without remeshing*, International Journal for Numerical Methods in Engineering, Vol. 46, 1999.

[2] M. Arfaoui, K. Mansouri, A. Rezgui, *An asymptotic finite plane deformation analysis of the elastostatic fields at a notch vertex of an incompressible hyperelastic material*, Comptes Rendus Mecanique, Vol. 336, 737-743, 2008.

[3] E. Béchet, H. Minnebo, N. Moës, B. Burgardt, *Improved implementation and robustness study of the X-FEM for stress analysis around cracks*, International Journal for Numerical Method in Engineering 2005; 64, 1033-1056.

[4] E. Chahine, P. Laborde, Y. Renard, *Crack tip enrichment in the XFEM using a cutoff function*, International Journal for Numerical Methods in Engineering 75, 6, 2008.

[5] P. Laborde, J. Pommier, Y. Renard, M. Salaün, *High order extended finite element method for cracked domains*. International Journal for Numerical Method in Engineering 2005; 64, 354-381.

[6] G. Legrain, N. Moes, E. Verron, *Robust and direct evaluation of J(2) in linear elastic fracture mechanics with the X-FEM*, International Journal for Numerical Methods in Engineering 76, 2008.

[7] G.A. Holzapfel, *On large strain viscoelasticity: Continuum formulation and finite element applications to elastomeric structures*, International Journal for Numerical Methods in Engineering, 26, 1996.
<figure><img src="image_5.png" /><figcaption>Fig. 4 – Stress relaxation of reference solution.</figcaption></figure>