M.N. Ichchou Sbastien Besset

L Jézéquel

A coupled BEM and energy flow method for mid-high frequency internal acoustic

Keywords: Energy flow, BEM, Energy method, Injected power density, Energy density, Intensity

This paper deals with the behaviour of acoustic cavities in the mid-high frequency range. The method proposed here is based on an energy flow method named Simplified Energy Method (MES). MES method is quite efficient in the mid-high frequency range but the directivity of the boundary sources is not well estimated. We propose a hybrid method which couples MES and the Boundary Element Method (BEM). Thus, the BEM method is used to estimate the direct field, considering a "correct" directivity. As a complete calculation is not adapted to BEM in mid-high frequency range because of the calculation costs, we only apply BEM on the domains including boundary sources. Other parts of the system and the reverberated field are estimated by the mean of MES method. This hybrid method leads to a consistent prediction of injected power densities. Numerical comparisons prove the efficiency of the proposed reformulation.

Introduction

Medium and high frequency range is a fundamental topic of research in the context of viroacoustics. In fact, well established and deeply studied numerical methods like finite element or boundary equation based formulations run into important limitations as frequency increases. The most important limitation is the size of the meshing that must be very fine and leads to high calculation costs. Modal methods base on finite element formulations have been develop [START_REF] Besset | Dynamic Substructuring based on a Double Modal Analysis[END_REF][START_REF] Besset | Vibroacoustical Analysis Based on a Multimodal Strategy: Triple Modal Synthesis[END_REF], but the calculation costs remain high and these methods are limited to medium frequency range.

Among possible alternatives, analysis methods which are energy variables based are often studied. The Statistical Energy Analysis (SEA) [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF] is the recognized mother and gives the mechanical energy of complex buildup structures. However SEA suffers from its axiomatic formulation and still requires improvements. Beyond these studies, a number of works attempts to enhance the Statistical Energy Analysis robustness and 1 predictivity. Among those tentatives, lets report the earlier work of Nefske and Sung [START_REF] Nefske | Power flow finite element analysis of dynamic systems: Basic theory and application to beams[END_REF] who proposed the use of an energy diffusion model to predict the space spread of energy density within subsystems. Basically, this way of thinking can be viewed as a local energy formalism whereas the SEA formalism is based on global energies of finite subsystems. This model has been improved in [START_REF] Wohlever | Mechanical energy flow models of rods and beams[END_REF][START_REF] Bouthier | Simple models of energy flow in vibrating membranes[END_REF][START_REF] Bouthier | Simple models of energy flow in vibrating plates[END_REF][START_REF] Lase | Energy analysis of bars an beams: Theoretical formulations[END_REF][START_REF] Bot | Equations énergétiques en mécanique vibratoire: application au domaine des moyennes et hautes fréquences[END_REF][START_REF] Ichchou | Formulations énergétiques pour l'étude moyennes et hautes fréquences des systèmes : théorie et applications[END_REF]. Among other contributors to this subject, let mention references [START_REF] Langley | On the vibrational conductivity approach to high frequency dynamics for two-dimensional structural components[END_REF][START_REF] Carcaterra | Energy density equations and power flow in structures[END_REF][START_REF] Carcaterra | Wavelength scale effects on energy propagation in structures[END_REF][START_REF] Xing | A power flow analysis based on continum dynamics[END_REF] where interesting discussions are given. In [START_REF] Lase | Energy analysis of bars an beams: Theoretical formulations[END_REF][START_REF] Bot | Equations énergétiques en mécanique vibratoire: application au domaine des moyennes et hautes fréquences[END_REF][START_REF] Ichchou | Formulations énergétiques pour l'étude moyennes et hautes fréquences des systèmes : théorie et applications[END_REF] energy models started with the so-called general energy method. The main goal of this method was to reformulate the classical displacement models using four energy variables: the total energy as well as the Lagrangian energy density, the active and the reactive energy flow. From the elimination of the Lagrangian energy density and the reactive energy flow, this general energy formulation leads to an interesting energy model well suited for medium and high frequency dynamics. In fact, from a propagative approach, it has been shown that the Lagrangian energy density and the reactive energy flow is mainly linked to wave interferences or singularities. So that, taking into account only the incoherent contribution of waves to the energy variables leads to the formulation of the local energy approach concerned here. References [START_REF] Bot | Equations énergétiques en mécanique vibratoire: application au domaine des moyennes et hautes fréquences[END_REF][START_REF] Ichchou | Formulations énergétiques pour l'étude moyennes et hautes fréquences des systèmes : théorie et applications[END_REF][START_REF] Moron | Mise en oeuvre de méthodes énergétiques dans létude du comportement vibratoire de structures complexes dans le domaine des moyennes et hautes fréquences[END_REF][START_REF] Houillon | Modélisation vibratoire des carrosseries automobiles en moyennes et hautes fréquences[END_REF][START_REF] Bocquillet | Méthodes énergétiques de caractérisations vibroacoustiques des réseaux complexes[END_REF] give in depth the formulations used in the context of local energy approach applications. The MES has already been applied to various domains: beams [START_REF] Lase | Energy analysis of bars an beams: Theoretical formulations[END_REF][START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF], membranes and plates [START_REF] Bouthier | Simple models of energy flow in vibrating membranes[END_REF][START_REF] Bouthier | Simple models of the energetics of transversely vibrating plates[END_REF][START_REF] Langley | On the vibrational conductivity approach to high frequency dynamics for two-dimensional structural components[END_REF][START_REF] Lase | Energy analysis of bars an beams: Theoretical formulations[END_REF][START_REF] Moron | Mise en oeuvre de méthodes énergétiques dans létude du comportement vibratoire de structures complexes dans le domaine des moyennes et hautes fréquences[END_REF] and acoustical radiation [START_REF] Cotoni | Modélisation de phénomènes vibroacoustiques en moyenne et hautes fréquences par méthode énergétique locale[END_REF][START_REF] Cotoni | High frequency radiation of l-shaped plates by a local energy flow approach[END_REF][START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF]. In this paper, the main concern is the improvement of MES prediction for internal acoustics [START_REF] Hardy | Absorption coefficient and energy flow path identification by means of inverse local energy method[END_REF].

The main addressed issue in this work is the injected power density for use in the MES formulation. The injected power density is usually described from the input acoustic mobility. The later is easily expressed for unbounded media. A major concern in the MES prediction quality is the input source directivity which needs to be fixed. So the major question this paper tries to solve is the concise estimation of injected power density. So far, the MES energy method needs an estimation of the entering power, and its directivity, which cannot easily be estimated and is usually considered to be lambertian. This is an assumption that is not always verified and the proposed hybrid method aims at estimating correctly this directivity. Particularly, the computation of the direct field is the key aspect. A coupled BEM and MES strategy is chosen. The algebraic coupled formulation is expressed and implemented. Numerical experiments was achieved comparing the coupled strategy with the conventional one. The interest and the feasibility of the finding is proved. It must be noticed here that authors already studied the coupling between deterministic and statistical subsystems. Among them, Shorter and Langley [START_REF] Shorter | Vibro-acoustic analysis of complex systemsr. s[END_REF] considered the coupling between several subsystems, each subsystem being described with deterministic or statistical methods. This work has lead to coupling methods that are currently used in softwares as VAOne. Nevertheless, the method we propose does not deal with coupling substructures, but allows to describe a single structure considering a mixed formulation.

First, the Simplified Energy Method is recalled. The hybrid formulation is then presented. Principle of the method is given and a matrix formulation is proposed. At last, numerical calculations are processed to show the efficiency of the method.

Simplified Energy Method (MES)

Assumptions

In view of describing the energy transfer inside the medium, two continuous energy fields are introduced. The first energy quantity is nothing else than the total energy density W ( s, t) defined as the sum of the potential energy density and the kinetic energy density. The second energy variable I( s, t) is the energy flow. These quantities are local in opposition with energies per subsystem involved in SEA. The energy balance equation from continuum mechanics obeying conservation principle, which govern the energy density in various vibrating system can be described by:

∂W ∂t = -π diss -∇ • I ( 1
)
where ∇ is the gradient operator, π diss is the power density dissipated. The damping model adopted here is the same as in SEA (power density being dissipated is proportional to the energy density). Hence:

π diss ( s) = ηωW ( s) (2)
where η is the damping loss factor and ω the circular frequency. The validity of this relationship has been discussed in the literature [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF] about SEA. Equation ( 1) is the energy balance relation for all elastoacoustic media and it is regarded as a valid method for steady or transient analysis. In order to derive the energy density equations, a wave description of vibrational-acoustical behaviour is definitely adopted. In the subsequent presentation planes and symmetrical propagating disturbances in a medium are considered.

Each propagating waves involve partial energy quantities defined as the energy variables associated to those waves. In order to establish a relationship between the partial energies W i ( s) and I i ( s) and the total ones W ( s), I( s), an additional assumption is introduced. In fact, a superposition principle is assumed and applied here for all quadratic variables:

£( s) = i £ i ( s) (3)
£, £ i being a quadratic variables, representing either a partial energy or a global energy associated to a particular wave field. This assumption is often used in statistical phenomenon in physics. We can find such an assumption in room acoustics, electromagnetism,. . . In the following developments an intrinsic energy law will be used. In fact, this superposition principle is valid if a term due to the phase is neglected. Considering two waves denoted v 1 and v 2 , we can write the following equation:

|u 1 + u 2 | 2 = |u 1 | 2 + |u 2 | 2 + u * 1 u 2 + u 1 u * 2 neglected terms (4)
These terms can be neglected because the quantities considered in MES are frequency average quantities. Therefore, we can write:

< u * 1 u 2 >= 0 and < u 1 u * 2 >= 0 (5)
where < u.v >= 1 ∆f ∆f u.v df , ∆f being the considered frequency range for the frequency average. This will mean that MES does not take into account modal information (linked to the phase of the waves).

The following equation [START_REF] Brillouin | Wave propagation in periodic structures[END_REF][START_REF] Morse | Methods of theoretical physics[END_REF] is often introduced in order to define the wave velocity c. In fact:

I i ( s) = c • W i ( s) n (6)
This expression was shown to be a valid law for non-damped waves field.

Elementary fields: planes waves

In this section, a wave approach very similar to the one given in [START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF], will be used. Hence, expression (6) will be first introduced in the energy flow balance (1), leading to:

c ∇ • W ( s) n + ηωW ( s) = 0 (7)
Equations 6 and 7 lead then to the following relation:

c ∇W ( s) n. n + ηω c I( s) = 0 (8)
Thus, it can be readily shown using planes waves characteristics that the total energy density and the active energy flow are linked by:

-c 2 ∇W ( s) = ηω I( s) (9)
∆W ( s) -η 2 ω 2 W ( s) = 0 (10)
where ∆ = ∇ • ∇. This is the local energy equation for an elementary planes waves field. This is an hyperbolic equation of telephone format which reduces to the simple wave equation when no dissipation mechanics is introduced. This equation needs energy boundary conditions and coupling conditions to be defined.

Elementary fields: symmetrical waves

Let us consider now a symmetrical waves field. Namely, we will consider cylindrical waves in dimension 2, and spherical waves in dimension 3. The notations used in this section are recalled on figure 1. I d , W d is the direct field whereas I r , W r is the reverberated field. n is the unit normal vector at boundary. In this cases, combining ( 6) and (1) leads to:

1 r n-1 ∂ ∂r r n-1 I + ηωW n = 0 (11)
using [START_REF] Langley | On the vibrational conductivity approach to high frequency dynamics for two-dimensional structural components[END_REF], one can readily show that:

-c 2 1 r n-1 ∂ ∂r r n-1 W n = ηω I (12)
Φ inp Ω ∂Ω I r I d I = I d + I r W = W d + W r n Figure 1: Notations used in section 2.3
Putting this intrinsic law in the energy flow balance (1) leads to:

1 r n-1 ∂ ∂r r n-1 W -η 2 ω 2 W = 0 (13)
This is the local energy equation for symmetrical waves.

Figure 2 shows two distinguish manner of energy reconstruction. In both techniques, energy variables are defined as the supperposition of a direct field and a reverberant field. In both techniques, the direct field is calculated assuming a symmetrical waves fields like described before. The difference comes from the reverberant field treatment.

W = W d + W r (14)
and

I = I d + I r (15)
Precisely, in figure 2(a), the reverberant field is assumed to be a set of planes waves. This approach is known as the hybrid energy flow method. Here we will designate it by the mixed wave solution. This solution was proposed by Smith [START_REF] Smith | A hybrid energy method for predicting vibrational response of point loaded plates[END_REF]. Further developments and considerations are given in [START_REF] Ichchou | Transient local energy approach for multi-dimensional, a review with extensions[END_REF]. The second approach considers the reverberant filed as the results of secondary sources located at boundaries as shown in figure 2(b). So, only symmetrical waves fields is considered. This solution was first proposed by Kuttruff [START_REF] Kuttruff | Simulierte nachhallkurven in rechteckräumen mit diffusem schallfeld[END_REF][START_REF] Kuttruff | A simple iteration scheme for the computation of decay constants in enclosures with diffusely reflecting boundaries[END_REF][START_REF] Kuttruff | Room acoustics[END_REF][START_REF] Kuttruff | Energetic sound propagation in rooms[END_REF]. It is also known as the radiosity method. Furthers contributions and details can be found in [START_REF] Carroll | Decay of reverberent sound in a spherical enclosure[END_REF][START_REF] Miles | Sound field in a rectangular enclosure with diffusely reflecting boundaries[END_REF] and [START_REF] Bot | A vibroacoustic model for high frequency analysis[END_REF]. It is this second approach that we will use in this paper.

Here, both the direct and the reverberant field will be represented by symmetrical waves. The elementary solutions in terms of energy density and active intensity will be denoted G and H. and r = --→ P M . Hence, using a superposition principle the energy solution is then expressed using the contribution of the primary source (direct field) and the contribution of the fictitious sources. The fictitious sources are energy unknown parameters representing the contribution of the reverberant filed. Precisely:

G( --→

P M ) = 1 γ 0 c e -ηωr r n-1 H( --→ P M ) = 1 γ 0 e -ηωr r n-1 u r (16)
W (M ) = Ω ρ(P ) G( --→ P M ) dP + ∂Ω σ(P ) u r . n G( --→ P M ) dP (17)
The power flow can be described by a similar expression:

I(M ) = Ω ρ(P ) H( --→ P M ) dP + ∂Ω σ(P ) u--→ P M . n H( --→ P M ) dP (18)
where ρ(P ) corresponds to volumic sources and σ(P ) corresponds to surfacic sources. Energy variables are given as a solution of a Fredholm equation, corresponding to an energy balance at the boundary of the domain [START_REF] Kuttruff | Simulierte nachhallkurven in rechteckräumen mit diffusem schallfeld[END_REF].

Discretized formulation

Let us consider an isolated system as shown in figure 3.

The boundaries of the cavity are discretized (N e elements). Properties are given for each element i of the meshing: • σ i (ω) is the reverberated power flow;

• Φ inj i (ω) is the input power flow; Φ inj Ω ∂Ω I rev I dir
• α i (ω) is the absorption coefficient.

The input power flow Φ inj i (ω) is injected in the cavity. In order to calculate I and W in the cavity, the directivity of these input power flows must be known. For high frequency, we assume the following directivity:

d(θ) = cos θ π ( 19
)
where θ is defined in figure 4. There is no way to evaluate the "correct" directivity through MES method, and the lambertian directivity d(θ) = cos θ π is not always a good approximation. The hybrid method proposed in the paper should aim at considering the correct directivity to calculate the direct field.

W dir (M ) = Ne i=1 Φ inp i c 0 Si d i (θ) r 2 dP ( 20
)
I dir (M ) = Ne i=1 Φ inp i Si d i (θ) r 2 u r dP (21)
where r = PM and u r = --→ MP --→ MP (see figure 1. S i corrsponds to the facet number i.

The reverberated field can be expressed the same way:

W rev (M ) = Ne i=1 σ inp i c 0 Si d i (θ) r 2 dP ( 22
)
I rev (M ) = Ne i=1 σ inp i S d i (θ) r 2 u r dP (23)
Considering equations 14 and 15, it is then possible to write W and I. In the following, we will consider I = I. n I where n I is a given vector, which corresponds to measurements of intensity among vector n I .

where P is a point situated on an element S i . d i (θ) has been defined in equation [START_REF] Bouthier | Simple models of energy flow in vibrating membranes[END_REF]. N e is the number of elements.

Considering one element i, the power balance can be expressed as follows:

σ i = 1 -α i S i    Ne j=1 j =i Φ inp j Si Sj cos θ i d j (θ j ) r 2 ij dP i dP j + Ne j=1 j =i σ j Si Sj cos θ i d j (θ j ) r 2 ij dP i dP j    (24)
This equation can be elaborated considering incoming and outcoming energy on the considered facet.

Matrix equations

Equation 24 can be written for i = {1, 2, . . . , N e } (N e equations can be written). These equations can be expressed using matrices as follows:

Mσ = B ( 25
)
with the following notations:

M = Id -(Id -A) T (26) B = (Id -A) TΦ inp (27)
T and A are defined as follows:

T ij = 1 S j Si Sj d i (θ i ) cos θ j r 2 ij dP i dP j (28)
A ij = δ ij α i (29)
Notice that T must verify two properties, due to the conservation of energy:

Ne j=1 T ij = 1 ∀i ( 30
) Ne i=1 S i S j T ij = 1 ∀j (31)
Equation 25 allows to write the expression of σ:

σ = M -1 B = M -1 (Id -A) TΦ inp (32)
Direct and reverberated fields expressed in equations 20, 21, 22 and 23 can now be written as follows:

W = Y W Φ inp + Y W σ (33) I = Y I Φ inp + Y I σ (34)
Considering equation 32, W and I can be expressed as follow:

W =          W 1 W 2 . . . W Ne          = S W Φ inp (35)
I =          I 1 I 2 . . . I Ne          = S I Φ inp (36)
where:

S W = Y W Id + M -1 (Id -A) T (37)
S I = Y I Id + M -1 (Id -A) T (38)
3 Hybrid Simplified Energy Method

Principle

The main issue of the Simplified Energy Method is the unknown directivity of the sources Φ inp . In the previous section, it is assumed to be d(θ) = cos θ π (see equation [START_REF] Bouthier | Simple models of energy flow in vibrating membranes[END_REF]. The aim of the present section is to improve Simplified Energy Method in order to take into account the genuine directivity. Thus, Hybrid Simplified Energy Method (HMES) is based on a BEM calculation for the direct field. The reverberated field is calculated as for Simplified Energy Method, as shown in figure 5.

The hybrid method allows to combine the advantages of MES and BEM. Indeed, calculating the entire problem through BEM would lead to very high calculation costs. That is why BEM will only be applied around the boundary sources. The vibrating structure is first studied using a finite element code. The pressure and celerity fields on the structure is then computed in order to obtain the pressure field for the acoustic domain using a BEM code. Figure 6 shows the various fields considered in HMES. HMES is the computed using the fields denoted I dir (BEM) and I rev (MES) (an equivalent figure can of course be plot for W fields):

FEM calculation

⇒ v on the structure BEM calculation (free field) ⇒ p, v MES calculation (σ ⇒ reverberated field) W dir , I dir W rev , I rev W , I v p, v ⇒ Φ inp
I tot (HMES) = I dir (BEM) + I rev (MES) (39) W tot (HMES) = W dir (BEM) + W rev (MES) (40) Φ inj Ω ∂Ω I rev (MES) I dir (BEM) I dir (MES)
Figure 6: Considered system

Direct field calculation

The input power Φ inp is due to a structure that is first studied using a FEM software. This calculation provides the celerity field on the structure. A BEM software is then used to obtain the pressure field on the structure, considering anechoic conditions. Thus, pressure and celerity are both known on the structure and it is possible to calculate the pressure field in the cavity, taking the directivity into account. Let (p I , v I ) be the pressure and celerity fields on the structure. These fields can be provided by any FEM software and a BEM calculation. A BEM calculation is then performed to obtain the pressure field p(M ) in the cavity:

cp(M ) = S p I (P ) ∂g ∂ n -g -→ MP ∂p I ∂ n dP, P ∈ S (41)
Notice that ∂pI ∂ n = -ρω v I . n, where ρ is the density of the fluid. n is the surface unit normal vector anc c is defind as follows:

c =    1, M ∈ Ω 1 2 , M is on a smooth boundary of Ω S 4π , M is on a nonsmooth boundary of Ω (42)
where S is the solid angle. g is the free-space green function defined as follows:

g -→ MP = e -ik -→ MP 4π -→ MP (43)
The pressure field in the cavity allows to calculate the direct fields W dir and I dir :

W dir (M ) = 1 2ρc 2 p(M )p * (M ) (44)
I dir (M ) = 1 2 ℜ [p(M )v * (M )] (45)
v = v. n p (and v * ) is obtained thanks to p using equation ∂p ∂ np = -ρω v. n p where n p is a freely chosen unit vector. Thus I corresponds to the projection of I among the direction n p .

Reverberated field calculation

The reverberated field is being calculated using classical MES theory described in section 2. The input power Φ inp is provided by the pressure and celerity fields on the structure:

Φ inp = 1 2 ℜ(p i v * i ) (46)
σ is then calculated, which allows to obtain W rev and I rev , as explained in equations 22 and 23.

Coupled matrix formulation

Using a matrix formulation, the pressure computed by BEM can be written as follows:

p I = H p p str + G p v str (47)
where H p and G p are BEM matrices, and p str and v str are the pressure and celerity fields on the structure (computed by a finite element code). An equivalent equation can be written for the celerity field:

v I = H v p str + G v v str (48)
Thanks to equations 44 and 45, equations 47 and 48 lead to the following equations concerning the direct field:

W dir n = 1 2ρc 2 p * str v * str Λ 1 p str v str (49)
I dir n = 1 2 ℜ p * str v * str Λ 2 p str v str (50)
where W dir n is the n th value of W dir and I dir n is the n th value of I dir . Λ 1 and Λ 2 are defined as follows:

Λ 1 = H * p H p H * p G p G * p H p G * p G p (51)
Λ 2 = H * p H v H * p G v G * p H v G * p G v (52)
Concerning the reverberated field, W rev and I rev are obtained using equations 35 and 36:

W rev n = 1 2 S W n ℜ [(pv * ) str ] (53)
I rev n = 1 2 S I n ℜ [(pv * ) str ] (54)
where S W n is the n th line of matrix S W n and S I n is the n th line of matrix S I n . W tot and I tot can then be written:

W tot n = p * str v * str Λ1 p str v str (55)
I tot n = ℜ p * str v * str Λ2 p str v str (56)
where Λ1 and Λ2 are defined as follows:

Λ1 = 1 2ρc 2 H * p H p H * p G p G * p H p G * p G p + 1 4 0 diag S W n diag S W n 0 (57) Λ2 = 1 2 H * p H v H * p G v G * p H v G * p G v + 1 4 0 diag S I n diag S I n 0 (58)

Results

The structure and the cavity used for the examples are shown on figure 7. Only one panel is able to move (the "structure"). The others are considered to be infinitely rigid (i.e.

-→ grad p = 0). The calculations have been made using MATLAB and OPENBEM [START_REF] Juhl | The boundary element method for sound field calculations[END_REF]. Of course, the meshing of the structure has to be quite fine in order to provide good FEM and BEM results. On the contrary, the MES calculation does not need fine meshes.

W is calculated using classical MES (figure 8) and the method we propose (figure 9). Concerning the BEM calculation, a frequency average has been done. The data used for this calculation are recalled in table 2.

Notice that the dimensions of the box allow to justify the mid-high frequency range situation. Indeed, the acoustic modes are situated at the following frequencies: 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000 0000000000000000000000000000000000000

f = c 2 n x L x 2 + n y L y 2 + n z L z 2 (59)
The modal density around f = 2000Hz can then be evaluated to d = 30 modes/Hz, which is enough to consider high frequency domain.

Figure 10 show a BEM calculation taking into account the impedance on the boundaries. This calculation is considered to be "exact". Considering these results, it appears that HMES gives much better results that MES: taking into account the right directivity of the source allows to provide more accurate results. MES remains quite efficient if only the W average in the cavity is considered, but locally results are far from "exact" results (BEM results), especially around the source point.

Concerning the calculation costs, table 1 shows the differences between MES, BEM and the hybrid method.

Method

Conclusion

The method proposed in this paper allows to predict the acoustical behaviour of a cavity. This energy method is more precise that existing methods like MES or SEA because it is able to take into account the directivity of the sources by FEM/BEM calculations without meshing the cavity. Moreover, it is possible to take the damping into account, even if it is not low. Compared to a full BEM or FEM calculation, the method we propose only requires finer meshing around the sources. The other parts of the cavity do not need finer meshing. The advantages and disadvantages of the method are The results show that the total field is improved, even in a simple case of a box as presented in figure 7. Nevertheless, the calculations are longer than for the MES method because of the FEM and BEM calculations. Moreover, unlike MES method, the BEM need quite fine meshes and several calculation because of the frequency average.

Acknowledgement

This work was made within the framework of the European project CREDO. The authors thank all the partners of this project.

Figure 2 :

 2 Figure 2: Energy fields reconstruction: (a) mixed wave assumption and (b) fictitious sources at boundaries.

Figure 3 :

 3 Figure 3: Considered system

Figure 4 :

 4 Figure 4: Notation for directivity It is now possible to define W dir (M ) and I dir (M ) in function of Φ inp . Considering that η << 1 in the fluid, equation 13 leads to the following equations:

Figure 5 :

 5 Figure 5: HMES method

Figure 7 :

 7 Figure 7: Structure used for the calculations

12 Figure 8 : 12 Figure 9 :

 128129 Figure 8: W computed by MES

12 Figure 10 :

 1210 Figure 10: W computed by BEM

Table 1 :

 1 Comparison between the methods -calculation costs

		Number of elements	Calculation cost	Error (% vs BEM method)
	MES		around 100	a few minutes	around 10% in front of the source
					50-200 % far from the source
	BEM	more than 5000	2 days	reference
	Hybrid method		around 500	around 15 minutes around 10% in front of the source
					10-40 % far from the source
	ρ	c	Dimensions	Frequency	α (absoption coefficient)
	1, 25 kgm -3 340 ms -1 2 × 2 × 4 m 3 1500 -2500 Hz	0, 3

Table 2 :

 2 Data used for calculations


Mohamed A Chabchoub

Sébastien Besset

M N Ichchou

Structural sources identification through an inverse mid-high frequency energy method

Keywords: Energy flow, MES, Forces identification, Forces localization, Injected power density, Energy density, Intensity

come

Introduction

The identification of force features that serve to excite structural components is a key concern in the field of structural engineering. The localization and quantification of input loads is of considerable interest; frequently, the direct measurement of sources proves quite complicated to perform. Indirect means for defining excitation often provides an effective alternative. A tremendous amount of work has been performed and published regarding this issue [START_REF] Liu | Input force estimation of a cantilever plate by using a system identification technique[END_REF][START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF][START_REF] Ji | A study on an estimation method for applied force on the rod[END_REF][START_REF] Pezerat | Identification of vibration sources[END_REF][START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF][START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF][START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF][START_REF] Noiseux | Measurement of power flow in uniform beams and plates[END_REF][START_REF] Pavic | Measurement of structure borne wave intenisty: Part I, formulation of the methods[END_REF][START_REF] Hamdi | Hybrid inverse boundary element method for the determination of the optimal spectral characteristics of complex radiating noise sources[END_REF][START_REF] Bai | Application of BEM-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries[END_REF]. The problem being considered herein may also have applications in structural health monitoring. Damage, whose features often remain unknown, leads to changes in the vibration energy distribution of structural members. Both the magnitude and position detection of damage thus become key parameters.

In general, two approaches are available in the literature relative to input force identification aspects. The first has more or less been inspired by identification approaches for the black box-like behavior of dynamic systems, coupled with the introduction of estimation theories. As an example, in references [START_REF] Liu | Input force estimation of a cantilever plate by using a system identification technique[END_REF][START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF], the author uses an online recursive inverse method in order to estimate the input forces of beam and plate structures. The inverse method has been based on the Kalman filter plus a recursive least-squares algorithm that includes just a small number of measurement responses within the structure. Several types of input forces have been successfully identified. This method performs effectively in estimating the input forces of beam structural systems from noisy measurements. It has been shown [START_REF] Ji | A study on an estimation method for applied force on the rod[END_REF] that a Kalman filter method may be used to estimate a timewise variation rod force. The second approach consists of defining the input force that allows the measured set of inputs to verify the equilibrium equations given in either continuous or discretized form. For example, a number of authors have developed a pertinent strategy, called the RIFF method, for identifying vibration sources [START_REF] Pezerat | Identification of vibration sources[END_REF][START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF][START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF]. The underlying principle entails both the discretization of motion equations via a finite difference scheme and the computation of force distribution from a set of displacements measured on the structure. Regularization is then achieved through a spatial window and a wave number filter. Feasibility studies were conducted for beam and plate elements, and an application to cylindrical shells has been proposed [START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF]. Recently, this approach was coupled to near-field acoustic holography in order to make a successful identification from acoustic measurements [START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF]. Perotin and Granger offered another approach, whereby a modal model of a structure along with a spatial orthonormal decomposition of the excitation field are used in order to derive the force identification technique [START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF]. The problem associated with the vibration input force localization has also been addressed in the structural intensity community; moreover, instead of applying an analytical model [START_REF] Pezerat | Identification of vibration sources[END_REF] or modal model [START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF] in the inversion process, discretized finite element (FE) or boundary element (BE) can be introduced into the identification [START_REF] Hamdi | Hybrid inverse boundary element method for the determination of the optimal spectral characteristics of complex radiating noise sources[END_REF][START_REF] Bai | Application of BEM-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries[END_REF]. Noiseux [START_REF] Noiseux | Measurement of power flow in uniform beams and plates[END_REF][START_REF] Pavic | Measurement of structure borne wave intenisty: Part I, formulation of the methods[END_REF] wound up using intensity vectors to indicate the energy flow and source localization.

It should be noted however that most of what has been proposed in the literature focuses on a limited frequency band, more precisely the so-called low-frequency region where modal overlap is weak and where resonant behavior of the structure still dominates. As frequency increases, existing methods lose their advantage. The first such drawback pertains to the robustness of strategies offered relative to noise and uncertainties. It is expected that approaches based on the spatial discretization of data will be very sensitive to discretization errors over short wavelengths. A costly and time-consuming treatment is thus expected in order to deal with high-frequency band features. Moreover, identification strategies that use finite element, boundary element or simplified analytical modeling of the targeted structure are naturally confronted with intrinsic limitations when modeling in the mid-high frequency regime. The need is therefore expressed to establish an identification tool that takes into account both medium and high-frequency particularities. The use of an energy-based method, which seeks to produce a robust estimation of the average vibration level, is the approach adopted in this paper to treat the problem.

Energy methods are often used as alternatives in the medium and high frequency ranges whenever numerical methods like finite element or boundary equation-based formulations encounter significant limitations. Among such energy methods, the most widespread remains the Statistical Energy Analysis (SEA) [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF], which yields the mechanical energy of complex built-up structures. Many energy methods have been based on SEA, in an attempt to enhance the Statistical Energy Analysis robustness and predictability. Let's cite the work of Nefske and Sung [START_REF] Nefske | Power Flow Finite Element Analysis of Dynamic Systems: Basic Theory and Application to Beams[END_REF], who derived a local energy formalism. Many authors have improved on this model [START_REF] Wohlever | Mechanical Energy Flow Models of Rods and Beams[END_REF][START_REF] Lase | Energy Analysis of bars an beams: Theoretical Formulations[END_REF], leading to the Simplified Energy Method (MES), which has been applied to various domains, including beams, plates and acoustic cavities [START_REF] Lase | Energy Analysis of bars an beams: Theoretical Formulations[END_REF][START_REF] Besset | A coupled BEM and energy flow method for mid-high frequency internal acoustic[END_REF] as a "direct" prediction tool.

In this paper, the primary concern revolves around the use of MES in predicting energy sources [START_REF] Hardy | Absorption coefficient and energy flow path identification by means of inverse local energy method[END_REF]. This issue necessitates inverting the direct formulation. The paper will be structured as follows. To ensure clarity, the direct MES formulation will be reviewed in Section II, before providing the inverse MES formulation. A twodimensional case will then be considered and the algebraic matrix format of the inversion problem offered in Section III. In order to validate the proposed method, numerical simulations will be run in Section IV. A step-by-step strategy will be employed during the validation process. An initial step introduces direct estimated energies into the inversion sequence, while a second attempt simulates the input energy using a finite element code. A 2D plate will be studied, in addition to identifying surface sources. Ultimately, a comparison of the inverse MES results with an inverse SEA formulation will serve to complete the numerical experiments.

Overview of the direct energy flow method MES

Assumptions

The Simplified Energy Method (MES) is based on a description of two local energy quantities. The first is the total energy density W defined as the sum of the potential energy density and the kinetic energy density. The second energy variable I is the energy flow. The energy balance can then be written as follows:

∇ • I + π diss = π inp ( 1
)
where ∇ is the gradient operator, π diss is the dissipated power density and π inp the input power density. MES adopts the same damping model as SEA, as previously discussed in the literature [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF], which can be expressed as:

π diss = ηωW ( 2
)
where η is the loss factor and ω the circular frequency. The propagating waves considered in this MES formulation are made of partial energy quantities corresponding to both direct and reverberated fields. Since the considered fields are quadratic variables, the superposition principle can thus be applied:

£ = i £ i (3)
where (£, £ i ) are quadratic variables corresponding to either global energy or a partial energy associated with wave fields (W, W i ) or I, I i . In the following, an intrinsic energy law will be introduced, and this often enables defining the wave velocity c:

I i = c • W i n (4)
This expression simply defines the energy velocity, which in some cases is the same as the group velocity, as the ratio of energy flow to convected energy density. This expression was demonstrated in a number of textbooks [START_REF] Ohayon | Structural Acoustics and Vibration: Mechanical Models, Variational Formulations and Discretization[END_REF] and applies to any kind of linear elastic propagating wave.

Direct formulation

In the following, the derivation of the MES formulation is offered by considering a nDimensional (nD) space. The following equations are valid for n = 2 or n = 3 in the general case. In the next sections, we will consider the n = 2 case, whereas the n = 3 case will be treated in a subsequent article. In considering a symmetrical wave field and Equations ( 4) and ( 1), the MES approach [START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF] leads to the following relation:

1 r n-1 ∂ ∂r r n-1 I + ηωW n = 0 (5)
which can also be written using Equation ( 4) as:

-c 2 1 r n-1 ∂ ∂r r n-1 W n = ηω I (6)
Moreover, this equation can now be rewritten by considering just the W field, thus providing a local energy equation for symmetrical wave energies, i.e.:

1 r n-1 ∂ ∂r r n-1 W + ηω c W = 0 (7)
Energy variables are defined as the superposition of a direct field and a reverberated field, both of which can be summed to obtain the total field, as explained in Equation [START_REF] Ji | A study on an estimation method for applied force on the rod[END_REF]. A representation of these fields is given in Figure [START_REF] Liu | Input force estimation of a cantilever plate by using a system identification technique[END_REF], where only the intensity field has been shown. I dir is the direct field, while I rev is the reverberated field.

In Figure 1, Φ inp is the input power, Ω is the system surface and ∂Ω is the system boundaries. The elementary solutions in terms of energy density and active intensity are denoted G and H, respectively.

G(S, M ) = 1 SM e -ηω c SM H(S, M ) = c SM e -ηω c SM u SM (8)
where SM = SM . The kernels are determined at each point M of the considered space and created by a source located in S. By implementing a superposition principle, the energy solution can then be expressed using the contribution of the primary source ρ (direct field) and fictitious sources σ. The fictitious sources are unknown energy parameters representing the contribution of the reverberant field. More specifically, from [START_REF] Besset | A coupled BEM and energy flow method for mid-high frequency internal acoustic[END_REF]:

Φ inp Ω ∂Ω I rev I dir
W (M ) = Ω ρ(S)G(S, M )dS + ∂Ω σ(P ) u SM n P G(P, M )dP (9)
The active energy flow can be described by a similar expression:

I(M ) = Ω ρ(S)H(S, M )dS + ∂Ω σ(P ) u SM n P H(P, M )dP (10)
where ρ corresponds to internal sources in the considered domain (in this case, a two dimensional) and σ corresponds to sources at the boundaries of the two dimensional space.

The inverse energy flow method IMES formulation

This section will focus on introducing a discretized format of the energy integral equation, beginning with the definition of T as the (N, N ) geometric interaction matrix between N two-dimensional edge elements. This operator is given below in reference to Figure [START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF], with P and P ′ corresponding to P i and P j (point S only concerns equation 12.

T ij = ∂Ω ∂Ω 1 2 1 P i P j e -ηω c PiPj u PiPj n Pi dP i dP j (11)
In reference to this same Figure [START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF], another (N, k) matrix is defined and called herein Q. This operator represents the geometric interaction taking place between N boundary elements and k sources. Figure (2) was drawn with k set equal to 1. When considering the inverse problem, matrices are built assuming k sources, which will comprise the unknowns of this problem. For example, setting k = m will lead to a square problem to be solved. Hence, even if it were known that just one source had to be found, the matrices would still be built with k > 1, where k sources were located on Figure 2: element/element interaction the plate. The identified point should be the nearest to the source, and the other (k -1) sources should be identified close to 0. Matrix Q is thus expressed as follows:

Q ij = ∂Ω 1 2 1 S j P i e -ηω c SjPi u SjPi n Pi dP i (12)
Two matrices, labeled Y W and Y I (with dimension (k, m) where m denotes the number of measurements) can then be defined in order to represent a geometric interaction between sources and measurement points. These matrices are determined respectively for the energy density W and the energy flow I such that:

Y W ij = ∂Ω 1 S j P i e -ηω c Sj Pi dP i (13)
Y I ij = 1 c ∂Ω 1 S j P i e -ηω c Sj Pi u SjPi dP i (14)
Ultimately, A and B can be defined as follows (where Id being the identity matrix and ρ the injected power vector):

A = Id -T (15)
B = Qρ (16)
The power vector σ can thus be readily expressed as a function of the injected power vector ρ as:

σ = A -1 B = A -1 Qρ (17)
By combining Equations ( 9), ( 10) and ( 17), the energy density vector W and energy flow vector I can be obtained from:

W = S W ρ = R W ρ + Y W σ (18)
I = S I ρ = R I ρ + Y I σ ( 19
)
where R W and R I are defined as follows:

R W ij = 1 S i M j e --ηω c SiMj (20)
R I ij = 1 S i M j e --ηω c SiMj u SiMj (21)
S W and S I are then given by:

S W = R W + Y W A -1 Q (22) S I = R I + Y I A -1 Q (23)
Lastly, both the W and I fields can be expressed thanks to a linear operator denoted F as:

W, I = F (ρ) (24)
The IMES formulation seeks to invert the operator F proposed in Equation ( 24). Sources are to be detected by means of measurements conducted on the structure. Since the number of measurements is not equal to the number of sources, this problem is not directly invertible, which prompts us to seek an operator G that yields σ and ρ as a function of W and I.

ρ = G W, I (25)
Now, let's introduce the quantity R (X ) defined as follows:

R (X ) = W, I -F oX W, I (26)
where the operator o is defined as f og(x) = f (g(x)). We now need to minimize the quantity R (X ) over the operator X . In the following discussion, the IMES will be applied to a 2D plate. Consequently, Ω corresponds to the plate surface, whereas ∂Ω corresponds to the boundaries. In this case, only surface sources are being treated. Since the access to structural intensities is a difficult task to perform experimentally, the following simulations have been run by only considering energy density W measurements at certain locations of the structure. As a result, we will replace S W by S for the sake of clarity. MES equation ( 9) can then be discretized to provide the following matrix formulation:

         W 1 W 2 . . . W m          =     
              ρ 1 ρ 2 . . . ρ k          ( 27
)
where m is the number of measurements and k is the number of sources. To invert this formulation, according to Equation (26), a matrix X must be found that minimizes the quantity R(X):

R(X) =          W 1 W 2 . . . W m          -    
          X 11 X 12 . . . X 1m X 21 X 22 . . . X 2m . . . . . . . . . . . . X k1 X k2 . . . X km               W 1 W 2 . . . W m          (28) = (I m -SX)          W 1 W 2 . . . W m          ( 29
)
Matrix X is in fact chosen as a generalized inverse of matrix S:

X = S T S -1 S T (30)
The inverse formulation can then be performed. From a set of input energy parameters, the inversion can be achieved in order to characterize input sources. Section 4 below will describe the efficiency and limitations of the proposed strategy through numerical experimentation.

Numerical simulations and Results

This section will perform a number of numerical tests on a plate in order to validate the formulation. The first test case considered will be very simple; it consists of implementing the inverse energy flow approach when using a set of energy predictions based on the direct energy flow method. This set-up constitutes an MES/IMES simulation. A more realistic test case will be examined afterwards. In this situation, the input variables are computed from a full finite element model. A comparison between the IMES formulation provided in this paper and an inverse SEA model will also be given. Lastly, the effect of spatial discretization applied to the input set of energy quantities on the inversion process will be studied. 3.

Case number Direct method

• The first one will allow verifying the potential for method inversion. If the matrix given in equation 30 were indeed invertible, then the results should be perfect.

• The second case is intended to apply IMES in data originating from a finite element code. First of all, the structure is modeled using the finite element software COMSOL. Data are then extracted and treated through the IMES code, as explained in the previous section.

• The third case is similar to the second, yet an Inverse SEA software is used instead of the IMES code. The goal here is to compare the efficiency of IMES and Inverse SEA.

MES/IMES simulation

An initial simulation, referred as the MES/IMES simulation, involved the processing of IMES calculations using energy densities W derived from the direct MES formulation.

The model under study in this case is a thin plate, whose characteristics are described in Table (2). This plate is excited by an input power of 1W/m 2 . In assuming N to be the number of edge plate elements, the IMES result found for N = m = 64 is presented in Figure ( 4) (let's recall that m denotes the number of measurements).

It appears that the source has been correctly identified in this simulation (for both localization and level). The inversion of S T S has thus been validated. To illustrate the effect of the number of measurements m, we repeated the MES/IMES simulation while decreasing m. The IMES results are shown in figure [START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF]. Next, let's then introduce the factor ξ = m N to quantify method efficiency (ξ = 1 corresponds to a square matrix to be inverted).

It also appears that by increasing the number of measurements m spread over the plate surface corresponding to ξ = 0.06, ξ = 0.14 and ξ = 0. 3) has been lowered. Nevertheless, the source input power is equal to 0.7 W/m 2 for ξ = 0.06, 0.75 W/m 2 for ξ = 0.14 and 0.82 W/m 2 for ξ = 0.5 whereas it should equal 1 W/m 2 . The correct value of input power is obtained for higher values of ξ. The IMES tends to accurately localize the input structural force yet underestimates the amount of injected power.

FEM/IMES simulation

The IMES/MES subsection primarily allowed validating the S T S inversion. This subsection will consider the validation of an approach introducing FEM-simulated enξ Error 0.06 30 % 0.14 25 % 0.5 18 % Table 3: Error on the input power ergies. The IMES calculations were in fact processed using energy densities W stemming from a finite element software (i.e. Comsol in this study). An FEM model of the plate (whose characteristics are provided below) has been processed with 1W/m 2 excitation and the same number of elements N . Given the input power ρ, the force F applied to a point on the plate is recalculated by: ρ = 1 2

F 2 8 √ Dρ s h ( 31
)
Where D is the plate stiffness, defined by

D = Eh 3 12(1-ν 2 )
. The FEM calculation is then processed on the plate, thus providing the energy density W field used as an input for IMES algorithm. The FEM calculations are first spatially averaged and then frequency averaged. Octave bands are considered along with various center frequencies f c . The values of W output by the FEM calculation and performed around f c , are input into the IMES algorithm. Results found for ξ = 0.06, ξ = 0.14, ξ = 0.5 and ξ = 1 are all presented in figure [START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF].

The influence of f c has also been assessed. Results found for ξ = 1 and various values of f c are given in figure [START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF]. This figure reveals an accurate source localization regardless of the f c value within the medium to high frequency range.

Figure [START_REF] Noiseux | Measurement of power flow in uniform beams and plates[END_REF] displays the evolution of a localization term P r L with respect to the center frequency. P r L is defined by the following equation:

P r L = 1 N N i=1 (ρ i -ρ exact i ) 2 (ρ exact i ) 2 (32)
where ρ i is the i component of the input power vector and ρ exact i is the i component of the exact input power vector. It would appear that a good identification (both localization and level) is derived for the source with values of f c = 950Hz and higher. For lower-frequency bands, the proposed method will only be able to localize the sources.

SEA / MES comparison

The statistical Energy Analysis (SEA) is widely recognized as the original class of energy methods. SEA provides the mechanical energy of complex built-up structures. The goal of this subsection is to compare identification results obtained for FEM / SEA and FEM / IMES simulations. Towards this end, Figures [START_REF] Pavic | Measurement of structure borne wave intenisty: Part I, formulation of the methods[END_REF] and [START_REF] Hamdi | Hybrid inverse boundary element method for the determination of the optimal spectral characteristics of complex radiating noise sources[END_REF] indicate the identification tendencies for various damping values when using FEM/IMES and FEM/SEA identification strategies, respectively.

In all cases, the source seems to be correctly localized, though the estimated value of the input power ρ deviates from the expected value ρ exact with increasing η. To identify the localization sensitivity in figure (11(a)), we introduce the same localization term P r L given in Equation (32); moreover, representation of the quantification sensitivity in Figure (11(b)) has led us to define another quality term P r Q obtained from the following equation:

P r Q = (ρ source -ρ exact source ) 2 (ρ exact source ) 2 (33)
It appears that for both the localization and quantification of sources, the FEM/IMES results underestimates the FEM/SEA tendencies, which reflects that IMES that IMES results are more efficient than ISEA results in identifying sources for the studied plate. This finding would also suggest that the capacity of FEM/IMES localization remains more stable when varying η than FEM/ISEA localization results, for which this precision term increases when increasing η. Moreover, from both the FEM/IMES and FEM/ISEA quantification trends, it can be remarked that good results are obtained with weak damping values.

Effect of spatial discretization (sensors spread away from the source)

In this subsection, the IMES validation will be considered in a realistic case, according to which sensors are removed from the source. The IMES simulations were performed using a set of uniformly-distributed sensors far from the source location. Results relative to these FEM/IMES simulation runs and sensors distribution in the plate for m = 55 and m = 39 are given in figure [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF].

The results obtained indicate that the source has been correctly localized, yet a discrepancy remains between the estimated input power ρ is and the expected power ρ exact due to the lack of sensors close to the source. To analyze the localization sensitivity of this method with respect to parameter ξ defined below, Figure [START_REF] Nefske | Power Flow Finite Element Analysis of Dynamic Systems: Basic Theory and Application to Beams[END_REF] depicts the evolution in quality index P r L evolution. This figure provides a presentation of localization effectiveness.

It has thus been demonstrated that by using the IMES formulation, source localization can be achieved while limiting the number of sensors and distributing them far from the source.

Conclusion

The aim of this paper has been to investigate the effectiveness of the Simplified Energy Method (MES) in identifying sources from a set of given measured energies within a plate. The MES formulation and related assumptions were first detailed. The inverse formulation for identifying sources through data energy recovered from the plate was then provided. Afterwards, the MES could be validated using FEM / MES simulations by inputting numerical data from the FEM software. A comparison with SEA was also conducted and revealed that IMES is indeed more accurate. Lastly, a robustness study was performed, displaying the MES performance when sensors were removed from the source. Validation of this strategy for 3D acoustic cavities is now underway. In subsequent work, it will be worthwhile to test this method by considering two or more sources. From a theoretical perspective, the matrices and formulation proposed herein should be capable of treating this type of problem, yet numerical considerations complicate the identification step. Moreover, it would be beneficial to validate the method through experimental testing since the noise induced by "real" measurements should lead to findings regarding the inversion. From a numerical standpoint, it would be interesting to introduce a random noise on the input data in order to estimate the robustness of the inverse method proposed in this paper.

Figure 1 :

 1 Figure 1: Considered system

Figure 3 :

 3 Figure 3: Numerical simulation methodology

Figure 4 :

 4 Figure 4: MES/IMES simulation for the plate

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 5678910 Figure 5: MES/IMES simulations for m = 4(ξ = 0.06) (a), m = 9(ξ = 0.14) (b), m = 32(ξ = 0.5) (c) and m = 64(ξ = 1) (d)

Figure 11 :Figure 12 :

 1112 Figure 11: location (a) and quantification (b) sensitivities to the variation of damping for FEM / IMES and FEM / SEA

Figure 13 :

 13 Figure 13: location sensitivity to the variation of ξ

  S 11 S 12 . . . S 1k S 21 S 22 . . . S 2k

	. . .	. . .	. . .	. . .

S m1 S m2 . . . S mk



  S 11 S 12 . . . S 1k S 21 S 22 . . . S 2k

	. . .	. . .	. . .	. . .
	S m1 S m2 . . . S mk

Table 1 :

 1 Numerical simulations Table 1 summarizes the three test cases being considered herein. The corresponding methodology is presented in Figure

			Inverse method
	1	MES	IMES
	2	FEM	IMES
	3	FEM	Inverse SEA

List of Figures


S Besset

M M Ichchou

/ Ltds

Acoustic absorbtion material optimisation in the mid-high frequency range

Keywords: Absorption coefficient, Mid-frequency range, Energy methods

Acoustic modelling, applicable to poroelastic materials in a wide frequency range, is time consuming. This paper offers an energy method for optimising absorption coefficients at the boundaries of an acoustic cavity. The influence of absorption coefficients on this proposed energy method will be examined first. The next step will seek to optimise the location of absorbing materials. Numerical results will prove the method's effectiveness. Improvements in the overall optimisation process will also be proposed. A model based on polynomial interpolations will be developed in order to further reduce time consumption.

Introduction

Acoustical insulation is a major issue in engineering applications such as transportation (aircraft, automotive industry) or civil engineering. Noise mitigation is indeed required for the sake of public comfort. The level of noise can often be reduced through absorbent materials. Poroelastic materials such as foams are among the noise treatment packets. The weight and location of such additional materials however become key concerns. The optimisation process has received considerable attention as part of this noise abatement strategy. Predictive numerical modelling serves as the basis for determining optimisation quality. Noise prediction focuses on the widest frequency band of interest. The entire audible frequency range (up to 20 kHz) should in fact be taken into account by the numerical modelling approach. Classically speaking, pressure noise displays three different types of behaviour as frequency increases (see Figure 1. The first band defines the low-frequency region, in which finite element methods and/or boundary element methods are well adapted. The second and third frequency regions define the mid-and high-frequency domains. In the open literature, many optimisation processes are based on finite element or boundary element methods, and only very few focus on the mid-high frequency domain. The main goal of this paper is to address noise optimisation issues in the mid-high frequency range.

HF BF MF

Pressure noise

Frequency Figure 1: LF, MF and HF behaviours

The finite element method (FEM) is one of the most widely employed tools available for assessing the dynamic behaviour of such optimisation problems [START_REF] Silva | An Evolutionary Optimization Method applied to Absorbing Poroelastic Systems[END_REF][START_REF] Dutilleux | Low-frequency assessment of the in situ acoustic absorption of materials in rooms: an inverse problem approach using evolutionary optimization[END_REF]. Nevertheless, the excessive computational cost associated with solving largesize models, as required for example in the frequency response calculation, constitutes a major limitation with this method. Another limitation of the FEM is related to model reduction aspects when convergence of the "classical" reduced modal basis is not necessarily satisfied a priori in the short wavelength domain. More precisely, Silva et al. [START_REF] Silva | An Evolutionary Optimization Method applied to Absorbing Poroelastic Systems[END_REF] proposed a method to maximise the absorption performance of insulating poroelastic systems using a coupled finite element model and evolutionary strategies. In reference [START_REF] Dutilleux | Low-frequency assessment of the in situ acoustic absorption of materials in rooms: an inverse problem approach using evolutionary optimization[END_REF], the main intention was to propose a numerical model of a room that allows computing the pressure field within the volume, in combination with a global optimisation algorithm. Analytical methods can also be used in the optimisation process for some situations. For instance, Dupont and Galland [START_REF] Dupont | Active absorption to reduce the noise transmitted out of an enclosure[END_REF] applied analytical methods to optimise absorbing material locations in order to reduce the sound power transmitted by a plate. However, such an approach is limited to well defined and classical geometries.

Many authors have explored the topic of absorbing materials. Among them, let's cite M.A. Biot [START_REF] Biot | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range, ASA[END_REF][START_REF] Biot | Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher-Frequency Range[END_REF], K. Attenborough [START_REF] Attenborough | Acoustical characteristics of rigid fibrous absorbents and granular media[END_REF] and J.F. Allard [START_REF] Allard | Acoustic properties of partially reticulated foams with high and medium ow resistance[END_REF]. Our aim here is not to deal with such complex formulations, but instead to apply an energy method in order to optimise absorption coefficients on the boundaries of an acoustic cavity. In this context, a number of authors have already studied acoustical behaviour in damped cavities. Let's mention the work conducted for years in the area of prediction in cavities [START_REF] Dowell | Acoustoelasticity: General theory, acoustic natural modes and forced response to sinusoidal excitation, including comparisons with experiment[END_REF][START_REF] Trochidis | Sound transmission through double partitions with cavity absorption[END_REF][START_REF] Suzuki | Boundary element analysis of cavity noise problems with complicated boundary conditions[END_REF][START_REF] Sakagami | Sound absorption of a cavity-backed membrane: A step towards design method for membrane-type absorbers[END_REF][START_REF] Craggs | Coupling of finite element acoustic absorption models[END_REF]. More recently, authors have been focusing on simplified models [START_REF] Besset | Optimization of characteristics of porous materials based on a modal synthesis method[END_REF] and the optimisation of structures coupled with acoustic cavities [START_REF] Yamamoto | Thickness optimization of a multilayered structure on the coupling surface between a structure and an acoustic cavity[END_REF][START_REF] Yamamoto | Topology design of multi-material soundproof structures including poroelastic media to minimize sound pressure levels[END_REF]. Energy methods are expected to offer many advantages. Such methods are well adapted to medium and high-frequency situations and yield smaller matrices and fast optimisation processes. Moreover, energy methods are often used as alternatives in the high and medium frequency ranges. Among these methods, the most widespread remains the Statistical Energy Analysis (SEA) [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF], which provides the mechanical energy of complex built structures. We have chosen to use a local energy formalism, as first proposed by Nefske and Sung [START_REF] Nefske | Power flow finite element analysis of dynamic systems: Basic theory and application to beams[END_REF] and improved by many subsequent authors [START_REF] Lase | Energy analysis of bars and beams: Theoretical formulations[END_REF][START_REF] Sui | Prediction of vibroacoustics energy using a discretized transient local energy approach and comparison with TSEA[END_REF][START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF]8,[START_REF] Hardy | A hybrid local energy formulation for plates mid-frequency flexural vibrations[END_REF][START_REF] Ichchou | A transient local energy approach as an alternative to transient sea: Wave and telegraph equations[END_REF][START_REF] Fs Sui | Prediction of vibroacoustics energy using a discretized transient local energy approach and comparison with TSEA[END_REF][START_REF] Sui | Prediction of Time-Varying Vibroacoustic Energy Using a New Energy Approach[END_REF]. In the following discussion, this energy method will be called the Simplified Energy Method (or MES). MES has already been evaluated and validated for various elastic media such as beams [START_REF] Lase | Energy analysis of bars and beams: Theoretical formulations[END_REF][START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF][START_REF] Ichchou | Guided waves group and energy velocities via finite elements[END_REF], membranes and plates [START_REF] Lase | Energy analysis of bars and beams: Theoretical formulations[END_REF][START_REF] Hardy | A hybrid local energy formulation for plates mid-frequency flexural vibrations[END_REF], and acoustic cavities [START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF][START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF]. This method has also been considered in both the transient and stationary cases [START_REF] Sui | Prediction of vibroacoustics energy using a discretized transient local energy approach and comparison with TSEA[END_REF][START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF]8].

The present paper is organised as follows. The first section will describe the optimisation problem. A definition of the inverse problem in terms of inputoutput data will then be presented. The input parameters for this optimisation process will also be recalled. The second section will summarise the Simplified Energy Method (MES) formulation. The direct MES problem will be raised in terms of calculating energy densities from both input power and dissipated power at the boundaries. The sensitivity of energy density to acoustic absorption will then be discussed. This study will show that a 200% change in the absorption coefficient α leads to only a 60% change in the pressure field inside the cavity: the results on α therefore will not be very robust. The intention here is for the outcome of these optimisation issues to yield an "on-off" response1 . The lack of robustness in results is thus not considered to be that important. In order to improve time consumption statistics, a meta-model of energy distribution inside the cavity will be proposed. In this paper, the term "meta-model" refers to an analytical model of system behaviour. More specifically, it is a mathematical model along the lines of MES. This meta-model avoids the need to recalculate energy matrices during the optimisation loop. Lastly, a series of numerical experiments will be introduced and discussed. The optimisation routine proposed in this paper allows for fast calculations thanks to the reduced model size. Nevertheless, a number of limitations remain: the material is not being explicitly modelled; the porous medium is being modelled using an absorption coefficient; and damping inside the cavity has been ignored.

Simplified Energy Method (MES) 2.1 MES formulation

The Simplified Energy Method is based on a description of two local energy quantities. The first one is the total energy density W , defined as the sum of the potential energy density and the kinetic energy density. The second energy quantity I is the energy flow. The energy balance of such a system can be written as follows:

π diss = -∇ • I (1)
where ∇ is the gradient operator. π diss is the dissipated power density and may be neglected because the MES damping model, which is the same as SEA [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF], can be written π diss = ηωW with the damping loss factor η << 1 in acoustic fluids.

In this paper, we are particularly interested with the energy density W since it is proportional to the square pressure p 2 in the cavity, which is directly correlated with the sound level inside the cavity. Since W is a quadratic variable made of partial energy quantities corresponding to both direct and reverberated fields, the superposition principle can be applied:

W = W dir + W rev (2)
In the following, an intrinsic energy law will be introduced (this law is often used to define the wave velocity c):

I = c • W . n (3)
where n is the normal outgoing unit vector. Equation 3 offers a basic definition of energy velocity; it clearly states that both energy flow and energy density serve to define energy velocity. This expression is well known and has already been demonstrated for general elastic media. A proof is available for instance in [START_REF] Brillouin | Tensors in mechanics and elasticity[END_REF]. As an intrinsic law, this expression is adapted to all kinds of waves under eventual consideration.

The MES formulation is derived by considering a 3D space. Let's assume a symmetrical wave field along with Equations 3 and 1; the MES approach [START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF] leads to the following relation:

1 r 2 ∂ ∂r r 2 I = 0 (4)
which can also be written using equation 3, i.e.:

-c 1 r 2 ∂ ∂r r 2 W n = 0 (5)
Equation 4 stems directly from the energy balance, which is expressed in Equation 1. For unloaded lossless media (in this case, the acoustic domain) in the steady-state regime, Equation 1 under axisymmetric conditions simply leads to Equation 4.

The solution to Equation 5 corresponds to the radiosity method first proposed by Kuttruff [START_REF] Kuttruff | Energetic sound propagation in rooms[END_REF].

A representation of both direct and reverberated fields is given in Figure 2.

The elementary solution of energy density is denoted G: where γ 0 is the solid angle. The energy density inside the cavity can then be expressed as a function of the primary sources and fictitious sources located on the boundaries: [START_REF] Wang | Development of an energy boundary element formulation for computing high-frequency sound radiation from incoherent intensity boundary conditions[END_REF] where M is the integration point,

G(r) = 1 γ 0 c 1 r 2 (6) Φ inj Ω ∂Ω W rev W dir
W (P ) = ∂Ω Φ(M ) u P M . n(M )G(M ) d∂Ω + ∂Ω σ(M ) u P M . n(M )G(M )d∂Ω
u P M = --→ P M --→ P M
and n(M ) is the unit normal vector on ∂Ω at point M . Φ(P 0 ) is the acoustic boundary source and σ(M ) are fictitious boundary sources. We will use the term "boundary source" to denote those sources located on the cavity boundary (which may be due, for example, to external excitations).

The power balance on a point situated on a boundary allows expressing σ(P ) as follows:

σ(P ) = [1 -α(P )] ∂Ω,M =P σ(M ) u P M . nG(M )d∂Ω + Ω Φ(M )G(M )dΩ (8)
where α(P ) is the absorption coefficient at point P . By discretizing Equation 8, σ can be expressed as a function of Φ.

Application

In this paper, MES is being applied through a matrix formulation. The equations from Section 2.1 are expressed as follows:

σ = ([Id] -[α]) [T ] σ + ([Id] -[α]) [Q] Φ (9)
where [Id] is the identity matrix and [α] the diagonal matrix of the absorption coefficients. [T ] and [Q] are matrices corresponding to the discretization of the integral formulations of Equation 8. Let's note that mesh size does not need to be fine for energy methods like MES, which only involve quadratic variables. The meshing criteria used in FEM or BEM would thus not be useful herein. A coarse mesh size is sufficient and allows for a rapid optimisation process. Equation 9 directly leads to σ, as shown below:

σ = ([Id] -[T ] + [α][T ]) -1 ([Id] -[α]) [Q]. Φ (10)
W can then be expressed by discretizing equation 7, i.e.:

W = [R] Φ + [S] σ (11)
Using equation 10, W can be rewritten as follows:

W = [R] + ([Id] -[T ] + [S][α][T ]) -1 ([Id] -[α]) [Q] [M(α)] (n × m matrix) Φ (12)
It is now possible to express the objective function F ( α):

F (α) = Φ T [M (α)] T [M (α)] Φ (13)

Optimisation problem

Let's consider the 3D acoustic cavity shown in Figure 3. Acoustic sources may be applied on the boundaries. Point sources inside the cavity will not be examined since the context of this paper is tied to the aircraft and automotive industries.

Although the work programme outlined herein is quite academic and not directly related to industry, only those excitations caused externally will receive our attention. Such structures have already been studied by the authors in [START_REF] Besset | A coupled BEM and energy flow method for mid-high frequency internal acoustic[END_REF] and comparisons have been drawn between the MES method and both BEM and FEM methods, including a hybrid HMES method. The cavity structure has been discretised; for each boundary element, an absorption coefficient (which may depend on frequency) has been assigned. Once again, the aim of this paper is to optimise the distribution of absorption on elements in order to control acoustic energy within the acoustic domain.

The Simplified Energy Method (MES) has been described in depth in Section 2.1 above, and its implementation leads to a matrix formulation of the problem that can be written as follows: n represents the number of points where W will be calculated, while m is the number of facets originating from the discretisation of the structure. W = {W 1 , W 2 , . . . , W n } T is the energy density vector inside the cavity, correlated with the square pressure. This vector is computed at n points, which will be referred to as test points in the following discussion. Φ = {Φ 1 , Φ 2 , . . . , Φ m } T is the vector of input power on the boundaries, and α = {α 1 , α 2 , . . . , α m } T is the vector of absorption coefficients on the boundaries. W can thus be computed for n points inside the cavity. Let's now define the quantity F to be minimised:

W (α) = [M (α)] Φ (14)
F ( α) = W = n i=1 W 2 i ( 15
)
Function F is computed for a single frequency; a frequency averaging technique may be employed for instance whenever a more realistic definition of the impedance condition is introduced. Note that this function F is in fact identical to the L-2 norm of vector W . Optimisation of noise level within the cavity can then be performed as follows:

• choosing n points in the cavity where the sound level is to be minimised;

• discretising the structure and computing the matrix [M (α)];

• applying an optimisation algorithm to minimise objective function F , with α i (1 < i < m) being the variables.

The goal of this paper is not to develop new optimisation algorithms, hence classical routines will be employed. In the next sections, we will develop the MES principle and its application to optimisation problems.

Results

Effects of α on W

The aim of this section is to analyse the influence of α on the W field in the cavity. Figure 4 shows the variation in W as a function of the variation in α.

Optimisation of α distribution

In this section, the position of absorbing materials will be optimised. The data used for this optimisation algorithm are given below in Table 1. It should be noted that these optimisation steps are performed at a single test point; consequently, results will depend on this particular point and should differ if the noise needs to be reduced at another point.

It is also important to notice that absorption coefficients for absorbing materials are always frequency-dependent. In this paper, we have not considered the frequency dependency aspect, but it could be included without any difficulty by solving the problem for each frequency or frequency band in the considered domain. Moreover, the results of this optimisation algorithm will demonstrate that the optimal value for is either 0 or 1. Including frequency-dependent coefficients would therefore yield the same results. Nevertheless, the energy method considered in this paper is only valid in the context of mid-high frequency prob-lems. The optimisation problem is thus to be solved in a frequency range where modal overlap is quite high. The optimisation algorithm is performed using a local gradient method. The aim of this paper has not been to develop new optimisation methods; consequently, an optimisation protocol subject to constraints has been adopted using the fmincon function of the matlab software. This optimisation problem can then be formulated as follows:

Minimising F ( α) = Φ T [M (α)] T [M (α)] Φ so that      0 < α < 1 Nα i=1 α i N α = 0.3 ( 16
)
where N α is the number of absorbing patches. Let's recall that every patch covers an equal area so that the mean value of α over these two quantities is the same. Keep in mind that the constraints on α prevent obtaining a trivial result (i.e. α = 1 on all boundaries). In fact, these constraints are equivalent to setting the number of absorbing patches.

The result of this optimisation algorithm is given in Figures 5 and6. Let's recall that the source is located at (1.33, 0.6, 0.6), hence on the side not shown in the figures.

It is apparent that this algorithm leads to values of α equalto α min or α max . The optimisation method is thus efficient for the purpose of optimising the position of absorbing materials, yet does not prove useful in determining the exact values of α. Regardless, we have shown in Section 4.1 that this kind of optimisation would not have been sufficiently robust.

The algorithm used in this optimisation problem is a local gradient algorithm, meaning that it requires initial values for α (see Table 1). Results depend on these initial values: the solution we derived may be a local solution and moreover it may be possible to obtain other solutions if other initial values were to be considered. For this reason, the next section will examine method robustness.

Method robustness

In order to study the robustness of this method, several calculations have been carried out by considering various initial values.

From Figure 7, it can be observed that the final positions of absorbing materials remain independent of the initial values assigned to α. Moreover, certain initial values lie outside the range 0 < α < 0.6 and lead to similar results, which demonstrates method robustness. In conclusion, the method proves to be quite robust and suitable for use with any initial condition, α init . Moreover, Figure 8 shows that noise level in the cavity after optimisation does not depend on init (at least for the case treated herein). Note that in Figure 8 the level of noise reduction can be defined as follows:

N r = W init W optim × 100 (17)
A good distribution of absorbing patches can thus serve to significantly reduce noise inside the cavity (at least at the considered test point). A number of calculations have also been performed by considering two test points and several initial values (see Table 2). Results are shown in Figure 9: these results are also independent of the initial values of α.

Notice that the MES calculation protocol assumes lambertian directivity (ie. ). These results may differ when considering an alternative directivity. We have sought to use a uniform directivity d(θ) = 1 4π . Figure 10 displays the pertinent results; this figure has been obtained from a single test point and should be compared to results presented in Figure 7, which are significantly different. It would appear that knowing the directivity of the source is quite important.

1 st point 2 nd point x = 2.3 x = 1.2 y = 1.5 y = 1.5 z = 1.2 z = 1.2
α init = 0.2 α init = 0.3 α init = 0.4 α init = 0.5 α init = 0.7 α init = 0.9

Introduction of the meta-model for fast optimization

This section presents a meta-model designed for the absorption optimisation problem. Two reasons offer motivation for using meta-modelling techniques.

The first relates to the type of predicted fields. The MES formulation normally predicts smooth energy behaviour, which can be very effectively represented through polynomial interpolation. MES also contains a number of drawbacks, including the high computational cost of energy resulting from the kinds of ma- The idea being introduced herein is a novel one compared to the authors' previous output. Moreover, calculations using the gradient method are costly because the gradients have to be calculated using discretisation methods. For these reasons, a meta-model has been proposed in this paper.

Since the problem solution is quite robust, we also suggest performing a polynomial interpolation in order to quickly derive the solution to the problem.

In considering Equation 14 applied to a single test point, it becomes possible to write the matrices and vectors corresponding to the problem solution as follows:

W sol = W (α sol ) = [M (α sol )] Φ (18)
α sol denotes the vector corresponding to the optimal solution of this problem. The dimension of α sol is thus equal to m, i.e. equal to the number of facets stemming from discretisation of the structure. In this section, we have rewritten Equation 18 as follows:

α(x, y, z) = [Z(x, y, z)] ξ (19)
where ξ = x, y, z, x 2 , y 2 , z 2 , xy, xz, yz T and the dimension of matrix Z is (n, N ), with N dependent on the degree of polynomial interpolation (N = 10 in the present case). Matrix Z can then be built by:

• performing p optimisation algorithms → [α 1 sol ], [α 2 sol ], . . . , [α p sol ] ; 1 st point 2 nd point 3 rd point 4 th point x = 2.3 x = 1.3 x = 3.1 x = 0.7 y = 1.5 y = 0.5 y = 1.2 y = 1.2 z = 1.2 z = 1.7 z = 1.4 z = 0.4
Table 3: Coordinates of the 4 test points used in these calculations

• linking these solutions to the test point positions: α j sol = Z j ξ j , where ξ j = x j , y j , z j , (x j ) 2 , (y j ) 2 , (z j ) 2 , x j y j , x j z j , y j z j T , j ∈ [1, p] and the dimension of matrix Z j equals (n, N );

• interpolating the system in order to be able to consider any test point;

• changing the values obtained of α in order to satisfy the constraints (0 ≤ α ≤ 0.6). The results reported in section 4.2 indicate that the values of α were either α = 0 or α = 0.6, which leads us to apply the following changes to the results:

-α < 0.3 ⇒ α = 0 -α > 0.3 ⇒ α = 0.6
The results of this optimisation method are illustrated in Figure 11. The noise level has been set equal to 1 in the case of a non-optimised configuration (α = 0.3 everywhere on the boundaries). These calculations were performed for 4 different test points, which have been summarised in Table 3.

Figure 11 shows that the polynomial interpolation approach is not as efficient as full optimisation; nevertheless, it provides for good results at a very low calculation cost. Moreover, these results yield an initial result that is applicable to the initial values for full optimisation. Other polynomial interpolations have been conducted by considering 4th and 6th-order interpolations. Results are given in Figures 12 and13. Since these results show no improvement, the second order proves to be sufficient (in the special case of a box-like cavity, as studied in this work). The reason why orders above 2 do not improve results should be related to the fact that approximations are inappropriate due to the 0-1 nature of absorption coefficient values. Moreover, the exact form of the solution should not match any of the higher polynomial forms, and the "right" solution seems to be quite regular (i.e. the solution should be an approximation when considering a simple function).

Conclusion

In this paper, we have demonstrated the efficiency of optimising the position of absorbing materials on boundaries of a cavity. Results prove to be satisfactory even though the optimisation algorithm requires initial values; these results do not depend on initial values in the special case of the box-like cavity studied herein. It has been found however that optimisation results do depend on the source directivity, which must be evaluated before conducting the MES calculations. In order to perform a rapid optimisation, we have also proposed a polynomial interpolation, obtained using just 8 calculation points; these calculations have been run using only 8 fully-optimised points. Results fall short of the full optimisation results, yet the method allows reducing noise level in the cavity.

Figure 2 :

 2 Figure 2: MES formulation: direct and reverberated fields

Figure 3 :

 3 Figure 3: The acoustic finite element domain

Figure 4 :Figure 4

 44 Figure 4: Variation in W vs. variation in α Figure 4 indicates that large variations in α lead to minor variations relative to W . For example, a 100% variation in α leads to a 40% change in W , which therefore makes it impossible to obtain accurate results on the post-optimisation values of α. Nevertheless, optimisation can allow identifying positions where α should be increased.

Figure 5 :Figure 6 :

 56 Figure 5: Optimisation result: distribution of α on the boundaries (top view)white: α = 0, black: α = 0.6

9 Figure 7 :Figure 8 :

 978 Figure 7: Study of method robustness

Figure 9 :

 9 Figure 9: Method robustness: 2 test points

9 Figure 10 :

 910 Figure 10: Method robustness: Uniform directivity

Figure 11 :

 11 Figure 11: Noise reduction: comparison between full optimisation and polynomial interpolation: (1) initial values, (2) full optimisation, (3) fast optimisation

Figure 12 :

 12 Figure 12: Noise reduction: comparison between full optimisation and 4 th degree polynomial interpolation: (1) initial values, (2) full optimisation, (3) fast optimisation

Figure 13 :

 13 Figure 13: Noise reduction: comparison between full optimisation and 6 th degree polynomial interpolation: (1) initial values, (2) full optimisation, (3) fast optimisation

Table 1 :

 1 Input data used for optimisation routine

	n	1
	m	360
	α min	0
	α max	0.6
	α (mean)	0.3
	First step	α = 0.3
	Box dimensions	4 × 2 × 2 m
	Position of the source (2.83, 1.17, 0)
	Position of test point (1.33, 0.6, 0.6)

Table 2 :

 2 Coordinates of the 2 test points used in these calculations

	0.4	0.4
	0.2	0.2
	0.4	0.4
	0.2	0.2
	0.4	0.4
	0.2	0.2

"On" means that the absorption coefficient is set at its maximum value, while "off" means otherwise