# Méthodes de dynamique singulière pour l'élastodynamique avec contact

Yves Renard

► To cite this version:

Yves Renard. Méthodes de dynamique singulière pour l'élastodynamique avec contact. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408470

HAL Id: hal-01408470

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

Submitted on 5 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
# Méthodes de dynamique singulière pour l'élastodynamique avec contact

Y. Renard¹

¹ Université de Lyon, CNRS,
INSA-Lyon, ICJ UMR5208, LaMCoS UMR5259,
F-69621, Villeurbanne, France.
Yves.Renard@insa-lyon.fr

**Résumé** — On présente une nouvelle méthode d'approximation par éléments finis des problèmes d'élastodynamique avec contact unilatéral qui conduit à des semi-discrétisations en espace stables et bien posées à la différence des approximations classiques. Il s'agit en quelque sorte d'une évolution de la méthode de redistribution de masse [7, 8]. Ces méthodes permettent aux schémas classiques d'intégration en temps (Newmark, HHT ...) d'être stables et de calculer une approximation exploitable de la contrainte de contact.

**Mots clés** — élastodynamique avec contact, méthodes d'éléments finis, inéquations variationnelles évolutives.

## 1 Introduction

Il est bien connu que les conditions de contact unilatéral dans les problèmes d'élastodynamique sont source d'instabilités pour les schémas numériques. Ainsi, la grande majorité des approximations classiques font apparaître des oscillations parasites sur le déplacement normal et sur les contraintes de contact ou font intervenir un amortissement excessif pour stabiliser ces instabilités (voir par exemple [7, 4, 5]). De plus, ces oscillations n'ont pas tendance à disparaître lorsque le pas de temps diminue. Pour beaucoup de schémas, elles ont au contraire tendance à augmenter. Ces difficultés, très spécifiques aux conditions de contact en dynamique, on donné lieu a un nombre important de travaux. Nous renvoyons par exemple à [15, 2, 6, 16, 10, 9, 4, 1, 5, 13]. Il s'agit en général, soit d'ajouter une dissipation suffisante pour stabiliser le schéma, soit au contraire, d'imposer la conservativité du schéma. La condition de contact sans frottement n'étant pas dissipative, il est clair qu'une dissipativité trop importante est contraire à la physique du problème. Les schémas conservatifs en énergie proposés dans ces études font quant à eux apparaître ou des non respects locaux de la contrainte de contact qui peuvent être importants ou bien des oscillations importantes sur la contrainte de contact (voire les deux). De notre point de vue, les solutions proposées dans ces études ne sont pas pleinement satisfaisantes. La raison principale est qu'elles sont centrées sur une adaptation du schéma temporel alors que le véritable problème se situe plutôt sur l'approximation en espace dont il est connu qu'elle est mal posée.

L'objectif de ce papier est de présenter une nouvelle classe de méthode dont le principe est de faire une semi-discrétisation en espace différente de la vitesse et du déplacement. Une condition de type inf-sup entre les deux espaces élément fini et la condition de contact discrète assure le caractère bien posé de cette semi-discrétisation. Les méthodes proposées sont du même type que la méthode de redistribution de masse présentée dans [8, 3] dans le sens où cela amène à une matrice de masse singulière. Cela justifie le nom proposé de méthodes de dynamique singulière.
<figure><img src="image_2.png" /><figcaption>Figure 1 – Solide élastique en contact avec une fondation rigide plane</figcaption></figure>

A noter que la présentation est faite en petites déformations, sans frottement et contre une fondation rigide. La condition de frottement de Coulomb, bien que posant d'autres types de problèmes, peut être assez facilement prise en compte dans des schémas de manière stable due à sa dissipativité. Les généralisations pour des problèmes à deux corps élastiques, ou en grandes déformations et/ou avec frottement ne pose pas de difficultés majeures de mise en œuvre.

## 2 Problème d'élastodynamique avec contact

On considère le problème élastodynamique d'un solide élastique $\Omega$ astreint à une condition de contact sans frottement contre une fondation rigide plane sur une partie $\Gamma_c$ de son bord (voir figure 1). Sur le reste du bord, on impose une condition de Dirichlet sur $\Gamma_D$ et une condition de Neumann sur $\Gamma_N$. Cela donne le problème suivant :

$$ \rho\ddot{u} - \operatorname{div} \boldsymbol{\sigma}(u) = f, \quad \text{dans } ]0, T] \times \Omega, \quad (1) $$

$$ \boldsymbol{\sigma}(u) = \mathcal{A}\boldsymbol{\varepsilon}(u), \quad \text{dans } ]0, T] \times \Omega, \quad (2) $$

$$ \boldsymbol{\sigma}(u)\mathbf{n} = g, \quad \text{on } ]0, T] \times \Gamma_N, \quad (3) $$

$$ u = 0, \quad \text{on } ]0, T] \times \Gamma_D, \quad (4) $$

$$ u_N \le 0, \quad \boldsymbol{\sigma}_N(u) \le 0, \quad u_N \boldsymbol{\sigma}_N(u) = 0, \quad \boldsymbol{\sigma}_T(u) = 0, \quad \text{on } ]0, T] \times \Gamma_C, \quad (5) $$

$$ u(0,x) = u_0(x), \quad \dot{u}(0,x) = u_1(x), \quad \text{dans } \Omega, \quad (6) $$

où $u$ est le déplacement, $\boldsymbol{\sigma}$ le tenseur des contraintes, $\mathcal{A}$ le tenseur des raideurs élastiques, $u_N = u \cdot n$ le déplacement normal sur $\Gamma_c$, $\boldsymbol{\sigma}_N = (\boldsymbol{\sigma}n) \cdot n$ la contrainte normale et $\boldsymbol{\sigma}_T$ la contrainte tangentielle sur $\Gamma_c$. La forme faible de ce problème peut s'écrire de la manière suivante :

$$ \left\{ \begin{array}{l} \text{Trouver } u : [0, T] \longrightarrow V \text{ et } \lambda_N : [0, T] \longrightarrow X'_N \text{ qui satisfont} \\ \langle \rho \ddot{u}, v \rangle_{V',V} + a(u,v) = l(v) + \langle \lambda_N, v_N \rangle_{X'_N \times N}, \quad \forall v \in V. \\ \\ \lambda_N \in \Lambda_N, \quad \langle \mu_N - \lambda_N, u_N \rangle_{X'_N \times N} \ge 0, \quad \forall \mu_N \in \Lambda_N, \\ \\ u(0) = u_0, \quad \dot{u}(0) = u_1. \end{array} \right. $$

où

$$ V = \{v \in H^1(\Omega; \mathbb{R}^3) : v = 0 \text{ sur } \Gamma_D\}, $$

$$ a(u,v) = \int_{\Omega} \mathcal{A}\varepsilon(u) : \varepsilon(v), \quad l(v) = \int_{\Omega} f.vdx + \int_{\Gamma_N} g.vd\Gamma, $$

$$ X_N = \{v_N |_{\Gamma_C} : v \in V\}, \Lambda_N = \{\mu_N \in X'_N : \langle \mu_N, v_N \rangle_{X'_N \times N} \ge 0, \forall v_N \in X_N, v_N \le 0\}. $$

Autrement dit, $V$ est l'espace des déplacements qui sont d'énergie élastique finie, $X_N$ est l'espace des traces normales sur $\Gamma_c$ et $\Lambda_N$ est l'ensemble des contraintes normales admissibles. Dans le problème, $\lambda_N$ représente la contrainte de contact. Les notations $\langle \ddot{u}, v \rangle_{V',V}$ et $\langle \lambda_N, v_N \rangle_{X'_N \times N}$ dénotent les produits de dualité entre les espaces respectifs qui coïncident avec $\int_{\Omega} \ddot{u}.vdx$ et $\int_{\Gamma_c} \lambda_N v_N d\Gamma$ dans le cas régulier.
### 3 Semi-discrétisation classique par une méthode d'éléments finis

Une discrétisation classique par éléments finis fait intervenir un espace élément fini $V^h$ qui est en général un sous-espace de dimension finie de $V$ et une approximation de l'ensemble des contraintes normales admissibles $\Lambda_N^h$ qui n'est pas nécessairement un sous-ensemble de $\Lambda_N$. La formulation du problème semi-discrétisé est alors donnée par

$$ egin{cases} 
\text{Trouver } u^h : [0, T] \to V^h, \text{ et } \lambda_N^h : [0, T] \to X_N^h \text{ qui satisfont} \\ 
\int_{\Omega} \ddot{u}^h \cdot v^h dx + a(u^h, v^h) = l(v^h) + \int_{\Gamma_c} \lambda_N^h v_N^h, \quad \forall v^h \in V^h, \forall t \in ]0, T], \\ 
\lambda_N^h \in \Lambda_N^h, \quad \int_{\Gamma_c} (\mu_N^h - \lambda_N^h) u_N^h \ge 0, \quad \forall \mu_N^h \in \Lambda_N^h, \forall t \in ]0, T], \\ 
u^h(0) = u_0^h, \quad \dot{u}^h(0) = u_1^h, 
\end{cases} 
\quad (7) $$

où $u_0^h, u_1^h \in V^h$ sont des approximations des conditions initiales. Si on note $\varphi_1, \varphi_2, \cdots \varphi_N$ les fonctions de forme de $V^h$ de telle manière que $u$ s'écrive $u(t,x) = \sum_{i=1}^N u^i(t)\varphi_i(x)$ on arrive à la formulation matricielle suivante du problème semi-discrétisé

$$ egin{cases} 
M\ddot{U}(t) + KU(t) = L + B_N^T \lambda_N(t), \quad \forall t \in ]0, T], \\ 
\lambda_N(t) \in \bar{\Lambda}_N^h, \quad (\lambda_N(t) - \mu_N)^T B_N U \ge 0, \quad \forall \mu_N \in \bar{\Lambda}_N^h, \forall t \in ]0, T], \\ 
U(0) = U_0, \quad V(0) = U_1. 
\end{cases} 
\quad (8) $$

où $M$ est la matrice de masse, $K$ la matrice de rigidité, $U(t) \in \mathbb{R}^{N_d}$ le vecteur des degrés de liberté en déplacement. L'expression de la matrice $B_N$ dépend de l'approximation $\bar{\Lambda}_N^h$ de l'ensemble des contraintes normales admissibles. Ce problème peut être vu comme étant une inclusion différentielle à solution mesure (voir [12, 13] par exemple). Malheureusement, ce type de problème est mal posé et possède en général une infinité de solutions. Cette non-unicité peut être illustrée en considérant le système très simple représenté sur la figure 2.

Équation d'évolution :

$$ egin{cases} 
m\ddot{u} + ku = -F_N, \\ 
F_N = (F_N - ru_N)_-, \quad \forall r > 0, \\ 
u(0) = u_0, \quad \dot{u}(0) = u_1. 
\end{cases} $$

Solutions pour $u_0 = 1$ et $u_1 = 0$ :

$$ u(t) = \cos(\sqrt{k/m} t), \quad \text{for } t < \frac{\pi}{2} \sqrt{m/k}, $$

$$ u(t) = a \cos(-\sqrt{k/m} t), \quad \text{for } \frac{\pi}{2} \sqrt{m/k} < t < \frac{3\pi}{2} \sqrt{m/k}. $$

<figure><img src="image_3.png" /><figcaption>Figure 2 - Système simple en contact unilatéral</figcaption></figure>

Ces solutions sont données pour tout $a \ge 0$. On a donc bien une infinité de solutions qui correspondent en fait à des vitesses différentes de rebond allant de $a = 0$ pour un choc inélastique à $a$ aussi grand que l'on veut avec un gain d'énergie cinétique lors de l'impact. Ce système très simple illustre le comportement du problème semi-discrétisé en chaque nœud de contact potentiel.

En conséquence, lorsque l'on approche ce système à l'aide d'un schéma en temps, on ne maîtrise pas la vitesse de rebond des nœuds en contact. La plupart des schémas (mise à part le schéma d'Euler implicite) donnent des vitesses après impact supérieures aux vitesses avant impact et ceci d'autant plus que le pas de temps est petit.

Une remarque importante est d'ailleurs que l'on ne peut pas sélectionner de solution acceptable du problème (8). Si on choisit des rebonds complètement inélastique, alors une énergie est perdue à
<figure><img src="image_4.png" /></figure>

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

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

<figure><img src="image_7.png" /><figcaption>Figure 3 – Énergie totale et contrainte de contact en un point de Γc pour un schéma de Newmark β = γ = 1/2 et une discrétisation par éléments finis classique.</figcaption></figure>

chaque impact (cette perte d'énergie est indépendante du pas de temps) ce qui n'est pas compatible avec le caractère non dissipatif de la condition de contact. Si on choisit des choc élastiques, les nœuds en contact potentiel ne peuvent entrer en contact persistant et rebondissent indéfiniment ce qui entraîne des oscillations parasites. Cette analyse montre bien que les semi-discrétisations classiques ne sont pas satisfaisantes. On peut voir sur la figure 3 un test fait avec un schéma de Newmark (β = γ = 1/2) qui paraît stable pour un pas de temps de 10⁻³ et devient instable pour un pas de temps de 10⁻⁴. A noter que le schéma de Newmark β = 1/4, γ = 1/2 qui est conservatif en énergie pour l'élastodynamique linéarisée est instable pour n'importe quel pas de temps en présence de contact.

4 La méthode de redistribution de masse

La méthode de redistribution de masse a été introduite par Khenous, Laborde et Renard dans
[8]. Le constat de départ est que le caractère mal posé du problème semi-discrétisé est principa-
lement dû au fait que les nœuds en contact potentiel possèdent une inertie dans les problèmes
semi-discrétisés alors que ce n'est pas le cas en continu. L'idée, très simple, est alors de faire en
sorte que ces nœuds ne possèdent plus d'inertie en redistribuant la masse de manière à ce que la
matrice de masse modifiée soit telle que pour U.Nᵢ, i ∈ Ic représentant les déplacements normaux
sur le bord de contact on ait

$$
N_i^T M N_j = 0, \quad \forall i, j \in I_c.
$$

Dans [8] un algorithme très simple est proposé qui permet de modifier la matrice de masse tout en
conservant ses principaux moments. Il est de plus montré que les propriété suivantes :

- le problème semi-discrétisé devient bien posé et peut-être interprété comme un système d'équations différentielles ordinaires lipschitziennes.

- L'unique solution conserve l'énergie totale.

La conséquence est que la plupart des schémas d'intégration en temps classiques qui sont
stables pour l'élastodynamique deviennent convergents pour ce problème. C'est le cas sur la fi-
gure 4 pour le schéma de Newmark (ici pour $\beta = 1/2, \gamma = 1/2$) à comparer avec la figure 3.
<figure><img src="image_8.png" /></figure>

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

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

<figure><img src="image_11.png" /><figcaption>Figure 4 – Énergie totale et contrainte de contact en un point de Γ<sub>c</sub> pour un schéma de Newmark β = γ = 1/2 et une redistribution de masse</figcaption></figure>

5 La nouvelle méthode

L'idée, développée dans [14], est de construire une méthode ayant les propriétés de la méthode
de redistribution de masse de manière plus “naturelle”, c’est à dire sans modification *a posteriori*
de la matrice de masse. On arrive à ce résultat en effectuant des approximations de la vitesse et
du déplacement par des méthodes d’éléments finis différentes. Autrement dit, on garde $V^h$ pour
l’approximation du déplacement et on considère un second espace élément fini $H^h \subset V$ destiné à
approcher la vitesse. Le problème semi-discrétisé s’écrit alors :

$$
\left\{
\begin{array}{ll}
\text{Trouver } u^h : [0, T] \to V^h, v^h : [0, T] \to H^h \text{ et } \lambda_N^h : [0, T] \to X_N^{h'} \text{ qui satisfont} & \\
\displaystyle\int_{\Omega} \rho \frac{\partial v^h}{\partial t} (w^h - u^h) dx + a(u^h, w^h - u^h) = l(w^h - u^h) + \langle \lambda_N^h, v^h \rangle_{X_N^{h', X_N}}, \quad \forall w^h \in V^h, \forall t \in ]0, T], & \\
\displaystyle\int_{\Omega} (v^h - \frac{\partial u^h}{\partial t}) q^h dx = 0 \quad \forall q^h \in H^h, \forall t \in ]0, T], & \\
\lambda_N^h \in \Lambda_N^h, \quad \langle \mu_N^h - \lambda_N^h, u_N^h \rangle_{X_N^{h', X_N}} \ge 0, \quad \forall \mu_N^h \in \Lambda_N^h, & \\
u^h(0) = u_0^h, \quad v^h(0) = v_0^h. & (9)
\end{array}
\right.
$$

Bien sûr quand $H^h = W^h$ on retrouve l'approximation classique. Voyons maintenant l'expression matricielle de ce problème :

$$
\left\{
\begin{array}{l}
B^T \dot{V}(t) + AU(t) = L + B_N^T \lambda_N(t), \quad \forall t \in ]0, T], \\
CV(t) = B \dot{U}(t), \quad \forall t \in ]0, T], \\
\lambda_N(t) \in \bar{\Lambda}_N^h, \quad (\lambda_N(t) - \mu_N)^T B_N U \ge 0, \quad \forall \mu_N \in \bar{\Lambda}_N^h, \quad \forall t \in ]0, T], \\
U(0) = U_0, \quad V(0) = V_0.
\end{array}
\right.
\qquad (10)
$$

où C est la matrice de masse obtenue à l'aide de l'élément fini $H^h$ et B est la matrice $B_{i,j} =$
$\int_{\Omega} \rho \varphi_i \psi_j dx$ où $\psi_j$ sont les fonctions de forme de $H^h$. Comme la matrice C est toujours inversible,
on a $V(t) = C^{-1}B\dot{U}(t)$ et donc en notant

$$M = B^T C^{-1} B,$$

L'inconnue $V$ peut donc être éliminée et le problème (10) peut alors se récrire sous la forme (8) à ceci près que la matrice de masse $M$ peut être singulière ici.

Dans [14], est donnée une condition de type inf-sup qui lie $V^h$, $H^h$ et $\bar{\Lambda}_N^h$ pour que le problème (10) soit bien posé. Cette condition n'est satisfaite que si $H^h$ est strictement moins riche que $V^h$ sur le bord de contact et dépend aussi de la manière de prendre en compte la condition de contact. On montre aussi dans [14] que la solution est alors d'énergie constante.

En particulier, une discrétisation de type quadratique (élément $P_2$) pour le déplacement $u^h$ et affine (élément $P_1$) pour la vitesse $v^h$ avec une condition de contact aux sommets des faces sur $\Gamma_c$ est une discrétisation qui satisfait expérimentalement la condition inf-sup. Les tests numériques donnent des résultats semblables dans ce cas à la méthode de redistribution de masse. Le point commun avec cette méthode est que lorsque la condition inf-sup est satisfaite la matrice de masse est aussi singulière. En revanche, une approximation en espace d'ordre quelconque est envisageable et une masse est présente sur le bord de contact.

# 6 Références bibliographiques

[1] M. Barboteu. *An efficient algorithm to solve large non linear elastodynamic problems with contact friction*. Bull. Math. Soc. Sc. Math. Roumanie, 48:96:2(2005), pp119-137.

[2] N. J. Carpenter. *Lagrange constraints for transient finite element surface contact*, Int. J. for Num. Meth. Eng., 32 (1991), pp 103-128.

[3] C. Hager, B. Wohlmuth. *Analysis of a modified mass lumping method for the stabilization of frictional contact problems*, to appear.

[4] P. Hauret. *Numerical methods for the dynamic analysis of two-scale incompressible nonlinear structures*, Thèse de Doctorat, Ecole Polytechnique, France, 2004.

[5] P. Hauret, P. Letallec. *Energy controlling time integration methods for nonlinear elastodynamics and low-velocity impact*, Comput. Meth. Appl. Mech. Engrg., 195 (2006), pp. 4890-4916.

[6] T. J.R. Hugues, R. L. Taylor, J. L. Sackman, A. Curnier, W. KanoKnukulchai. *A finite method for a class of contact-impact problems*, 8 (1976), pp 249-276.

[7] H.B. Khenous. *Problèmes de contact unilatéral avec frottement de Coulomb en élastostatique et élastodynamique. Etude math. et résolution numérique*, thèse de l'INSA de Toulouse, France, 2005.

[8] H.B. Khenous, P. Laborde, Y. Renard. *Mass redistribution method for finite element contact problems in elastodynamics*. to appear in Eur. J. Mech., A/Solids.

[9] T.A. Laursen & V. Chawla. *Design of energy conserving algorithms for frictionless dynamic contact problems*. Int. J. Num. Meth. Engrg, 40 (1997), 863-886.

[10] T.A. Laursen, G.R. Love. *Improved implicit integrators for transient impact problems-geometric admissibility within the conserving framework*. Int. J. Num. Meth. Eng, 53 (2002), 245-274.

[11] J.J. Moreau. *Liaisons unilatérales sans frottement et chocs inélastiques*. C.R.A.S. série II, 296 pp 1473-1476, 1983.

[12] J.J. Moreau, *Numerical aspects of the sweeping process*. Comp. Meth. Appl. Mech. Engrg., 177 (1999), 329-349.

[13] L. Paoli. *Time discretization of vibro-impact*. Phil. Trans. R. Soc. Lond. A., 359 (2001), 2405-2428.

[14] Y. Renard. The singular dynamic method for constrained second order hyperbolic equations, submitted.

[15] K. Schweizerhof, J.O. Hallquist, D. Stillman. *Efficiency Refinements of Contact Strategies and Algorithms in Explicit Finite Element Programming*. In Computational Plasticity, eds. Owen, Onate, Hinton, Pineridge Press, pp 457-482, 1992.

[16] R. L. Taylor, P. Papadopoulos. *On a finite element method for dynamic contact-impact problems*, Int. J. for Num. Meth. Eng., Vol 36 (1993), pp 2123-2140.