# Simulation temporelle d'un modèle d'impédance de liner en utilisant la représentation diffusive d'opérateurs

Florian Monteghetti, Denis Matignon, Estelle Piot, Lucas Pascal

► To cite this version:

Florian Monteghetti, Denis Matignon, Estelle Piot, Lucas Pascal. Simulation temporelle d'un modèle d'impédance de liner en utilisant la représentation diffusive d'opérateurs. 13e Congrès Français d'Acoustique (CFA 2016), Société Française d'Acoustique (SFA), Apr 2016, Le Mans, France. pp. 2549-255. hal-01342492

HAL Id: hal-01342492

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

Submitted on 6 Jul 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
CFA/VISHNO 2016

Simulation temporelle d'un modèle d'impédance de liner
en utilisant la représentation diffusive d'opérateurs

F. Monteghetti<sup>a</sup>, D. Matignon<sup>b</sup>, E. Piot<sup>a</sup> et L. Pascal<sup>a</sup>

<sup>a</sup>ONERA – The French Aerospace Lab, 2 Avenue Edouard Belin, 31055 Toulouse, France

<sup>b</sup>ISAE-SUPAERO, 4 Avenue Edouard Belin, 31400 Toulouse, France

florian.monteghetti@onera.fr

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

LE MANS
Une nouvelle condition aux limites temporelle de type impédance (TDIBC), adaptée à la simulation large-bande d'un liner, est présentée. Elle se base sur des modèles physiques existants, qui séparent impédance de cavité et de perforation. Pour la cavité, une approximation du nombre d'onde aboutit à un modèle EHR. Pour la perforation, une analyse temporelle théoriquement exacte du terme en dérivée fractionnaire mène à une formulation équivalente par équations différentielles, sans convolution récursive, qui se couple à tout schéma numérique dès le niveau semi-discret. Une étude numérique dans une méthode Galerkin Discontinu, qui résout les équations d'Euler linéarisées monodimensionnelles, met en évidence la neutralité de la TDIBC vis à vis du CFL.

# 1 Introduction

Afin de réduire le bruit de fan, les parois internes des turboréacteurs sont tapissées d'un matériau absorbant acoustique (liner). De par leur simplicité d'exploitation, les liners à un degré de liberté (SDOF) sont d'un emploi courant. L'objet de la présente étude est de développer une condition aux limites temporelle de type impédance (TDIBC) qui permet de réaliser une simulation temporelle large bande de la réponse acoustique d'un liner SDOF. Malgré les avancées considérables accomplies depuis le début des années 90, une telle simulation reste un défi. En effet, les modèles temporels disponibles ne représentent pas toute la physique en jeu : le terme fractionnaire $\sqrt{j\omega}$, qui traduit la diffusion visco-thermique dans la perforation, est systématiquement négligé (voire considéré comme mal posé [1]), et la modélisation de la cavité n'est que rarement large bande.

La démarche suivie pour la conception de cette TDIBC prend comme point de départ une analyse temporelle des modèles physiques existants (qui donnent déjà des résultats satisfaisants en fréquentiel). Une telle analyse permet de comprendre la structure temporelle du modèle, et d'en déduire une stratégie de discrétisation *ad hoc*.

Cet article présente un résumé de l'application de cette démarche à l'impédance de perforation. Le cas de l'impédance de cavité sera présenté lors de la conférence; ici, elle se réduit à un terme analogue à celui présent dans le modèle *Extended Helmholtz Resonator* (EHR) [1], qui est satisfaisant pour l'application présentée.

Le présent article est organisé comme suit : la section 2 présente le type de liner considéré et la modélisation retenue. Les deux sections suivantes se concentrent sur l'impédance de perforation du point de vue théorique (Sec. 3) et numérique (Sec. 4), où une nouvelle TDIBC est présentée. L'impédance de cavité est brièvement traitée Sec. 5. Enfin, la section 6 présente l'intégration de cette TDIBC dans une méthode Galerkin Discontinu (DG) pour les équations d'Euler linéarisées monodimensionnelles. L'impact de la TDIBC sur la limite de stabilité linéaire, d'une importance primordiale pour les applications pratiques, est étudiée.

# 2 Modélisation physique d'un liner SDOF

## 2.1 Généralités

Le principe des liners SDOF (voir Fig. 1) est basé sur celui du résonateur de Helmholtz : une plaque perforée, siège de pertes visco-thermiques, est superposée à un nid d'abeilles, qui fixe la fréquence de résonance. Lorsqu'ils sont considérés comme *localement réactifs*, leur caractérisation se fait naturellement en régime harmonique établi, par la

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Constitution d'un liner SDOF : superposition d'une plaque perforée et d'une cavité.</figcaption></figure>

<table><caption>TABLE 1 – Dimensions du liner SDOF considéré.</caption><thead><tr><th></th><th>σ</th><th>5 %</th><th>Porosité</th></tr></thead><tbody><tr><td rowspan="2">Plaque perforée</td><td>l</td><td>0.8 mm</td><td>Longueur</td></tr><tr><td>d</td><td>0.3 mm</td><td>Diamètre</td></tr><tr><td rowspan="2">Cavité</td><td>l<sub>c</sub></td><td>20 mm</td><td>Longueur</td></tr><tr><td>d<sub>c</sub></td><td>9 mm</td><td>Diamètre</td></tr></tbody></table>

définition d'une impédance de surface $\hat{z}_s$ liant pression $\hat{p}$ et vitesse acoustique normale (entrante) $\hat{u} \cdot \underline{n}$ par :

$$ \forall s \in \mathbb{C} : \Re[s] > 0, \quad \hat{p}(s) = \hat{z}_s(s) \hat{u}(s) \cdot \underline{n}. \quad (1) $$

Relation extensible aux cas présentant un écoulement de base et/ou de forts niveaux sonores. Le matériau est alors modélisé comme un filtre linéaire invariant en temps. Dans ce document, on utilise la convention $e^{+st} (e^{+j\omega t})$ à la synthèse pour la transformée de Laplace (Fourier).

## 2.2 Modèle de référence

Le matériau considéré ici est un liner micro-perforé, tel que représenté Fig. 1, dont les dimensions sont données Tab. 1. La longueur $l$ de la perforation étant très faible devant les longueurs d'ondes d'intérêt en acoustique, l'écoulement en son sein est supposé incompressible. Cette hypothèse autorise la décomposition de l'impédance de surface $\hat{z}_s$ en une impédance de perforation $\hat{z}_p$ et une impédance de cavité $\hat{z}_c$ selon

$$ \begin{aligned} \hat{z}_s &= \frac{\hat{p}}{\hat{u} \cdot \underline{n}} = \frac{\hat{p} - \hat{p}_{cav}}{\hat{u} \cdot \underline{n}} + \frac{\hat{p}_{cav}}{\hat{u} \cdot \underline{n}} \\ &= \frac{\hat{z}_p}{\sigma} + \hat{z}_c, \end{aligned} $$

où l'impédance de cavité s'exprime

$$ \hat{z}_c = z_0 \coth [jk_c(s) l_c], $$

en notant l'impédance caractéristique de l'air $z_0 = \rho_0 c_0$.

Un modèle physique satisfaisant peut s'obtenir en combinant les résultats de Crandall (perforation) [2, Ann. A] et de Bruneau (cavité) [3, Sec. 3.7] :

$$ \begin{aligned} \hat{z}_p^{\text{ref}}(s) &= l\rho_0 [1 - \Lambda(k_v d/2)]^{-1} s \\ jk_c^{\text{ref}}(s) &= \frac{s}{c_0} \left( \frac{1 + (\gamma - 1) \Lambda (k_\alpha d_c/2)}{1 - \Lambda (k_v d_c/2)} \right)^{1/2}, \end{aligned} $$
où $\Lambda(s) := \frac{2}{s} \frac{I_1}{I_0}(s)$ ($I_1$ et $I_0$ désignent les fonctions de Bessel modifiées d'ordre 1 et 0), $k_v := \sqrt{s/v}$ et $k_\alpha := \sqrt{s/\alpha}$ sont les nombres d'onde associés à la diffusion visqueuse et thermique dans la cavité et $\gamma$ est le rapport des capacités thermiques isobare et isochore. Le groupement sans dimension $|k_v d/2|$ est appelé le nombre de Stokes et est proportionnel au rapport entre le diamètre $d$ et l'épaisseur de la couche limite acoustique.

## 2.3 Modèle approximé

Le modèle d'impédance cible pour la simulation, $\hat{z}_s^{\text{phys}}$, est une approximation du modèle de référence $\hat{z}_s^{\text{ref}}$

$$ \hat{z}_s^{\text{phys}}(s) = 3 \frac{\hat{z}_p^{\text{phys}}(s)}{\underbrace{\overbrace{\sigma(d/2)^2}^{a_0}}_{a_1/2} + \underbrace{\overbrace{\sigma d/2}^{a_1}}_{a_1} \sqrt{s} + \underbrace{\overbrace{\sigma d/2}^{a_1} s + \dots} + \dots + \underbrace{\overbrace{z_0 \coth[(b_0 + b_1 s) l_c]}^{\hat{z}_c^{\text{phys}}(s)}}_{\hat{z}_c^{\text{phys}}(s)}. \quad (2) $$

L'impédance de perforation $\hat{z}_p^{\text{phys}}$ est l'approximation à haut nombre de Stokes en $O(|k_v d/2|^{-1})$ de $\hat{z}_p^{\text{ref}}$, et l'impédance de cavité $\hat{z}_c^{\text{phys}}$ consiste à approximer le nombre d'onde $jk_c^{\text{ref}}$. Les coefficients positifs $b_1$ et $b_0$ sont calculés via

$$ b_0 = \mathfrak{R}[jk(j\omega_b)], \quad b_1 = \frac{\mathfrak{I}[jk(j\omega_b)]}{\omega_b}, $$

où $\omega_b$ est une pulsation choisie (par exemple à la résonance). Une comparaison entre $\hat{z}_s^{\text{phys}}$ et des mesures effectuées en tube à impédance est proposée Fig. 2, où l'on observe un accord satisfaisant, une fois quelques corrections effectuées sur les paramètres géométriques.

Dans l'expression (2) on observe la présence d'un terme fractionnaire $\sqrt{s}$, qui est la manifestation de la diffusion visqueuse dans la perforation ; on a en effet

$$ \sqrt{\nu\omega} \propto \nu\delta^{-1}, $$

en notant $\delta$ l'épaisseur de la couche limite acoustique. Si l'on sépare résistance et réactance, on a, pour $\omega \in \mathbb{R}$ :

$$ \mathfrak{R}[\hat{z}_p^{\text{phys}}(j\omega)] = a_0 + \frac{1}{\sqrt{2}} a_1 \sqrt{\omega} \\ \mathfrak{I}[\hat{z}_p^{\text{phys}}(j\omega)] = a_1 \omega + \frac{1}{\sqrt{2}} a_1 \sqrt{\omega}. $$

On constate alors que l'effet du terme fractionnaire se manifeste principalement sur la résistance, puisque sa contribution relative sur la réactance est en $O(|k_v d/2|^{-1})$ :

$$ \frac{\mathfrak{I}[a_1 \sqrt{s}]}{\mathfrak{I}[a_1 s]} = \frac{1}{\sqrt{2}} \frac{a_1}{a_1 \sqrt{\omega}} = \frac{\sqrt{2}}{|k_v d/2|}. $$

# 3 Analyse temporelle de l'impédance de perforation

On a vu à la section 2 que, physiquement, l'impédance de perforation se ramenait à un polynôme fractionnaire

$$ \hat{z}_p^{\text{phys}}(s) = a_0 + a_1 \sqrt{s} + a_1 s \quad (\mathfrak{R}[s] > 0), \quad (3) $$

<figure><img src="image_4.png" /><figcaption>(a) Impédance normalisée.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Coefficient d'absorption $\alpha_b := 1 - |\frac{z/z_0 - 1}{z/z_0 + 1}|^2$.</figcaption></figure>

FIGURE 2 – Comparaison entre mesures et modèles. (—○—) : mesures en tube à impédance pour le liner SDOF décrit Tab. 1. (—---) : modèle physique $\hat{z}_s^{\text{phys}}$ (2) avec $f_b = 2\pi\omega_b = 2.6$ kHz (résonance), $l$ corrigé de +3 %, $l_c$ corrigé de -3 % et $\sigma$ à 6.5 %. (—) : modèle physique $\hat{z}_s^{\text{phys}}$ (2) avec $f_b = 2\pi\omega_b = 2.6$ kHz (résonance), sans corrections.

où les coefficients $a_0, a_1$ et $a_1/2$ sont positifs. Cette section se concentre sur le terme fractionnaire $\sqrt{s}$, et est organisée en deux temps : 3.1 montre l'admissibilité de $\sqrt{s}$ (récemment remise en question dans la littérature [1]); 3.2 présente une représentation temporelle de $\sqrt{s}$, qui permet de formuler une TDIBC continue, qui sera discrétisée dans la section 4.

## 3.1 Admissibilité

En suivant Rienstra [1], un modèle d'impédance $\hat{z}(s)$ est dit admissible s'il respecte les conditions de passivité, réalité et causalité.

**Passivité** Le modèle doit dissiper l'énergie acoustique : pour $\omega \in \mathbb{R}$, $\mathfrak{R}[\hat{z}(j\omega)] \ge 0$.

**Réalité** La réponse impulsionnelle $z_p$ doit être à valeurs réelles.

Il suffit de vérifier la symétrie hermitienne du modèle, i.e. $\hat{z}(s^*) = \hat{z}(s)^*$. Remarquons que cette condition fait défaut à certains modèles physiques ; par exemple, tous les modèles qui incluent une correction visqueuse en $\sqrt{\omega}$ (au lieu de $\sqrt{j\omega}$).

**Causalité** La réponse impulsionnelle doit être nulle aux instants négatifs, i.e. $z_p(t) = 0$ pour $t < 0$.

La vérification rigoureuse de cette condition demande un certain soin, selon l'espace auquel appartient $\hat{z}$; d'autant plus que l'utilisation purement formelle de la transformée de Fourier (i.e. « $s = j\omega$ ») peut mener
à des contresens pour certains modèles. On réfère le lecteur au chapitre 1 de [4], qui traite du cas où $\hat{z}$ est dans l'espace de Hardy $H^2$ (théorème 1.6.1) et celui où $\hat{z}$ est à croissance polynomiale (théorème 1.8.2) : c'est de ce dernier théorème que découle la causalité de $z_p^{\text{phys}}$.

En résumé, l'étude de la fonction $\hat{z}_p^{\text{phys}}$ nous permet de conclure à sa validité. Une autre vision nous est offerte par une théorie du calcul fractionnaire [5]. En réécrivant (3) sous la forme

$$ \hat{z}_p^{\text{phys}}(s) = a_0 + a_1 \hat{h}(s) s + a_2 s \quad (\mathbb{R}[s] > 0), \quad (4) $$

où $\hat{h}(s) = 1/\sqrt{s}$, on peut montrer que le modèle (3) mène, en temporel, à l'équation différentielle fractionnaire suivante

$$ p = a_0 u_n + a_1 \frac{1}{2} d^{\frac{1}{2}} u_n + a_2 \dot{u}_n, \quad (5) $$

où $u_n$ désigne $u \cdot n$, et $\dot{u}_n$ sa dérivée temporelle. La dérivée fractionnaire présente dans (5) est la dérivée au sens de Caputo [6], qui s'exprime

$$ d^{\frac{1}{2}} u_n = I^{\frac{1}{2}} \dot{u}_n = \int_0^t \frac{1}{\sqrt{\pi \tau}} \dot{u}_n (t - \tau) \, d\tau, \quad (6) $$

où l'on reconnait l'intégrale de Riemann-Liouville appliquée à $\dot{u}_n$. Le calcul numérique de la convolution (6) pourrait se faire par des techniques de convolutions récursives, bien documentées dans la littérature sur les TDIBCs. Cependant, il est possible de ramener le calcul de (6) à celui d'équations différentielles ordinaires (EDO) du premier ordre, et ce de façon théoriquement exacte. C'est l'objet de la prochaine section, 3.2.

## 3.2 Représentation diffusive

On a vu en 3.1 que l'intégrateur fractionnaire $I^{\frac{1}{2}}$ est un opérateur de convolution, dont le noyau $h$ est

$$ \forall t > 0, \quad h(t) = \frac{1}{\sqrt{\pi t}}, \quad (7) $$

qui admet pour transformée de Laplace :

$$ \forall s \in \mathbb{C} : \mathbb{R}[s] > 0, \quad \hat{h}(s) = \frac{1}{\sqrt{s}}. $$

Cette section fournit une autre expression du noyau (7), plus adaptée à la simulation numérique : sa représentation diffusive (ou représentation par coupure [7]). Pour obtenir cette représentation, on (re)calcule la transformée de Laplace inverse de $\hat{h}$ :

$$ \forall t > 0, \quad h(t) = \frac{1}{2j\pi} \int_{c+j\mathbb{R}} \hat{h}(s) e^{st} ds. \quad (8) $$

Le calcul de (8) se mène avec une technique bien connue en analyse complexe [8] : l'utilisation conjointe du théorème des résidus et du lemme de Jordan sur un contour fermé. La subtilité du calcul réside dans le fait que $\hat{h}(s)$ n'est défini que dans le demi-plan droit : toute extension au demi-plan gauche mène à une fonction multivaluée. Le choix d'une coupure sur $\mathbb{R}_-$ permet de définir une fonction étendue $\hat{h}_e$, qui coïncide

<figure><img src="image_6.png" /><figcaption>FIGURE 3 – Contour de Bromwich pour le calcul de (8). Le saut de $\hat{h}_e$ à travers la coupure est à l'origine de la représentation diffusive (9).</figcaption></figure>

avec $\hat{h}(s)$ dans le demi-plan droit, et préserve la symétrie hermitienne :

$$ \forall s \in \mathbb{C} \setminus \mathbb{R}_-, \quad \hat{h}_e(s) = \frac{1}{\sqrt{s}}. $$

Sur le contour de Bromwich représenté Fig. 3, le théorème des résidus donne, pour tout $\epsilon > 0$ et $R > \epsilon$ :

$$ \int_{\Gamma_R \cup \Gamma_e \cup C_\epsilon} \hat{h}_e(s) e^{st} ds + \int_{C-jR}^{c+jR} \hat{h}_e(s) e^{st} ds = 0. $$

Après passage à la limite sur $R$ (lemme de Jordan) puis sur $\epsilon$, seule reste la contribution de la coupure, qui se regroupe en une seule intégrale :

$$ h(t) = \int_{\mathbb{R}_+} \mu(\xi) e^{-\xi t} d\xi, \text{ soit } \hat{h}(s) = \int_{\mathbb{R}_+} \frac{\mu(\xi)}{s+\xi} d\xi. \quad (9) $$

L'expression (9) est la représentation diffusive (ou par coupure) de $h(t)$. Le poids diffusif $\mu(\xi)$, qui s'exprime en fonction du saut de $\hat{h}_e$ à travers la coupure [7], vaut ici :

$$ \forall \xi > 0, \quad \mu(\xi) = \frac{1}{2j\pi} [\hat{h}_e(-\xi^{-}) - \hat{h}_e(-\xi^{+})] = \frac{1}{\pi\sqrt{\xi}}. \quad (10) $$

La représentation (9) nous renseigne sur la dynamique de l'intégration fractionnaire, qui consiste en un continuum de dynamiques du premier ordre, s'étendant des hautes fréquences ($\xi \gg 1$, mémoire courte) aux basses fréquences ($\xi \ll 1$, mémoire longue).

## 3.3 TDIBC continue

La représentation diffusive de l'intégrateur fractionnaire (9) permet de ré-écrire l'impédance de perforation (4) sous la forme théoriquement exacte :

$$ \hat{z}_p^{\text{phys}}(s) = a_0 + a_1 \left[ \int_{\mathbb{R}_+} \frac{\mu(\xi)}{s+\xi} d\xi \right] s + a_2 s, \quad (11) $$

Cette expression permet de ramener l'équation différentielle fractionnaire (5) à

$$ p(t) = a_0 u_n + a_1 \frac{1}{2} \int_{\mathbb{R}_+} \mu(\xi) (e^{-\xi \cdot \star \dot{u}_n}(t)) d\xi + a_1 \dot{u}_n(t). \quad (12) $$

L'intérêt de se ramener aux convolutions $e^{-\xi \cdot \star \dot{u}_n}$ est que celles-ci peuvent se réduire à des EDO du premier ordre :
après quelques calculs élémentaires, l'expression (12) peut
se ré-écrire

$$
p(t) = a_0 u_n + a_{\frac{1}{2}} \int_{\mathbb{R}_+} \mu(\xi) [-\xi \varphi_{\xi}(t) + u_n(t)] d\xi + a_1 \dot{u}_n(t), \quad (13)
$$

où les variables diffusives $\varphi_{\xi}$ suivent une dynamique du premier ordre (filtre passe-bas de fréquence de coupure $\xi$) :

$$
\forall \xi > 0, \quad \left\{
\begin{array}{l}
\dot{\varphi}_{\xi}(t) = -\xi\varphi_{\xi}(t) + u_n(t) \\
\varphi_{\xi}(0) = \frac{1}{\xi}u_n(0)
\end{array}
\right.
\qquad (14)
$$

Mathématiquement, la dérivation fractionnaire (6) est ainsi
vue comme un observateur sur l'état (de dimension *infinie*)
($\varphi_{\xi}$)$_{\xi>0}$. La condition initiale non-nulle marque la différence
entre la dérivation fractionnaire au sens de Caputo et de
Riemann-Liouville. Les équations (13,14) constituent une
TDIBC *continue*, dont la discrétisation est traitée en Sec. 4.

4 Modèle numérique pour l'impédance
de perforation

4.1 TDIBC discrète

Le modèle numérique pour l'impédance de perforation, $\hat{z}_p^{\text{num}}$, se déduit simplement de l'analyse de $\hat{z}_p^{\text{phys}}$ effectuée à la section 3 ; la discrétisation de la représentation diffusive (11) donne :

$$
\hat{z}_p^{\text{num}}(s) = a_0 + a_{\frac{1}{2}} \left[ \sum_{k=1}^{N_{\varepsilon}} \frac{\tilde{\mu}_k}{s + \xi_k} \right] s + a_1 s, \quad (15)
$$

où les poids $\tilde{\mu}_k$ et pôles $\xi_k$ sont à choisir (voir Sec. 4.2).
L'expression temporelle correspondante est

$$
p = a_0 u_n + a_{\frac{1}{2}} \sum_{k=1}^{N_{\varepsilon}} \tilde{\mu}_k [-\xi_k \varphi_k(t) + u_n(t)] + a_1 \dot{u}_n, \quad (16)
$$

où les $N_\xi$ variables diffusives $\varphi_k$ suivent chacune une EDO
d'ordre 1 :

$$
\forall k \in [1, N_{\xi}], \quad \left\{
\begin{array}{l}
\dot{\varphi}_k(t) = -\xi_k \varphi_k(t) + u_n(t) \\
\varphi_k(0) = \frac{1}{\xi_k} u_n(0)
\end{array}
\right.
. \quad (17)
$$

Les équations (16,17) constituent la TDIBC discrète approximant l'impédance de perforation. L'avantage d'une formulation par EDO sur une méthode de convolution récursive a récemment été mise en évidence par Dragna et al. [9].

4.2 Calcul des poids par optimisation

Les degrés de liberté de la TDIBC (15) sont les $N_\xi$ pôles $\xi_k$ et poids $\tilde{\mu}_k$. Numériquement, l'on souhaite réaliser l'approximation

$$
\sqrt{\bar{j}\omega} \approx \left[ \sum_{k=1}^{N_{\epsilon}} \frac{\tilde{\mu}_k}{j\omega + \xi_k} \right] j\omega \quad (18)
$$

sur une bande d'intérêt [$\omega_{\min}$, $\omega_{\max}$], pour $N_\xi$ le plus faible
possible. Pour des applications aéroacoustiques, une TDIBC

<figure><img src="image_7.png" /><figcaption>FIGURE 4 – Approximation de $\sqrt{j\bar{\omega}}$, optimisation avec $w(\omega) = \omega$ et placement des pôles linéaires entre 10 Hz et 10 kHz. Tracé : $\mathbb{R}[a_{\frac{1}{2}}\sqrt{j\bar{\omega}}/z_0]$ (—), approximation avec 2 pôles (—–) et (○), 6 pôles (—) et (x).</figcaption></figure>

peut être considéré large-bande (broadband) si la largeur de
bande est de l'ordre de 10 kHz.

L'utilisation de l'expression analytique (10) avec des quadratures numériques s'avère peu intéressante : les approximations réalisées ne sont pas parcimonieuses. On recourt alors à une méthode d'optimisation, décrite par exemple dans [7]. La fonction coût *J* proposée est un écart quadratique dans le domaine fréquentiel :

$$
J(\tilde{\mu}, \xi) = \int_{\omega_{\min}}^{\omega_{\max}} w(\omega) \left| \sum_{k=1}^{N_{\xi}} \frac{\tilde{\mu}_k}{j\omega + \xi_k} j\omega - \sqrt{j\bar{\omega}} \right|^2 d\omega, \quad (19)
$$

où la fonction w est un poids à choisir. Pour approximer $\sqrt{j\bar{\omega}}$, il n'est pas indispensable d'optimiser simultanément les pôles et les poids : les premiers peuvent être choisis. La procédure d'optimisation est alors :

1. Choix d'une distribution de $N_\xi$ pôles $(\xi_k^\circ)_k$. Typiquement, une répartition linéaire ou logarithmique sur $[\omega_{\min}, \omega_{\max}]$ convient.

2. Calcul des poids optimaux $\mu_{\text{opt}}$ par minimisation de $\tilde{\mu} \mapsto J(\tilde{\mu}, \xi^\circ)$. Il s'agit d'un problème d'optimisation quadratique (sur-déterminé), dont la solution s'obtient instantanément par pseudo-inverse.

Une illustration de l'approximation de $\sqrt{j\bar{\omega}}$ est proposée Fig. 4, pour respectivement 2 et 6 pôles répartis linéairement entre 10 Hz et 10 kHz. On observe que quelques pôles suffisent à réaliser une approximation large-bande satisfaisante.

5 Modèle numérique pour l'impédance
de cavité

Comme indiqué dans l'introduction, le modèle physique
de cavité

$$
\hat{z}_c^{\text{phys}}(s) = z_0 \coth \left[ (b_0 + b_1 s) l_c \right] \quad (\mathfrak{R}[s] > 0) \quad (20)
$$

correspond au terme de cavité du modèle numérique EHR
[1], on a donc $\hat{z}_c^{\text{num}} = \hat{z}_c^{\text{phys}}$. En utilisant temporairement les
notations de Rienstra, on a en effet

$$
\hat{z}_{\text{EHR}}(s) = \beta \coth \left[ \frac{\epsilon}{2} + \frac{\nu\Delta t}{2}s \right].
$$

L'identification entre $(\beta, \epsilon, \nu\Delta t)$ et $(z_0, b_0 l_c, b_1 l_c)$ fournit
alors un moyen physique de régler l'EHR, sans passer par un
processus d'optimisation.
Pour la simulation numérique de (20), plusieurs méthodes sont documentées dans la littérature. Pascal et al. [10] ont récemment proposé une méthode permettant d'éviter l'emploi de convolutions récursives, au prix de la résolution des équations d'Euler monodimensionnelles dans la cavité du liner.

Comme indiqué dans l'introduction, on se concentre ici sur la simulation numérique de l'impédance de perforation : c'est l'objet de la section 6.

# 6 Étude numérique dans une méthode Galerkin Discontinu

Cette section présente une intégration de la TDIBC (16,17) dans une méthode de type DG, pour les équations d'Euler linéarisées monodimensionnelles. Après un bref rappel de la formulation DG en Sec. 6.1, ce cas élémentaire nous permet, en Sec. 6.2, de valider la fonction de flux numérique proposée, puis d'étudier l'impact de la TDIBC sur le CFL en Sec. 6.3.

On considère ici le cas du tube à impédance monodimensionnel : une source est placé en $x = 0$, le liner étant en $x = L$. En notant le vecteur des variables acoustiques $\underline{q} = [u, p/z_0]^T$, les équations du problème sont

$$ (\text{LEE}) \qquad \partial_t \underline{q}(x,t) + A_x \cdot \partial_x \underline{q}(x,t) = \underline{0} \qquad (21) $$

$$ (\text{Source}) \qquad \underline{q}_{\text{in}}(t) \qquad (22) $$

$$ (\text{Impédance}) \qquad \frac{p}{z_0}(L,t) = \left[ \frac{z_p \text{num}}{z_0} \star u(L, \cdot) \right](t), \qquad (23) $$

où la jacobienne de la fonction de flux est simplement :

$$ A_x = c_0 \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}. $$

## 6.1 Formulation Galerkin Discontinu

Un exposé concis et pratique de la méthode Galerkin Discontinu est disponible dans l'ouvrage de Hesthaven et Warburton [11]. Il s'agit d'une méthode à stencil compact, d'ordre élevé, dont les propriétés de dissipation et de dispersion numériques la rendent adaptée à l'acoustique.

On discrétise l'espace $[0, L]$ en $N_K$ éléments ; à chaque élément est associé $N$ noeuds, $N$ étant l'ordre spatial visé. Le point de départ de la méthode est la formulation faible de (21), écrite sur chacun des $N_K$ éléments ; en choisissant pour fonctions tests les polynômes $(\psi_i)_i$ de degré au plus $N-1$, cette dernière s'écrit :

$$ \begin{gathered} \forall k \in [1, N_K], \quad \forall i \in [1, N], \\ \langle \psi_i, \partial_t \underline{q} \rangle_{\Omega_k} - \left\langle \frac{d\psi_i}{dx}, A_x \cdot \underline{q} \right\rangle_{\Omega_k} = - \int_{\partial\Omega_k} \psi_i n_x A_x \cdot \underline{q} \, dx, \end{gathered} \qquad (24) $$

où la normale unitaire se réduit ici à $\pm 1$. Puisque $\underline{q}$ n'est pas défini de façon univoque à la frontière entre éléments, la connexion entre éléments se fait via une fonction de flux numérique, notée $\underline{f}^*$ ; (24) est alors reformulé en :

$$ \begin{gathered} \forall k \in [1, N_K], \quad \forall i \in [1, N], \\ \langle \psi_i, \partial_t \underline{q} \rangle_{\Omega_k} - \left\langle \frac{d\psi_i}{dx}, A_x \cdot \underline{q} \right\rangle_{\Omega_k} = - \int_{\partial\Omega_k} \psi_i \underline{f}^* \, dx. \end{gathered} \qquad (25) $$

Après projection de $\underline{q}$ sur la base polynomiale $(\psi_i)_i$, et concaténation de (25) pour tous les éléments, on aboutit, in fine, à la formulation globale suivante :

$$ \underline{M} \cdot \dot{\underline{Q}} + \underline{K} \cdot \underline{Q} = - \left[ \left[ \int_{\partial \Omega_k} \psi_i \underline{f}^* dx \right]_{i \in [1, N]} \right]_{k \in [1, N_K]} \qquad (26) $$

où $\underline{Q}$ est le vecteur colonne regroupant les $N_K \times N \times 2$ degrés de libertés $(u$ et $p/z_0$ sur chacun des $N_K \times N$ noeuds).

À l'intérieur du domaine de calcul et pour la source en $x = 0$, on utilise pour la fonction de flux numérique $\underline{f}^*$ une décomposition de flux (Flux Vector Splitting), qui résout (ici) de façon exacte le problème de Riemann. L'enjeu réside dans la définition d'une fonction de flux numérique pour la condition d'impédance (23) en $x = L$ : c'est l'objet de la prochaine section.

## 6.2 Flux numérique pour la TDIBC

À la frontière du domaine, on choisit un flux numérique centré pour imposer la TDIBC ; ce flux se formule :

$$ \underline{f}^*_{\underline{z}}(\underline{q}, \underline{n}) := n_x A_x \left( \frac{\underline{q} + \underline{q}_z}{2} \right), \qquad (27) $$

où $\underline{q}$ est l'état à la frontière du domaine (ici, en $x = L$), et $\underline{q}_z$ est un état fictif « dans » le liner, tel que schématisé Fig. 5. L'utilisation de (27) revient à reconstruire l'état à la frontière par une moyenne des états gauche et droite.

<figure><img src="image_8.png" /><figcaption>FIGURE 5 – Schéma illustrant l'état fictif $\underline{q}_z$ pour la définition du flux centré (27).</figcaption></figure>

L'état fictif $\underline{q}_z$ est construit de sorte que, si $\underline{q}$ vérifie la condition d'impédance, alors $\underline{q}_z = \underline{q}$. Pour représenter la condition (5), on propose l'état fictif suivant :

$$ \begin{align} u_z &= u \nonumber \\ p_z &+ \left[ (1 - \beta_0)(1 - \beta_{\frac{1}{2}})(1 - \beta_1) - 1 \right] p = \nonumber \\ &\quad (1 + \beta_0)(1 - \beta_{\frac{1}{2}})(1 - \beta_1)u_n \nonumber \\ &\quad + (1 - \beta_0)(1 + \beta_{\frac{1}{2}})(1 - \beta_1)\frac{1}{2}u_n(t) \nonumber \\ &\quad + (1 - \beta_0)(1 - \beta_{\frac{1}{2}})(1 + \beta_1)\dot{u}_n. \tag{28} \end{align} $$

Pour des raisons de conditionnement numérique, il est préférable d'employer les coefficients de réflexion

$$ \beta_i := \frac{a_i - 1}{a_i + 1} $$

plutôt que les $a_i$, puisque $|\beta_i| < 1$ pour $a_i > 0$. Avec (27) et (28), la formulation globale (26) devient finalement

$$ \underline{M} \cdot \dot{\underline{Q}} + \underline{K} \cdot \underline{Q} = \underline{F}_{\underline{0}} \cdot \underline{Q} + \underline{F}_{\underline{\frac{1}{2}}} \mathrm{d}^{\underline{1/2}}[\underline{c} \cdot \underline{Q}] + \underline{F}_{\underline{1}} \cdot \dot{\underline{Q}} + \underline{F}_{\underline{\text{in}}} \cdot \underline{q}_{\text{in}}(t), \qquad (29) $$

où $\underline{c}$ est nul sauf pour la vitesse acoustique au noeud situé sur le liner, en $x = L$. L'utilisation de la forme harmonique de (29) a permis de valider la fonction de flux numérique proposée jusqu'à l'ordre spatial $N = 8$.
Pour mener le calcul en temporel, comme exposé en Sec. 4, on calcule la dérivée fractionnaire à l'aide des $N_\xi$ variables diffusives $\varphi_k$ ; en notant $\underline{\underline{\Lambda}} := \text{diag}(\xi_k)_k$, on a :

$$d_t^{1/2} [\underline{\underline{c}} \cdot \underline{\underline{q}}] (t) \approx -\underline{\underline{\mu}} \cdot \underline{\underline{\Lambda}} \cdot \underline{\underline{\varphi}} (t) + \left( \sum_{i=1}^{N_\xi} \underline{\underline{\mu}}_i \right) \underline{\underline{c}} \cdot \underline{\underline{q}} (t) . \quad (30)$$

L'équation effectivement intégrée pour la simulation temporelle s'obtient en injectant (30) dans (29). On aboutit à une extension de l'état sur la formulation globale, i.e. l'intégration temporelle s'effectue non plus seulement sur l'état $\underline{\underline{Q}}$, mais sur l'état étendu $[\underline{\underline{Q}}, \underline{\underline{\varphi}}]^T$ :

$$\left[ \begin{array}{c} \dot{\underline{\underline{Q}}} \\ \underline{\underline{\varphi}} \end{array} \right] = \left[ \begin{array}{c|c} \underline{\underline{A}} + (\sum_i \underline{\underline{\mu}}_i) \underline{\underline{F}}_i' \otimes \underline{\underline{c}} & -\underline{\underline{F}}_i' \otimes (\underline{\underline{\mu}} \cdot \underline{\underline{\Lambda}}) \\ \hline \underline{\underline{1}}_{N_\xi} \otimes \underline{\underline{c}} & -\underline{\underline{\Lambda}} \end{array} \right] \cdot \left[ \begin{array}{c} \underline{\underline{Q}} \\ \underline{\underline{\varphi}} \end{array} \right], \quad (31)$$

où, par comparaison à (29), on a omis le terme source lié à $\underline{\underline{q}}_{in}$, et noté $A := (M - F_1)^{-1} (F_0 - K)$ et $\underline{\underline{F}}_i' := (M - F_1)^{-1} \underline{\underline{F}}_i$. Ici encore, la forme harmonique permet de valider (31) et de mettre en évidence un des intérêts de l'état étendu : l'ordre spatial sur $\underline{\underline{Q}}$ reste l'ordre visé, N. La formulation diffusive (30) est donc satisfaisante sur le plan de la précision : son coût pour une simulation temporelle est étudiée en Sec. 6.3.

## 6.3 Impact sur le CFL

Pour mesurer le coût d'une simulation temporelle de (31), on considère le nombre CFL maximum

$$\text{CFL}_{\text{max}} := \frac{c_0}{\Delta x} \Delta t_{\text{max}}, \quad (32)$$

où $\Delta t_{\text{max}}$ est le pas de temps correspondant à la limite de stabilité du schéma utilisé pour l'intégration en temps. Dans cette étude, on a utilisé un Low-Storage Runge-Kutta 4 explicite, dont la limite de stabilité linéaire peut être estimée de façon conservative via

$$\Delta t_{\text{max}} := \frac{4.4}{\max_i |\lambda_i|}, \quad (33)$$

en notant $\lambda_i$ les valeurs propres de la matrice (31). Ces valeurs propres dépendent des poids diffusifs $\underline{\underline{\mu}}$, du pôle de plus basse fréquence $\xi_{\text{min}}$ et du pôle de plus haute fréquence $\xi_{\text{max}}$. Numériquement, on vérifie que c'est $\xi_{\text{max}}$ qui contrôle l'impact de la TDIBC sur le CFL : l'augmentation de $\xi_{\text{max}}$ tend à réduire $\text{CFL}_{\text{max}}$.

Cependant, sur le cas du tube à impédance étudié ici, une étude paramétrique a mis en évidence que cet impact n'était pas sensible. En effet, la valeur de $\xi_{\text{max}}$ a partir de laquelle $\text{CFL}_{\text{max}}$ diminue est bien supérieure à la fréquence maximale $f_{\text{max}}$ résolue par le schéma DG : $\xi_{\text{max}}^{\text{lim}} \gg f_{\text{max}}$. En pratique, le choix de $\xi_{\text{max}} \approx f_{\text{max}}$ permet de s'assurer que l'impact de la condition d'impédance sur le pas de temps soit neutre.

## 7 Conclusion

Dans cet article, la démarche suivie pour concevoir la TDIBC repose sur l'analyse des modèles physiques. La représentation diffusive de l'impédance de perforation (11), présentée en Sec. 3, a permis la formulation, en Sec. 4, d'une TDIBC (16,17) basée uniquement sur des EDO du premier

ordre. L'étude présentée en Sec. 6 a exposé l'intérêt de la formulation retenue et fournit un critère simple pour le choix de la fréquence maximale dans le calcul de la représentation diffusive discrète (15).

Une extension de la présente analyse à d'autres modèles physiques fait l'objet d'un article prochainement soumis ; une telle démarche renseigne sur l'applicabilité des TDIBC actuellement disponibles dans la littérature. Des pistes envisageables pour la poursuite des travaux sont : l'application à un conduit bi-dimensionnel ; l'extension aux cas avec écoulement et/ou forts niveaux sonores (ces derniers nécessitent une TDIBC non-linéaire).

## Remerciements

Ces travaux ont été initialement financés par l'ISAE-SUPAERO, et sont actuellement financés conjointement par l'ONERA et la DGA. Le premier auteur remercie l'ISAE-SUPAERO pour le financement de la mission.

## Références

[1] Sjoerd W. RIENSTRA. "Impedance Models in Time Domain, Including the Extended Helmholtz Resonator Model". In: 12th AIAA/CEAS Aeroacoustics Conference. AIAA Paper 2006-2686. Cambridge, MA, USA, 8-10 mai 2006.

[2] Irving Bardshar CRANDALL. Theory of vibrating systems and sound. D. Van Nostrand Company, 1926.

[3] Michel BRUNEAU. Fundamentals of acoustics. John Wiley & Sons, 2013.

[4] Herch Moyses NUSSENZVEIG. Causality and dispersion relations. Academic Press New York, 1972.

[5] Denis MATIGNON. "An introduction to fractional calculus". In: Scaling, Fractals and Wavelets. Sous la dir. de Patrice ABRY, Paulo GONÇALVÊS et Jacques LEVY-VEHEL. ISTE - Wiley, 2009, p. 237-277.

[6] Michele CAPUTO. "Vibrations of an infinite plate with a frequency independent Q". In: J. Acoust. Soc. Am. 60.3 (1976), p. 634-639.

[7] Thomas HÉLIE et Denis MATIGNON. "Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses". In: Math. Mod. Meth. Appl. S. 16.04 (2006), p. 503-536.

[8] Dean G. DUFFY. Transform methods for solving partial differential equations. 2e éd. CRC press, 2004.

[9] D. DRAGNA, P. PINEAU et P. BLANC-BENON. "A generalized recursive convolution method for time-domain propagation in porous media". In: J. Acoust. Soc. Am. 138.2 (2015), p. 1030-1042.

[10] L. PASCAL, E. PIOT et G. CASALIS. "A New Implementation of the Extended Helmholtz Resonator Acoustic Liner Impedance Model in Time Domain CAA". In: J. Comput. Acoust. (2015).

[11] Jan S. HESTHAVEN et Tim WARBURTON. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Texts in Applied Mathematics. Springer-Verlag New York, 2008.