# On the use of a tailored fluid-fluid Green's function to predict scattering from two-phase fluid interfaces

Louise Pacaut, Gilles Serre, Jean-François Mercier, Stéphanie Chaillat

► To cite this version:

Louise Pacaut, Gilles Serre, Jean-François Mercier, Stéphanie Chaillat. On the use of a tailored fluid-fluid Green's function to predict scattering from two-phase fluid interfaces. INTER-NOISE 2023 - 52nd International Congress and Exposition on Noise Control Engineering, Nov 2023, Chiba, Japan. pp.434-443, 10.3397/IN_2023_0076. hal-04875581

HAL Id: hal-04875581

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

Submitted on 10 Jan 2025

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 Attribution 4.0 International License
# On the use of a tailored fluid-fluid Green's function to predict scattering from two-phase fluid interfaces

Louise Pacaut<fnref n="1" />

Gilles Serre<fnref n="2" />

Naval Group

199 Avenue Pierre-Gilles de Gennes, 83190 Ollioules, France

Jean-François Mercier<fnref n="3" />

Stéphanie Chaillat<fnref n="4" />

POEMS, ENSTA Paris, CNRS, INRIA, Institut Polytechnique de Paris

828 Boulevard des Maréchaux, 91120 Palaiseau, France

## ABSTRACT

In naval applications, the precise knowledge of the acoustical behaviour of two-phase fluids is of interest. Indeed, they appear in many configurations as bubbles curtains, cavitation (along propellers and pipes) or two-phase turbulent boundary layers and wakes (gas exhaust). To model this problem, we focus on a Boundary Element/ Boundary Element coupling for low to moderate frequencies. We propose two approaches: (a) solving the Helmholtz equations in each phase by introducing the free field Green's function and (b) determining a tailored Green's function, taking into account the presence of the two-phase fluid. The determination of a tailored Green's function has two main advantages: (i) it allows a reduction of the numerical model, since it contains all the information on the fluid-fluid coupling, in particular the transmission conditions and resonances and (ii) it is source-independent and thus it gives directly the answer to any source. Numerical tests on the evaluation of radiated noise are performed to determine the efficiency of the approach based on a tailored Green's function compared to the free-field based one.

## 1. INTRODUCTION

Two-phase fluids are often encountered in the naval field, as bubble curtains, gas exhaust or two-phase turbulent boundary layers and wakes. A good knowledge of the scattered pressure by a two-phase fluid is of interest, to optimize its characteristics. Various numerical methods to compute fluid-structure or fluid-fluid interactions have been proposed. In the literature BEM-BEM couplings are based on the solution of the Helmholtz equation with the free field Green's function ([1], [2]). In this work, we focus on BEM-BEM coupling for fluid-fluid interaction and we determine a Green's function tailored to this type of interaction. Then, we validate this method by comparison to analytical results and to a classical method using two free field Green's functions. The validation case is a bubble of air in water. Finally, we focus on the case study of an air-bubble curtain in water.

<footnote n="1">¹louise.pacaut@ensta-paris.fr</footnote>

<footnote n="2">²gilles.serre@naval-group.com</footnote>

<footnote n="3">³jean-francois.mercier@ensta-paris.fr</footnote>

<footnote n="4">⁴stephanie.chaillat@ensta-paris.fr</footnote>
## 2. NUMERICAL DETERMINATION OF A TAILORED GREEN'S FUNCTION IN A FLUID-FLUID COUPLED PROBLEM

Let a fluid of volume $\Omega_2$, density $\rho_2$ and boundary $\Gamma$ be in fluid of volume $\Omega_1$ and density $\rho_1$.

Noting $p_1$ the exterior pressure in $\Omega_1$ and $p_2$ the interior pressure in $\Omega_2$, we assume the exterior pressure is decomposed as

$$p_1 = p_i + p_{1s}, \qquad (1)$$

where $p_i$ is the incident pressure wave and $p_{1s}$ is the scattered pressure. $p_{1s}$ and $p_2$ satisfy

$$ \left\{ \begin{array}{l} (\Delta_z + k_1^2)p_{1s}(\mathbf{x}) = 0 \quad \text{for } \mathbf{x} \in \Omega_1, \\ (\Delta_y + k_2^2)p_2(\mathbf{x}) = 0 \quad \text{for } \mathbf{x} \in \Omega_2, \\ p_{1s} - p_i = p_2 \quad \text{on } \Gamma, \\ \displaystyle\frac{1}{\rho_1}\frac{\partial(p_{1s} - p_i)}{\partial n} = \frac{1}{\rho_2}\frac{\partial p_2}{\partial n} \quad \text{on } \Gamma. \end{array} \right. \qquad (2) $$

### 2.1. Boundary integral representation in each domain

The solution of the Helmholtz equation can be obtained by introducing an arbitrary Green's function $\tilde{G}_i$ defined, by

$$(\Delta + k_i^2)\tilde{G}_i(\mathbf{x}, \mathbf{y}) = -\delta(\mathbf{x} - \mathbf{y}), \quad \forall (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^3 \times \mathbb{R}^3. \qquad (3)$$

In the fluid domain $\Omega_1$, the scattered pressure is solution of the Helmholtz equation

$$(\Delta + k_1^2)p_{1s}(\mathbf{x}) = 0, \quad \forall \mathbf{x} \in \Omega_1. \qquad (4)$$

Multiplying (4) by $\tilde{G}_1$ and integrating over $\Omega_1$, we get

$$ \int_{\Gamma} \left( p_{1s}(\mathbf{Y}) \frac{\partial \tilde{G}_1}{\partial n_Y}(\mathbf{x}, \mathbf{Y}) - \frac{\partial p_{1s}}{\partial n_Y}(\mathbf{Y}) \tilde{G}_1(\mathbf{x}, \mathbf{Y}) \right) d\mathbf{Y} = \begin{cases} p_{1s}(\mathbf{x}) & \forall \mathbf{x} \in \Omega_1, \\ 0 & \forall \mathbf{x} \in \Omega_2, \\ \frac{p_{1s}(\mathbf{x})}{2} & \forall \mathbf{x} \in \Gamma. \end{cases} \qquad (5) $$

In the fluid domain $\Omega_2$, the interior pressure is solution of the Helmholtz equation

$$(\Delta + k_2^2)p_2(\mathbf{x}) = 0, \quad \forall \mathbf{x} \in \Omega_2. \qquad (6)$$

In the same way, multiplying (6) by $\tilde{G}_2$ and integrating over $\Omega_2$, we get

$$ \int_{\Gamma} \left( \frac{\partial p_2}{\partial n_Y}(\mathbf{Y}) \tilde{G}_2(\mathbf{x}, \mathbf{Y}) - p_2(\mathbf{Y}) \frac{\partial \tilde{G}_2}{\partial n_Y}(\mathbf{x}, \mathbf{Y}) \right) d\mathbf{Y} = \begin{cases} 0 & \forall \mathbf{x} \in \Omega_1, \\ p_2(\mathbf{x}) & \forall \mathbf{x} \in \Omega_2, \\ \frac{p_2(\mathbf{x})}{2} & \forall \mathbf{x} \in \Gamma. \end{cases} \qquad (7) $$

### 2.2. Formulation based on the free field Green's function

A simple choice for $\tilde{G}_i$ is to consider the free field Green's functions $G_i^{(0)}$ in each domain $\Omega_i$. These functions are defined by

$$ G_i^{(0)}(\mathbf{x}, \mathbf{y}) = \frac{e^{ik_i |\mathbf{x} - \mathbf{y}|}}{4\pi |\mathbf{x} - \mathbf{y}|}, \quad \forall (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^3 \times \mathbb{R}^3, \mathbf{x} \neq \mathbf{y}. \qquad (8) $$

It is more convenient to express (5) in terms of the total field. To do so, we remember that the incoming pressure wave $p_i$ is solution of the Helmholtz equation

$$(\Delta + k_1^2)p_i(\mathbf{x}) = 0, \quad \forall \mathbf{x} \in \Omega_2, \qquad (9)$$
such that multiplying (9) by $G_1^{(0)}$ and integrating over $\Omega_2$ leads to

$$ \int_{\Gamma} \left( \frac{\partial p_i}{\partial n_y}(\mathbf{Y}) G_1^{(0)}(\mathbf{x}, \mathbf{Y}) - p_i(\mathbf{Y}) \frac{\partial G_1^{(0)}}{\partial n_y}(\mathbf{x}, \mathbf{Y}) \right) d\mathbf{Y} = \begin{cases} 0 & \forall \mathbf{x} \in \Omega_1, \\ \frac{p_i(\mathbf{x})}{2} & \forall \mathbf{x} \in \Omega_2, \\ \frac{p_i(\mathbf{x})}{2} & \forall \mathbf{x} \in \Gamma. \end{cases} \quad (10) $$

(5) - (10), and considering (1), lead to

$$ \int_{\Gamma} \left( p_1(\mathbf{Y}) \frac{\partial G_1^{(0)}}{\partial n_y}(\mathbf{x}, \mathbf{Y}) - \frac{\partial p_1}{\partial n_y}(\mathbf{Y}) G_1^{(0)}(\mathbf{x}, \mathbf{Y}) \right) d\mathbf{Y} = \frac{p_{1s}(\mathbf{x}) - p_i(\mathbf{x})}{2} = \frac{p_1(\mathbf{x})}{2} - p_i(\mathbf{x}), \quad \forall \mathbf{x} \in \Gamma. \quad (11) $$

Introducing the single and double layer potentials:

$$ (S_i^{(0)} p)(\mathbf{x}) = \int_{\Gamma} G_i^{(0)}(\mathbf{x}, \mathbf{Y}) p(\mathbf{Y}) d\mathbf{Y}, \quad (D_i^{(0)} p)(\mathbf{x}) = \int_{\Gamma} \frac{\partial G_i^{(0)}(\mathbf{x}, \mathbf{Y})}{\partial n} p(\mathbf{Y}) d\mathbf{Y}, $$

(11) and (7) can be written as

$$ \left\{ \begin{aligned} \frac{p_1}{2} &= p_i + D_1^{(0)} p_1 - S_1^{(0)} \left( \frac{\partial p_1}{\partial n} \right) \quad \text{in } \Gamma, \\ \frac{p_2}{2} &= -D_2^{(0)} p_2 + S_2^{(0)} \left( \frac{\partial p_2}{\partial n} \right) \quad \text{in } \Gamma, \end{aligned} \right. \quad (12) $$

where $S_i$ and $D_i$ are the Dirichlet traces of $S_i$ and $D_i$. We obtain the system

$$ \left\{ \begin{aligned} \left( \frac{I}{2} - D_1^{(0)} \right) p_1 &= p_i - S_1^{(0)} \left( \frac{\partial p_1}{\partial n} \right) & \text{on } \Gamma, \\ \left( \frac{I}{2} + D_2^{(0)} \right) p_2 &= S_2^{(0)} \left( \frac{\partial p_2}{\partial n} \right) & \text{on } \Gamma, \end{aligned} \right. \quad (13) $$

with the transmission conditions at the interface $\Gamma$, (13) can be given using only the unknowns $p_1$ and $\frac{\partial p_1}{\partial n}$

$$ \left\{ \begin{aligned} \left( \frac{I}{2} - D_1^{(0)} \right) p_1 + S_1^{(0)} \left( \frac{\partial p_1}{\partial n} \right) &= p_i \quad \text{on } \Gamma, \\ \left( \frac{I}{2} + D_2^{(0)} \right) p_2 - \frac{\rho_2}{\rho_1} S_2^{(0)} \left( \frac{\partial p_1}{\partial n} \right) &= 0 \quad \text{on } \Gamma. \end{aligned} \right. \quad (14) $$

## 2.3. Tailored Green's function for the case of the fluid-fluid interaction

Another possible choice for $\tilde{G}_i$ is to consider a tailored Green's function for the fluid-fluid coupled problem [3]. The definition of a tailored Green's function is useful and physically meaningful since it contains all the information on the fluid-fluid coupling, in particular it includes the transmission conditions and internal resonances and, as we will show, it gives directly the scattered pressure for any source (plane wave, turbulent boundary layer, shock wave).

We introduce the tailored Green's function $G$, defined in $\mathbb{R}^3 \times \mathbb{R}^3$. From now, $\mathbf{x} \in \Omega_1$ and we note $G_1 = G|_{\Omega_1}$ and $G_2 = G|_{\Omega_2}$ the restrictions of $G$. $G$ satisfies

$$ \left\{ \begin{aligned} (\Delta_z + k_1^2) G_1(\mathbf{x}, \mathbf{y}) &= -\delta(\mathbf{x} - \mathbf{y}) & \text{for } \mathbf{y} \in \Omega_1, \\ (\Delta_y + k_2^2) G_2(\mathbf{x}, \mathbf{y}) &= 0 & \text{for } \mathbf{y} \in \Omega_2, \\ G_1 &= G_2 & \text{on } \Gamma, \\ \frac{1}{\rho_1} \frac{\partial G_1}{\partial n} &= \frac{1}{\rho_2} \frac{\partial G_2}{\partial n} & \text{on } \Gamma. \end{aligned} \right. \quad (15) $$
(5) and (7) give

$$ egin{cases} p_{1s}(\mathbf{x}) = \int_{\Gamma} \frac{\partial G_1(\mathbf{x}, \mathbf{Y})}{\partial n} p_{1s}(\mathbf{Y}) - G_1(\mathbf{x}, \mathbf{Y}) \frac{\partial p_{1s}(\mathbf{Y})}{\partial n} d\mathbf{Y}, \\ 0 = \int_{\Gamma} -\frac{\partial G_2(\mathbf{x}, \mathbf{Y})}{\partial n} p_2(\mathbf{Y}) + G_2(\mathbf{x}, \mathbf{Y}) \frac{\partial p_2(\mathbf{Y})}{\partial n} d\mathbf{Y}. \end{cases} \quad (16) $$

Moreover, using the transmission conditions,

$$ egin{aligned} \int_{\Gamma} \frac{\partial G_1}{\partial n} p_{1s} - G_1 \frac{\partial p_{1s}}{\partial n} &= \int_{\Gamma} \frac{\partial G_1}{\partial n} (p_1 - p_i) - G_1 \frac{\partial (p_1 - p_i)}{\partial n}, \\ &= - \int_{\Gamma} \left( \frac{\partial G_1}{\partial n} p_i - G_1 \frac{\partial p_i}{\partial n} \right) + \frac{\rho_1}{\rho_2} \int_{\Gamma} \left( \frac{\partial G_2}{\partial n} p_2 - G_2 \frac{\partial p_2}{\partial n} \right), \\ &= \int_{\Gamma} G_1 \frac{\partial p_i}{\partial n} - \frac{\partial G_1}{\partial n} p_i. \end{aligned} $$

Therefore, (16) can be written $\forall \mathbf{x} \in \Omega_1$

$$ p_1^s(\mathbf{x}) = \int_{\Gamma} -\frac{\partial G_1(\mathbf{x}, \mathbf{Y})}{\partial n} p_i(\mathbf{Y}) + G_1(\mathbf{x}, \mathbf{Y}) \frac{\partial p_i(\mathbf{Y})}{\partial n} d\mathbf{Y}. \quad (17) $$

Therefore, the scattered pressure is directly (contrary to (14)) given by the knowledge of the tailored Green's function and the incident pressure wave. The only remaining step is now moved to the determination of the Green's function $G$.

Now, we will show that to determine $G$ we have to solve

$$ egin{cases} \left(\frac{I}{2} - D_1^{(0)}\right) G_1 = G_1^{(0)} - S_1^{(0)} \left(\frac{\partial G_1}{\partial n}\right) & \text{on } \Gamma, \\ \left(\frac{I}{2} + D_2^{(0)}\right) G_1 = S_2^{(0)} \frac{\rho_2}{\rho_1} \left(\frac{\partial G_1}{\partial n}\right) & \text{on } \Gamma. \end{cases} \quad (18) $$

Equation satisfied by $G_1$ We know that, for $(\mathbf{x}, \mathbf{y}, \mathbf{z}) \in \Omega_1^3$:

$$ egin{cases} (\Delta_z + k_1^2) G_1(\mathbf{x}, \mathbf{z}) = -\delta(\mathbf{x} - \mathbf{z}), \\ (\Delta_z + k_1^2) G_1^{(0)}(\mathbf{y}, \mathbf{z}) = -\delta(\mathbf{y} - \mathbf{z}). \end{cases} \quad (19) $$

Multiplying the first equation by $G_1^{(0)}(\mathbf{x}, \mathbf{z})$ and the second by $G_1(\mathbf{x}, \mathbf{z})$, combining both equations, then integrating over $\Omega_1$ lead to:

$$ \int_{\Omega_1} G_1^{(0)}(\mathbf{y}, \mathbf{z})(\Delta_z + k_1^2)G_1(\mathbf{x}, \mathbf{z})d\mathbf{z} = -\left[\int_{\Gamma} G_1^{(0)}(\mathbf{y}, \mathbf{Z})\frac{\partial G_1(\mathbf{x}, \mathbf{Z})}{\partial n_z} - G_1(\mathbf{x}, \mathbf{Z})\frac{\partial G_1^{(0)}(\mathbf{y}, \mathbf{Z})}{\partial n_z}d\mathbf{Z}\right] + \int_{\Omega_1} G_1(\mathbf{x}, \mathbf{z})(\Delta_z + k_1^2)G_1^{(0)}(\mathbf{x}, \mathbf{z})d\mathbf{z}, \quad (20) $$

which leads to

$$ \int_{\Omega_1} G_1^{(0)}(\mathbf{y}, \mathbf{z})\delta(\mathbf{x} - \mathbf{z})d\mathbf{z} = \int_{\Gamma} G_1^{(0)}(\mathbf{y}, \mathbf{Z})\frac{\partial G_1(\mathbf{x}, \mathbf{Z})}{\partial n_z} - G_1(\mathbf{x}, \mathbf{Z})\frac{\partial G_1^{(0)}(\mathbf{y}, \mathbf{Z})}{\partial n_z}d\mathbf{Z} + \int_{\Omega_1} G_1(\mathbf{x}, \mathbf{z})\delta(\mathbf{y} - \mathbf{z})d\mathbf{z}. \quad (21) $$
Since $(\mathbf{x}, \mathbf{y}) \in \Omega_1^2$, we get

$$G_1(\mathbf{x}, \mathbf{y}) - G_1^{(0)}(\mathbf{x}, \mathbf{y}) = - \left[ \int_{\Gamma} G_1^{(0)}(\mathbf{y}, \mathbf{Z}) \frac{\partial G_1(\mathbf{x}, \mathbf{Z})}{\partial n_{\mathbf{z}}} - G_1(\mathbf{x}, \mathbf{Z}) \frac{\partial G_1^{(0)}(\mathbf{y}, \mathbf{Z})}{\partial n_{\mathbf{z}}} d\mathbf{Z} \right], \quad (22)$$

or, in a more compact form,

$$G_1 = G_1^{(0)} + \mathcal{D}_1^{(0)} G_1 - S_1^{(0)} \left( \frac{\partial G_1}{\partial n} \right) \quad \text{in } \Omega_1. \quad (23)$$

**Equation satisfied by $G_2$.** We know that, for $\mathbf{x} \in \Omega_1$, $(\mathbf{y}, \mathbf{z}) \in \Omega_2^2$:

$$\left\{ \begin{array}{l} (\Delta_z + k_2^2) G_2(\mathbf{x}, \mathbf{z}) = 0, \\ (\Delta_z + k_2^2) G_2^{(0)}(\mathbf{y}, \mathbf{z}) = -\delta(\mathbf{y} - \mathbf{z}). \end{array} \right. \quad (24)$$

Multiplying the first equation by $G_2^{(0)}$ and the second by $G_2$, combining both equations, then integrating over $\Omega_2$ lead to:

$$\int_{\Omega_2} G_2^{(0)}(\mathbf{y}, \mathbf{z})(\Delta_z + k_2^2)G_2(\mathbf{x}, \mathbf{z})d\mathbf{z} = \int_{\Gamma} G_2^{(0)}(\mathbf{y}, \mathbf{Z}) \frac{\partial G_2(\mathbf{x}, \mathbf{Z})}{\partial n_{\mathbf{z}}} - G_2(\mathbf{x}, \mathbf{Z}) \frac{\partial G_2^{(0)}(\mathbf{y}, \mathbf{Z})}{\partial n_{\mathbf{z}}} d\mathbf{Z} + \int_{\Omega_2} G_2(\mathbf{x}, \mathbf{z})(\Delta_z + k_2^2)G_2^{(0)}(\mathbf{y}, \mathbf{z})d\mathbf{z}, \quad (25)$$

which leads to, for $\mathbf{x} \in \Omega_1$, $\mathbf{y} \in \Omega_2$

$$G_2(\mathbf{x}, \mathbf{y}) = \int_{\Gamma} G_2^{(0)}(\mathbf{y}, \mathbf{Z}) \frac{\partial G_2(\mathbf{x}, \mathbf{Z})}{\partial n_{\mathbf{z}}} - G_2(\mathbf{x}, \mathbf{Z}) \frac{\partial G_2^{(0)}(\mathbf{y}, \mathbf{Z})}{\partial n_{\mathbf{z}}} d\mathbf{Z}, \quad (26)$$

or in a more compact form

$$G_2 = -\mathcal{D}_2^{(0)} G_2 + S_2^{(0)} \left( \frac{\partial G_2}{\partial n} \right) \quad \text{in } \Omega_2. \quad (27)$$

From the boundary integral representations (23) and (27) and the transmission conditions at the interface, we deduce the system of boundary integral equations:

$$\left\{ \begin{array}{l} \left(\frac{I}{2} - D_1^{(0)}\right) G_1 = G_1^{(0)} - S_1^{(0)} \left(\frac{\partial G_1}{\partial n}\right) \quad \text{on } \Gamma, \\ \left(\frac{I}{2} + D_2^{(0)}\right) G_1 = S_2^{(0)} \frac{\rho_2}{\rho_1} \left(\frac{\partial G_1}{\partial n}\right) \quad \text{on } \Gamma. \end{array} \right.$$

### 3. NUMERICAL VALIDATIONS WITH ANALYTIC SOLUTIONS

To validate our new formulation, we consider a bubble of air in water, of radius $a$. The source is a plane wave, incoming towards the sphere along increasing $\mathbf{z}$ (Figure 1).

**Analytic solution.** Considering the observer point $\mathbf{x}_{\text{obs}}$ is placed in the $(\mathbf{x}, \mathbf{z})$ plane (Figure 1, with $\varphi = 0$), $\mathbf{x}_{\text{obs}}$ can be written: $\mathbf{x}_{\text{obs}}=r(\sin \theta, 0, \cos \theta)$. Thus, we can write the pressures as ([4], [5]):

$$p_i(\theta, r) = A \sum_{n=0}^{\infty} (i)^n (2n + 1) j_n(k_1 r) P_n(\cos \theta) \quad \text{for } \mathbf{x} \text{ in } \Omega_1, \quad (28)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Definition of the coordinates used in this study.</figcaption></figure>

$$p_{1s}(\theta, r) = \sum_{n=0}^{\infty} A_n h_n(k_1 r) P_n(\cos \theta) \text{ for } x \text{ in } \Omega_1, \qquad (29)$$

$$p_2(\theta, r) = \sum_{n=0}^{\infty} B_n j_n(k_2 r) P_n(\cos \theta) \text{ for } x \text{ in } \Omega_2, \qquad (30)$$

where A is the amplitude of the incoming plane wave, $j_n$ and $h_n$ are spherical Bessel's functions, $P_n$ are Legendre polynomials, $A_n$ and $B_n$ are constants to be determined.

Since the Legendre polynomials are linearly independent, the boundary conditions in (2) lead to

$$B_n = A(i)^n (2n + 1) \frac{j_n(k_1 a) h'_n(k_1 a) - j'_n(k_1 a) h_n(k_1 a)}{j_n(k_2 a) h'_n(k_1 a) - \frac{\rho_1 c_1}{\rho_2 c_2} j'_n(k_2 a) h_n(k_1 a)}, \qquad (31)$$

and

$$A_n = -A(i)^n (2n + 1) \frac{-\frac{\rho_1 c_1}{\rho_2 c_2} j_n(k_1 a) j'_n(k_2 a) + j'_n(k_1 a) j_n(k_2 a)}{j_n(k_2 a) h'_n(k_1 a) - \frac{\rho_1 c_1}{\rho_2 c_2} j'_n(k_2 a) h_n(k_1 a)}. \qquad (32)$$

**Evaluation of the boundary integral representation.** We set the radius of the sphere at $a = 1$ m. We use a mesh with 5659 nodes and 11314 elements, designed to have at least 10 points per wavelength for all the frequencies considered. The relative error over the scattered pressures is defined by

$$err = \frac{|p_s^{num} - p_{1s}|}{|p_{1s}|}, \qquad (33)$$

where $p_s^{num}$ is the computed scattered pressure and $p_{1s}$ is the analytical scattered pressure (29).

Figures 2 to 3 show a good agreement between numerical and analytical results. The relative error is, for most frequencies, smaller than 2% which are acceptable errors over pressure. On Figure 4, various peaks at typical resonance frequencies are visible. They correspond to the cancellation of the denominator of the pressure's analytical expressions (31) and (32). Figure 5 shows the cancellations denominator of $A_n$ (32) for the first 3 modes.
<figure><img src="image_3.png" /><figcaption>Figure 2: Comparison between formulation (17) with $G_1$, formulation (14) with $G_1^{(0)}$ and the analytical solution (29) at frequency $f = 600$ Hz. The source is a plane wave incoming along increasing $z$. The observers are moving on a circle of radius $r = 10$ m with $\theta$ varying.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Computation of the exterior scattered pressure, for an observer at $r = 10$ m, $\varphi = 0$ and $\theta = 0$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Relative error between formulation (17) with $G_1$ and the analytical solution (29) averaged on 60 random points at 10 m from the sphere. A similar plot is found for the formulation with $G_1^{(0)}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Denominator of $A_n$ (32) for the first 3 modes of the analytical pressure.</figcaption></figure>

## 4. CASE STUDY: CURTAIN OF BUBBLES OF AIR IN WATER

A curtain of bubbles of air in water is used in various applications. For instance, it is used to reduce pile diving noise [6], to cover tracks of military ships or to protect submerged structures against explosive shock waves [7]. A good modeling of the curtain allows to optimize its characteristics, depending on the targeted application. Since the sources that can excite an air-bubble curtain are various, it is of interest to be able to compute the scattered pressure for any source and this can be done using the formulation (17) with $G_1$.

We consider an idealized a curtain of bubbles of air consisting of a finite plate containing bubbles of equal radii and distributed homogeneously in the curtain ([7], [8]). We consider the following characteristics:

- Radius of the bubbles: $a = 0.5$ mm

- Width: $D_x = 12.5$ cm

- Height and length: $D_y = D_z = 1$ m

- Fractional air volume in the curtain: $s = 0.06\%$

Such configuration can be modeled by an homogeneous plate [7] with an acoustic wave velocity varying with frequency as shown Figure 6. The plate is excited by a plane wave incoming with $\varphi_{PW} = \pi$ and $\theta_{PW} = \frac{\pi}{2}$. The mesh of the plate has 5 312 nodes and 10 620 elements and is designed to have at least 6 points per wavelength at each frequency.

<figure><img src="image_7.png" /><figcaption>Figure 6: Acoustic wave velocity of an bubbles of air / water mixture, with $a = 0.5$ mm and $s = 0.06\%$.</figcaption></figure>
Figure 7 shows a good agreement between the formulation (17) with $G_1$ and reference formulation (14) with $G_1^{(0)}$. A finer mesh considering at least 10 points per wavelength could improve the results, especially around 500 Hz.

<figure><img src="image_8.png" /><figcaption>Figure 7: Magnitude ($20 \log_{10}(|p_1|)$) of the total pressure $p_1$ computed with formulation (17) with $G_1$ and formulation (14) with $G_1^{(0)}$. The observer is located $r = 0.125$ m from the plate, $\varphi = \pi$ and $\theta = \frac{\pi}{2}$.</figcaption></figure>

# 5. CONCLUSION

We have proposed a BEM-BEM coupling method to determine a Green's function tailored to a fluid obstacle in a fluid domain. The advantage of considering a tailored Green's function is that its computation does not rely on the source (or excitation) type (e.g. shock wave, plane wave, turbulent boundary layer) and only depends on the geometry of the obstacle and on the fluids characteristics. The formulation with the tailored Green's function (17) is compared to formulation (14) with free field Green's functions and an analytical solution (29) for a sphere and applied to the case study of a simplified air-bubble curtain. The comparisons show satisfactory results.

# ACKNOWLEDGEMENTS

The authors would like to thank Naval Group and Association Nationale de la Recherche et de la Technologie (ANRT) for supporting this work.

# REFERENCES

1. Daniel R Wilkes and Alec J Duncan. Acoustic coupled fluid-structure interactions using a unified fast multipole boundary element method. *The Journal of the Acoustical Society of America*, 137(4):2158-2167, 2015.

2. Rudyard Jerez Boudesseul and Elwin van't Wout. Modeling frequency shifts of collective bubble resonances with the boundary element method. *The Journal of the Acoustical Society of America*, 153(3):1898-1911, 2023.

3. Stéphanie Chaillat, Benjamin Cotté, Jean-François Mercier, Gilles Serre, and Nicolas Trafny. Efficient evaluation of three-dimensional helmholtz green's functions tailored to arbitrary rigid geometries for flow noise simulations. *Journal of Computational Physics*, 452:110915, 2022.

4. Elizabeth A Skelton, Jimmy H James, and David Feit. *Theoretical acoustics of underwater structures*. Acoustical Society of America, 2000.
5. Victor C Anderson. Sound scattering from a fluid sphere. *The Journal of the Acoustical Society of America*, 22(4):426-431, 1950.

6. Mathias H Andersson, Brodd Leif Andersson, Jörgen Pihl, Leif KG Persson, Peter Sigray, Sandra Andersson, Andreas Wikström, Jimmy Ahlsén, and Jonatan Hammar. *A framework for regulating underwater noise during pile driving*. Naturvårdsverket, 2017.

7. SN Domenico. Acoustic wave propagation in air-bubble curtains in water; part i, history and theory. *Geophysics*, 47(3):345-353, 1982.

8. Hervé Grandjean. *Propagation d'une onde de choc dans un liquide aéré: modélisation et application aux rideaux de bulles*. PhD thesis, Université de Bretagne occidentale-Brest, 2012.