[FEM-BEM applications in vibro-acoustics using G YPSILAB]
Matthieu AUSSAL(1) , François ALOUGES(2) , Marc BAKRY(3) and Gilles SERRE(4)(1) Ecole Polytechnique, France, matthieu.aussal@polytechnique.edu(2) Ecole Polytechnique, France, francois.alouges@polytechnique.edu(3) Ecole Polytechnique, France, marc.bakry@polytechnique.edu(4) Naval Group Research, France, gilles.serre@naval-group.com

Abstract
The G YPSILAB framework is an open-source M ATLAB library which allows the user to solve complex problems
using the well-known Finite Element Method (FEM) or the Boundary Element Method (BEM). Cutting edge
algorithms are included like a complete H -matrix algebra for the compression of the matrices arising in the
BEM. In this paper, we study the strong coupling between an incoming sound wave and a vibrating structure offinite dimension. We discuss its implementation with G YPSILAB. The elastic part is solved using the FEM as
it is the most suited for such computations because the medium may be non-homogeneous or anisotropic. The
acoustic propagation is computed using the BEM because it is well-suited for infinite media and very accurate.
We show that the matrices arising in both methods are assembled and manipulated just like the native M ATLAB
matrices thanks to operator overloads and only require a minimal implementation effort.
Keywords: FEM-BEM coupling, vibro-acoustics, open-source, G YPSILAB

1INTRODUCTION

The strong interaction between an acoustic wave and an incoming acoustic wave and an elastic structure isa problem which may arise in multiple fields of engineering. A first example may be the interaction of asound wave with some reflective or absorbing material in the design-process of concert halls or anti-noise
walls. Another example is the propagation of underwater acoustic wave produced by submarine machinery.
These waves may go through the hull and may be heard by other submarines, giving away their location.
We study this interaction in the frequency domain. In our context, we suppose that the acoustic medium ishomogeneous and isotropic. This is the case in water or air on reasonably long distances (as long as themedium properties do not change because of the temperature for example). On the other side, the elastic
medium may be fully inhomogeneous and anisotropic, absorbing, dispersive, etc.
There are many possibilities for the numerical computation of such problems. For complex (multi-)physics problem, the most popular method is probably the Finite Element Method (FEM). It is based on the discretizationof the variational formulations associated to the phenomenological equations. The matrices are sparse and canbe easily manipulated and inverted thus enabling the resolution of problems with hundreds of millions of unknowns for the most complex cases. The FEM is widely used in the engineering community and may be foundin many free or commercial solvers. Unfortunately, the FEM has some major drawbacks. A first one is thehigh numerical dispersion resulting from the discretization of the domain and the equation. A second drawback,arising in the case of wave propagation in an infinite medium, is the need to mesh a large domain in orderto mimic an infinite medium by pushing away the exterior boundary. A third drawback inheriting the second
one is the kind of the boundary condition which must be applied on the exterior boundary of the mesh. Onceagain there are multiple possibilities like an approximation of the Dirichlet-to-Neumann (DtN) operator or the
Perfectly-Matched-Layers (PML) which are additional elements with absorbing properties designed to fully absorb the outbound wave without reflection. As a consequence of the assumed complexity of the elastic medium,
the FEM is the best choice for the computation of the propagation of an elastic wave.
As mentioned earlier, for FEM, the propagation of a discretized wave induces numerical dispersion increasing7526

with the distance to the source [10]. Fortunately, this issue may be circumvented by using the aforementioned
BEM [3, 5]. The BEM is based on the discretization of the Boundary Integral Equations related (BIE) to theconsidered wave equation under the assumption that the propagation medium is homogeneous, isotropic andinfinite. They only require the discretization of the boundary of the scatterers and the radiation condition is,by construction, taken into account, thus featuring a minimal numerical dispersion. Of course, the BEM also
have major drawbacks. For example, they require the numerical computation of weakly singular integrals. The
most limiting one is that the arising matrices are dense and non-symetric as the BIEs are inherently convolution
products. In order to compute problems featuring more than a few hundred thousands degrees of freedom, it
is mandatory to compress these matrices. Existing methods are the Fast Multipole Method [7], the hierarchicalmatrices (H -matrices) [8], and more recently the Sparse Cardinal Sine Decomposition (SCSD) [6]. Thesedifficulties are one of the reasons why the BEM is much less popular. The coupling between FEM and BEM
is even less popular.
This paper is organized as follows. First we recall the well-known elasticity problem. In the second part, werecall the BEM and we justify the choice of our formulation. In the third part, we introduce the G YPSILAB
framework, we give the coupling condition between the elastic medium and the acoustic medium and we deducethe variational formulations which will be discretized. We describe the key points for the implementation using
G YPSILAB and the resolution procedure. Finally, we give a validation example.

2

ELASTIC WAVEIn the following, we assume that the elastic medium is homogeneous and isotropic as all the following development remain valid with little changes.
Let Ω ⊂ Rd , d = 2, 3 the elastic medium and Γ its boundary with outbound normal n, we aim at computing thedisplacement u = [u1 , . . . , ud ]T .
Let ω the angular frequency of the incoming acoustic wave, the corresponding elastic-wave equation in thefrequency domain in the absence of source terms is, following [9],div(σ (u)) + ρS ω 2 u = 0,
where ρS is the density of the material. The constraint tensor σ (u) reads asσ (u) = µ ∇u + ∇T u + λ (∇ · u) I,

(1)

(2)

where (λ , µ) are the Lamé coefficients and I is the identity matrix. Please note that this equation may be fully
re-written by making the longitudinal (cL ) and shear wave (cT ) celerities appear by settingµ = ρS · c2T

λ = ρS · (c2L − 2c2T ).

and

(3)

In the context of absorbing media, cT and cL are complex.

3

BOUNDARY INTEGRAL EQUATIONS FOR THE HELMHOLTZ EQUATION

We aim at computing the pressure field p in some propagation domain. We re-use the notations introduced in
section 2. The pressure field may be computed by solving the Helmholtz equation, e.g. [5],−∆p − k02 p = 0+

boundary conditions on Γ,(4)

in the whole propagation domain, where k0 = ω/c0 is the wave number and c0 is the celerity of sound wave inthe propagation medium. The field p must additionaly obey the Sommerfeld radiation condition
lim |r|(d−1)/2 (∂r p − ik0 p) = 0,

r→∞

(5)

7527

where r is the distance to the origin and d the dimension of the propagation medium.
There exists particular solutions of (4) called fundamental solutions or Green kernels G, see [3, 5]. They arei (1)G(x, y) = H0 (k|x − y|)4(1)where H0

for d = 2

andG(x, y) =1 eik|x−y|
4π |x − y|for d = 3,

(6)

is the Hankel function of the first kind of order 0.
(1)

Remark 1. Close to r → 0, H0 (kr) behaves like −1/(2π) log(r).
By multiplying (4) by G, then integrating by parts, using the Cauchy formula and the radiation condition (5),one can show that there exist λ and µ such that we have the following integral representation formulap(x) =Z
ΓG(x, y)µ(y) dγy −∂G(x, y)λ (y) dγy ,Γ ∂ ny

Z

∀x ∈
/ Γ.

(7)

Remark 2. For given λ and µ, p computed using (7) is a solution of (4) obeying (5).
Boundary integral equations are obtained by applying the trace operator to (7). The first step is to decompose
the total field ptot = pi + ps where pi is the known incident field and ps is the unknown scattered field. It is~
not required for pi to follow (5); this is the case when it is a plane wave pi (x) = eik0 ·x .
For example, the scattering by a sound-hard obstacle (homogeneous Neumann boundary condition on Γ) can becomputed by solving the following equation, see for example [5],
∂ pi
(x) =
∂n

∂ 2G
(x, y) λ (y) dγy ,Γ ∂ nx ∂ nyZx ∈ Γ.

(8)

One way to solve equation (8) can be solved by using a Galerkin formulation which reads: find λ such that for
any complex-valued test function wZ∂ pi(x) dγx =k02w̄(x)G(x, y)λ (y)(nx · ny ) dγx dγy −
∂n
Γ×Γ
Z

w̄(x)
Γ

Z
Γ×Γ

(nx × ∇Γ w̄(x)) · (ny × ∇Γ λ (y)) G(x, y) dγx dγy .

(9)
(10)

Some advantages of the BEM appear immediately. In particular, one only need a discretization of Γ and there
is no need to enforce the radiation condition as it is intrinsinc !
In general, the systems arising in the BEM are very ill-conditioned and ill-posed for certain frequencies. Fortunately, some formulations circumvent these limitations at the cost of an increased assembling time. In thefollowing, we will use the Brakhage & Werner [4] representation formula for the scattered field which reads

Z 
∂G
s
p (x) =G(x, y) + iβ
(x, y) α(y) dγy ,
x∈
/ Γ,
(11)
∂ ny
Γ

Z 
∂ ps
∂G
∂ 2G
(x) =
(x, y) + iβ
(x, y) α(y) dγy ,(12)∂n∂ nx ∂ nyΓ ∂ nxor for x ∈ Γ
Z 


∂G
α(x)
p (x) =G(x, y) + iβ
(x, y) α(y) dγy + iβ,
x ∈ Γ,
∂ ny
2
Γ

Z 
∂ ps
α(x)
∂G
∂ 2G
(x) = −+
(x, y) + iβ
(x, y) α(y) dγy ,∂n2∂ nx ∂ nyΓ ∂ nx
s

7528

(13)
(14)

where β is a parameter whose optimal value is β = min(1/k0 , 2) (see [4]) and α some unknown. For convenience, we write
ps (x) = BWα

and

∂ ps
∂ BW
(x) =
α,
∂n
∂n

x ∈ Γ.

(15)

We insist on the fact that α does not need to hold a physical signification as the value of interest is ps which
can be computed in the whole domain using the Brakhage & Werner representation formula (11).

4

COUPLING, VARIATIONAL FORMULATIONS AND RESOLUTION

4.1 A few words about G YPSILABG YPSILAB is an open-source M ATLAB prototyping framework available at [1] under the GPL3.0 licence. It
aims at simplifying the development of numerical methods for multiphysics problems. It focuses mainly on the
FEM and BEM formulations.
The matrices arising in such methods can be computed thanks to multiple overloads of the base integral()function of M ATLAB whose return type depends on the input parameters. In the case of the FEM, the matrices
are automatically sparse. When dealing with BEM, the matrices may be full (in that case the user should bevery careful if the number of unknowns is too high) or in the hierarchical format (H -matrices). In fact, H matrices in G YPSILAB can be considered as an extension to the M ATLAB language since a complete algebra
has been implemented: matrix-vector product, multiplication by a scalar, matrix-matrix multiplication, addition,
addition with a sparse matrix, concatenation, LU-decomposition, etc. Consequently, the user does not need toknow which type he is currently manipulating and the coding effort is minimal. Other possibilities offered by
the software are basic mesh manipulations, mesh refinement with MMG [11], and ray-tracing. Many examples areavailable in the ./nonRegressionTests/ subdirectory of the main G YPSILAB directory. For more details, werefer to [2].
4.2 Coupling
The coupling between the elastic problem and the acoustic problem is done easily by enforcing two continuity
condition at the interface Γ between the two domains. First we enforce the continuity of the displacement at
the boundary. which reads, extrapolating from [5] (2.1.1) p. 9,
ρ0 ω 2 (u · n) =∂ ptot,∂n(16)

where ρ0 is the density of the acoustic medium. Second we enforce the continuity of the constraint which reads
σ (u) · n + ptot n = 0.

(17)

Next, we introduce the test function v for the elasticity problem and w the test function for the acoustic problem.
By testing (16) with w and (17) and (1) with v, we obtain the following coupled system of equations in theGalerkin frameworkZZZσ (u) : ∇v̄ − ρs ω 2 u · v̄ dΩ + (v̄ · n)(BWα) dγ = − (v̄ · n)pi dγ,
(18)
Ω
Γ
Γ


Z
Z
Z
∂ BW
∂ pi
w̄(u · n) dγ − w̄
α dγ = w̄
dγ.
(19)
∂n
∂n
Γ
Γ
Γ
This is a 2 × 2 block-system of equations such that    
A11 A12ub·= 1 ,A21 A22
α
b2

(20)

7529

where the unknowns are u (displacement in Ω, vector) and α (defined on Γ only, no physical signification,scalar). A11 is the pure FEM block which admits a sparse discretization. A22 is the pure BEM block which
admits a full or hierarchical representation. The coupling term A21 is sparse while A12 can be either full or in
the H -matrix format.
We choose to solve the block-system of equations (20) using a Schur-complement. The acoustic-unknown α issolution of

−1
(21)
A22 − A21 · A−1
11 · A12 α = b2 − A21 · A11 · b1 ,which is solved iteratively with the GMRES algorithm using the LU-decomposition of A22 .
4.3 Some aspects of the implementation with G YPSILAB
In this part, we describe some aspects of the implementation of the coupled FEM-BEM problem within the
G YPSILAB framework. Unless specified otherwise, the names of the variables in the code reflect the name given
previously and we assume that the material properties have already been defined. For the sake of simplicity, we
solve the two-dimensional scattering problem.
Pre-processing The first step is to obtain the mesh (in our case it can be generated using G YPSILAB) andextract its boundary, define the integration domain and the finite element spaces. This is done with% meshmesh = mshSquare (N , [L , e ]) ; % N : number of nodes , [L , e ]: dimensions
bnd = mesh . bnd ;% boundary of mesh
% integration domain
Omega = dom ( mesh , 3) ; % integration over mesh , 3 quadrature nodes
Gamma = dom ( bnd , 3) ; % integration over bnd , 3 quad . nodes% finite element spaces
U = fem ( mesh , ' P1 ') ; % linear functions
V = fem ( mesh , ' P1 ') ;a = fem ( bnd , ' P1 ') ;w = fem ( bnd , ' P1 ') ;Assuming the incident wave is a plane wave, it is defined as
Ui = @( X ) exp (1 i * k0 * X * X0 ') ;% plane wave
gxUi {1} = @( X ) 1 i * k0 * X0 (1) .* Ui ( X ) ; % 1 st component of gradient
gxUi {2} = @( X ) 1 i * k0 * X0 (2) .* Ui ( X ) ;where X0 is the direction of propagation.
Assembling operators In order to build the BEM matrices, the user needs to define the Green kernels. Fortunately, some are already implemented.
% Green kernel in 2 D
Gxy = @( X , Y ) femGreenKernel (X ,Y , '[ H0 ( kr ) ] ' , k0 ) ;% gradient components of the Green kernel
gyGxy {1} = @( X , Y ) femGreenKernel (X ,Y , ' grady [ H0 ( kr ) ]1 ' , k0 ) ;gyGxy {2} = @( X , Y ) femGreenKernel (X ,Y , ' grady [ H0 ( kr ) ]2 ' , k0 ) ;One important thing to know is that G YPSILAB does not know wether an unknown or test function is a vector
or a scalar field. Consequently, integrals which do not involve uniquely thescalar product of twoR vector fields
Rshouldbeexpandedcomponentbycomponent.
For
example,theintegral(v̄· n)pi dγ becomes Γ v̄1 n1 pi dγ +ΓRiΓ v̄2 n2 p dγ

7530

rhsFEM = -1 i /4*( integral ( Gamma , ntimes (V ,1) , Ui ) + ...
integral ( Gamma , ntimes (V ,2) , Ui ) ) ;R

In the computation of the full-FEM block, we need to compute integrals of the form µ · Ω ∂1 v2 · ∂2 u1 which
can be implemented easily like
mug1g2 = mu * integral ( Omega , grad (V ,1) , grad (U ,2) ) ;For the BEM matrices, we use once again the integral() function. By passing a tolerance parameter tol, thematrix returned is in the compressed H -matrix format. When assembling the pure acoustic block containing(14), one needs to assemble at some point the variational formulation (9). Since the integrals are singular andintegral() integrates using a Gauss quadrature, we must regularize it by removing the singularity integrated
by quadrature and replacing it with a semi-analytic formula. This is done by calling the regularize() function
which returns a sparse matrix. In our example, it reads
A22 = 1 i /4* ...
( k0 ^2* integral ( Gamma , Gamma , ntimes ( w ) ,Gxy , ntimes ( a ) , tol ) - ...
integral ( Gamma , Gamma , nxgrad ( w ) ,Gxy , nxgrad ( a ) , tol ) ) ;A22 = A22 - 1/(2* pi ) * ...
( k0 ^2* regularize ( Gamma , Gamma , ntimes ( w ) , '[ log ( r ) ] ' , ntimes ( a ) ) - ...
regularize ( Gamma , Gamma , nxgrad ( w ) , '[ log ( r ) ] ' , nxgrad ( a ) ) ) ;where A22 is a H -matrix. When tol is removed, the result is a full matrix.
The source code for the complete construction of the different matrices and right hand sides is available inthe corresponding folder in ./nonRegressionTest/ within the G YPSILAB directory. It covers many of thepossibilities covered by the framework.
Solving the problem The Schur complement (21) is implemented as follows
A11 = decomposition ( A11 ) ; % decomposition of A11
rhs = b2 - A21 *( A11 \ b1 ) ;[ LL , UU ] = lu ( A22 ) ;% lu of H - matrix A22
prec = @( V ) UU \( LL \ V ) ;% definition of preconditioner
lhs = @( V ) A22 * V - A21 *( A11 \( A12 * V ) ) ; % matrix - vector productalpha = gmres ( lhs , rhs ,[] , tol ,100 , prec ) ; % compute solution
The Brakhage-Werner formulation being well-conditioned, the use of the preconditioner may not be mandatory.
Post-processing In the previous step, we computed the unknown α. We assume that the post-processing
consists in computing the total field at a given point xs in the propagation domain sufficiently far from theboundary so that there is no need for ragularizing. Other kinds of post-processing could be carried out inthe same fashion. We will use the representation formula (11) implemented by first computing the "radiation
matrix" rad, then using it to compute the scattered field ps. Finally, we add the incident field to obtain the
total field ptot. The implementation is as follows
rad= 1 i /4*( integral ( Xs , Gamma , Gxy , a ) + ...
1 i * beta * integral ( Xs , Gamma , gyGxy , ntimes ( a ) ) ) ;ps= rad * alpha ;% scattered field at Xs
ptot = ps + Ui ( Xs ) ; % total field at Xs

7531

~y
~x

Teiky
500 mmeiky + Re−iky

Figure 1. Left: Setting for the numerical computation of the transmision and reflexion coefficients T and R.
Right: cut-off function for the real part of the incident plane wave associated to the measurement points (bluefor reflexion and right for transmission).

5

NUMERICAL EXAMPLE

We present here a validation example available with G YPSILAB in 2D and 3D. The goal is to compute dereflexion and transmission coefficients R and T for a three-layered medium: the incoming wave propagates
through three successive media. A part is reflected and transmitted at each interface. By definition, for anormalized incident wave, the transmission coefficient is the amplitude of the total field in the third medium
and the reflexion coefficient is the amplitude of the backward-propagating wave in the first medium. In our
setting, we suppose moreover that the first and the third medium have the same physical properties.
Of course, the theoretical problem is purely one-dimensional and the purpose of our study is to show that one
can recover the analytical values by solving a two-dimensional problem for a slab with the correct hypotheses.
Steel - d=500 mm

Reflexion coefficient

0Transmission coefficient

0
Analytic
Gypsilab

Analytic
Gypsilab
−5

−4

−10

R (dB)T (dB)

−2

−6

−15

−8

−20

−10

0

2000

4000
6000
Frequency (Hz)

8000

10000

−25

0

2000

4000
6000
Frequency (Hz)

8000

10000

Figure 2. Top left: scattered field at 10000 Hz (used to compute the reflexion coefficient). Top right: Total fieldat 10000 Hz (used to compute the transmission coefficient). Bottom: Reflexion and transmission coefficients fora steel plate with thicknes d = 500 mm in the range 100 − 10000 Hz.

7532

First we set the thickness of the slab and a frequency. The length of the plate is a fixed multiple of the
wavelength. In order to avoid acoustic scattering at the corners of the slab, we localize the plane wave by
multiplying it with a cut-off function centered on the middle of the slab in the propagating direction whose
bandwidth is significantly smaller than the width of the slab but is also a fixed multiple of the wavelength. We
also set the shear-wave celerity cT = 0 (see (3)) since we forbid any displacement along the longest dimensionof the slab (direction ~x on Fig. 1). The exterior medium is water (c0 = 1500 m.s−1 and ρ0 = 1000 kg) and theelastic medium is steel (cL = 6000 m.s−1 and ρs = 8000 kg).
The coefficients T and R are computed in the range f ∈ [100, 10000] Hz by probing the total field at a distance
one wavelength above the slab, resp. one wavelength below. The results are represented on Fig. 2. We observethat the analytical computation and the G YPSILAB computation overlap thus validating our setting.

6

CONCLUSIONSIn this paper, we showed how one can efficiently implement strongly coupled FEM-BEM problems by combining the powerful high-level interfaces provided by both the native M ATLAB API and the G YPSILAB framework.
We demonstrated that the variational formulations can be easily implemented and the development time is short
making it a prime framework for fast prototyping. The three-dimensional case has also been validated but requires much more computational time as the number of degrees of freedom increases much faster with thefrequency than in the two-dimensional case. It is not presented in the context of this paper as we lack space
but the implementation is done similarly.
In the future, this work will be applied to more complex structures of research or industrial dimension. Some ofour topics of interest are the design of music instruments, acoustics of underwater structure or non-homogeneous
anisotropic media.

REFERENCES[1] https://github.com/matthieuaussal/gypsilab[2] Alouges, F.; Aussal, M. FEM and BEM simulations with the GYPSILAB framework, SMAI-JCM, Vol. 4,
2018, pp. 297–318[3] Colton, D. ; Kress, R. Integral Equation Methods in Scattering Theory, SIAM Classiccs in Applied Mathematics, 2013[4] Kress, R Minimizing the condition number of boundary integral operators in acoustic and electromagnetic
scattering, Q. J. Mech. Appl. Math., Vol. 38 (2), 1985, pp. 323–341[5] Nédélec, J.-C. Acoustic and Electromagnetic Equations, Integral Representations for Harmonic Problems,
Springer Applied Mathematical Sciences 144, 2001
[6] Alouges, F.; Aussal, M. The Sparse Cardinal Sine Decomposition and its Application to fast numericalconvolution, Numerical Algorithms, Vol. 70 (2), 2015, pp. 427–448[7] Greengard, L. The rapid evaluation of potential fields in particle systems, MIT Press, 1988
[8] Börm, S.; Grasedyck, L.; Hackbusch, W. Hierarchical Matrices, technical. report, 2015[9] Skelton, E. A.; James, J. H. Theoretical Acoustics od Underwater Structures, Imperial College Press, 1997[10] Terrasse,
I.;
Abboud,T.
Modélisationdesphènomènes
https://www.math.u-bordeaux.fr/ abachelo/Abboud-Terrasse.pdf

depropagationd’onde,[11] https://www.mmgtools.org/7533