# Éléments finis étendus en temps : un autre regard sur les méthodes de Galerkin discontinues en temps

Julien Réthoré, Anthony Gravouil, Alain Combescure

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

## To cite this version:

Julien Réthoré, Anthony Gravouil, Alain Combescure. Éléments finis étendus en temps : un autre regard sur les méthodes de Galerkin discontinues en temps. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01814816

HAL Id: hal-01814816

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

Submitted on 13 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Éléments finis étendus en temps :
un autre regard sur les méthodes de Galerkin discontinues en temps.

J. Réthoré — A. Gravouil — A. Combescure

LaMCoS UMR-5514
INSA Lyon
20, Av. Albert Einstein
69621 Villeurbanne Cedex
Julien.Rethore@insa-lyon.fr

ABSTRACT. The Newmark method [NEW 59] has been considered as a reference time integrator for second order equations. Zienkiewicz [ZIE 77] showed that it is possible to recover particu- lar Newmark schemes using a Finite Element formulation of the time problem. The FE method presents a lot of advantages but the integration of discontinuities is still a problem. Concerning space discontinuous problems, the Partition of Unity method [BAB 97] provides a theoretical framework for many developments like eXtended FEM [BEL 01]. We propose here the use of an enriched basis of time shape functions. We also show that the obtained method is equivalent to some time discontinuous Galerkin methods [HUL 92].

MOTS-CLÉS : dynamique, X-FEM, méthode de Galerkin discontinue

KEYWORDS: dynamics, X-FEM, discontinuous Galerkin method

<footnote>RÉSUMÉ. Les schémas de Newmark [NEW 59] sont considérés comme une référence pour l'intégration des équations du second ordre. Zienkiewicz [ZIE 77] a montré qu'il était possible de retrouver certains des schémas de la famille de Newmark en adoptant un vision Eléments Finis (EF) du problème en temps. La méthode des EF présente de nombreux avantages mais le traitement des discontinuités reste un problème. Concernant les problèmes discontinus en espace, les méthodes de Partition de l'unité (PU) [BAB 97] ont fourni un cadre théorique à de nombreux développements comme les EF Étendus (X-FEM) [BEL 01]. On propose ici d'utiliser une base de fonctions de forme en temps enrichie. On montre aussi que le schéma ainsi obtenu est équivalent à certaines méthodes de Galerkin discontinues (TDGM) [HUL 92].</footnote>
# 1. Méthode des Eléments Finis Etendus en temps (TX-FEM)

Dans cette partie, on se propose d'utiliser la méthode de PU présentée dans [BAB 97] afin d'enrichir l'interpolation de la vitesse. Comme l'ensemble des fonctions linéaires $\{\lambda_i\}_{i=0..N}$ constitue une partition de l'unité dans l'intervalle $[0; T]$, on peut enrichir cette base de fonction de forme linéaire :

$$v(t) = \sum_{i=0}^{N} \lambda_i(t)v_i^c + \sum_{i=0}^{N} \sum_{j=0}^{M} \lambda_i(t)\phi_j(t)v_{i,j}^e \quad [1]$$

où $v_i^c$ sont les degrés de liberté classiques, $\phi_j$ les fonctions de forme enrichies et $v_{i,j}^e$ les degrés de liberté additionnels correspondant à ces fonctions enrichies. Les $\phi_j$ peuvent alors être choisis de manière à capturer précisément la discontinuité. On considère ici l'utilisation de la fonction de Heaviside $H$ comme fonction d'enrichissement. On place ces fonctions à chaque piquet de temps (cf. figure 1). Considérons

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Fonctions de forme dans $I_{n+1}$</figcaption></figure>

maintenant l'intervalle de temps $I_{n+1} = ]t_n; t_{n+1}]$. Dans $I_{n+1}$, la vitesse et le déplacement sont décomposés en une contribution continue et une contribution discontinue :

$$ \begin{aligned} v(t) &= v^c(t) + v^e(t) \\ u(t) &= u^c(t) + u^e(t) \end{aligned} \quad [2] $$

avec

$$ \begin{aligned} v^c(t) &= v_n^c \lambda_n(t) + v_{n+1}^c \lambda_{n+1}(t) \\ v^e(t) &= v_{n+1}^e \lambda_n(t) H(t - t_n) \end{aligned} \quad [3] $$

On propose ensuite d'imposer de façon forte la continuité du déplacement aux piquets de temps $(u(t_n^+)=u(t_n^-))$ et la relation cinématique entre le déplacement et la vitesse $(\dot{u}=v)$. On peut alors écrire les parties continue et discontinue du déplacement comme suit :

$$ \begin{aligned} u^c(t) &= u_n^c + \int_{t_n}^{t} v^c(\tau)d\tau \\ u^e(t) &= u_n^e + \int_{t_n^+}^{t} v^e(\tau)d\tau \end{aligned} \quad [4] $$
En utilisant la méthode des résidus pondérés dans $I_{n+1}$, la continuité de la vitesse étant imposée de façon faible à l'instant $t_n$, on a :

$$ \int_{t_n^+}^{t_{n+1}^-} W(t) [m\dot{v}(t) + ku(t)] dt + W(t_n^+)mv_{n+1}^e = \int_{t_n^+}^{t_{n+1}^-} W(t)f(t)dt \quad [5] $$

Comme des degrés de liberté ont été ajoutés, il est nécessaire de choisir deux fonctions poids indépendantes $\{W_i\}_{i=1,2}$ pour résoudre le problème qui peut être écrit de la façon suivante : connaissant $q_n = [u_n^c, u_n^e, v_n^c, v_n^e]$, trouver $q_{n+1} = [u_{n+1}^c, u_{n+1}^e, v_{n+1}^c, v_{n+1}^e]$ tels que :

$$ \begin{align*} u_{n+1}^c &= u_n^c + \frac{\Delta t}{2} (v_{n+1}^c + v_n^c) \\ u_{n+1}^e &= u_n^e + \frac{\Delta t}{2} v_{n+1}^e \\ [m + \beta_1 \Delta t^2 k] v_{n+1}^c + [(\delta_1 - 1)m + (\gamma_1 - \beta_1) \Delta t^2 k] v_{n+1}^e &= -\Delta t k u_n^c - \Delta t k u_n^e \\ & \quad + [m + (\beta_1 - \gamma_1) \Delta t^2 k] v_n^c + \Delta t [(1 - \gamma_1) f_n + \gamma_1 f_{n+1}] \\ [m + \beta_2 \Delta t^2 k] v_{n+1}^c + [(\delta_2 - 1)m + (\gamma_2 - \beta_2) \Delta t^2 k] v_{n+1}^e &= -\Delta t k u_n^c - \Delta t k u_n^e \\ & \quad + [m + (\beta_2 - \gamma_2) \Delta t^2 k] v_n^c + \Delta t [(1 - \gamma_2) f_n + \gamma_2 f_{n+1}] \end{align*} $$

où $\{\beta_i, \gamma_i, \delta_i\}_{i=1,2}$ sont des paramètres qui dépendent des fonctions de poids $\{W_i\}_{i=1,2}$. Pour résumer, cette formulation a les propriétés suivantes : le déplacement est continu, la vitesse est linéaire dans l'intervalle de temps et peut présenter des discontinuités aux piquets de temps (car la continuité de la vitesse est imposée de façon faible), la contrainte cinématique est satisfaite. Il reste à choisir les fonctions poids $W_1, W_2$ pour donner à ce schéma d'intégration en temps les meilleures propriétés.

REMARQUE. — L'accélération peut être définie comme une fonction constante par morceaux dont la valeur dans $I_{n+1}$ est obtenue par dérivation de l'équation (2) :

$$ a(t) = a_{n+1} = \frac{1}{\Delta t} (v_{n+1}^c - v_{n+1}^e - v_n^c) \quad [6] $$

## 2. Une approche alternative des méthodes de Galerkin discontinues

Comme il est possible de retrouver les formules de la méthode de Newmark en utilisant une vision EF du problème en temps [ZIE 77], le concept TX-FEM (Time eXtended Finite Element) permet d'approcher les TDGM sous un angle nouveau. Considérons le changement de variables suivant :

$$ \begin{align*} v_n^+ &= v_n^c + v_{n+1}^e \\ v_{n+1}^- &= v_{n+1}^c \end{align*} \quad [7] $$
et choisissons les fonctions poids comme étant les fonctions de forme linéaires :

$$ 
\begin{aligned}
W_1 &= \lambda_n \\
W_2 &= \lambda_{n+1}
\end{aligned}
\qquad [8]
$$

Le système à résoudre est alors le même que celui obtenu avec une approximation P3-P1 TDGM [HUL 92, EKE 02, LI 03] ou avec une formulation en vitesse comme dans [MIC 03] :

$$ 
\begin{gathered}
u_{n+1} = u_n + \frac{\Delta t}{2} (v_{n+1}^{-} + v_n^{+}) \\
\left[ \begin{array}{cc} m + \frac{\Delta t^2}{12} k & -\frac{\Delta t^2}{12} k \\ \frac{\Delta t^2}{3} k & m + \frac{\Delta t^2}{6} k \end{array} \right] \left[ \begin{array}{c} v_n^{+} \\ v_{n+1}^{-} \end{array} \right] = \left[ \begin{array}{c} m v_n^{-} + \frac{\Delta t^2}{6} (f_n - f_{n+1}) \\ m v_n^{-} - \Delta t k u_n + \frac{\Delta t^2}{2} (f_n + f_{n+1}) \end{array} \right]
\end{gathered}
$$

### 3. Stabilité et précision

Dans les TDGM la notion de vecteur d'état n'a pas de sens et il n'est donc pas possible d'en étudier la stabilité en utilisant les méthodes classiques. Par contre, l'approche TX-FEM (lié par une changement de variable (7) aux méthodes TDGM) permet de définir le vecteur d'état $q_n$. Nous étudions donc la stabilité et la précision de la méthode TX-FEM avec les fonctions poids $W_1 = \lambda_n, W_2 = \lambda_{n+1}$. Les résultats sont donc valables pour ce choix particulier de $W_i$ et ils le sont par conséquent pour les TDGM mentionnées plus haut. Le rayon spectral $\rho$ de la matrice d'amplification du

<figure><img src="image_4.png" /><figcaption>**Figure 2. Rayon spectral**</figcaption></figure>

schéma est comparé à celui obtenu dans l'étude de la méthode de Newmark de l'accélération moyenne ($\gamma_N = \frac{1}{2}, \beta_N = \frac{1}{4}$) et de la méthode HHT [HIL 77] avec $\alpha = -0.3$. La figure 2 montre l'évolution du rayon spectral pour ces trois méthodes en fonction du paramètre fréquentiel $\omega\Delta t$ ($\omega = \sqrt{k/m}$). La première observation est que la méthode est inconditionnellement stable pour ce choix de fonctions poids. Contrairement au cas de la méthode de l'accélération moyenne de Newmark, le rayon n'est pas
toujours égal à 1 et la méthode présentera de l'amortissement numérique comme la méthode HHT. On distingue deux régions de $\omega\Delta t \approx \pi$ et $\omega\Delta t \ge 10$ pour lesquelles le comportement de la méthode diffère de celui de Newmark et HHT. En effet, dans ces zones on montre que la réponse du schéma est apériodique. Une conséquence de cette propriété est que la réponse du schéma sera exempte d'oscillations numériques pour les valeurs correspondantes du paramètre fréquentiel.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Erreurs en norme $L_2$</figcaption></figure>

De plus, contrairement à la méthode de Newmark dont l'ordre de convergence est inférieure ou égale à 2, on obtient ici une schéma du troisième ordre comme le montre la figure 3. Ce résultat concorde avec celui obtenu pour les TDGM.

## 4. Rupture d'un barreau élastique

On considère un barreau dont la longueur est 10m et la section $0.5m^2$. Il est constitué d'un matériau homogène de densité $\rho = 8000kgm^{-3}$ dont le comportement est considéré élastique linéaire, homogène et isotrope. Le paramètre caractérisant ce comportement est le module d'Young $E = 210GPa$. On soumet cette structure à une contrainte de traction $\sigma = 500MPa$ à une de ces extrémités, l'autre étant fixée. On impose la rupture du barreau en son milieu au temps $t/t_c = 2.5$ ($t_c$ valant la moitié du temps de parcours des ondes de compression dans le barreau). La figure 4 compare les résultats obtenus avec la méthode de Newmark et la méthode TX-FEM. Pour cette dernière la taille du pas de temps a été choisie de façon à ce que les résultats concernant le déplacement aient la même précision que ceux obtenus avec la méthode de Newmark. La taille de ce pas de temps est trois fois celle utilisée avec la méthode de Newmark. Les oscillations numérique du schéma de l'accélération moyenne sont grandement atténuées avec la méthode proposée ici. Cela donne à la méthode des Éléments Finis Etendus en temps une aptitude remarquable pour traiter les discontinuité en temps et ce, d'autant plus que ces discontinuités sont dues à l'apparition soudaine de discontinuité en espace.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Vitesse et accélération normalisées à l'interface créé</figcaption></figure>

**5. Bibliographie**

[BAB 97] BABUSKA I., MELENK J., « The Partition of unity method », *International Journal for Numerical Methods in Engineering*, vol. 40, 1997, p. 727-758.

[BEL 01] BELYTSCHKO T., MOÈS N., USUI S., PARIMI C., « Arbitrary discontinuities in finite elements », *International Journal for Numerical Methods in Engineering*, vol. 50, n° 4, 2001, p. 993-1013.

[EKE 02] EKEVID T., WIBERG N., « A comparison of parallel implementation of explicit DG and central difference method », *Communications in Numerical Methods in Engineering*, vol. 18, 2002, p. 585-597.

[HIL 77] HILBER H., HUGHES T., TAYLOR R., « Improved numerical dissipation for the time integration algorithms in structural dynamics », *Earthquake Engineering and Structural Dynamics*, vol. 5, 1977, p. 283-292.

[HUL 92] HULBERT G., « Discontinuity-capturing operators for elastodynamics », *Computer Methods in Applied Mechanics and Engineering*, vol. 96, 1992, p. 409-426.

[LI 03] LI X., YAO D., LEWIS R., « A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non-linear solids and saturated porous media », *International Journal for Numerical Methods in Engineering*, vol. 57, 2003, p. 1775-1800.

[MIC 03] MICHLER C., VAN BRUMMELEN E., HULSHOFF S., DE BORST R., « The relevance of conservation for stability and accuracy of numerical methods for fluid-structure interaction », *Computer Methods in Applied Mechanics and Engineering*, vol. 192, 2003, p. 4195-4215.

[NEW 59] NEWMARK N., « A method of computation for structural dynamics », *Proc. A.S.C.E.*, vol. 85, 1959, p. 67-94.

[ZIE 77] ZIENKIEWICZ O., « A new look at Newmark, Houbolt and other time stepping formulae. A weighted residual approach », *Earthquake Engineering and Structural Dynamics*, vol. 5, 1977, p. 413-418.