# Calcul du rayonnement acoustique d'onde d'instabilité utilisant une PML (Perfectly Matched Layer)

Thomas Emmert, Christophe Bailly

► To cite this version:

Thomas Emmert, Christophe Bailly. Calcul du rayonnement acoustique d'onde d'instabilité utilisant une PML (Perfectly Matched Layer). 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-00089619

HAL Id: hal-00089619

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

Submitted on 25 Apr 2023

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 du rayonnement acoustique d'ondes d'instabilité utilisant une PML
(Perfectly Matched Layer)

Thomas Emmert, Christophe Bailly

École Centrale de Lyon

Laboratoire de Mécanique des Fluides et d'Acoustique, LMFA UMR CNRS 5509

36 avenue Guy de Collongue

69134 Ecully, France

thomas.emmert@edf.fr

## Abstract :

In this paper a Perfectly Matched Layer (PML) formulation was tested as boundary condition in linearized Euler's Equations. Counted among absorbing boundary conditions, it was recently introduced in Computational Aeroacoustics (CAA). A formulation is presented that makes the implementation of the PML very easy and computationally efficient. This technique was evaluated with typical problems in CAA by using the linearized Euler equations in 2-D such as a convected incompressible vortex and a Kelvin-Helmholtz instability wave. The PML was compared directly to a standard sponge zone and shows an excellent performance.

## Mots-clefs :

### aéroacoustique numerique, Perfectly Matched Layer (PML), instabilité de Kelvin-Helmholtz

#### 1 Introduction

The importance of non-reflecting boundary conditions in numerical simulations with an open domain, such as those that occur in many practical problems in CAA and computational fluid dynamics, is crucial for time and cost accurate solutions. The most widely spread non reflecting boundary condition is the characteristic based one of Thompson [1]. Tam and Webb [2] followed another strategy and developed a boundary condition starting from the asymptotic solution of the linearized Euler equations in 2-D. As radiation boundary condition it performs very well for all angle of incidences and was therefore favored over the characteristic boundary condition.

However for CAA problems radiation boundary conditions are not efficient enough to exit turbulent structures. Once a vortex has entered the boundary, acoustic waves are generated and contaminate the solution. Therefore vorticity waves that are upwinded by the flow have to be dissipated before meeting the boundary condition. A very simple and common method is to apply a Laplacian filter on the fluctuations of the flow in a damping zone or sponge zone. Following Bogey & Bailly [3] after each time step the fluctuations of $\mathbf{u} = [\rho, u, v, p]^T$ are updated such a way that

$$ \mathbf{u}_i = \mathbf{u}_i - \sigma_{\text{LF}} \left[ \frac{1}{2} \mathbf{u}_i - \frac{1}{4} (\mathbf{u}_{i+1} + \mathbf{u}_{i-1}) \right] \quad \sigma_{\text{LF}} = \sigma_{\text{LF,max}} \left| \frac{x - x_b}{x_e - x_b} \right|^s \quad (1) $$

where $x_b$ and $x_e$ represent the beginning and the end of the damping zone respectively. $\sigma_{\text{LF}}$ is the maximum absorption rate reached at $x_e$ and $s$ determines the smoothness of the filter. However the efficiency of this damping zone with reasonable widths is not satisfying for many aero-acoustic applications and for even smooth damping, reflections at the beginning of the
zone can be observed. Berenger [4] introduced a absorbing technique called Perfectly Matched Layer (PML) that was builded up to be theoretically transparent for all fluctuations that entering in it.

In this paper the PML will be applied to the linearized Euler equations only, because in many aero-acoustic problems viscous and nonlinear effects can be neglected [5]. The linearized Euler equations for non uniform flows can be written as :

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{A} \frac{\partial \mathbf{u}}{\partial x} + \mathbf{B} \frac{\partial \mathbf{u}}{\partial y} + \mathbf{C}_x \mathbf{u} + \mathbf{C}_y \mathbf{u} = 0 \quad (2) $$

where

$$ \mathbf{A} = \begin{pmatrix} \bar{u} & \bar{\rho} & 0 & 0 \\ 0 & \bar{u} & 0 & 1/\bar{\rho} \\ 0 & 0 & \bar{u} & 0 \\ 0 & \gamma\bar{\rho} & 0 & \bar{u} \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} \bar{v} & 0 & \bar{\rho} & 0 \\ 0 & \bar{v} & 0 & 0 \\ 0 & 0 & \bar{v} & 1/\bar{\rho} \\ 0 & 0 & \gamma\bar{\rho} & \bar{v} \end{pmatrix}, $$

$$ c_x = \begin{pmatrix} \bar{u}_x & \bar{\rho}_x & 0 & 0 \\ \bar{u}\bar{u}_x & \bar{u}_x & 0 & 0 \\ \bar{u}\bar{v}_x & \bar{v}_x & 0 & 0 \\ 0 & \bar{p}_x & 0 & \gamma\bar{u}_x \end{pmatrix}, c_y = \begin{pmatrix} \bar{v}_y & 0 & \bar{\rho}_y & 0 \\ \bar{v}\bar{u}_y & 0 & \bar{u}_y & 0 \\ \bar{v}\bar{v}_y & 0 & \bar{v}_y & 0 \\ 0 & 0 & \bar{p}_y & \gamma\bar{v}_y \end{pmatrix}. $$

Note, that the terms $C_x$, $C_y$ containing the derivatives of the mean properties (denoted with an overbar) become zero for a uniform mean flow field.

## 2 Perfectly Matched Layer

The PML equations can be derived by applying a complex change of variable to the Fourier transformed Euler equations (2) such as

$$ x' = x + \frac{i}{\omega} \int_{x_b}^{x} \sigma_{PML} dx, $$

where $\sigma_{PML}$ presents the absorption rate and $x_b$ indicates the position of the interface between the Euler and the PML domain. One can demonstrate easily by impedance theory that the PML is theoretically transparent for disturbances that enters in it. The absorption rate can vary within the absorbing layer along the coordinate direction vertical to the absorbing layer interface. The distribution of $\sigma_{PML}$ was chosen according to (1).

Unfortunately the derivation of the two and more dimensional case leads to a system of equations that support exponentially growing unstable solutions [6]. Hu [7] proofed that unstable solutions occur when the component of the group velocity and the phase velocity in the corresponding direction doesn't have the same sign. Note, that this is only possible for convec- ted acoustic modes whereas the group and phase velocity of the vortical and entropic modes are always of the same sign.

To get the PML equations stable, Hu proposed to apply a proper time space transformation before performing the complex change of variable in order to get a propagation of acoustic waves in a medium at rest. This transformation is similar to the Prandtl-Glauert transformation in aerodynamics and can be written as

$$ \bar{t} = t + \beta x, $$
<figure><img src="image_2.png" /><figcaption>FIG. 1 – (a) Implementation of the PML as absorbing boundary layer. The vector **q** that contains the auxiliary variables only have to be calculated and stored in the layer. (b) Time evolution of $L_p$ for incompressible vortex — $D = 0, \dots, D = 7, \dots, D = 14$.</figcaption></figure>

where $\beta$ can be found for the case of a uniform mean flow in $x$-direction analytically as $\beta = \bar{u}/(1 - \bar{u}^2)$. Solving the compressible Rayleigh equation, Hu [8] found out that even for the nonuniform case such as subsonic shear layers, jet flows or Poiseuille flows a $\beta$ can be found that is constant. For the velocity and density mean profile of a jet, he found $\beta = 1/0.25$.

Starting from the formulation of Hu [7], the PML equations can be rewritten in a quasi conservative form, that is very simple to implement. They can be written in respect to the matrix form (2) as follows :

$$ \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{A}^* \mathbf{u} \mathbf{q}_x}{\partial x} + \frac{\partial \mathbf{B}^* \mathbf{u} \mathbf{q}_y}{\partial y} + \mathbf{H} \mathbf{q} + \mathbf{G} = 0, \quad (3) $$

$$ \frac{\partial \mathbf{q}}{\partial t} = \mathbf{u}, $$

where

$$ \mathbf{U} = [\rho, \bar{\rho}u, \bar{\rho}v, p], \quad \mathbf{u} \mathbf{q}_x = \mathbf{u} + \sigma_y \mathbf{q}, \quad \mathbf{u} \mathbf{q}_y = \mathbf{u} + \sigma_x \mathbf{q}, $$

$$ \mathbf{H} \mathbf{q} = \left( \mathbf{C}_x^* - \frac{\partial \mathbf{A}^*}{\partial x} \right) \mathbf{u} \mathbf{q}_x + \left( \mathbf{C}_y^* - \frac{\partial \mathbf{B}^*}{\partial y} \right) \mathbf{u} \mathbf{q}_y, \quad \mathbf{G} = \left[ (\sigma_x + \sigma_y)\mathbf{u} + \sigma_x \sigma_y \mathbf{q} + \sigma_x \beta \mathbf{A} \mathbf{u} \mathbf{q}_x \right]^* $$

and * indicates that the velocity components are multiplied with the mean density $\bar{\rho}$. The auxiliary variable **q** and its derivatives are only needed to be calculated and stored in the absorbing layer as Fig. 1(a) illustrates, what reduces the computational and storage effort. Note, that for $\sigma_x, \sigma_y = 0$ the PML equations become the ordinary linearized Euler in the conservative form.

### 3 Performance of the PML

The 2-D PML equations were implemented in a linearized Euler 2-D code that uses an optimized 11-pts finite difference scheme for the spatial discretization (4 points/λ) and an optimized 6-steps Runge-Kutta algorithm [9] to advance in time with CFL = 0.9. Also selective filtering was used ($D_s = 0.3$) in order to suppress grid-to-grid oscillations.

For testing the efficiency of the PML a benchmark test of [10] were used, consisting of an incompressible vortex convected in a uniform mean flow. The tests were done on a uniform
grid. The physical domain defined by a $N^2$ uniform grid with $N = 101$ was created such that $-50 \le x, y \le 50$. To limit the computational domain the boundary conditions of Tam and Webb [2] were applied on three additional points along all boundaries.

To estimate the magnitude of acoustic waves reflected back into the computational domain after the exit of disturbances, the time evolution of the residual fluctuating pressure $L_p$ is recorded. This residual pressure is based on the following $L_2$ norm :

$$L_p = \left[ \frac{1}{N^2} \sum_{i,j}^{N^2} p'_{i,j} \right]^{1/2}.$$

## 3.1 Vortex in a uniform flow

The vortex is upwinded in a uniform flow with $\bar{u} = (M, 0, 0)$, constant mean density $\bar{\rho} = 1$ and pressure $\bar{p} = 1/\gamma$ were tested and $\gamma = 1.4$. The initial conditions are

$$\begin{cases} p, \rho = & 0, \\ u = & x_1 \epsilon \exp[-\alpha(x_1^2 + x_2^2)], \\ v = & -x_2 \epsilon \exp[-\alpha(x_1^2 + x_2^2)], \end{cases}$$

where $\alpha = (\ln 2)/b^2$, the flow Mach number $M = 0.5$, the Gaussian half-width $b = 3$ and the amplitude $\epsilon = 10^{-3}$. Tests with different PML width $D$ were done. The parameters of the PML are $\sigma_{PML} = 2$ and $s = 2$.

Figure 1 (b) illustrates the time evolution of the residual pressure. Once the vortex has impacted the boundary at about $t = 70$, $L_p$ strongly increases if only the radiation boundary conditions are applied. Even with a layer of 7 points the reflections reduces 2 orders of magnitude. A layer of 14 points makes the reflections neglectable. Furthermore only neglectable reflections can be observed when the vortex enters the PML. In addition with grid stretching, the reflections can be reduced by another order of magnitude using a 7 point PML.

## 3.2 PML vs. Laplacian filter

In the following the performance of the PML is compared with the performance of a simple Laplacian filter according to (1). The maximum absorption rate was chosen to be $\sigma_{max} = 0.6$ and the smoothness to be $s = 1.5$ for both absorbing layers. The damping zones were applied over 30 additional points that were added at the outflow boundary condition. Fig. 2 (a) displays the results for the residual pressure that was recorded only in the physical domain $N^2$. The solid line represents the residual pressure when no absorbing layer ($\sigma = 0$) is applied by using the same computational domain. The oscillations at $t = 180$ are due to purely numerical disturbances that are generated when the vortex leaves the computational domain. Propagated by the numerical schemes many times faster than physical acoustic waves they reach and contaminate the physical domain $N^2$ earlier than acoustical waves.

Figure 2 (a) underlines how much better the PML works compared to the Laplacian filter. Constructed to be transparent the PML generates only small reflections, when the vortex impacts in the damping zone indicated by the first peak, whereas the classical damping zone generates acoustic waves that are about one order higher in magnitude. Moreover the Laplacian filter can't reduce the amplitude of the vortex sufficiently Fig. 2 (b) leading to a high residual pressure. In practice the Laplacian filter only works well in addition with strong grid stretching.
<figure><img src="image_3.png" /></figure>

FIG. 2 – (a) Residual pressure without absorbing layer — σ = 0. , ........ σ = 0.3, — — σ = 0.6 ; upper curve Laplacian filter, lower curve PML. (b) Vorticity along y = 0 — — — initial condition; — PML, — — — Laplacian filter at t = 126 ; ........ beginning of damping zone.

4 Test cases for linearized Euler equations with non uniform mean flow

The results of a simulation of a more complex test case with the mean flow of 2-D subsonic jet is presented now. Having an inflexing point in the velocity profile, the considered mean flow support solutions that are unstable. In the real physical problem, these Helmholtz instabilities are limited by non-linear and viscous effects and are the reason for the creation of vortices and turbulent structures. In the linear case they rise exponentially and without limitation making it very difficult to exit them properly by means of ordinary sponge zones such as (1). The test that was made, is the benchmark problem proposed by Morris et al. [11] that tries to model the radiation and refraction of sound waves through a two-dimensional subsonic and sub critical jet by means of the linearized Euler equations. The instability wave is excited by an acoustical source and absorbed with the PML technique.

The mean flow variables are given by

$$
\bar{u} = \begin{cases} u_j \exp[-(\ln 2)(\frac{y}{b} - \frac{h}{b})^2] & 0 \le y \le h \\ u_j & y \ge h \end{cases}
$$

The mean density $\bar{\rho}$ is derived by using the Crocco-Busemann relation

$$
\frac{1}{\bar{\rho}} = -\frac{1}{2} \frac{\gamma - 1}{\gamma \bar{\rho}} (\bar{u} - u_j) \bar{u} + \frac{1}{\rho_j} \frac{\bar{u}}{u_j} + \frac{1}{\rho_\infty} \frac{u_j - \bar{u}}{u_j},
$$

where the index *j* designs the other parameters of the jet and $\infty$ the ambient properties. The jet Mach number is calculated by $M_j = u_j/c_j = 0.756$. The speed of sound $c_j$ is calculated by using the ideal gas law $c_j = \sqrt{\gamma RT_j}$. The mean pressure was assumed to be constant in the whole computational domain. The parameters of the problem can be found in [11]. The Strouhal number based on the jet velocity $u_j$ and an estimation of the jet diameter $2b$ defines $\omega_0$ such as $St = 4\pi\omega_0 b/u_j = 0.85$.

The Simulation was made on a uniform grid 547 × 443 with −55.25 m ≤ x ≤ 159.2 m, −59.2 m ≤ y ≤ 59.2 m and carried out over 6000 iterations. As the jet is considered to be symmetric, a symmetry boundary condition was applied along y = 0. The PML is applied over 20 points at the inflow and outflow boundaries and over 17 points at the upper and lower
boundaries. The maximum absorption rates of the PML are taken to be $\sigma_{x,max} = 1/\Delta$ and $\sigma_{y,max} = 2/\Delta$, where $\Delta$ is the length of a grid cell.

Figure 3 presents a screen shot of the pressure field $p/\epsilon$ that was taken after 3600 iterations. As expected an instability develops and is convected downstream. The PML has no difficulty to dissipate and exit this instability wave without any spurious reflections.

<figure><img src="image_4.png" /><figcaption>FIG. 3 – Screen shot of the pressure contours $p/\epsilon$ between $\pm0.003$ and $\pm0.1025$ following a geometrical ratio of 2.</figcaption></figure>

## 5 Conclusion

The PML is a very efficient absorbing boundary condition for linearized Euler equations for subsonic problems in CAA. In combination with the radiation boundary condition this technique is able to exit exponentially growing solutions very efficiently. This is very difficult and computational costly, using the classical Laplacian filter for example. The instability of the PML could always suppressed by tuning the crucial value $\beta$ and the maximum absorption $\sigma_{max}$. Because of the efficiency it will be possible to calculate very difficult problems of CAA, where very exact methods are required as for numerical integration of the Green function. Finally nonlinear cases should be considered in future.

## Références

[1] THOMPSON, K.W. *Time Dependent Boundary Conditions for Hyperbolic Systems*. J. of Computational Physics, **68** (1978), 1-24.

[2] TAM, C.K.W & WEBB, J.C. *Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics*. J. of Computational Acoustics Dynamics, **107** (1993), 262-281.

[3] BOGEY, C & BAILLY, C. *Three-dimensional non-reflective boundary conditions for acoustic simulation: far field formulation and validation test cases*. Acta Acustica United With Acustica, **88** (2002), 463-471.

[4] BÉRENGER. *A perfectly matched layer for the absorption of electromagnetic waves*. J. of Computational Physics, **114** (1994), 185-200.

[5] BAILLY, C. *Computation of generation and radiation of acoustic waves in supersonic flows*. Third CAA Workshop on Benchmark Problems, NASA CP-2000-209790, (1999), 283-294.

[6] TAM, C.K.W, AURIAULT, L. & CAMBULI F. *Perfectly Matched Layer as an Absorbing Boundary Condition for the Linearized Euler Equations in Open and Ducted Domains*. J. of Computational Physics, **144** (1998), 213-234.

[7] HU, F. Q. *A stable Perfectly Matched Layer for Linearized Euler Equations in Unsplit Physical Variables*. J. of Computational Physics, **173** (2001), 455-480.

[8] HU, F. Q. *On Constructing Stable Perfectly Matched Layers as an Absorbing Boundary Condition for Euler Equations*. 40th AIAA Aerospace Sciences Meeting and Exhibit, (2004).

[9] BOGEY, C, BAILLY, C. & JUVÉ D. *A family of low dispersive and low dissipative schemes for Large Eddy simulations and for sound propagation*. J. of Computational Physics, **194** (2002), 194-214.

[10] J.C. HARDIN, J.R. RISTORCELLI AND C.K.W. TAM. *Workshop on benchmark problems in CAA*. NASA-CP-3300, (1995).

[11] MORRIS, P.J., AGARWAL, A. *Problem 1 - Radiation and refraction of sound waves through a two-dimensional shear layer*. NASA, ICASE, (2003), 21-22.