# Schémas numériques conservatifs pour des problèmes de vibro-impacts de poutres et de plaques

Cédric Pozzolini, Yves Renard, Michel Salaün

► To cite this version:

Cédric Pozzolini, Yves Renard, Michel Salaün. Schémas numériques conservatifs pour des problèmes de vibro-impacts de poutres et de plaques. 10e colloque national en calcul des structures, May 2011, Giens, France. hal-00592794

HAL Id: hal-00592794

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

Submitted on 3 May 2011

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
# SCHEMAS NUMERIQUES CONSERVATIFS POUR DES PRO- BLEMES DE VIBRO-IMPACTS DE POUTRES ET DE PLAQUES

C. POZZOLINI<fnref n="1" />,<fnref n="2" />, Y. RENARD<fnref n="3" />, M. SALAUN<fnref n="4" />

**Résumé** — Nous proposons une adaptation de la méthode de Dynamique Singulière pour certains schémas numériques modélisant les vibro-impacts d'une poutre ou d'une plaque sur des obstacles rigides. Elle est appliquée sur deux problèmes pour lesquels on a démontré l'existence d'une solution, mais sans avoir d'information sur le bilan énergétique. Par cette nouvelle méthode d'approximation par éléments finis les semi-discrétisations en espace sont stables et bien posées à la différence des approximations classiques.

**Mots clés** — Inéquations Variationnelles, Méthodes Eléments Finis, Plaques et poutres élastiques, Dynamique sous contraintes unilatérales, Coefficient de restitution.

## 1 Introduction

Dans le but de mieux comprendre le comportement non-linéaire des panneaux solaires d'un satellite durant le lancement et éviter l'endommagement de la structure, nous considérons la dynamique d'une poutre ou d'une plaque élastiques vibrant entre deux obstacles rigides. Nous modélisons le contact par des conditions de type Signorini entre le déplacement et la réaction de contact. Pour cette classe de problèmes, nous proposons différentes familles de schémas totalement discrétisés et nous explorons l'in- fluence d'une loi d'impact donnée par un coefficient de restitution entre les vitesses avant et après choc. De façon classique les schémas totalement discrétisés pour ce type de problèmes combinent une approxi- mation éléments finis en espace, et un schéma en temps. Mais, dans ce cadre la plupart des méthodes sont instables, ou présentent des oscillations numériques parasites. Nous avons adapté la méthode de Dynamique Singulière proposée par Y. Renard [10] (inspirée de méthode de redistribution de masse [5]). Il s'agit de faire une formulation semi-discrète en espace qui résout ces difficultés en changeant la matrice de masse dans l'inéquation variationnelle. Ceci empêche que les forces inertielles créent des oscillations numériques parasites et fournit des solutions conservatives quand le coefficient de restit- tion est égal à 1, ce qui assure la stabilité. Nous décrivons et comparons dans cet exposé les solutions de ces familles de schémas totalement discrétisés et leur convergence vers le problème continu pour les impacts d'une poutre ou d'une plaque. De plus, nous avons étendu au cas des plaques de Kirchhoff- Love [8] un théorème d'existence [3], de convergence et de stabilité relatif à la dynamique d'une poutre d'Euler-Bernouilli contrainte par des obstacles rigides.

## 2 Problèmes d'élasto-vibro-impacts

Pour le problème de vibro-impact entre une poutre d'Euler-Bernoulli et un obstacle rigide, nous disposons du résultat de [3], où la convergence d'une solution du problème complètement discrétisé vers celle du modèle continu est établie, mais sans que rien ne soit dit sur la conservation de l'énergie à la limite. Une des difficultés provient de la discontinuité du champ de vitesses aux instants d'impact. De plus des contre-exemples sont connus à l'unicité d'une solution. Dans [8] nous avons exhibé un schéma convergeant de type Newmark pour le cas de la dynamique sans frottement du modèle de plaque de Kirchhoff-Love (K-L) sous les conditions de non-pénétration de Signorini.

<footnote n="1">¹ Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, Villeurbanne, France.</footnote>

<footnote n="2">² Centre National d'Etudes Spatiales, 18 Avenue Edouard Belin, 31401 Toulouse, France.</footnote>

<footnote n="3">³ Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, LaMCoS UMR5259, F-69621, Villeurbanne, France.</footnote>

<footnote n="4">⁴ Université de Toulouse, ICA, ISAE, 10 Avenue Edouard Belin, 31055 Toulouse, France.
email: cedric.pozzolini@insa-lyon.fr; Yves.Renard@insa-lyon.fr; michel.salaun@isae.fr</footnote>
Considérons une plaque mince d'épaisseur $2\epsilon$ et de plan moyen $\Omega \subset \mathbb{R}^2$ un ouvert, borné, connexe du plan, ayant pour bord $\partial\Omega$. Nous la supposons libre sur $\Gamma_f = \partial\Omega \setminus \Gamma_c$, et encastrée sur un ensemble de mesure non-nulle $\Gamma_c \subset \partial\Omega$

$$ \frac{\partial u}{\partial x}(x,t) = 0, \quad u(x,t) = 0, \quad \forall x \in \Gamma_c. \qquad (1) $$

Le bord $\partial\Omega$ est initialement fixé et plat et sa normale extérieure est notée $n$. Par ailleurs, le matériau constituant la plaque est supposé homogène et isotrope, de masse volumique $\rho$, de module d'Young $E$ et de coefficient de Poisson $\nu$. Dans la suite, on notera $D = \frac{2E\epsilon^3}{3(1-\nu^2)}$ la rigidité à la flexion. La plaque est soumise à un chargement vertical $f: \Omega \times [0, T] \to \mathbb{R}$. Le déplacement vertical $u: \Omega \times [0, T] \to \mathbb{R}$ est limité par des obstacles rigides représentés par deux applications $g_1$ et $g_2$ de $\Omega$ dans $\mathbb{R}$ telles qu'il existe $g > 0$ :

$$ g_1(x) \leq -g < 0 < g \leq g_2(x), \quad \forall x \in \Omega. $$

$u$ appartient donc au convexe admissible :

$$ \mathbb{K} = \{v \in \mathcal{V} / g_1(x) \leq v(x) \leq g_2(x), \forall x \in \Omega\}, $$

avec

$$ \mathbb{H} = L^2(\Omega), \quad \mathcal{V} = \{w \in H^2(\Omega) / w(x) = 0 = \partial_n w(x), \forall x \in \Gamma_c\}, $$

puisque l'opérateur associé au modèle de K-L est un bilaplacien. Alors, sous forme faible, la flèche $u$ est solution de l'inéquation variationnelle :

$$ \left\{ \begin{array}{l} \text{Trouver } u: \Omega \times [0, T] \to \mathbb{K} \text{ tel que} \\ \displaystyle \int_0^T \int_{\Omega} \frac{\partial^2 u}{\partial t^2} (w-u) \, dxdt + \int_0^T a(u(\cdot,t), w(\cdot)-u(\cdot)) \, dt \ge \int_0^T \int_{\Omega} f(w-u) \, dxdt, \quad (2) \\ \forall w \in L^2(0,T; \mathbb{K}), \quad u(\cdot,0) = u_0 \in \mathbb{K}, \quad \dot{u}(\cdot,0) = v_0. \end{array} \right. $$

La résultante des efforts extérieurs est la forme linéaire :

$$ l(w) = \int_{\Omega} f(x,t) \cdot w(x,t) \, dx \qquad (3) $$

Le travail virtuel des efforts intérieurs est donné par la forme bilinaire :

$$ a(u,v) = \int_{\Omega} \frac{D}{2\rho\epsilon} \left[ (1-v) \partial_{\alpha\beta}^2 u + v \Delta u \delta_{\alpha\beta} \right] \partial_{\alpha\beta}^2 v \, dx. \qquad (4) $$

Les solutions de (2) représentent la flèche de la plaque depuis sa position initiale, avec $\rho$ sa densité, $D = \frac{2E\epsilon^3}{3(1-\nu^2)}$ la rigidité à la flexion, $E$ est le module de Young, et $0 < \nu < 0.5$ est le coefficient de Poisson.

Dans (2) les solutions ont une régularité faible en espace et en temps, en particulier les vitesses $\frac{\partial u}{\partial t} = \dot{u}$ peuvent être discontinues en temps, y compris pour des données régulières. De plus la discrétisation de la dynamique continue (2) ne décrit pas complètement le mouvement [6], cela oblige à introduire (au moins dans le schéma numérique) une loi d'impact. Par exemple, s'il y a impact en $(x_0, t_0)$, une loi peut être donnée (au moins formellement) par une relation entre la vitesse avant et après impact :

$$ \dot{u}(x_0, t_0^+) = -e\dot{u}(x_0, t_0^-), \quad \text{quand } u(\cdot, t_0) \in \partial\mathbb{K}. \qquad (5) $$

Dans le cas d'une poutre ou d'une plaque les équations sont obtenues par passage à la limite dans l'épaisseur, ce qui élimine de fait l'élasticité transversale. Du coup l'introduction du coefficient de restitution pallie cette absence au moment des impacts.
# 3 Discrétisations

## 3.1 Schémas de type Newmark-Dumont-Paoli

Dans cette section nous présentons une discrétisation spatiale et temporelle du déplacement. Les schémas numériques de type β-Newmark sont stables et convergeants quand la dynamique est linéaire, i.e. quand les contraintes ne sont pas actives. Du fait que l'opérateur différentiel de K-L gouvernant le problème d'impact pour la plaque est du quatrième ordre, il est nécessaire de faire appel à une méthode élémentaire finie de classe $\mathcal{C}^1$, si on souhaite utiliser une méthode conforme. Pour les tests numériques nous utiliserons les éléments de type HCT ou FVS [2]. Nous notons $\Delta t = \frac{T}{N}$ le pas de temps. Nous supposons $h > 0$ et notons $\mathbb{V}^h$ le sous-espace de dimension finie de $\mathbb{V}$ tel que, pour tout $w \in \mathbb{V}$, il existe une suite $(w^h)_{h>0}$ vérifiant : $\|w^h - w\|_{\mathbb{V}} \to 0, w^h \in \mathbb{V}^h, \forall h > 0, h \to 0$. Nous initialisons la discrétisation de (2), par $u_0^h$ et $u_1^h$ dans $\mathbb{K}^h := \mathbb{K} \cap \mathbb{V}^h$, tels que $(u_1^h)_{h>0}$ reste bornée, et $\lim_{h \to 0, \Delta t \to 0} \|u_0^h - u_0\|_{\mathbb{V}} + |\frac{u_1^h - u_0^h}{\Delta t} - v_0|_{\mathbb{H}} = 0$. La force est donnée par :

$$f^n = \frac{1}{\Delta t} \int_{n\Delta t}^{(n+1)\Delta t} f(\cdot, s)ds. \qquad (6)$$

Alors, pour tout $n \in \{1, \dots, N-1\}$ nous résolvons l'inclusion suivante :

$$\text{Trouver } u_{n+1}^h \in \mathbb{K}^h \ / \ (\mathbf{M} + \beta \Delta t^2 \mathbf{K}) u_{n+1}^h + \Delta t^2 \partial_{\mathbb{K}^h} \left( \frac{u_{n+1}^h + e u_{n-1}^h}{1+e} \right) \ni f_n^h \qquad (7)$$

avec $f_n^h = (2\mathbf{M} - \Delta t^2(1-2\beta)\mathbf{K})u_n^h - \beta(\mathbf{M} + \Delta t^2\mathbf{K})u_{n-1}^h + \Delta t^2 F^{n,\beta}$, et

$$F^{n,\beta} = \beta F^{n+1} + (1-2\beta)F^n + \beta F^{n-1}, \qquad (8)$$

avec $F^k$ de composantes $F_i^k = \int_{\Omega} f(x, k\Delta t) \cdot \psi_i(x) d\Omega$. où les $\psi_i$ désignent les $N_{\mathbb{V}}$ fonctions de bases du schéma élémentaire choisis.

La matrice de masse $\mathbf{M} = (\mathbf{M}_{i,j})_{i,j}$ est donnée par

$$\mathbf{M}_{i,j} = \langle \psi_i, \psi_j \rangle_{\mathbb{H}}. \qquad (9)$$

et la matrice de rigidité $\mathbf{K} = (\mathbf{K}_{i,j})_{i,j}$ est donnée par

$$\mathbf{K}_{i,j} = a(\psi_i, \psi_j). \qquad (10)$$

### Théorème 1 (Résultat d'existence [8])

Soient $e = 0, \beta = 1/2$, pour tout $h > 0$ et $N \ge 1$ pour $n \in \{1, \dots, N-1\}$, on note $(u_{n+1}^h)_{1 \le n \le N-1}$ les solutions de (7). La suite des approximations $(u_{h,N})_{h>0,N \ge 1}$ obtenues par interpolation linéaire :

$$u_{h,N}(x,t) = u_n^h(x) \frac{(n+1)\Delta t - t}{\Delta t} + u_{n+1}^h(x) \frac{t - n\Delta t}{\Delta t}, \text{ pour } n\Delta t \le t \le (n+1)\Delta t$$

admet une sous-suite qui converge faible-\* dans $\{w \in L^\infty(0,T; \mathbb{K}), \dot{w} \in L^\infty((0,T); \mathbb{H})\}$ vers une solution de (2).

De par l'hypothèse $e = 0$ (de chocs inelastique), ce schéma de type Newmark est dissipatif.

# 4 La méthode de dynamique singulière

L'idée, développée dans [10], est de construire une méthode ayant les propriétés de la méthode de redistribution de masse 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 $\mathbb{V}^h$ pour l'approximation du déplacement et on considère un second espace élément fini $\mathbb{H}^h$ destiné à approcher la vitesse. Les multiplicateurs $\lambda^h$
appartiennent à l'espace $\Lambda^h$, ils approximent les réactions de contact $\lambda \in \mathbb{V}$. Le problème semi-discrétisé s'écrit alors :

$$
\left\{
\begin{array}{l}
\text{Trouver } u^h : [0, T] \to \mathbb{V}^h, v^h : [0, T] \to \mathbb{H}^h \text{ et } \lambda^h : [0, T] \to \mathbb{V}^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^h, v^h \rangle, \quad \forall w^h \in \mathbb{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 \mathbb{H}^h, \forall t \in ]0, T], \\
\lambda^h \in \Lambda^h, \quad \langle \mu^h - \lambda^h, u^h \rangle \ge 0, \quad \forall \mu^h \in \Lambda^h, \\
u^h(0) = u_0^h, \quad v^h(0) = v_0^h.
\end{array}
\right.
\tag{11}
$$

Bien sûr quand $\mathbb{V}^h = \mathbb{H}^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(t) \in \bar{\Lambda}_N^h, \quad (\lambda(t) - \mu)^T B_N U \ge 0, \quad \forall \mu \in \bar{\Lambda}_N^h, \quad \forall t \in ]0, T], \\
U(0) = U_0, \quad V(0) = V_0.
\end{array}
\right.
\tag{12}
$$

où les vecteurs L, U et V ont pour coordonnées

$$
L_i = l(\psi_i), \quad u^h = \sum_{i=1}^{N_{\text{V}}} U_i \psi_i, \quad v^h = \sum_{i=1}^{N_{\text{H}}} V_i \phi_i.
$$

La matrice **K** a pour composante $\mathbf{K}_{i,j} = a(\psi_i, \psi_j)$, $C$ est la matrice de masse obtenue à l'aide des éléments finis $\mathbb{H}^h$ et $B$ est la matrice $B_{i,j} = \int_{\Omega} \rho \varphi_i \psi_j dx$ où $\varphi_i$ sont les $N_{\mathbb{H}}$ fonctions de forme de $\mathbb{H}^h$. Comme la matrice $C$ est toujours inversible, on a $V(t) = C^{-1}B\dot{U}(t)$ et donc en notant $\mathbf{M}_s = B^T C^{-1}B$, l'inconnue $V$ peut être éliminée et le problème (12) se réécrire sous la forme (7) à ceci près que la matrice de masse $\mathbf{M}_s$ est singulière ici.

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

En particulier, une discrétisation de type HCT, Argyris, ou FVS pour le déplacement $u^h$ et constante par éléments pour la vitesse $v^h$ avec une condition de contact aux noeuds, est une discrétisation qui satisfait expérimentalement la condition inf-sup. Les tests numériques donnent des résultats conservatifs comme dans le cas de la méthode de redistribution de masse. Le point commun avec cette dernière méthode est que lorsque la condition inf-sup est satisfaite la matrice de masse est aussi singulière.

4.1 Implémentation

Ici nous notons pour simplifier $u_n^h = U^n$. Pour expliquer comment l'algorithme fonctionne dans le
cas général, commençons par réécrire le problème avec $e = 0$ comme suit :

$$
\left\{
\begin{array}{l}
\text{Trouver } U^{n+1} \in \mathbb{K}^h \text{ tel que pour tout } W \in \mathbb{K}^h \\
\quad (W - U^{n+1})^T \left( \frac{\mathbf{M}}{\Delta t^2} U^{n+1} + \beta \mathbf{K} U^{n+1} \right) \ge (W - U^{n+1})^T \tilde{\mathbf{F}}^n
\end{array}
\right.
\tag{13}
$$

avec

$$
\tilde{F}^n = F^{n,\beta} - \mathbf{M} \frac{-2U^n + U^{n-1}}{\Delta t^2} - \mathbf{K} \left( (1-2\beta)U^n + \beta U^{n-1} \right),
$$

$F^{n,\beta}$ donné par (8). Puisque la matrice $\mathbf{A} \equiv \frac{1}{\Delta t^2} \mathbf{M} + \beta \mathbf{K}$ est définie positive et symétrique, l'inéquation (13) est équivalente au problème d'optimisation

$$
U^{n+1} = \operatorname{Argmin}_{W \in \mathbb{K}^h} \left( \frac{1}{2} W^T A W - (\tilde{F}^n)^T W \right) .
$$
Remarquons qu'il s'agit d'un problème de minimisation sous contraintes convexes $\mathbb{K}^h$ qui peut être classiquement résolu de différentes façons. Comme dans [3], nous utilisons la fonction Matlab "quadprog", qui est basées sur une méthode de multiplicateurs de Lagrange.

Pour le cas $e \neq 0$, nous avons la formulation suivante :

$$ egin{cases} 
\text{Trouver } U^{n+1,e} = \frac{U^{n+1} + eU^{n-1}}{1+e} \in \mathbb{K}^h \text{ tel que pour tout } W \in \mathbb{K}^h \\ 
(W - U^{n+1,e})^T \left( \mathbf{M} \frac{U^{n+1} - 2U^n + U^{n-1}}{\Delta t^2} + \mathbf{K} (\beta U^{n+1} + (1-2\beta)U^n + \beta U^{n-1}) \right) \\ 
\geq (W - U^{n+1,e})^T F^{n,\beta} 
\end{cases} $$

Alors, il est possible d'exprimer $U^{n+1}$ comme $U^{n+1} = (1+e) U^{n+1,e} - e U^{n-1}$ et de le remplacer dans l'inégalité précédante

$$ egin{cases} 
(W - U^{n+1,e})^T \left( \mathbf{M} \frac{(1+e)U^{n+1,e} - 2U^n + (1-e)U^{n-1}}{\Delta t^2} \right) \\ 
+ \mathbf{K} (\beta(1+e)U^{n+1,e} + (1-2\beta)U^n + \beta(1-e)U^{n-1}) \geq (W - U^{n+1,e})^T F^{n,\beta} 
\end{cases} $$

ce qui nous conduit à devoir résoudre

$$ egin{cases} 
\text{Trouver } U^{n+1,e} \in \mathbb{K}^h \text{ tel que pour tout } W \in \mathbb{K}^h \\ 
(W - U^{n+1,e})^T \left( \frac{\mathbf{M}}{\Delta t^2} U^{n+1,e} + \beta \mathbf{K} U^{n+1,e} \right) \geq (W - U^{n+1,e})^T G^n 
\end{cases} $$

avec

$$ G^n = \frac{1}{1+e} \left( \tilde{F}^n + e \left( \frac{1}{\Delta t^2} \mathbf{M} + \beta \mathbf{K} \right) U^{n-1} \right) . $$

Il s'apparente à (13) et ce résoud de la même manière.

En conséquence la procédure de résolution pour (7) suit celle de Paoli-Schatzman [7]. Nous commençons par calculer $Q^{n+1} = \mathbf{A}^{-1} \tilde{F}^n$. Puis nous arrivons à l'alternative suivante :

* Si $\frac{Q^{n+1} + eU^{n-1}}{1+e}$ appartient à $\mathbb{K}^h$, alors l'équation (14) est vérifiée, et nous posons $U^{n+1} = Q^{n+1}$;

* Si $\frac{Q^{n+1} + eU^{n-1}}{1+e}$ n'appartient pas à $\mathbb{K}^h$, nous résolvons (14) par la fonction "quadprog", et nous avons $U^{n+1} = (1+e) U^{n+1,e} - e U^{n-1}$.

Bien entendu le fait que $U^{n+1,e}$ soit dans $\mathbb{K}^h$ ne garanti pas que $U^{n+1}$ soit aussi dans $\mathbb{K}^h$, mais seulement proche. Ceci est aussi souligné par Paoli-Schatzman [7].

**Remark 1** Le problème avec la masse singulière se traite exactement de la même façon car la matrice $\frac{1}{\Delta t^2} \mathbf{M}_s + \beta \mathbf{K}$ a les même propriétés que la matrice $\mathbf{A}$.

# 5 Résultats numériques

Nous considérons une plaque rectangulaire $\Omega$ en acier de longueur 120cm, de largeur 40cm, et d'épaisseur $\varepsilon = 1cm$ avec $E = 210 \text{ GPa}$, $\rho = 7.77 \times 10^3 \text{ kg/m}^3$, et $\nu = 0.3$. Supposons que la plaque est au dessus d'un obstacles rigide plan i.e. $g_2 = +\infty$ et :

$$ g_1(x) = -0.1, \quad \forall x = (x_1, x_2) \in \Omega. $$

Tous les tests numériques ont été réalisé avec GETFEM++ [4] via l'interface Matlab. Nous souhaitons mettre en évidence la conservation de l'énergie par nos schémas numériques. Pour cela nous imposons
<figure><img src="image_2.png" /><figcaption>FIG. 1 – Etat initial de la plaque sous le chargement $f_0$, maillage quadrangulaire (FVS).</figcaption></figure>

un chargement initial $f(x,t=0) = f_0 = 8600 \text{ N}$ (puis $f(x,t) = 0$ pour $(x,t) \in \Omega \times ]0, T]$), nous déduisons
le déplacement $u_0$ associé à l'instant $t=0$, et la vitesse à l'instant $t=0$, $v_0 = 0$.

L'énergie totale associée s'écrit

$$
\mathcal{E}(w, \dot{w}) := \int_{\Omega} \left[ \frac{1}{2}(\dot{w})^2(x,t) + \frac{1}{2}a(w,w) - f(x,t)w(x,t) \right] d\Omega,
$$

la dissipation signifie que $\mathcal{E}(u^{n+1}, v^{n+1}) \le \mathcal{E}(u^n, v^n)$. Rappelons que l'énergie cinétique est conservée si et seulement si $e=1$. Nous observons ce que fournissent les calculs quand $\Delta x$ et $\Delta t$ tendent vers zero.

5.1 Oscillations parasites

Nous traçons ci-dessous les déplacements des coins libres et celui du milieu des coins libres. Nous obtenons qualitativement les même résultats que dans [9] pour des poutres impactant des obstacles rígides.

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Déplacements des coins libres et celui du milieu des coins libres, méthode NDP-FVS masse régulière, $e = 0$, $\Delta t = 10^{-5} \text{ sec.}$, 140 éléments.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Zoom sur les premiers impacts des coins libres et celui du milieu des coins libres, méthode NDP-FVS masse régulière, $e = 0$, $\Delta t = 10^{-5}$ sec., 140 éléments.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4 – Zoom sur les premiers impacts des coins libres et celui du milieu des coins libres, méthode NDP-FVS masse régulière, $e = 0$, $\Delta t = 5 \times 10^{-6}$ sec., 140 éléments.</figcaption></figure>

## 5.2 Stabilité énergétique

Quand $e = 0$ nous obtenons aussi les même résultats qualitatifs que [1] et [3] pour une poutre.

D'autres résultats numériques avec différentes valeurs de $e$ et avec différentes masses singulières seront présentés lors du colloque.
<figure><img src="image_6.png" /><figcaption>FIG. 5 – Evolution de l’énergie totale au cours du temps (3 sec.), méthode NDP-FVS masse régulière, $e = 0$, 140 éléments.</figcaption></figure>

# 6 Remerciements

Nous remercions le CNES et l'Institut Camille Jordan pour leur soutien financier.

## Références

[1] Ahn, J. and Stewart, D.E. An Euler-Bernoulli beam with dynamic contact : discretisation, convergence and numerical results, SIAMJ. Numer. Anal., Vol. 43, No. 4, pp.1455-1480, (2005).

[2] Ciarlet P.G., The finite element method for elliptic problems, North-Holland, (1978).

[3] Dumont, Y. and Paoli, L., Vibrations of a beam between stops : convergence of a fully discretised approximation, ESAIM : Math. Modelling Num. Ana. (M2AN), Vol. 40, No. 4, pp.705-734, (2006).

[4] Y. Renard, J. Pommier, An open source generic C++ library for finite element methods.
http://home.gna.org/getfem/

[5] H.B. Khenous, P. Laborde, Y. Renard, Mass redistribution method for finite element contact problems in elastodynamics, Eur. J. Mech., A/Solids, 27(5), 918-932, (2008).

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

[7] L. Paoli, M. Schatzman, Numerical simulation of the dynamics of an impacting bar, Comput. Methods Appl. Mech. Eng., 196, 2839-2851, 2007.

[8] Pozzolini C., Renard Y., and Salaun M., Vibro-Impact of a plate on rigid obstacles : existence theorem, convergence of scheme and numerical simulations, submitted

[9] Pozzolini C. and Salaun M., Some energy conservative schemes for vibro-impacts of a beam on rigid obstacles, submitted

[10] Renard Y., The singular dynamic method for constrained second order hyperbolic equations. application to dynamic contact problems, J. Comput. Appl. Math., 234(3):906-923, 2010.