# Calcul de spectre sous contraintes. Cas de l'opérateur grad(div)

Majdi Azaiez, Etienne Ahusborde

► To cite this version:

Majdi Azaiez, Etienne Ahusborde. Calcul de spectre sous contraintes. Cas de l'opérateur grad(div). CFM 2009 - 19e Congrès Français de Mécanique, Aug 2009, Marseille, France. hal-03378317

HAL Id: hal-03378317

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

Submitted on 14 Oct 2021

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
# Calcul de spectre sous contrainte : Cas de l'opérateur grad(div)

ETIENNE AHUSBORDEª ET MEJDI AZAÏEZª

a. Institut Polytechnique de Bordeaux. Laboratoire TREFLE UMFR 8508.

## Résumé :

Le calcul de spectre d'opérateurs sous contrainte est connu pour être difficile problème numérique difficile. La présence d'une ou de plusieurs contraintes conduit à la mise en place de schémas spécifique. Selon la nature de la contrainte on a recourt à des techniques de stabilisation appropriées. Le but de cette contribution et de présenter une technique d'approximation du type hp ayant l'avantage de s'adapter à différents type de contraintes. Nous décrivons l'approche et nous apportons quelques résultats numériques.

## Abstract :

This paper describes a new hp method for the numerical solution of the two-dimensional -grad(div) eigenvalue problem in primal variational formulation. In standard methods with triangular or quadrangular finite elements, or with spectral elements, the spectrum contains spurious modes if the mesh is non-Cartesian. With the new hp method described in this paper we show that, although it delivers a $\mathbb{P}_p - \mathbb{P}_p$ approximation on the same grid for both velocity components, the spectrum is represented without any spurious mode whether a Cartesian or a non-Cartesian quadrangular grid is chosen. Spectra computed with standard high-order methods and with the new element are presented and compared with the analytic solutions.

**Mots clefs :** modes parasites, modes poluants, valeurs propres, stable, grad(div), rot rot

## 1 Introduction

Depuis leur introduction en mécanique des fluides numérique, les méthodes des éléments spectraux ont été utilisées avec succès dans de nombreuses applications. A travers les années, les communautés des mécaniciens et des mathématiciens ont apporté de grandes contributions dans le domaine de l'approximation d'opérateurs divers et variés. On peut trouver dans les ouvrages parus depuis quelques décennies (voir [1], [2], [3], [4], [5]) des idées de schémas analysés et numériquement validés pour approcher des opérateurs tels que le Laplacien, Stokes, Darcy et Maxwell.

A notre connaissance de la littérature, l'opérateur **grad**(div) souffre d'une absence notoire d'études dans le cadre des méthodes spectrales et des éléments spectraux, et ce malgré sa présence dans de multiples modèles mathématiques. Le but de cet article est d'apporter une contribution dans ce domaine.

Plus exactement, nous allons proposer une technique d'approximation pour calculer les valeurs propres et les vecteurs propres de cet opérateur.

### 1.1 État de l'art

L'opérateur **grad**(div) a déjà fait l'objet de nombreuses investigations numériques. Plusieurs travaux ont déjà été effectués notamment en éléments finis afin de trouver des éléments stables pour résoudre le problème aux valeurs propres : Trouver $\lambda^2 \in \mathbb{R}^+$ et **u** solutions de :

$$ -\mathrm{grad}(\mathrm{div} \, \mathbf{u}) = \lambda^2 \mathbf{u}, \quad \text{dans } \Omega, \quad (1) $$

$$ \mathbf{u} \cdot \mathbf{n} = 0, \quad \text{sur } \partial\Omega, \quad (2) $$

où $\mathbb{R}^+$ désigne l'ensemble des réels positifs ou nuls.

En plus de la formulation primale exprimée par (1)-(2), ce problème aux valeurs propres peut être exprimé sous une formulation mixte qui s'écrit : Trouver $\lambda^2 \in \mathbb{R}_+^*$, **u** et $\varphi$ tels que :

$$ \mathbf{u} - \mathrm{grad} \varphi = 0, \quad \text{dans } \Omega, $$

$$ \mathrm{div} \, \mathbf{u} = -\lambda^2 \varphi, \quad \text{dans } \Omega, $$

$$ \mathbf{u} \cdot \mathbf{n} = 0, \quad \text{sur } \partial\Omega, $$
où $\mathbb{R}_*^+$ désigne l'ensemble des réels strictement positifs.

Chacune de ces deux formulations est considérée et étudiée dans la littérature sous sa forme variationnelle. Ce qui consiste pour la primale à : Trouver $\lambda^2 \in \mathbb{R}^+$ et $\mathbf{u} \in H_0(\text{div}; \Omega)$ tels que :

$$ (\text{div } \mathbf{u}, \text{div } \mathbf{v}) = \lambda^2 (\mathbf{u}, \mathbf{v}), \quad \forall \mathbf{v} \in H_0(\text{div}; \Omega). $$

(..) désigne le produit scalaire dans $L^2(\Omega)$ ou bien dans $(L^2(\Omega))^2$. La formulation variationnelle mixte s'écrit : Trouver $\lambda^2 \in \mathbb{R}_*^+$, $\mathbf{u} \in H_0(\text{div}; \Omega)$ et $\varphi \in L_0^2(\Omega)$ tels que :

$$ (\mathbf{u}, \mathbf{v}) + b(\mathbf{v}, \varphi) = 0, \quad \forall \mathbf{v} \in H_0(\text{div}; \Omega), \quad (3) $$

$$ b(\mathbf{u}, q) = -\lambda^2(\varphi, q), \quad \forall q \in L_0^2(\Omega), \quad (4) $$

où $b(\mathbf{v}, q) = (\text{div } \mathbf{v}, q)$.

C'est cette dernière formulation qui a été la plus abordée dans la littérature pour traiter le problème (1)-(2).

La seconde méthode permettant de traiter la contrainte $\text{rot } \mathbf{u} = 0$ consiste à régulariser le problème en considérant une formulation de pénalisation. Soit $\alpha$ un nombre réel positif, la formulation pénalisée du problème (1)-(2) s'écrit : Trouver $\lambda^2 \in \mathbb{R}_*^+$ et $\mathbf{u}$ solutions de :

$$ -\text{grad}(\text{div } \mathbf{u}) + \alpha \text{rot}(\text{rot } \mathbf{u}) = \lambda^2 \mathbf{u}, \quad \text{dans } \Omega, \quad (5) $$

$$ \mathbf{u} \cdot \mathbf{n} = 0, \quad \text{sur } \partial\Omega. \quad (6) $$

Ce problème est alors associé à un opérateur elliptique et l'avantage de cette démarche est de pouvoir utiliser des éléments finis nodaux.

L'élément spectral que nous allons présenter a l'avantage d'offrir une discrétisation stable pour l'opérateur **grad**(div) dans des géométries "complexes" et être facilement utilisable pour d'autres types d'opérateurs tels que Stokes et rot rot( rot). L'élément clé est la définition d'une nouvelle famille de polynômes orthogonaux

## 2 Element stable pour grad(div)

Pour l'approximation par méthode spectrale nous allons introduire quelques outils et notations.

- $\mathbb{P}_p(\Lambda)$ est l'ensemble des polynômes de degré $\le p$ définis sur $\Lambda = ]-1, 1[$.

- Formule de quadrature de Gauss-Lobatto-Legendre (GLL) : $\xi_1 = -1$ et $\xi_{p+1} = 1$. Il existe un unique ensemble de $p-1$ noeuds $\xi_j$, $2 \le j \le p$, dans $]-1, 1[$ et $p+1$ réels, positifs (poids) $\rho_j$, $1 \le j \le p+1$ tels que :

$$ \forall \Phi \in \mathbb{P}_{2p-1}(\Lambda), \quad \int_{-1}^{1} \Phi(\xi) \, d\xi = \sum_{j=1}^{p+1} \Phi(\xi_j) \rho_j. $$

Les $(p+1)$ point $\xi_j$ sont les solutions de l'équation $(1-x^2)L'_p(x) = 0$ avec $L_p$ le polynôme de Legendre de degré $p$.

- Formule de quadrature de Gauss-Legendre (GL) : Il existe un unique ensemble de $p$ noeuds $\zeta_j$, $1 \le j \le p$, dans $]-1, 1[$ et $p$ réels, positifs (poids) $\omega_j$, $1 \le j \le p$ tels que :

$$ \forall \Phi \in \mathbb{P}_{2p-1}(\Lambda), \quad \int_{-1}^{1} \Phi(\xi) \, d\xi = \sum_{j=1}^{p} \Phi(\zeta_j) \omega_j. $$

Les $p$ point $\zeta_j$ sont les solutions de l'équation $L_p(x) = 0$

Une base canonique pour $\mathbb{P}_p(\Lambda = ]-1, +1[$) construite sur les points de GLL est donnée par les $(p+1)$ polynômes de Lagrange $h_i(x) \in \mathbb{P}_p(\Lambda)$ : Pour tout $-1 \le x \le +1$ et $1 \le i \le p+1$,

$$ h_i(x) = - \frac{1}{p(p+1)} \frac{1}{L_p(\xi_i)} \frac{(1-x^2) L'_p(x)}{(x-\xi_i)}, \quad (7) $$

Ces polynômes de Lagrange vérifient : $h_i(\xi_j) = \delta_{ij}$, $1 \le i, j \le p+1$.

On introduit une famille originale de polynômes $g_i$
## 2.1 Une nouvelle base pour $\mathbb{P}_{p-1}(\Lambda)$

On introduit les $p+1$ polynômes $g_i$ de degré $p-1$ définis par :

$$ \int_{-1}^{+1} (g_i(x) - h_i(x)) x^j dx = 0, \quad 1 \le i \le p+1, \ 0 \le j \le (p-1). $$

Cette relation montre que les différences $(g_i - h_i)_{i=1,...,p+1}$ sont proportionnelles à $L_p$. Plus exactement, on peut montrer que ces nouveaux polynômes sont donnés par la relation :

$$ \begin{aligned} g_i(x) &= h_i(x) - \beta_i L_p(x), \\ \beta_i &= \frac{(h_i, L_p)}{(L_p, L_p)} = \frac{1}{(p+1) L_p(\xi_i)}, \quad 1 \le i \le p+1. \end{aligned} $$

Ils vérifient : $\forall i = 1, \dots, p+1$

$$ g_i(\zeta_j) = h_i(\zeta_j), \quad 1 \le j \le p, \tag{8} $$

$$ g'_i(\xi_j) = h'_i(\xi_j), \quad 2 \le j \le p. $$

## 2.2 Approximation polynomiale

On considère pour toute fonction donnée $f$ d'une part le polynôme d'interpolation de Lagrange de degré $p$ défini par :

$$ f_p(x) = \sum_{j=1}^{p+1} f(\xi_j) h_j(x) $$

et d'autre part le polynôme de degré $p-1$, $f_{p-1}^{(0)}$ :

$$ f_{p-1}^{(0)}(x) = \sum_{j=1}^{p} f_j^{(0)} g_j(x). $$

$f^{(0)}$ s'obtient à partir de $f_p$ en écrivant que :

$$ (f_p - f_{p-1}^{(0)}, g_i) = \int_{-1}^{1} (f_p - f_{p-1}^{(0)}) (x) g_i(x) dx. = 0, \quad \forall i = 1, \dots, p. \tag{9} $$

Les $p$ égalités (9) s'écrivent :

$$ \sum_{j=1}^{p} (g_j, g_i) f_j^{(0)} = \sum_{j=1}^{p} (h_j, g_i) f(\xi_j) + (h_{p+1}, g_i) f(\xi_{p+1}), \quad \forall i = 1, \dots, p. \tag{10} $$

**Remarque 1 :** La propriété (8) permet facilement de voir que :

$$ (g_j, g_i) = (h_j, g_i) = (g_j, h_i) = (h_j, h_i)_{GL}. $$

((, )$_{GL}$ désigne le produit scalaire discret construit grâce à la formule de Gauss-Legendre). Ce qui veut dire que nous pouvons remplacer (10) par :

$$ \sum_{j=1}^{p} (g_j, g_i) f_j^{(0)} = \sum_{j=1}^{p} (g_j, g_i) f(\xi_j) + (g_{p+1}, g_i) f(\xi_{p+1}), \quad \forall i = 1, \dots, p. \tag{11} $$

Pour la suite on pose $\mathbf{F}^{(0)}$ le vecteur de taille $p$ contenant les $(f_1^{(0)}, \dots, f_p^{(0)})$ et $\mathbf{F}$ le vecteur de taille $p+1$ contenant les $(f(\xi_1), \dots, f(\xi_{p+1}))$. On introduit les matrices $\mathbf{G}$ et $\mathbf{H}$ de taille respective $(p, p)$ et $(p, p+1)$ définies par :

$$ \mathbf{G}_{i,j} = (g_j, g_i), \quad 1 \le i, j \le p, \tag{12} $$

$$ \mathbf{H}_{i,j} = (g_j, g_i), \quad 1 \le i \le p, \ 1 \le j \le p+1. \tag{13} $$
Sous forme matricielle, (11) devient :

$$ \mathbf{GF}^{(0)} = \mathbf{HF}. \qquad (14) $$

On peut montrer facilement que $\mathbf{G}$ est inversible et que $\mathbf{G}^{-1}\mathbf{H}$ possède la forme suivante (p=4) :

$$ \begin{pmatrix} 1 & 0 & 0 & 0 & \alpha_1 \\ 0 & 1 & 0 & 0 & \alpha_2 \\ 0 & 0 & 1 & 0 & \alpha_3 \\ 0 & 0 & 0 & 1 & \alpha_4 \end{pmatrix} $$

Ce qui nous permet de conclure que :

$$ \forall i = 1, \dots, p \quad f_i^{(0)} = f(\xi_i) + \alpha_i f(\xi_{p+1}). \qquad (15) $$

**Remarque 2** : Si la fonction $f$ est nulle au bord, (cas des problème du type Dirichlet), alors les relations (15) montrent que :

$$ \forall i = 1, \dots, p \quad f_i^{(0)} = f(\xi_i) \quad \text{et} \quad f_{p-1}^{(0)}(x) = \sum_{j=1}^{p} f(\xi_j)g_j(x). \qquad (16) $$

Quand les conditions aux limites sont du type Neumann ou périodiques alors on se contente de la relation (15). Les relations (10) conjuguées avec (15) montrent que les coefficient $\alpha_i$ sont solutions du système :

$$ \forall i = 1, \dots, p \quad \sum_{j=1}^{p} (g_j, g_i) \alpha_j = (g_{p+1}, g_i). \qquad (17) $$

Les coefficients $\alpha_i$ sont issus de la résolution du système linéaire (17) et ne dépendent que de la base et peuvent
dont être calculés en prétraitement.

## 2.3 Cas multidimensionnel

L'extension au cas bidimensionnel se fait de manière naturelle. Le domaine $\Omega$ est le carré $]-1, +1]^2$ et soit une fonction $f$ continue sur $\bar{\Omega}$. Celle-ci peut être approchée par le biais de quatre polynômes de natures différentes

$$ \begin{align*} f_p(r, s) &= \sum_{k=1}^{p+1} \sum_{\ell=1}^{p+1} f_{k\ell} h_k(r) h_\ell(s), \quad f_p \in \mathbb{P}_p \otimes \mathbb{P}_p, \\ f_p^{(0)}(r, s) &= \sum_{k=1}^{p} \sum_{\ell=1}^{p} f_{k\ell}^{(0)} g_k(r) g_\ell(s), \quad f_p^{(0)} \in \mathbb{P}_{p-1} \otimes \mathbb{P}_{p-1} \\ f_p^{(1)}(r, s) &= \sum_{k=1}^{p+1} \sum_{\ell=1}^{p} f_{k\ell}^{(1)} h_k(r) g_\ell(s), \quad f_p^{(1)} \in \mathbb{P}_p \otimes \mathbb{P}_{p-1}, \\ f_p^{(2)}(r, s) &= \sum_{k=1}^{p} \sum_{\ell=1}^{p+1} f_{k\ell}^{(2)} g_k(r) h_\ell(s), \quad f_p^{(2)} \in \mathbb{P}_{p-1} \otimes \mathbb{P}_p \end{align*} $$

En reprenant les notations et les calculs menés dans le cas unidimensionnel, on peut montrer que :

- Pour $f^{(0)} \forall i, j = 1, \dots, p$

$$ f_{ij}^{(0)} = f(\xi_i, \xi_j) + \alpha_i f(\xi_{p+1}, \xi_j) + \alpha_j f(\xi_i, \xi_{p+1}) + \alpha_i \alpha_j f(\xi_{p+1}, \xi_{p+1}). \qquad (18) $$

- Pour $f^{(1)} \forall i = 1, \dots, p+1 \text{ et } \forall j = 1, \dots, p$

$$ f_{ij}^{(1)} = f(\xi_i, \xi_j) + \alpha_j f(\xi_i, \xi_{p+1}) \qquad (19) $$

- Pour $f^{(2)} \forall i = 1, \dots, p \text{ et } \forall j = 1, \dots, p+1$

$$ f_{ij}^{(2)} = f(\xi_i, \xi_j) + \alpha_i f(\xi_{p+1}, \xi_j) \qquad (20) $$

**Remarque 3** : De nouveau, on remarque que si la fonction, en question est nulle au bord, les trois polynômes bien que complètement différents ont les mêmes coefficients. On peut selon la régularité voulue choisir l'une ou l'autre des fonctions $f^{(j)}$
## 3 Spectre d'opérateurs sous contrainte

Pour illustrer l'intérêt de cet élément nous allons considérer pour $\Omega$ le carré de référence tourné d'un angle $\pi/4$ et divisé en 16 éléments Transformé sur l'élément de référence nous sommes amenés à approcher le problème

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

suivant : Trouver $\mathbf{u}_p$ et $\lambda^2 \in \mathbb{R}^+$ tels que $\forall \mathbf{v}_p$ :

$$ 
\begin{aligned}
& \int_{\hat{\Omega}} \left( \frac{\partial u_{px}}{\partial r} - \frac{\partial u_{px}}{\partial s} + \frac{\partial u_{py}}{\partial r} + \frac{\partial u_{py}}{\partial s} \right) \left( \frac{\partial v_{px}}{\partial r} - \frac{\partial v_{px}}{\partial s} + \frac{\partial v_{py}}{\partial r} + \frac{\partial v_{py}}{\partial s} \right) dr ds \\
& = \frac{\lambda^2 h^2}{2} \int_{\hat{\Omega}} (u_{px}v_{px} + u_{py}v_{py}) dr ds.
\end{aligned}
$$

Un choix naturel consiste à calculer $\mathbf{u}_p$ sur chaque quadrangle comme un polynôme de degré inférieur ou égal à $p$ pour les deux composantes et dans les deux directions. Le figure 3 illustre l'instabilité de cet élément. On constate sur la partie gauche l'apparition de valeurs propres non nulles. Ce phénomène est connu sous le nom de *pollution*.

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

On constate que le terme de gauche fait intervenir les dérivés des deux composantes et dans les deux directions, ce qui a pour conséquence de mettre en échec l'élément fini de Raviart-Thomas [6]. Pour contourner cette difficulté et éliminer la pollution on propose la version suivante :

$$ 
\begin{aligned}
& \int_{\hat{\Omega}} \left( \frac{\partial u_{px}^{(1)}}{\partial r} - \frac{\partial u_{px}^{(2)}}{\partial s} + \frac{\partial u_{py}^{(1)}}{\partial r} + \frac{\partial u_{py}^{(2)}}{\partial s} \right) \left( \frac{\partial v_{px}^{(1)}}{\partial r} - \frac{\partial v_{px}^{(2)}}{\partial s} + \frac{\partial v_{py}^{(1)}}{\partial r} + \frac{\partial v_{py}^{(2)}}{\partial s} \right) dr ds \\
& = \frac{\lambda^2 h^2}{2} \int_{\hat{\Omega}} \left( u_{px}^{(0)} v_{px}^{(0)} + u_{py}^{(0)} v_{py}^{(0)} \right) dr ds.
\end{aligned}
$$

Une fois la matrice issue de la discrétisation est assemblée, elle est diagonalisée grâce à une procédure de la bibliothèque Lapack¹

### 3.1 Résultats numériques

Le but de cette sous section est de prouver sur un exemple, le caractère stabilisant et convergent de notre approche. Pour cela on choisit de mesurer l'approximation de la valeur propre $\lambda_{55}$ en sachant que les valeurs propres non nulles sont données par $\lambda_{k,l}^2 = \pi^2/4(k^2 + l^2)$.
<figure><img src="image_4.png" /><figcaption>FIG. 1 – Erreur relative sur la valeur propre $\lambda_{5,5}^2$ en fonction de p.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 2 – Erreur relative sur la valeur propre $\lambda_{5,5}^2$ en fonction de N.</figcaption></figure>

N désigne le nombre d'éléments dans chaque direction et p le degré polynomiale par élément. La figure 1 montre que pour N=1, 2, 4 et 8, la convergence est exponentielle en p.

La figure 2 montre que pour p, degré polynômial, fixé, la convergence est algébrique en N.

On constate qu'en plus de la disparition de la pollution, la convergence est du type hp. Cette conclusion est confirmée par d'autres expériences numériques

Références

[1] Deville M., Fischer P., and Mund E. High-Order Methods for Incompressible Fluid Flow . Cambridge University Press, Cambridge, 2002.

[2] Canuto C., Hussaini M. Y., Quarteroni A., and Zang T. A. Spectral Methods in Fluid Dynamics. Springer-Verlag, Berlin, 1988.

[3] Bernardi C. and Maday Y. Spectral Methods in Handbook of Numerical Analysis . P.G. Ciarlet and J.L. Lions, North-Holland, 1997.

[4] Peyret R. Spectral Methods for Incompressible Viscous Flow. Springer, Berlin, 2002.

[5] Gottlieb D. and Orszag S. Numerical Analysis of Spectral Methods : Theory and Applications. CBMS-NSF Regional Conference Series in Applied Mathematics, 1977.

[6] Raviart P. and Thomas J. A mixed finite element method for second order elliptic problems. Mathematical Aspects of the Finite Element Methods, 606 of Lecture Notes in Mathematics, 292–315, 1977.

<footnote n="1">¹Linear Algebra Package</footnote>