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

HAL
open science

# Formulations and numerical techniques for computation of acoustic pressure and its gradient by integral method

Gilles Rahier, Jean Prieur

► To cite this version:

Gilles Rahier, Jean Prieur. Formulations and numerical techniques for computation of acoustic pressure and its gradient by integral method. Journal of Sound and Vibration, 2025, 597 Part B, pp.118801. 10.1016/j.jsv.2024.118801. hal-04799062

HAL Id: hal-04799062

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

Submitted on 22 Nov 2024

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 Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

# Formulations and numerical techniques for computation of acoustic pressure and its gradient by integral method

Gilles Rahier<fnref n="1" />,*, Jean Prieur<fnref n="2" />

DAAA, ONERA, Institut Polytechnique de Paris, Châtillon, F-92322, France

## ARTICLE INFO

**Keywords:**

Computational aeroacoustics

Acoustic integral methods

Ffowcs Williams-Hawkings formulation

Pressure gradient

## ABSTRACT

Starting from the Ffowcs Williams-Hawkings surface integral formulation for a moving medium, the article proposes rather simple expressions for the radiated pressure and its gradient in the time domain, that are valid for solid or porous, fixed or moving integration surfaces. Moreover, these original expressions allow calculations with integration surfaces in supersonic motion (such as rotating surfaces around propellers or rotors). Versions dedicated to fixed integration surfaces are also proposed. The usual locally compact and the fully non-compact integration techniques are recalled, with, for both, a detailed description of efficient calculation algorithms. Particular attention is devoted to the order of precision of the calculations. The time derivation techniques and integration schemes used in this study lead to a theoretical second order that can easily be increased for the locally compact integration method. The results of these expressions and integration methods are compared to the analytical solution for the case of a fixed monopole and for that of a rotating monopole. They clearly show the benefit of the direct gradient calculation over a calculation by finite differences of the pressure around the observation point, particularly for broadband signals.

## 1. Introduction

The analysis of the noise generated by aerodynamic disturbances most often consists in calculating the acoustic pressure they radiate at various observation points, in particular in the far field. The Ffowcs Williams-Hawkings (FW-H) formulation [1] is a particularly well-suited tool for this purpose. This formulation, based on Lighthill's acoustic analogy [2] and extended to permeable surfaces [3], provides an explicit expression of the radiated pressure in the form of an integral over known quantities on a control surface encompassing the noise sources. The relationship between this surface formulation with the Kirchhoff method has been shown, for example, in [4,5].

In practice, the assumption of uniform flow and of free-field acoustic propagation outside the control surface is generally made in the Green's function used for the radiation calculation. However, in many cases, free-field radiation calculation is unsuitable, as the presence of reflecting bodies (centerbody of a propeller, aircraft fuselage or wing...) outside the control surface significantly alters the acoustic field.

A solution for taking these reflecting bodies into account is to couple the FW-H calculation with BEM methods [6–8] or with the

<footnote>* Corresponding author.
E-mail address: gilles.rahier@onera.fr (G. Rahier).</footnote>

<footnote n="1">¹ Research Scientist at Aerodynamics Aeroelasticity Acoustics Department.</footnote>

<footnote n="2">² Previously Research Scientist at CFD and Aeroacoustics Department.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 1. Schematic of the coordinate systems and integration surface S.</figcaption></figure>

Kirchhoff method [9]. In both cases, knowledge of the acoustic velocity or pressure gradient on the reflecting wall is required (the relationship between the acoustic velocity and the pressure gradient is described in [10]). Independently of couplings, these vector quantities are interesting in their own right, as they enable more comprehensive analyses of the acoustic field than the pressure field alone.

Expressions resulting from the FW-H formulation have been developed for the pressure gradient, in the time domain [11–14] as well as in the frequency domain [15,16]. These expressions are sometimes rather complicated [14], most are only valid for a propagation medium at rest [11–13,15] and none are applicable to control surfaces in supersonic motion (e.g. rotating control surfaces around rotors or propellers).

This article proposes, in the time domain, relatively simple expressions - from our point of view - for the radiated pressure and its gradient, for a moving medium, as well as efficient algorithms for fixed surfaces, for surfaces in subsonic motion and for surfaces in supersonic motion. A particular attention is devoted to the order of precision of the calculations.

The article is organized as follows. In Section 2, starting from the time-domain expression of pressure, original formulas for the pressure gradient are developed, valid whatever the kinematics of the integration surface. Formulas more suited to fixed surfaces are also proposed. Section 3 then describes in detail efficient integration techniques for each kinematics. The results of these expressions and these integration techniques are analyzed and compared with analytical solutions in Section 4. Finally, conclusions are drawn in Section 5.

## 2. Formulation

### 2.1. Expression of the radiated pressure

We define a reference frame $R_Y$ in which the unperturbed fluid is at rest, and a reference frame $R_Y$ in translation at the subsonic velocity $-U_0Y_1$ ($Y_1$: unit vector of the $y_1$ axis) with respect to $R_Y$ (Fig. 1). Reference frame $R_Y$ corresponds to wind tunnel or flight aeroacoustic simulations. In this reference frame, the unperturbed fluid is thus in uniform, steady translation at velocity $U_0Y_1$.

In the following, $S$ is a fictitious closed porous surface, fixed or moving in $R_Y$, and enclosing noise sources (Fig. 1). The observation points are assumed to be fixed in $R_Y$. The space coordinates $X_i$ of an observation point $X$ and $Y_i$ of an emission point $Y$ are defined in $R_Y$. The velocities $\mathbf{u}$ of the fluid and $\mathbf{v}$ of the surface $S$ are relative to $R_Y$.

Starting from the extension of the Fowcs Williams and Hawkings formulation to porous surfaces [3], the pressure radiated to observer at location $X$ at time $t$ by the acoustic sources inside surface $S$ can be written as [17,18]:

$$p'(\mathbf{X}, t) = -\frac{\partial}{\partial X_i} \int_{\tau} \int_{S} A_i G dS d\tau + \frac{d}{dt} \int_{\tau} \int_{S} BG dS d\tau \quad (1)$$

where $A_i = p'n_i + \rho u_i(u_n - v_n)$ (viscous term assumed to be negligible), $B = \rho_0 v_n + \rho(u_n - v_n)$, $G = \frac{\delta_1 g}{4\pi d}$ the free-space Green's function in $R_3$ for the convective wave equation, $\tau$ an emission time, $t$ a reception time and $\frac{d}{dt} = (\frac{\partial}{\partial t})_X + U_0 \frac{\partial}{\partial X_i}$, with $p' = p - p_0$ the pressure perturbation, $\rho_0$ the pressure in the unperturbed fluid, $\rho$ the density, $\rho_0$ the density in the unperturbed fluid, $u_n = \mathbf{u.n}$, $v_n = \mathbf{v.n}$, $\mathbf{n}$ the normal unit vector on the surface $S$, pointing outward, $\delta(f)$ the Dirac distribution, $g = \tau - t + \sigma/a_0$, $\sigma = (d - M_0(X_1 - Y_1))/\beta^2$, $d = \sqrt{C_i(X_i - Y_i)^2}$, $C_i = (1 - \beta^2)\delta_{1i} + \beta^2$, $\delta_{ij}$ the Kronecker symbol, $\beta^2 = 1 - M_0^2$, $M_0 = U_0/a_0$, $a_0$ the speed of sound in the unperturbed medium.

Moving the space derivatives inside the integrals leads to:

$$p'(\mathbf{X}, t) = -\int_{\tau} \int_{S} A_i \frac{\partial G}{\partial X_i} dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} BG dS d\tau + \int_{\tau} \int_{S} U_0 B \frac{\partial G}{\partial X_1} dS d\tau \quad (2)$$

Remarks:
- The time derivative is not moved inside the integral because the internal time derivative reinforces the numerical singularity present in the integral for a surface S in supersonic motion (see Section 3.3), and the problem due to this singularity is then no longer resolved by the fully non-compact integration technique described in Section 3.3.

- It is worth noting that Eq. (2) remains valid for a solid integration surface (e.g. the surface of a rotor or propeller blade). For such a surface, $u_n - v_n = 0$ and we find the well-known loading and thickness noises [19], respectively in the terms linked to $A_i = p^n_i$ and $B = \rho_0 v_n$. This remark applies to all the expressions which follow, for the radiated pressure as well as for its gradient. However, the thickness and loading terms do not have any real physical meaning in the more general framework of a porous surface, which is the focus of this study (Fig. 1). Therefore, they will not be discussed further.

Knowing that $\partial g / \partial t = -1$, the spatial derivative $\partial G / \partial X_i$ can write:

$$ \frac{\partial G}{\partial X_i} = -\frac{1}{4\pi d} \left( \frac{\partial(\delta(g))}{\partial t} \frac{\partial g}{\partial X_i} + \frac{\delta(g)}{d} \frac{\partial d}{\partial X_i} \right) $$

Therefore Eq. (2) becomes:

$$ \begin{aligned} p'(\mathbf{X}, t) = & \int_{\tau} \int_{S} \frac{A_i}{4\pi d} \frac{\partial g}{\partial X_i} \frac{\partial(\delta(g))}{\partial t} dS d\tau + \int_{\tau} \int_{S} \frac{A_i}{d} \frac{\partial d}{\partial X_i} \frac{\delta(g)}{4\pi d} dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} \frac{B}{4\pi d} \frac{\delta(g)}{d} dS d\tau \\ & - \int_{\tau} \int_{S} \frac{U_0 B}{4\pi d} \frac{\partial g}{\partial X_1} \frac{\partial(\delta(g))}{\partial t} dS d\tau - \int_{\tau} \int_{S} \frac{U_0 B}{d} \frac{\partial d}{\partial X_1} \frac{\delta(g)}{4\pi d} dS d\tau \end{aligned} \quad (3) $$

For fixed observation points, d and $\partial g / \partial X_i$ do not depend on the reception time t. Putting the time derivative outside the integrals gives:

$$ p'(\mathbf{X}, t) = \int_{\tau} \int_{S} \frac{1}{d} \left( A_i \frac{\partial d}{\partial X_i} - U_0 B \frac{\partial d}{\partial X_1} \right) G dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} \left( A_i \frac{\partial g}{\partial X_i} + B \left( 1 - U_0 \frac{\partial g}{\partial X_1} \right) \right) G dS d\tau \quad (4) $$

By developing $\partial g / \partial X_i = \frac{1}{a_0 \beta^2} \left( \frac{C_i(X_i - Y_i)}{d} - M_0 \delta_{1i} \right)$ and $\partial d / \partial X_i = \frac{C_i(X_i - Y_i)}{d}$, Eq. (4) transforms into:

$$ \begin{aligned} p'(\mathbf{X}, t) = & \int_{\tau} \int_{S} \left( \frac{A_i C_i (X_i - Y_i) - U_0 B (X_1 - Y_1)}{d^2} \right) G dS d\tau \\ & + \frac{\partial}{\partial t} \int_{\tau} \int_{S} \left( \frac{A_i}{a_0 \beta^2} \left( \frac{C_i (X_i - Y_i)}{d} - M_0 \delta_{1i} \right) + B \left( 1 - \frac{U_0}{a_0 \beta^2} \left( \frac{X_1 - Y_1}{d} - M_0 \right) \right) \right) G dS d\tau \end{aligned} \quad (5) $$

Knowing that $A_i C_i (X_i - Y_i) = \beta^2 A_i (X_i - Y_i) + M_0^2 A_1 (X_1 - Y_1)$ and $1 + U_0 M_0 / a_0 \beta^2 = 1 / \beta^2$, Eq. (5) writes:

$$ p'(\mathbf{X}, t) = \int_{\tau} \int_{S} F G dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} F G dS d\tau \quad (6) $$

where $F = \frac{1}{d^2} (\beta^2 A_i (X_i - Y_i) + (M_0^2 A_1 - U_0 B) (X_1 - Y_1))$

$$ F' = \frac{1}{a_0} \left( A_i \frac{(X_i - Y_i)}{d} + \frac{M_0 A_1}{\beta^2} \left( M_0 \frac{(X_1 - Y_1)}{d} - 1 \right) \right) + \frac{B}{\beta^2} \left( 1 - M_0 \frac{(X_1 - Y_1)}{d} \right) $$

We finally arrive to the following expressions of the terms F and F':

$$ F = \frac{\Sigma_1^2}{d^2} \Sigma_1, F' = \frac{\Sigma_1}{a_0 d} + \frac{\Sigma_2}{\beta^2} \text{ with } \Sigma_1 = A_i(X_i - Y_i) + \frac{(M_0^2 A_1 - U_0 B)}{\beta^2} (X_1 - Y_1) \text{ and } \Sigma_2 = B - \frac{M_0 A_1}{a_0} $$

## 2.2. Expression of the radiated pressure gradient

Starting from Eq. (6), the pressure gradient can write:

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{\tau} \int_{S} \frac{\partial F}{\partial X_j} G dS d\tau + \int_{\tau} \int_{S} F \frac{\partial G}{\partial X_j} dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} \frac{\partial F}{\partial X_j} G dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} F \frac{\partial G}{\partial X_j} dS d\tau \quad (7) $$

As done for the terms $A_i \frac{\partial G}{\partial X_i}$ and $U_0 B \frac{\partial G}{\partial X_1}$ in Eq. (2) of the pressure:

$$ \int_{\tau} \int_{S} F \frac{\partial G}{\partial X_j} dS d\tau = - \int_{\tau} \int_{S} \frac{1}{d} \frac{\partial d}{\partial X_j} F G dS d\tau - \frac{\partial}{\partial t} \int_{\tau} \int_{S} \frac{\partial d}{\partial X_j} F G dS d\tau \text{ and } \int_{\tau} \int_{S} F' \frac{\partial G}{\partial X_j} dS d\tau = - \frac{\partial}{\partial t} \int_{\tau} \int_{S} \frac{1}{d} \frac{\partial d}{\partial X_j} F' G dS d\tau - \frac{\partial^2}{\partial t^2} \int_{\tau} \int_{S} \frac{\partial d}{\partial X_j} F' G dS d\tau $$

Therefore, the pressure gradient writes:
$$ \frac{\partial p'(\mathbf{X}, t)}{\partial \mathbf{X}_j} = \int_{\tau} \int_{S} \left( \frac{\partial F}{\partial \mathbf{X}_j} - \frac{1}{d} \frac{\partial d}{\partial \mathbf{X}_j} F \right) G dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} \left( \frac{\partial F}{\partial \mathbf{X}_j} - \frac{1}{d} \frac{\partial d}{\partial \mathbf{X}_j} F' - \frac{\partial g}{\partial \mathbf{X}_j} F \right) G dS d\tau - \frac{\partial^2}{\partial t^2} \int_{\tau} \int_{S} \frac{\partial g}{\partial \mathbf{X}_j} F' G dS d\tau \quad (8) $$

Noting $\hat{r}_i = \frac{x_i - y_i}{d}$:

$$ \frac{\partial F}{\partial \mathbf{X}_j} = \frac{\beta^2}{d^2} \left( \frac{\partial \Sigma_1}{\partial \mathbf{X}_j} - \frac{2\Sigma_1}{d} C_j \hat{r}_j \right) $$

$$ \frac{\partial d}{\partial \mathbf{X}_j} \frac{F}{d} = C_j \hat{r}_j \frac{\beta^2}{d^3} \Sigma_1 $$

$$ \frac{\partial F'}{\partial \mathbf{X}_j} = \frac{1}{a_0 d} \frac{\partial \Sigma_1}{\partial \mathbf{X}_j} - \frac{\Sigma_1}{a_0 d^2} C_j \hat{r}_j + \frac{1}{\beta^2} \frac{\partial \Sigma_2}{\partial \mathbf{X}_j} $$

$$ \frac{\partial d}{\partial \mathbf{X}_j} \frac{F}{d} = C_j \hat{r}_j \left( \frac{\Sigma_1}{a_0 d^2} + \frac{\Sigma_2}{\beta^2 d} \right) $$

$$ \frac{\partial g}{\partial \mathbf{X}_j} F = \left( \frac{C_j \hat{r}_j - M_0 \delta_{1j}}{a_0 d^2} \right) \Sigma_1 $$

$$ \frac{\partial g}{\partial \mathbf{X}_j} F' = \left( \frac{C_j \hat{r}_j - M_0 \delta_{1j}}{a_0 \beta^2} \right) \left( \frac{\Sigma_1}{a_0 d} + \frac{\Sigma_2}{\beta^2} \right) $$

Since $A_i$ and $B$ do not depend on $\mathbf{X}_j$, $\frac{\partial \Sigma_2}{\partial \mathbf{X}_j} = 0$.

Grouping terms and noting $\Sigma'_j = \frac{\partial \Sigma_1}{\partial \mathbf{X}_j} - 3C_j \hat{r}_j \frac{\Sigma_3}{d} = A_j + \frac{(M_0^2 A_1 - U_0 B)}{\beta^2} \delta_{1j} - 3C_j \hat{r}_j \frac{\Sigma_1}{d}$, the pressure gradient finally writes:

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial \mathbf{X}_j} = \int_{\tau} \int_{S} F_{1j} G dS d\tau + \frac{\partial}{\partial t} \int_{\tau} \int_{S} F_{2j} G dS d\tau + \frac{\partial^2}{\partial t^2} \int_{\tau} \int_{S} F_{3j} G dS d\tau \quad (9) $$

$$ F_{1j} = \beta^2 \frac{\Sigma'_j}{d^2} $$

$$ F_{2j} = \frac{\Sigma'_j}{a_0 d} + M_0 \delta_{1j} \frac{\Sigma_1}{a_0 d^2} - \frac{C_j \hat{r}_j \Sigma_2}{\beta^2 d} $$

$$ F_{3j} = - \frac{(C_j \hat{r}_j - M_0 \delta_{1j})}{a_0 \beta^2} \left( \frac{\Sigma_1}{a_0 d} + \frac{\Sigma_2}{\beta^2} \right) $$

## 2.3. Other expressions for fixed integration surfaces

Moving the time derivatives outside the integrals is a necessary condition for the space-time integration method designed for integration surfaces in supersonic motion (Section 3.3). On the other hand, it has the disadvantage of amplifying high-frequency oscillations due to numerical calculation inaccuracies. For fixed integration surfaces in the reference frame $R_Y$ (Fig. 1), using expressions with inner time derivatives does not lead to more complicated expressions nor numerical singularities, as shown in the following. Applications have shown that they lead to more regular signals than those obtained with outer time derivatives.

Moving the time derivative inside the integral, Eq. (6) becomes:

$$ p'(\mathbf{X}, t) = \int_{\tau} \int_{S} F G dS d\tau + \int_{\tau} \int_{S} \frac{\partial}{\partial t} (F' G) dS d\tau \quad (10) $$

For fixed surface $S$ and fixed observation points, $F'$ does not depend on the reception time $t$ and $\frac{\partial G}{\partial t} = - \frac{\partial G}{\partial \tau}$.

$$ \text{So: } \int_{\tau} \int_{S} \frac{\partial}{\partial t} (F' G) dS d\tau = - \int_{\tau} \int_{S} F' \frac{\partial G}{\partial \tau} dS d\tau = \int_{\tau} \int_{S} \frac{\partial F}{\partial \tau} G dS d\tau $$

Therefore, the pressure can be written in the following form:

$$ p'(\mathbf{X}, t) = \int_{\tau} \int_{S} \left( F + \frac{\partial F}{\partial \tau} \right) G dS d\tau \quad (11) $$

In the same way, from Eq. (9), the pressure gradient can be written in the following form:
$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{\tau} \int_{S} \left( F_{1j} + \frac{\partial F_{2j}}{\partial \tau} + \frac{\partial^2 F_{3j}}{\partial \tau^2} \right) G dS d\tau \quad (12) $$

### 3. Integration techniques

#### 3.1. Integration technique well-suited to fixed integration surfaces

For a fixed integration surface, $\partial g / \partial \tau = 1$ and changing the integration variable $\tau$ in g in Eqs. (11) and (12) leads to the following expressions for the pressure and its gradient:

$$ p'(\mathbf{X}, t) = \int_{g} \int_{S} \left( F + \frac{\partial F}{\partial \tau} \right) \frac{\partial(g)}{4\pi d} dS dg \quad (13) $$

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{g} \int_{S} \left( F_{1j} + \frac{\partial F_{2j}}{\partial \tau} + \frac{\partial^2 F_{3j}}{\partial \tau^2} \right) \frac{\partial(g)}{4\pi d} dS dg \quad (14) $$

That is, finally:

$$ p'(\mathbf{X}, t) = \int_{S} \frac{1}{4\pi d} \left[ F + \frac{\partial F}{\partial \tau} \right]_{\text{ret}} dS \quad (15) $$

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{S} \frac{1}{4\pi d} \left[ F_{1j} + \frac{\partial F_{2j}}{\partial \tau} + \frac{\partial^2 F_{3j}}{\partial \tau^2} \right]_{\text{ret}} dS \quad (16) $$

where the quantities inside the brackets have to be evaluated for $g = 0$, that is at the emission time $\tau = t - \sigma / a_0$ (called retarded time).

Usually, the integration surface S is discretized in source elements $\delta S$ treated as point sources. This locally compact approach, referred to as compact grid cells (CGC) approach in the following, considerably simplifies the spatial integration. It then comes down to a sum over the cells of the mesh, of quantities calculated at their center. A fast method is to start from the emission time and proceed forward in time. There is no need then to calculate any retarded quantities. The values at times $t_i$ used to describe the received signals can be determined by linear interpolation from the values obtained for the reception times $\tau_j + \sigma / a_0$ resulting of the emission at times $\tau_j$. A particularly efficient algorithm consists of describing the received signals with the same time step as that of the emission and distributing the quantity $Q = [\cdot] \frac{\delta S}{4\pi d}$ calculated for each emission time $\tau_j$ over the two reception time $t_k$ and $t_{k+1}$ on either side of the reception time $t = \tau_j + \sigma / a_0$, using the following simple operations:

$$ p'(\mathbf{X}, t_k) \leftarrow p'(\mathbf{X}, t_k) + \frac{t_{k+1} - t}{t_{k+1} - t_k} Q $$

$$ p'(\mathbf{X}, t_{k+1}) \leftarrow p'(\mathbf{X}, t_{k+1}) + \frac{t - t_k}{t_{k+1} - t_k} Q $$

This time integration scheme, based on linear interpolations, is thus of second order. In this study, the first and second derivatives in time are calculated by second-order centered finite differences. Therefore, time derivation and integration schemes are consistent with each other, and the complete time scheme is of second order, as confirmed by the results presented in section 4.1. The space integration then consists of summing the contributions of each cell of surface S to the pressure and its gradient at the observation point.

#### 3.2. Integration technique well-suited to subsonically moving integration surfaces

Changing the integration variable $\tau$ in g in Eqs. (6) and (9), leads to the following expressions:

$$ p'(\mathbf{X}, t) = \int_{g} \int_{S} \frac{F}{|\partial g / \partial \tau|} \frac{\partial(g)}{4\pi d} dS dg + \frac{\partial}{\partial t} \int_{g} \int_{S} \frac{F}{|\partial g / \partial \tau|} \frac{\partial(g)}{4\pi d} dS dg \quad (17) $$

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{g} \int_{S} \frac{F_{1j}}{|\partial g / \partial \tau|} \frac{\partial(g)}{4\pi d} dS dg + \frac{\partial}{\partial t} \int_{g} \int_{S} \frac{F_{2j}}{|\partial g / \partial \tau|} \frac{\partial(g)}{4\pi d} dS dg + \frac{\partial^2}{\partial t^2} \int_{g} \int_{S} \frac{F_{3j}}{|\partial g / \partial \tau|} \frac{\partial(g)}{4\pi d} dS dg \quad (18) $$

That is, finally:

$$ p'(\mathbf{X}, t) = \int_{S} \left[ \frac{F}{4\pi d |\partial g / \partial \tau|} \right]_{\text{ret}} dS + \frac{\partial}{\partial t} \int_{S} \left[ \frac{F'}{4\pi d |\partial g / \partial \tau|} \right]_{\text{ret}} dS \quad (19) $$
<figure><img src="image_6.png" /><figcaption>Fig. 2. Schematic of the time integration process for a subsonically moving point source.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Schematic of the time integration process for a non-compact source element $\delta S$, for two pressure rate shape functions: (a) rectangular and (b) triangular.</figcaption></figure>

$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_S \left[ \frac{F_{1j}}{4\pi d|\partial g/\partial \tau|} \right]_{\text{ret}} dS + \frac{\partial}{\partial t} \int_S \left[ \frac{F_{2j}}{4\pi d|\partial g/\partial \tau|} \right]_{\text{ret}} dS + \frac{\partial^2}{\partial t^2} \int_S \left[ \frac{F_{3j}}{4\pi d|\partial g/\partial \tau|} \right]_{\text{ret}} dS \quad (20) $$

$$ \text{with } \frac{\partial g}{\partial \tau} = 1 + \frac{1}{a_0 \beta^2} \left( M_0 v_1 - \frac{C_i(X_i - Y_i)}{d} v_i \right). $$

$$ \frac{\partial g}{\partial \tau} \text{ can be expressed in the form } 1 - M_R \text{ with } M_R = \frac{1}{a_0 \beta^2} \left( \frac{C_i(X_i - Y_i)}{d} v_i - M_0 v_1 \right). $$

$M_R$ can be interpreted as the 'acoustic Mach number' of the center Y of the cell in the direction of the observation point X, taking into account the velocity $U_0$ of the flow. For a surface S in subsonic motion relative to the observer location X, $\partial g / \partial \tau$ is never null ($M_R < 1$). The CGC approach described in the previous section still makes it possible to reduce spatial integration to a simple summation over the cells of the mesh of quantities evaluated at their center.

For each of the two pressure terms and three terms of its gradient, the contributions $Q_k$ received at a reception time $t_k$ can still be determined by linear interpolation from the values $Q_j$ and $Q_{j+1}$ obtained at the reception times $t_j = \tau_j + \sigma_j / a_0$ and $t_{j+1} = \tau_{j+1} +$
<figure><img src="image_8.png" /><figcaption>Fig. 4. Example of evolutions of the reception time interval between two emission times $\tau_j$ and $\tau_{j+1}$ and example of integration bounds $\tau_a$ and $\tau_b$, for subsonic, transonic and supersonic kinematics of the emission cell.</figcaption></figure>

$\sigma_{j+1}/a_0$ on both sides of $t_k$ (Fig. 2). The calculation is, however, a little longer than for a fixed integration surface because the propagation times (hence the interpolation coefficients) linked to an emission point are no longer constant because of its movement. Note that the particular summation algorithm described in Section 3.1 does not work for moving integration surfaces, because the reception times $t_j$ are no longer equidistant and several reception times $t_k$ can be included in a reception interval $[t_j, t_{j+1}]$.

Once all the spatial summations performed, the time derivatives can be calculated by finite differences (second-order centered differences in the present study) and the contributions of the different terms are added to provide the final result.

### 3.3. Integration technique well-suited to supersonically moving integration surfaces

For surfaces in supersonic motion, the integrals in Eqs. (19) and (20) become singular because the term $\partial g/\partial \tau$ can equal zero for certain observation points. They can no longer be calculated numerically using the usual CGC approach (cf. Sections 3.1 and 3.2). A very interesting study has been carried out by Ianniello to solve the problem raised by this singularity [20,21], but the calculation cost, due to the complexity of the method, seems excessive for industrial applications. This problem, together with an original alternative solution, has already been presented in previous articles and a paper [5,18,22]. To sum up, the usual CGC approach does not reflect the continuity of the integration surface – especially the fact that at a given emission time, the different points of a cell are not all at the same distance from the observer. In fact, because of these different distances, the emission of a cell at time $\tau$ is not received at a single instant $t$, but during a time interval $[t_1, t_2]$ of which boundaries depend on the size, shape and orientation of the cell. Therefore, the CGC approach leads to focusing the radiation of the whole cell onto the same reception time. This phenomenon of acoustic focusing is accentuated when the velocity of the source point towards the observation point tends towards the speed of sound. The radiation then "piles up" even more on the same reception time, which explains the divergence in the results for the transonic moving cells. Rather than modeling the emission of the cell at the time $\tau$ by a point source, a relatively simple way of taking into account the continuity of the integration surface consists in modeling the reception of the radiated signal at the observation point by a pressure rate $\Pi$ (in Pa s⁻¹) over the time interval $[t_1, t_2]$. This is what is performed in the so-called fully non-compact (FNC) integration technique, which does not use the analytical time integration leading to Eqs. (19) and (20), but restarts from Eqs. (6) and (9).

In practice, the reception intervals $[t_1, t_2]$ are determined for each cell $\Delta_i$ of the integration surface and for each emission time $\tau_j$, starting from the position of the observation point and the position of the nodes of the cell and of its center, at time $\tau_j$. The associated radiated pressure rate $\Pi_{ij}$ is calculated as a function of the considered source term $F$ (which can be evaluated at the center of the cell), of the surface $|\Delta_i|$ of the cell, and of the reception time duration $t_2-t_1$. For instance, if a constant pressure rate is considered, the pressure rate $\Pi_{ij}$ is then equal to $F|\Delta_i|/(t_2-t_1)$ over the time interval $[t_1, t_2]$ and to zero outside. The pressure received at the observation point at a given instant $t$ is then determined by summing this pressure rate over the emission times (see Fig. 3). This time integration can be quite simple (depending on the pressure rate model), since it consists in calculating the area of the hatched surface in Fig. 3 by linear interpolations starting from the reception time intervals and the pressure rates determined for the emission times $\tau_j$ and $\tau_{j+1}$. Examples of integration bounds $\tau_a$ and $\tau_b$ are given in Fig. 4. They show the effect of the cell’s kinematics in the summation of the pressure rate received at the observation point (integration range $[\tau_a, \tau_b]$ greatly increased for $M_R \approx 1$). The space integration is then carried out by summing the contributions of each cell of the mesh, for each integral of Eq. (6) and Eq. (9). The rest of the calculation process (time derivations and terms summation) is the same as for the previous method (section 3.2). It is worth pointing out that this integration technique gives the same result as the usual CGC method when the mesh is fixed. The FNC approach has been thoroughly validated [22] and is used by ONERA and by its industrial partners for rotors and propellers acoustic studies (for example [23,24]).

From a numerical point of view, calculating the area of the hatched surface in Fig. 3 requires logic tests to precisely determine its
geometry as a function of the reception time $t$. This process is time-consuming, as shown by the comparison of calculation times in Table 1. To facilitate this calculation, a pressure rate that is constant over the reception interval, and therefore discontinuous at the boundaries of this interval, was initially used (Fig. 3a). This very simple model has been fully satisfactory in the specific case of rotating meshes describing axisymmetric integration surfaces ([5,22], FNC integration technique applied to formulations with inner time derivatives, see Section 3.4). A continuous increasing/decreasing model (Fig. 3b) was then designed for any moving mesh. It is twice as computationally expensive because the hatched surface is more complicated to determine, but it leads to much regular results for formulations with time derivatives outside the integral [18]. In this study, all acoustic calculations are carried out using this second model.

For a cell, this time integration scheme, which is based on linear interpolations, is therefore of the second order, as the previous schemes, but the modeling of the radiated pressure rate introduces inaccuracies which reduce the actual order of complete calculations (see section 4.3). After summation on the cells, none of the two models presented leads to signals that are regular enough to prevent the external time derivatives – calculated by second-order centered finite differences – from being slightly noisy at the sampling frequency and a few subharmonics. These spurious high frequencies remain at an acceptable level, even for surfaces in supersonic motion, however they reduce or stop the decrease of the error – in $L_\infty$ norm – with that of the time step.

### 3.4. Particular case of fixed integration surfaces described by moving meshes

It can be shown that the expressions established for fixed surfaces (Section 2.3) remain valid in the special case of a moving mesh describing, over time, a geometrically identical surface [5]. This situation occurs, for example, with a rotating mesh describing a surface of revolution around the axis of rotation of a rotor or propeller. The FNC integration technique can then be applied to expressions with inner time derivatives (Eqs. (11) and (12)), using the simplest pressure rate model (rectangular shape function, Fig. 3a). This approach has been used to account for volume sources in the noise radiated by a transonic rotor without numerical difficulties, using rotating cylindrical surfaces of large radii to include these sources, therefore with supersonically moving meshes [5].

To conclude this part of the article, a summary of the formulas and integration techniques is provided in the Appendix. This Appendix indicates the formulas and integration methods best suited to each configuration, but they are not exclusive. For example, Eq. (A.5) and Eq. (A.6) with the FNC integration method can be applied to surfaces in subsonic motion or to fixed surfaces, although the calculations are unnecessarily more expensive. Table 1 provides an example of calculation times as a function of the surface kinematics and integration technique.

**Table 1**

Calculation time per observation point, on a processor running at 3 GHz, for an integration surface described by 50,000 cells and for 1000 emission times.

<table><thead><tr><th>Integration tech.</th><th>CGC (S fixed)</th><th>CGC (S moving)</th><th>FNC (rectangular pressure rate model)</th><th>FNC (triangular pressure rate model)</th></tr></thead><tbody><tr><td>Pressure</td><td>2.5 s</td><td>3.5 s</td><td>7 s</td><td>14 s</td></tr><tr><td>Pressure gradient</td><td>3 s</td><td>5 s</td><td>10 s</td><td>20 s</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Instantaneous pressure field on the integration surface. $M_0 = 0.5$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 6. Pressure time histories at the observation point (0, 10, 0). $M_0 = 0.5$, $dt = T/97$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Pressure error ($L_\infty$ norm) as a function of the time step $dt$, for the observation point (10, 10, 10).</figcaption></figure>

# 4. Validations

## 4.1. Formulas with inner time derivatives and CGC integration method

The expressions and integration methods described in Section 3 are first validated for formulas with inner time derivatives and the CGC integration method (Eqs. (15) and (16), Section 3.1), which are well-suited for fixed integration surfaces.

The configuration is as follows. The integration surface is a fixed sphere of radius $R = 1$ m represented by a uniform unstructured mesh composed of 71,060 triangles. The acoustic field is generated by a monopolar point source of amplitude $A$ and frequency $f$ in a uniform flow at velocity $U_0$ in the $X_1$ direction, placed at the sphere's center. The complex potential for the monopole, at a node $X$ of the integration surface and at an instant $t$, is [25,26]:

$$ \varphi(\mathbf{X}, \mathbf{Y}, t) = \frac{A}{4\pi d} \exp^{(2\pi f)(t-\sigma/a_0)} \text{ where } A \text{ is its strength and } \mathbf{Y} = (0, 0, 0). \text{ (See notations in Section 2.1)} $$

The inputs to the FW-H calculation are the real parts of the velocities $u'_i = \partial\varphi / \partial X_i$, pressure $p' = -\rho_0 (\frac{\partial\varphi}{\partial t} + U_0 \frac{\partial\varphi}{\partial X_1})$ and density $\rho' = p' / a_0^2$ disturbances. Their developed expressions are:

$$ p' = \frac{A\rho_0}{4\pi d} \left[ 2\pi f \left( 1 - \frac{U_0}{a_0} \frac{\partial\sigma}{\partial X_1} \right) \sin(2\pi f(t - \sigma / a_0)) + \frac{U_0}{d} \frac{\partial d}{\partial X_1} \cos(2\pi f(t - \sigma / a_0)) \right] \quad (21) $$
<figure><img src="image_12.png" /><figcaption>Fig. 8. Pressure gradient time histories at the observation point (0, 10, 0). M₀ = 0.5, dt = T/97.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 9. Pressure gradient error (L∞ norm) as a function of h in calculation by finite differences, and indication of the error level of the direct calculation (dashed lines). Observation point (10, 10, 10), dt = T/97.</figcaption></figure>

$$ \rho' = \frac{p'}{a_0^2} \qquad (22) $$

$$ u'_i = \frac{A}{4\pi d} \left[ \frac{2\pi f}{a_0} \frac{\partial \sigma}{\partial X_i} \sin(2\pi f(t - \sigma / a_0)) - \frac{1}{d} \frac{\partial d}{\partial X_i} \cos(2\pi f(t - \sigma / a_0)) \right] \qquad (23) $$

Calculations are performed for a₀ = 1 m s⁻¹, p₀ = 1 Pa, f = 1 Hz. For a fluid at rest (M₀ = 0), the wave length λ₀ is thus equal to unity. A is taken equal to 1/7 10⁻² m³ s⁻¹ so that, for a flow at rest, the amplitude of the pressure disturbance p' is equal to 10⁻⁴ Pa at 10 m from the point source. An example of instantaneous pressure field on the integration surface is given in Fig. 5. For the pressure, the results at an observation point (X₁, X₂, X₃) are compared to the analytical solution Eq. (21). For the pressure gradient, the reference solution is the pressure gradient obtained by centered finite differences of the analytical solution Eq. (21) for 8 points around the observation point with a very small space step h = 0.002 m (λ₀/500), identical in the three directions: (X₁ ± h/2, X₂ ± h/2, X₃ ± h/2). The deviation from the reference solution is quantified by L∞ norm divided by the maximum absolute value of the reference solution.

**Results for the pressure** - A comparison of FW-H calculation (Eq. (15)) with the analytical solution is given in Fig. 6 for the observation point (0, 10, 0), a flow Mach number M₀ = 0.5 and a time step dt = T/97 (T = 1/f). The difference between the signals is not
<figure><img src="image_14.png" /><figcaption>Fig. 10. Rotating point monopole, integration surfaces S<sub>0.7</sub>, S<sub>0.9</sub>, S<sub>1.2</sub>, S<sub>1.6</sub>, S<sub>2.0</sub> and observation point.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. Instantaneous pressure fields on the integration surface S<sub>0.7</sub>. (a) M<sub>0</sub> = 0 and (b) M<sub>0</sub> = 0.5.</figcaption></figure>

discernible on this plot. The evolution of the deviation from the analytical solution as a function of the time step *dt* (from *T*/5 to *T*/*97*) of the FW-H calculation is plotted Fig. 7 for the more general case of an observation point out of any plane of symmetry. This evolution provides more precise information on the precision of the calculations. In particular, it confirms that the time scheme is of second order, with or without flow.

**Results for the pressure gradient** - For the observation point (0, 10, 0) and identical parameters M₀ = 0.5, dt = T/97, the direct calculation using Eq. (16) provides the same signals as those of the reference solution, for all components of the pressure gradient (Fig. 8). For the observation point (10, 10, 10), the error (L∞ norm on the three components) of the gradient calculation by finite differences starting from the pressures provided by Eq. (15), is plotted as a function of the spatial step h, in Fig. 9. We find fairly well the second-order precision of the calculation by centered finite differences. The error level of the direct calculation by Eq. (16) clearly highlights the advantage of this approach for fixed integration surfaces. Note that deviation of FW-H direct calculation with respect to the reference solution depends on the time step *dt*, as for the pressure (Fig. 7). This deviation would therefore be lower with a time step smaller than T/97.

## 4.2. Formulas with outer time derivatives and CGC integration method

The expressions and integration methods described in Section 3 are now validated for formulas with outer time derivatives and the CGC integration method (Eqs. (19) and (20), Section 3.2), which are well-suited for subsonically moving integration surfaces.

The acoustic field is created by a point monopole placed in the plane X₁ = 0 and rotating around the X₁ axis at the frequency *f* on a circle of radius *r* = 1 m (Fig. 10). Starting from the general expressions established by Golstein for moving sources [27], the velocities, pressure and density disturbances (*u*<sub>i</sub>' = *u*<sub>i</sub> - *U*<sub>0</sub>*δ*<sub>1b</sub>, *p*' = *p* - *p*<sub>0</sub> and *ρ*' = *ρ* - *ρ*<sub>0</sub>, respectively) radiated by this rotating monopolar point source in a uniform flow at velocity *U*<sub>0</sub> in the X₁ direction, can be written:

$$u'_{1} = \frac{A}{4\pi d^2} \left[ (1-c) \frac{(X_1 - Y_1) - aM_0/a_0}{d} + M_0 c \right] \quad (24)$$
<figure><img src="image_16.png" /><figcaption>Fig. 12. Pressure time histories (CGC integration method). $dt = T/120$. (a) $M_0 = 0$ and (b) $M_0 = 0.5$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Pressure error ($L_\infty$ norm) as a function of the time step $dt$. Surface $S_{0.7}$, $M_0 = 0$.</figcaption></figure>

$$ u'_{2} = \frac{A}{4\pi d^2} \left[ (1-c) \frac{\beta^2(X_2 - Y_2)}{d'} - \frac{V_{e2}}{a_0} \right] \qquad (25) $$

$$ u'_{3} = \frac{A}{4\pi d^2} \left[ (1-c) \frac{\beta^2(X_3 - Y_3)}{d'} - \frac{V_{e3}}{a_0} \right] \qquad (26) $$

$$ p' = \frac{A\rho_0}{4\pi d^3} (\beta^2 a - db) - \rho_0 U_0 u'_{1} \qquad (27) $$

$$ \rho' = \frac{p'}{a_0^2} \qquad (28) $$

where $d = d - \frac{v_e \times}{a_0}$, $a = X_2 V_{e2} + X_3 V_{e3}$, $b = \frac{2\pi f}{a_0} (X_2 V_{e3} - X_3 V_{e2})$, $c = \frac{b}{\beta^2 a_0}$ and $V_e = \Omega \times Y = \begin{pmatrix} 2\pi f \\ 0 \\ 0 \end{pmatrix} \times \begin{pmatrix} Y_1 \\ Y_2 \\ Y_3 \end{pmatrix}$

Like in the previous section, the sound velocity $a_0$ and the pressure $p_0$ are respectively equal to 1 m s$^{-1}$ and 1 Pa. The rotation frequency $f$ of the monopole is taken equal to $1/4\pi$ so that its rotation Mach number $M_{\text{rot}} = 2\pi f / a_0$ is equal to 0.5. For a fluid at rest ($M_0 = 0$), the wave length $\lambda_0$ is thus equal to $4\pi$. A is taken equal to 0.0127 m$^3$ s$^{-1}$ so that, for a flow at rest, the amplitude of the
<figure><img src="image_18.png" /><figcaption>Fig. 14. Pressure gradient time histories (CGC integration method). Surface $S_{0.7}$, $M_0 = 0$, $dt = T/120$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15. Pressure gradient time histories (CGC integration method). $dt = T/120$. (a) $M_0 = 0$ and (b) $M_0 = 0.5$.</figcaption></figure>

pressure disturbance $p'$ is equal to $10^{-4}$ Pa at 10 m from the sphere center, in the plane of rotation of the monopole.

The integration surface is a sphere of radius $R'$, obtained by scaling the mesh of the sphere from the previous section. It rotates around the $X_1$ axis at the same speed as the monopole. Therefore the rotation Mach number depends on its radius. It is supersonic for radii $R'$ greater than 2. It can be noted that the perturbation fields are stationary in the rotating frame of the sphere. Five scaling factors from 1.4 to 4 are considered - i.e. 5 rotation Mach numbers from 0.7 to 2 - for the integration surface (Fig. 10) and two flow Mach numbers are considered: $M_0 = 0$ and $M_0 = 0.5$. Examples of instantaneous pressure fields on an integration surface are given in Fig. 11.

All the computations are performed for the observation point (0, 10, 0) in the rotation plan of the monopole (Fig. 10). Results for the pressure are compared to the analytical solution Eq. (27). For the pressure gradient, as done in the previous section, the reference solution is the pressure gradient obtained by centered finite differences of the analytical solution Eq. (27) for 8 points around the observation point with a very small space step $h = 0.02$ m ($\lambda_0/628$), identical in the three directions.

For the observation point considered, the 'acoustic Mach number' $M_R$ in Eqs. (19) and (20) remains subsonic only for the integration surfaces $S_{0.7}$ and $S_{0.9}$ for $M_0 = 0$, and only for the surface $S_{0.7}$ for $M_0 = 0.5$. Therefore, the results presented below concern only these three cases.

**Results for the pressure** – Calculations for the three cases mentioned above are first performed with a time step equal to $T/120$. Despite this moderate number of time steps to describe one monopole revolution, the pressure time histories provided by Eq. (19) compared very well to the analytical solution Eq. (27) (Fig. 12). The evolution of the error as a function of the time step (from $T/90$ to
<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Pressure gradient error ($L_\infty$ norm) as a function of $h$ in calculation by finite differences, and error level of the direct calculation. Surface $S_{0.7}$, $M_0 = 0$, $dt = T/360$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Pressure time histories as a function of the integration surface (FNC integration method). $dt = T/120$. (a) $M_0 = 0$ and (b) $M_0 = 0.5$.</figcaption></figure>

$T/720$) of the FW-H calculation is plotted in Fig. 13, for the integration surface $S_{0.7}$ and a fluid at rest. The order of precision is slightly lower than the second order expected for the complete temporal scheme (derivation and integration). This is due to the integration scheme. Indeed, because of the source motion, the interpolation described in section 3.2 is applied to non-equidistant reception times (Fig. 2), which tends to reduce the order of precision of the calculations. The plot also shows that the error presents a minimum value (close to $10^{-3}$ in the present case) linked to the space discretization of the integration surface.

**Results for the pressure gradient** – The three components of the pressure gradient are compared to the reference solution in Fig. 14, for the surface $S_{0.7}$ and a fluid at rest. This comparison is also made in Fig. 15 for the surface $S_{0.9}$ and for $M_0 = 0.5$, for the predominant component ($dp'/dx_2$) for the observation point considered. The results of the FW-H direct calculation by Eq. (20) are in excellent agreement with the reference solution. For the surface $S_{0.7}$ and a fluid at rest, the error ($L_\infty$ norm on the three components) of the gradient calculation by finite differences, starting from the pressures provided by Eq. (19), is plotted as a function of the spatial step $h$ in Fig. 16. The second-order precision expected from the space-centered finite difference calculations is not completely achieved. As with pressure, the slight deviation from the theoretical second-order precision can be attributed to the time integration scheme (inaccuracies due to non-equidistant reception times that depend on the observation point, hence on $h$). This influence of the time scheme on the precision of a gradient is not surprising for an acoustic propagation calculation. Furthermore, we can notice that the error of direct calculation can be slightly higher than that of calculation by finite differences. This can be explained by the presence of an outer second derivative in time in the direct calculation which tends to amplify errors in high frequencies. Direct gradient
<figure><img src="image_22.png" /><figcaption>Fig. 18. Pressure time histories as a function of the time step $dt$ (FNC integration method). $M_0 = 0$. (a) Surface $S_{0.7}$ and (b) Surface $S_{2.0}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 19. Zooms on pressure time histories. $M_0 = 0$. (a) Surface $S_{0.7}$ and (b) Surface $S_{2.0}$.</figcaption></figure>

calculation remains nonetheless, on the whole, more reliable.

### 4.3. Formulas with outer time derivatives and FNC integration method

For configurations with supersonically moving grids, calculations can no longer be performed using Eqs. (19) and (20) and the CGC integration method. Instead, Eqs. (6) and (9) and the FNC integration method must be used, as explained in Section 3.3.

The calculations carried out in the previous section are repeated, this time using the FNC integration method and are supplemented by calculations for surfaces for which the “acoustic Mach number” $M_R$ in the direction of the observation point is locally supersonic. This includes surfaces $S_{1.2}$, $S_{1.6}$ and $S_{2.0}$, for $M_0 = 0$, and surfaces $S_{0.9}$, $S_{1.2}$, $S_{1.6}$ and $S_{2.0}$, for $M_0 = 0.5$.

**Results for the pressure** – Like for the CGC integration method, the calculations are first carried out with a time step equal to $T/120$. For both Mach numbers $M_0$, the signals obtained appear identical whatever the integration surface (Fig. 17). The evolution of the signal as a function of the time step is plotted in Fig. 18 for a flow at rest ($M_0 = 0$), for the two extreme surfaces $S_{0.7}$ and $S_{2.0}$. This evolution highlights small oscillations at high frequencies for the supersonically rotating surface ($S_{2.0}$). Zooming in on the most affected time interval (Fig. 19) shows that these oscillations appear for small time steps, but have a relative amplitude that remains lower than 0.1. These oscillations are due to the modeling of the pressure rate upon which the FNC integration method is founded (section 3.3, second
<figure><img src="image_24.png" /><figcaption>**Fig. 20.** Pressure error ($L_\infty$ norm) as a function of the time step $dt$, for compact grid cells (CGC) and fully non-compact (FNC) integration methods. Surface $S_{0.7}$, $M_0 = 0$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 21.** Pressure gradient time histories (FNC integration method). $M_0 = 0$, $dt = T/120$. (a) Surface $S_{0.7}$ and (b) Surface $S_{2.0}$.</figcaption></figure>

model). Although this pressure rate is more progressive than that of the first model, its temporal derivative remains discontinuous (see Fig. 3). This leads to oscillations in the contribution of cells in supersonic motion for which the outer time derivative in Eq. (6) has a predominant contribution. A more elaborate pressure rate model could eliminate these small oscillations if it were necessary for certain applications. The results for the surface in subsonic motion ($S_{0.7}$) show, for their part, a monotonic convergence of the signal with the time step. For this integration surface, the evolution of the error as a function of the time step is plotted in Fig. 20. The error decreases slightly slower with the FNC integration method than with the CGC one. This slight degradation in accuracy is another consequence of the rather simple modeling of the radiated pressure rate. With the present discretization of the integration surface, the minimum error of the CGC method is not reached by the FNC approach, within the range of time steps explored.

**Results for the pressure gradient** – The three components of the pressure gradient are compared to the reference solution in Fig. 21, for the surfaces $S_{0.7}$ and $S_{2.0}$, and for a fluid at rest. This comparison is also made in Fig. 22 for all the integration surfaces and for $M_0 = 0.5$, for the predominant component ($dp'/dx_2$) for the observation point considered. The results of the FW-H direct calculation by Eq. (9) are in good agreement with the reference solution in all the cases. For surfaces in supersonic motion, the amplitude of the spurious oscillations are slightly superior to those observed for pressure. It is a logical consequence of the presence of a second time derivative in the expression of the gradient.

For the surface $S_{0.7}$ and a fluid at rest, the error ($L_\infty$ norm on the three components) of the gradient calculation by finite differences, starting from the pressures provided by Eq. (6), is plotted as a function of the spatial step $h$ in Fig. 23. For spatial steps $h$ greater than or
<figure><img src="image_26.png" /><figcaption>Fig. 22. Pressure gradient time histories (FNC integration method). $dt = T/120$. (a) $M_0 = 0$ and (b) $M_0 = 0.5$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 23. Pressure gradient error ($L_\infty$ norm) as a function of $h$ in calculation by finite differences, for CGC and FNC integration methods, and error levels of direct calculations. Surface $S_{0.7}$, $M_0 = 0$, $dt = T/360$.</figcaption></figure>

equal to $\lambda_0/100$, the precision of the gradient calculation by finite differences is identical for both integration methods (CGC, FNC). For smaller spatial steps, because of inaccuracies introduced by the radiated pressure rate model, the error with the FNC integration method no longer decreases, unlike with the CGC method. On the other hand, with both methods, the error of the direct gradient calculation is close to the minimum error of the calculation by finite differences, which consolidates the consistency of the results. As with the CGC method, direct calculation thus seems more reliable, especially for signals rich in harmonics (rotor or propeller impulsive noise, for example) for which $h/\lambda$ can take both small and large values in a gradient calculation by finite differences.

## 5. Conclusion

Starting from the Ffowcs Williams-Hawkings surface integral formulation for a moving medium, time-domain expressions were developed for the radiated pressure and its gradient. These expressions, with outer time derivatives, have the benefit of being relatively simple and valid for all integration surfaces (solid or porous, fixed or moving, rigid or deformable). As they do not comprise the $1/(1-M_R)$ singularity, these expressions have the considerable advantage of allowing calculations with integration surfaces in supersonic motion (case of rotating propeller or rotor meshes, for instance). An expression with inner time derivatives, likely to provide more regular signals but limited to fixed surfaces, was also proposed.
Three integration methods, based on the compact grid cells (CGC) usual simplification or a fully non-compact (FNC) approach, were presented, each well-suited to fixed integration surfaces, subsonically moving surfaces and supersonically moving surfaces, respectively. Calculation algorithms were proposed for each integration method, and their computational efficiency was illustrated by examples of calculation times.

The results of these expressions and integration methods were compared to analytical solutions for the case of a fixed monopole and for that of a rotating monopole. For fixed or subsonically moving surfaces, the computational error can be considered negligible for both the pressure and its gradient, for both CGC and FNC approaches. They clearly show the benefit of the direct gradient calculation over a calculation by finite differences of the pressure around the observation point, in particular for broadband signals for which the spatial step in the calculation by finite differences cannot be optimal for all wavelengths.

Particular attention was devoted to the order of precision of the calculations, which, to our knowledge, has never been done before. In the present case, the choice of time derivatives by centered finite differences and time integration schemes based on linear interpolations leads to a theoretical second order.

For the pressure, on the whole, this theoretical second order is found quite well. With the fully non-compact method, however, the pressure rate model, as it stands, slightly degrades the accuracy of the calculations. For surfaces in supersonic motion, which require this particular integration method, the error linked to this model is dominant and does not really depend on the time step. Nevertheless, the error remains below 10 % in the $L_\infty$ norm and is essentially due to spurious fluctuations at the sampling frequency and a few subharmonics that can be filtered, if necessary. It therefore seems acceptable for most applications.

For fixed or subsonically moving surfaces, the order of precision of the CGC integration method can be increased without difficulty. Indeed, this can be achieved simply by increasing the order of the finite-difference time derivatives and the order of the interpolations in the time integration process. However, it should be noted that the order of precision of the time integration scheme and that of the time derivatives must be identical for consistency (the final order of precision of the calculations being determined by the smaller of the two orders). Improving accuracy for surfaces in supersonic motion is also possible, but requires the design of a more elaborate model for the pressure rate on which the FNC integration method is based.

### CRediT authorship contribution statement

**Gilles Rahier:** Writing – review & editing, Writing – original draft, Validation, Software, Formal analysis. **Jean Prieur:** Validation, Formal analysis.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

Many thanks to Julie Rahier for her invaluable help in proofreading this article.

### Appendix. Summary of formulas and integration techniques

This appendix summarizes the best suited integration formulas and techniques, based on the kinematics of the integration surface (Table A.1). It should be noted that the similarity of the expressions simplifies the design of programming capable of handling all configurations with relative ease.

**Table A.1**

<table><caption>Best suited formulas and integration techniques based on the kinematics of the integration surface.</caption><thead><tr><th>Case</th><th>Formulas</th><th>Integration technique</th></tr></thead><tbody><tr><td>Fixed integration surface</td><td>p'(X, t) = ∫<sub>S</sub> (1 / (4πd)) [F + (∂F / ∂τ)]<sub>ret</sub> dS (A.1)<br>∂p'(X, t) / ∂X<sub>j</sub> = ∫<sub>S</sub> (1 / (4πd)) [F<sub>ij</sub> + (∂F<sub>2j</sub> / ∂τ) + (∂<sup>2</sup>F<sub>3j</sub> / ∂τ<sup>2</sup>)]<sub>ret</sub> dS (A.2)</td><td>CGC (locally compact) Constant interpolation coefficients See Section 3.1</td></tr><tr><td>Subsonically moving integration surface</td><td>p'(X, t) = ∫<sub>S</sub> [F / (4πd[1 - M<sub>R</sub>])]<sub>ret</sub> dS<br>+ ∂ / ∂t ∫<sub>S</sub> [F / (4πd[1 - M<sub>R</sub>])]<sub>ret</sub> dS (A.3)<br>∂p'(X, t) / ∂X<sub>j</sub> = ∫<sub>S</sub> [F<sub>ij</sub> / (4πd[1 - M<sub>R</sub>])]<sub>ret</sub> dS</td><td>CGC (locally compact) Non constant interpolation coefficients See Section 3.2</td></tr></tbody></table>

(continued on next page)
<table><caption>**Table A.1 (continued)**</caption><thead><tr><th>Case</th><th>Formulas</th><th>Integration technique</th></tr></thead><tbody><tr><td>Supersonically moving integration surface</td><td>$$ +\frac{\partial}{\partial t}\int_{S}\left[\frac{F_{2j}}{4\pi d|1-M_{R}|}\right]_{\text{ret}} dS +\frac{\partial^2}{\partial t^2}\int_{S}\left[\frac{F_{3j}}{4\pi d|1-M_{R}|}\right]_{\text{ret}} dS \quad (A.4) $$<br>$$ p'(\mathbf{X}, t) = \int_{\tau} \frac{F}{4\pi d} \delta(g) dS d\tau +\frac{\partial}{\partial t} \int_{S} \frac{F}{4\pi d} \delta(g) dS d\tau \quad (A.5) $$<br>$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{\tau} \frac{F_{1j}}{4\pi d} \delta(g) dS d\tau +\frac{\partial}{\partial t} \int_{S} \frac{F_{2j}}{4\pi d} \delta(g) dS d\tau +\frac{\partial^2}{\partial t^2} \int_{S} \frac{F_{3j}}{4\pi d} \delta(g) dS d\tau \quad (A.6) $$</td><td>FNC (fully non compact, triangular pressure rate shape function) See Section 3.3</td></tr><tr><td>Particular case of fixed integration surface described by a subsonically or supersonically moving mesh</td><td>$$ p'(\mathbf{X}, t) = \int_{\tau} \frac{1}{4\pi d} \left(F + \frac{\partial F}{\partial \tau}\right) \delta(g) dS d\tau \quad (A.7) $$<br>$$ \frac{\partial p'(\mathbf{X}, t)}{\partial X_j} = \int_{\tau} \frac{1}{4\pi d} \left(F_{1j} + \frac{\partial F_{2j}}{\partial \tau} + \frac{\partial^2 F_{3j}}{\partial \tau^2}\right) \delta(g) dS d\tau \quad (A.8) $$</td><td>FNC (fully non compact, rectangular pressure rate shape function) See Section 3.4</td></tr></tbody></table>

The terms used in the various expressions are developed below:

$$F = \frac{\beta^2}{d^2} \Sigma_1$$

$$F = \frac{\Sigma_1}{a_0 d} + \frac{\Sigma_2}{\beta^2}$$

$$F_{1j} = \beta^2 \frac{\Sigma'_j}{d^2}$$

$$F_{2j} = \frac{\Sigma'_j}{a_0 d} + M_0 \delta_{1j} \frac{\Sigma_1}{a_0 d^2} - \frac{C_j \hat{r}_j \Sigma_2}{\beta^2 d}$$

$$F_{3j} = -\frac{(C_j \hat{r}_j - M_0 \delta_{1j})}{a_0 \beta^2} \left( \frac{\Sigma_1}{a_0 d} + \frac{\Sigma_2}{\beta^2} \right)$$

$$M_R = \frac{1}{a_0 \beta^2} \left( C_i \left( \frac{X_i - Y_i}{d} \right) v_i - M_0 v_1 \right)$$

$$\Sigma_1 = A_i(X_i - Y_i) + \frac{(M_0^2 A_1 - U_0 B)}{\beta^2} (X_1 - Y_1)$$

$$\Sigma_2 = B - \frac{M_0 A_1}{a_0}$$

$$\Sigma'_j = A_j + \frac{(M_0^2 A_1 - U_0 B)}{\beta^2} \delta_{1j} - 3C_j \hat{r}_j \frac{\Sigma_1}{d}$$

$$A_i = \hat{p}' n_i + \rho u_i (u_n - v_n)$$

$$B = \rho_0 v_n + \rho (u_n - v_n)$$

See other notations and reference frame in Section 2.

## Data availability

Data will be made available on request.
References

[1] J.E. Frowcs Williams, D.L. Hawkings, Sound Generation by turbulence and surfaces in arbitrary motion, Ph. Trans. R. Soc. Lond. 264 (1151) (1969) 321–342,
https://doi.org/10.1098/rsta.1969.0031.

[2] M.J. Lighthill, On sound generated aerodynamically, I: general theory, P. Roy. Soc. Lond. A Mat. 211 (1107) (1952) 564–587, https://doi.org/10.1098/rspa.1952.0060.

[3] P. Di Francesantonio, A new boundary integral formulation for prediction of sound radiation, J. Sound Vib. 202 (4) (1997) 491–509, https://doi.org/10.1006/jsvi.1996.0843.

[4] K.S. Brentner, F. Farassat, An analytical comparison of the acoustic analogy and Kirchhoff formulation for moving surfaces, AIAA J 36 (8) (1998) 1379–1386,
https://doi.org/10.2514/2.558.

[5] J. Prieur, G. Rahier, Aeroacoustic integral methods, formulation and efficient numerical implementation, Aerosp. Sci. Technol. 5 (7) (2001) 457–468, https://
doi.org/10.1016/S1270-9638(01)01123-3.

[6] D. Casalino, M. Barbarino, A. Visingardi, Simulation of helicopter community noise in complex urban geometry, AIAA J 49 (8) (2011) 1614–1624, https://doi.
org/10.2514/1.J050774.

[7] J. Drofelnik, M. Andrejasic, B. Mocan, T. Kosel, J. Christophe, et al., Measurement and Modelling of Aero-Acoustic Installation Effects in Tractor and Pusher Propeller Architectures, AIAA Aviation Forum, 2005, p. 2301, https://doi.org/10.2514/6.2021-2301.

[8] Z. Wang, J. Huang, M. Yi, Acoustic scattering effect prediction of helicopter fuselage based on BEM and convective FW-H equation, Acta Acoustica 6 (24) (2022), https://doi.org/10.1051/aacus/2022019.

[9] M. Huet, F. Gand, G. Rahier, Simulation of isolated and installed jet noise at Mach=0.9: influence of numerical mesh and physical insights, Flow Turbul. Combust. (2023) 1–43, https://doi.org/10.1007/s10494-023-00461-y.

[10] S. Lee, K.S. Brentner, Comment on, Acoustic Velocity Formulation for Sources in Arbitrary Motion, AIAA J 54 (5) (2016) 1810–1811, https://doi.org/10.2514/1.J054845.

[11] F. Farassat, K.S. Brentner, The derivation of the gradient of the acoustic pressure on a moving surface for application to the fast scattering code (FSC), NASA/TM-2005-213777. https://core.ac.uk/download/pdf/10514412.pdf.

[12] S. Lee, K.S. Brentner, F. Farassat, P.J. Morris, Analytic formulation and numerical implementation of an acoustic pressure gradient prediction, J. Sound Vib. 319 (3–5) (2009) 1200–1221, https://doi.org/10.1016/j.jsv.2008.06.028.

[13] G. Ghorbaniasl, M. Carley, C. Lacor, Acoustic velocity formulation for sources in arbitrary motion, AIAA J 51 (3) (2013) 632–642, https://doi.org/10.2514/1.J051958.

[14] C.-X. Bi, Z.-H. Wang, X.-Z. Zhang, Analytic time-domain formulation for acoustic pressure gradient prediction in a moving medium, AIAA J 55 (8) (2017) 2607–2616, https://doi.org/10.2514/1.J055630.

[15] Y. Mao, Q. Zhang, C. Xu, D. Qi, Two types of frequency domain acoustic-velocity formulations for rotating thickness and loading sources, AIAA J 53 (3) (2015) 713–722, https://doi.org/10.2514/1.J053230.

[16] G. Ghorbaniasl, Z. Huang, L. Siozos-Rousoulis, C. Lacor, Analytical acoustic pressure gradient prediction for moving medium problems, Proc. R. Soc. Lond. Ser. A, Contain Pap. Math. Phys. Eng. Sci. 471 (2184) (2015) 1–14, https://doi.org/10.1098/rspa.2015.0342.

[17] G. Rahier, M. Huet, J. Prieur, Additional terms for the use of Frowcs Williams and Hawkings surface integrals in turbulent flows, Comput. Fluids 120 (2015) 158–172, https://doi.org/10.1016/j.compfluid.2015.07.014.

[18] G. Rahier, Comparison of surface and volume integral methods for transonic propeller acoustic predictions, Comput. Fluids 179 (2018) 178–193, https://doi.
org/10.1016/j.compfluid.2018.10.015.

[19] A. Najafi-Yazdi A, G.A. Brès, L. Mongeau L, An acoustic analogy formulation for moving sources in uniformly moving media, Proc. R. Soc. A 467 (2125) (2011) 144–165, https://doi.org/10.1098/rspa.2010.0172.

[20] S. Ianniello, Algorithm to integrate the Frowcs Williams-Hawkings equation on supersonic rotating domain, AIAA J 37 (9) (1999) 1040–1047, https://doi.org/
10.2514/2.831.

[21] S. Ianniello, New perspectives in the use of the Frowcs Williams-Hawkings equation for aeroacoustic analysis of rotating blades, J. Fluid Mech. 570 (2007) 79–127, https://doi.org/10.1017/S002211200600293X.

[22] G. Rahier, J. Prieur, An efficient Kirchhoff integration method for rotor noise prediction starting indifferently from subsonically or supersonically rotating meshes, in: 53rd AHS forum, 1997.

[23] Y. Colin, B. Caruelle, A.B. Parry, Computational strategy for predicting CROR noise at low-speed. Part III: investigation of the noise radiation with the Frowcs Williams-Hawkings analogy, in: 18th AIAA/CEAS Aeroacoustics Conference, 2012, p. 2223, https://doi.org/10.2514/6.2012-2223.

[24] A. Giauque, B. Ortun, B. Rodriguez, B. Caruelle, Numerical error analysis with application to transonic propeller aeroacoustics, Comput. Fluids 69 (2012) 20–34,
https://doi.org/10.1016/j.compfluid.2012.07.022.

[25] D. Lockard, A comparison of Frowcs Williams-Hawkings solvers for airframe noise applications, in: 8th AIAA/CEAS Aeroacoustics Conference & Exhibit, 2002, p. 2580, https://doi.org/10.2514/6.2002-2580.

[26] A.P. Dowling, J.E. Frowcs Williams, Sound and Sources of Sound, Chap, 9, Horwood Publishing, Westergate, 1983, pp. 207–208.

[27] M.E. Goldstein, Aeroacoustics, Mc Graw-Hill International Book Co., New York, 1976.