# Stabilisation interne optimale de l'équation des ondes par une méthode de lignes de niveau

Arnaud Münch

► To cite this version:

Arnaud Münch. Stabilisation interne optimale de l'équation des ondes par une méthode de lignes de niveau. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813075

HAL Id: hal-01813075

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

Submitted on 12 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_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Stabilisation interne optimale de l'équation des ondes par une méthode de level set

## Arnaud Münch

Laboratoire de Mathématique, UMR CNRS 6623
Université de Franche-Comté,
16, route de Gray 25030 - Besançon
arnaud.munch@math.univ-fcomte.fr

**RÉSUMÉ.** On considère une équation des ondes linéaire définie sur un domaine $Ω$ régulier du plan, et amortie sur un sous-domaine interne $ω ⊂ Ω$. On considère le problème de la position et de la forme optimale de $ω$ minimisant l'énergie du système à un instant $T > 0$. La méthode de dérivation de forme conduit à la variation de l'énergie vis-à-vis de $ω$ exprimée comme une intégrale curviligne le long de $∂ω$. La méthode des "Level Set" ramène alors le problème à la résolution d'une équation non linéaire d'Hamilton-Jacobi dont le terme d'advection est l'intégrand de la dérivée de forme. L'efficacité de la méthode est numériquement confirmée.

**ABSTRACT.** We consider a damped linear wave equation defined on a bi-dimensional domain $Ω$. Due to the damping term defined on the sub-domain $ω ⊂ Ω$, the system is dissipative. We address the problem of the optimal position and shape design of the support $ω$ in order to minimize the energy of the system at a given time $T$. Introducing an adjoint problem we first obtain explicitly the (shape) derivative of the energy at time $T$ with respect to the variation of $ω$. Expressed as a boundary integral on $∂ω$, this derivative is then used as an advection velocity in an Hamilton-Jacobi equation for changing the shape. We use the level-set methodology on a fixed working Eulerian mesh. We also consider the optimization with respect to the value of the damping. The numerical approximation is presented in detail and leads to several numerical experiments that indicate the efficiency of the method.

**MOTS-CLÉS:** Optimisation de forme, équation des ondes, level set, schemas aux differences finis hyperbolique.

**KEYWORDS:** Optimal shape design, wave equation, level set, hyperbolic schemes.
# 1. Introduction

On considère le problème évolutif suivant défini sur le cylindre $\Omega \times (0, T)$, $\Omega$ étant un domaine borné de classe $C^2(\mathbb{R}^2)$ :

$$ egin{cases} y''_{\omega,a}(\mathbf{x},t) - \Delta y_{\omega,a}(\mathbf{x},t) + a(\mathbf{x})y'_{\omega,a}(\mathbf{x},t) = 0 & \text{in } \Omega \times (0,T), \\ y_{\omega,a}(\mathbf{x},t) = 0 & \text{on } \partial\Omega \times (0,T), \\ y_{\omega,a}(\mathbf{x},0) = y_0(\mathbf{x}), \quad y'_{\omega,a}(\mathbf{x},0) = y_1(\mathbf{x}) & \text{in } \Omega. \end{cases} \quad (1) $$

Le symbole ' désigne la dérivée partielle par rapport au temps. On suppose qu'il existe une constante non négative $a$ telle que le potentiel dissipatif $a(\mathbf{x}) \in L^{\infty}(\Omega, \mathbb{R}^+)$ vérifie

$$ a(\mathbf{x}) \geq a \quad \text{a.e. } \mathbf{x} \in \omega, \quad \omega \subset \Omega, \omega \neq \emptyset. \quad (2) $$

Le sous-domaine $\omega$ indépendent du temps peut être composé d'un nombre fini de composantes disjointes. Lorsque les données initiales ($y_0, y_1$) indépendantes de $\omega$ et $a$ soient dans $H^2(\Omega) \cap H_0^1(\Omega)$, le système (1) est bien posé. De plus, l'énergie du système (1) donnée par

$$ E(\omega, a, t) = \frac{1}{2} \int_{\Omega} \{ |y'_{\omega,a}(\mathbf{x},t)|^2 + |\nabla y_{\omega,a}(\mathbf{x},t)|^2 \} dx, \quad \forall t \ge 0, \quad (3) $$

vérifie la loi de dissipation

$$ E'(\omega, a, t) = - \int_{\Omega} a(\mathbf{x}) |y'_{\omega,a}(\mathbf{x}, t)|^2 dx \le 0, \quad \forall t \ge 0. \quad (4) $$

Sous certaines conditions dites "d'optiques géométriques" sur $\omega$, le système (1) modélise la stabilisation d'une membrane par un actuateur interne. $y$ désigne le déplacement transversal de la membrane tandis que $y_0$ et $y_1$ désigne la position et la vitesse initiale. Précisément, lorsque $\omega$ est dans un voisinage de $\Gamma_0 = \{\mathbf{x} \in \partial\Omega; \mathbf{x} \cdot \mathbf{v} > 0\}$ où $\mathbf{v}$ désigne le vecteur normal unité extérieur à $\Omega$, la décroissance de $E$ est exponentielle en temps. Que $\omega$ satisfasse ou non une propriété d'optique géométrique, la question du placement et de la forme optimal de l'actuateur afin de minimiser l'énergie au temps $T$ apparaît. Sans restrictions sur $\omega$, la solution triviale est $\omega = \Omega$. On se place alors sur le sous-ensemble $V_\beta = \{\omega \subset \Omega, \text{aire}(\omega) = \beta \text{aire}(\Omega)\}$, $\beta \in (0, 1)$. Instons sur le fait que l'énergie au temps $T$ - appelée fonction cout dans la suite - n'est pas monotone par rapport à l'aire de $\omega$. On peut au plus conjecturer que l'inclusion $\omega_1 \subset \omega_2 \subset \Omega$ entraîne l'inégalité $E(\omega_2, a, T) \le E(\omega_1, a, T)$ pour $T$ suffisamment grand. Par ailleurs, le domaine $\omega$ étant fixé, on peut s'intéresser également à la minimization de la fonction coût vis-à-vis de $a$. A nouveau la fonction coût n'est pas monotone par rapport à $a$. Lorsque que $a$ tend vers l'infini, la solution $y_{\omega,a}$ tend vers une constante tandis $E(\omega, a, T)$ tend vers $E(\omega, a, 0) = E(\omega, 0, T)$ par valeur inférieure si $y_1 \equiv 0$ dans $\omega$. C'est le phénomène d'over-damping où $a$ se conduit comme un paramètre de pénalisation. Il en résulte qu'aucune restriction sur $a$ est nécessaire et les problèmes considérés sont

$$ (P_\omega) : \inf_{\omega \in V_\beta} E(\omega, a, T) \quad ; \quad (P_a) : \inf_{a \in L^\infty(\Omega, \mathbb{R}^+)} E(\omega, a, T). \quad (5) $$
Le problème ($\mathcal{P}_\omega$) est illustré sur la figure 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Optimisation de la forme de $\omega$, support de $a$ afin de minimiser $E(\omega, a, T)$.</figcaption></figure>

## 2. Dérivée de l'énergie E vis-à-vis de $\omega$ et de $a$ - Méthode de level set

Pour simplifier, on note simplement $y$ pour $y_{\omega,a}$ et on prend sans perte de généralité $a(\mathbf{x}) = a \mathbb{1}_{\{\mathbf{x} \in \omega\}}$. Afin de relaxer la contrainte $\omega \in V_\beta$, on considère la fonctionnelle penalisée

$$E^{\epsilon}(\omega, a, T) = E(\omega, a, T) + \frac{1}{2} \epsilon^{-1} (\text{aire}(\omega) - \beta \text{aire}(\Omega))^2. \quad (6)$$

La résolution des problèmes par une méthode de gradient recquiert la connaissance des dérivées de E par rapport à $\omega$ et $a$. On introduit pour cela le champs $\theta \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ avec $\theta|_{\partial\Omega} = 0$ non nul dans un voisinage.

**Théorème 1** La dérivée de $E^\epsilon$ par rapport à $\omega$ est donnée par :

$$\frac{\partial E^{\epsilon}(\omega, a, T)}{\partial \omega} \cdot \theta = \int_{\partial \omega} \left[ \underbrace{\epsilon^{-1}(\text{aire}(\omega) - \beta \text{aire}(\Omega)) + a \int_{0}^{T} y'(\mathbf{x}, t) p(\mathbf{x}, t) dt}_{j^{\epsilon}(y, p, a)} \right] \theta \cdot \nu d\sigma \quad (7)$$

où $\nu$ désigne la normale extérieure à $\omega$ et $p$ la solution du problème adjoint suivant :

$$\left\{
\begin{array}{ll}
p''(\mathbf{x}, t) - \Delta p(\mathbf{x}, t) - a(\mathbf{x})p'(\mathbf{x}, t) = 0 & \text{dans } \Omega \times (0, T), \\
p(\mathbf{x}, t) = 0 & \text{sur } (\partial\Omega \setminus \partial\omega) \times (0, T), \\
p(\mathbf{x}, T) = -y'(\mathbf{x}, T) & \text{in } \Omega, \\
p'(\mathbf{x}, T) = -a(\mathbf{x})y'(\mathbf{x}, T) - \Delta y(\mathbf{x}, T) & \text{in } \Omega.
\end{array}
\right.
\quad (8)$$

<figure><img src="image_3.png" /></figure>
**Théorème 2** La dérivée de $E(\omega, a, T)$ vis-à-vis de $a$ est donnée par :

$$ \frac{\partial E(\omega, a, T)}{\partial a} \cdot a^1 = \int_0^T \int_\omega^T a^1 y'(x, t) p(x, t) dx dt \quad (9) $$

où $p$ est solution de (8).

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

Compte tenu de ces résultats, l'algorithmé (éventuellement couplé) de descente prend
la forme suivante, étant donné une initialisation $(a_0, \omega_0)$ :

$$ \begin{pmatrix} \omega_{k+1} \\ a_{k+1} \end{pmatrix} = \begin{pmatrix} \omega_k \\ a_k \end{pmatrix} - \begin{pmatrix} \eta_1 \nu_k \\ \eta_2 \end{pmatrix} \cdot \left( \varepsilon^{-1}(\text{aire}(\omega_k) - \beta \text{aire}(\Omega)) + \int_0^T a_k y_k' p_k dt \right) \int_0^T \int_{\omega_k} y_k' p_k dx dt \quad (10) $$

$\eta = (\eta_1, \eta_2)\mathbb{R}^2$ suffisamment petit de façon à assurer la décroissance de la fonctionnelle. Afin de conserver un maillage fixe au cours des iterations, on découple la description de la solution $y$ de celle de $\omega$ en introduisant une fonction $\psi$ telle que $\partial\omega = \{x, \psi(x) = 0\}$ solution du système hyperbolique non linéaire suivant :

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial \psi}{\partial \tau} - j^{\varepsilon}(y_{\omega,a}, p_{\omega,a}, T) |\nabla \psi| = 0 & \text{in } \Omega \times (0, \infty), \\ \\ \psi(\cdot, \tau = 0) = \psi_0 & \text{in } \Omega, \\ \\ \psi = \psi_0 > 0, \quad \nabla \psi \cdot \nu = 0 & \text{on } \partial \Omega \times (0, \infty). \end{array} \right. \quad (11) $$

**3. Schéma numérique de résolution des systèmes hyperboliques en $y$, $p$ et $\psi$**

Les systèmes dont $y$, $p$ et $\psi$ sont solutions sont de nature hyperbolique. Les schémas numériques différences finis ou élément finis usuels n’assurent pas en général la convergence de la solution $y_h$ ou $p_h$ dans la norme de l’énergie. On montre en particulier que si les conditions initiales sont discontinues, le taux décroissance exponentiel associé à l’énergie discrète n’est pas minoré uniformément. Cela est dû à la dispersion numérique pour les hautes fréquences des schémas. Parmi les rémèdes permettant de restaurer la convergence uniforme de la solution, signalons celle qui consiste à rajouter un terme de régularisation à l’équation d’état :

$$ y'' - \Delta y + a(x)y' - h^2 \Delta y' = 0 \quad (12) $$

facteur de $h^2 - h$ désignant le paramètre de discrétisation en espace. Cette technique
analogue à une régularisation de Tychonov permet de redresser le spectre du schéma
pour les hautes fréquences. Une discrétisation usuelle centrée de l'équation précédente
permet de restaurer la stabilisation numérique uniforme sous la condition de stabilité
$\Delta t \le h/\sqrt{2}$ - $\Delta t$ désignant le paramètre de discrétisation en temps. La solution $p$ du
problème adjoint étant moins régulière que $y$, cette modification est à plus forte raison
nécessaire. Enfin, une alternative à cette technique pour résoudre le système en $\psi$ (où
une pathologie numérique apparaît également) est d'employer un schéma "upwind"
reposant sur la technique des caractéristiques.
## 4. Une application numérique

Sur le carré unité $\Omega = (0, 1)^2$, on considère les conditions initiales suivantes :

$$y_0(x) = 100 \sin(2\pi x_1) \sin(\pi x_2) \quad (13)$$

et $a(x) = 101_{(x \in \omega)}$, $\beta = 0.1$, $\varepsilon = 3.E + 5$, $T = 1$. Dans le cas conservatif, la solution est $y(x, t) = 100 \cos(\sqrt{5\pi t})y_0(x)$ de sorte que $(y'(x, t))^2$ admet pour tout temps $t$ deux maxima aux points $(1/4, 1/2)$ et $(3/4, 1/2)$. De la relation $E'(t) = -10 \int_{\omega} (y'_{\omega,a})^2 dx$, il est alors raisonable de conjecturer que le $\omega$ optimal est composé de deux parties centrées aux points de maximum. Les essais numériques confirment cette prédiction. La figure 2 donne l'évolution de $\{x \in \Omega, \psi_k(x) = 0\}$ en fonction de l'itération $k$ (ou de façon équivalente, en fonction du pseudo-temps $\tau$). Le domaine $\omega_0$ initial retenu est le disque de centre $(0.35, 0.35)$ et de rayon $\sqrt{0.1/\pi}$. Les différentes simulations numériques effectuées indiquent que la conjecture précédente est en général fausse et dépend de la valeur du coefficient d'amortissement. Conformément à la littérature récente, la méthode bien que permettant une réelle diminution de la fonction coût, ne conduit pas à des résultats intuitifs. Une analyse spectrale du système montre que cela est dù au comportement non linéaire de l'énergie vis-à-vis de $a$.

## 5. Bibliographie

* Allaire G., Jouve F., Toader A.M., *Structural optimization using sensitivity analysis and level-set methods*, J. Comp. Phys. **194**(1), 363-393 (2004).

* Banks H.T., Ito K., Wang B., *Exponentially stable approximations of weakly damped wave equations*, Ser. Num. Math. **100**, 1-33, Birkhäuser (1990).

* Bardos C., Lebeau G., Rauch J., *Sharp sufficient conditions for the observation, control and stabilization from the boundary*, SIAM J. Control and Opt., **30**, 1024-1065 (1992).

* Barles G., *Viscosity solutions of Hamilton-Jacobi equations*, Mathematics and Applications, Springer-Verlag, Paris, (1994).

* Burger M., Osher S.J., *A survey on level set methods for inverse problems and optimal design*, Preprint, (2004).

* Henrot A., Maillot H., *Optimization of the shape and the location of the actuators in an internal control problem*, Boll. Unione Mat. Ital. Sez. B Artic. Ric. Mat., **3**, 737-757 (2001).

* Münch A., *Family of implicit schemes uniformly controllable for the 1-D wave equation*, C.R.Acad. Sci. Paris, Série I **339**(10), 733-738, (2004).

* Münch A., *A uniformly controllable and implicit scheme for the 1-D wave equation*, To appear in Mathematical modelling and numerical analysis (2005).

* Münch A., *Optimal internal stabilization of a damped wave by a level set approach*, To appear in *Interface and Free boundaries* (2005).

* Münch A., Pazoto A., *Uniform stabilization of a numerical approximation of the locally damped wave equation*, To appear in *Control, Optimization and Calculus of Variation*, (2005).

* Münch A., Ousset Y., *Numerical simulation of delamination growth in curved interfaces*, *Comput. Methods Appl. Mech. Engrg* **191**, 2045-2067, (2002).
6 Nom de la revue ou conférence (à définir par \submitted ou \toappear)

<figure><img src="image_5.png" /><figcaption>**Figure 2.** $(y_0, y_1) = (100 \sin(2\pi x_1) \sin(\pi x_2), 0)$, $T = 1$, $a = 10$. - *Evolution de {x ∈ Ω, ψ_k(x) = 0} par rapport à k - h = 1/151.*</figcaption></figure>