Nicolas Hermant
email: nicolas.hermant@gipsa-lab.fr

Franz Chouly
email: franz.chouly@univ-fcomte.fr

Fabrice Silva
email: silva@lma.cnrs-mrs.fr

Paul Luizard
email: paul.luizard@gipsa-lab.fr

Numerical study of the vibrations of an elastic container filled with an inviscid fluid

Keywords: vibrations, fluid-structure interaction, finite element method, added-mass formulation, vocal folds replica

We investigate numerically the vibrational behavior of a simple finite element model that stands for an elastic container filled with an inviscid fluid. The underlying mathematical model is detailed and its spectra is characterized. The finite element method relies upon the added-mass formulation of Morand and Ohayon. A parametric study allows to characterize the system's response to dimensionless parameters in terms of eigenfrequencies. Then an insight into the mode shapes is provided, with a discussion on the presence and the behavior of singular modes caused by specific boundary conditions in the fluid.

Introduction

We consider the small vibrations of an elastic container filled with a liquid such as water. This situation appears frequently in engineering. It can be modelled conveniently in terms of fluid displacement and with a linearized Euler equation in the fluid, that avoids solving Navier-Stokes equations in a moving domain. Since the first pioneering works (see, e.g., [START_REF] Morand | Interactions fluides-structures[END_REF][START_REF] Schotté | Various modelling levels to represent internal liquid behaviour in the vibration analysis of complex structures[END_REF] and references therein) to derive models that correspond to such a situation, many numerical methods based on finite element discretization have been designed and analyzed.

First Kiefling & Feng [START_REF] Kiefling | Fluid-structure finite element vibrational analysis[END_REF] exhibited numerical instabilities, i.e., spurious modes occurring with almost zero eigenfrequency and showing circulation patterns of the fluid with no motion of the container. This has been confirmed by Hamdi et al. [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF]: they observe that the spurious modes can also appear away from zero, making it difficult to discriminate real modes from spurious ones. Bermúdez et al. [START_REF] Bermúdez | Finite element vibration analysis of fluid-solid systems without spurious modes[END_REF] provides a clear picture of the phenomenon, illustrating that the null eigenspace (with infinite dimension) contains exclusively pure rotational fluid motions. This justifies the variant with a penalty for the irrotationality of the fluid motion (see, e.g., [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF][START_REF] Olson | A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems[END_REF]), introducing the fluid vorticity as an additional variable. This method has two drawbacks: it increases the computational cost and spurious modes are not eliminated but pushed to higher frequencies. This could be overcome by an appropriate choice of the discretization as, e.g., the one suggested in Bermúdez & Rodriguez [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] where the fluid dispaclement variable is approximated with Raviart-Thomas finite element. This requires special consideration for coupling with the solid. For an elastoacoustic problem, Wang & Bathe [START_REF] Wang | Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems[END_REF] consider that the spurious modes can only be eliminated properly by means of a mixed displacement/pressure/vorticity formulation for the fluid. In Bermúdez et al. [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF], and following Morand & Ohayon [START_REF] Morand | Interactions fluides-structures[END_REF], an added-mass formulation is described and analyzed, as a simple way to get rid of the spurious modes. Extension to a fluid-solid problem involving a plate is proposed in [START_REF] Bermúdez | Finite element computation of sloshing modes in containers with elastic baffle plates[END_REF]. Barrientos et al. [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF] start from the same added-mass formulation, but discretize the fluid with boundary elements, an efficient technique that avoids to mesh the fluid domain.

In this paper we focus on a fluid model with a specificity: the liquid fills the container through a small hole in the boundary, that allows for pressure alimentation. This point is motivated by the physical study of an in-vitro vocal fold replica that is described in Sec. 6 as an example. Conversely, most of the previous works investigate either the problem of a fluid that fills completely a closed container [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF] or the problem of partial filling with a free surface where gravity waves and sloshing could damage the structure [START_REF] Berger | On a spectral problem in vibration mechanics: computation of elastic tanks partially filled with liquids[END_REF][START_REF] Morand | Interactions fluides-structures[END_REF]. In the first case, the whole fluid boundary is the interface with the elastic solid. The second case is usually modelled with an equation relating the pressure on the free surface to its vertical displacement (see, e.g., Boujot [START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF]). This case usually focuses on incompressible fluids, which induce a constraint accounting for the fluid volume conservation: the motion of the free surface has to compensate the normal displacement of the fluid-solid interface. At the opposite, in our work, we consider the case of a prescribed pressure on the «free» surface. This would relate to the usual free surface condition setting the gravity constant g to zero, but this has several important consequences. First there is now a Dirichlet condition on the portion of the fluid boundary that corresponds to the hole. This removes the need for the vertical displacement variable. Secondly, it releases the volume conservation constraint that is usually handled through a Lagrange multiplier [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF]. Finally, the Dirichet condition serves to provide a reference value for the fluid potential, that would otherwise require the introduction of a single static pressure value per fluid domain [START_REF] Olson | Eigenproblems from finite element analysis of fluid-structure interactions[END_REF].

Among the great amount of possibilities for the finite element discretization of such a problem, we selected the added-mass formulation introduced by Morand & Ohayon [START_REF] Morand | Interactions fluides-structures[END_REF] (see also, e.g., [START_REF] Felippa | Mixed variational formulation of finite element analysis of acoustoelastic/slosh fluid-structure interaction[END_REF][START_REF] Schotté | Various modelling levels to represent internal liquid behaviour in the vibration analysis of complex structures[END_REF]) and analyzed in Bermúdez et al. [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF] for a system with a free surface. This approach relies on the Helmholtz decomposition: the rotational motion being ignored, the fluid flow is potential, and the problem can be stated in terms of the solid displacement and the fluid potential. This has the advantage to replace the vectorial displacement variable by a scalar one and can not generate pure rotational flows. This may lead to non-symmetric mass and stiffness matrices for the solid displacement and fluid potential variables. This is circumvented by eliminating the fluid variable using an added-mass operator at the continuous level: the formulation then takes into account the fluid loading on the structure by means of a Neumannto-Dirichlet operator. At the algebraic level this means that we recover a generalized eigenvalue problem similar to what is obtained for a pure elastic problem: the contribution of the fluid is taken into account through an added-mass matrix that perturbs the solid mass matrix. The choice of the added-mass method is motivated by the following:

• It has a low computational cost, because the fluid variable is a scalar one, and the final eigenvalue problem to solve only involves the degrees of freedom associated to the solid displacement.

• It is simple to implement in a standard finite element library such as FreeFem++: it only makes use of Lagrange finite elements (for instance) and involves simple algebraic manipulations. No care has to be taken for the discretization of the fluid-solid interface (non-conforming discretizations are allowed). The symmetric eigenvalue problem that results can be solved with any standard algorithm.

• It has a firm mathematical background: the analysis detailed in [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF] ensures that no spurious modes can be computed, and optimal convergence to the exact eigenfrequencies and eigenmodes occurs when the size of finite elements gets smaller.

Our first contribution is to present a mathematical analysis of the added-mass formulation for the hydro-elastic spectral problem with prescribed pressure on part of the fluid boundary, that does not exist to the best of our knowledge. Then we present a detailed numerical parametric study, based on a dimensional analysis, that allows to understand the influence of each main physical parameter, on a configuration inspired from a real one: an in-vitro vocal fold replica. Finally we focus on eigenmodes and investigate the correlation between pure elastic and hydro-elastic modes, as well as the occurence of singular modes due to the specific boundary conditions in the fluid, associated to the pressure vent.

The paper is organized as follows: in Sec. 2 the mathematical model is presented as well as the underlying physical assumptions. Then, the added-mass formulation and the finite element discretization are detailed in Sec. 3. Sec. 4 is focused on the mathematical study of the yielding spectral problem. A parametric analysis, that allows to identify dimensionless parameters, follows in Sec. 5. The numerical results are described in Sec. 6, for an idealized vocal fold replica. Some conclusive remarks are stated in Sec. 7.

Problem statement

We consider a fluid enclosed in a deformable container (for an example made of latex, see Sec. 6). Notations are provided in Fig. 1. In what follows n denotes the unit normal to the interface Σ, outward to the solid domain Ω s . Modelling of each part as well as boundary and interface conditions are described below. While 3D brings no further problems (d = 3), a plane strain is assumed here for the sake of simplicity (d = 2). Within the vibrational analysis, small amplitude displacements in the vicinity of an equilibrium state are investigated.

Fluid Ω f Solid Ω s Σ Γ ext Γ l Γ m Γ m
Figure 1: Geometry of the coupled domains (with notations)

Solid domain: the elastic problem

The solid is assumed to be a linear elastic material with density ρ s . It is submitted to the fluid pressure fluctuation p on the interface Σ, and clamped on Γ s = ∂Ω s ∩ Γ m . All other volume forces or traction forces on Γ ext are set to zero, and harmonic solutions are sought (with convention e jωt ). This leads the following problem: find a pulsation ω ≥ 0 and a displacement u s :

Ω s → R d such that ρ s ω 2 u s + div σ(u s ) = 0 on Ω s , (1a)
u s = 0 on Γ s , (1b)
σ(u s )n = 0 on Γ ext , (1c)
σ(u s )n = p n on Σ. (1d)
Denoting σ(u s ) the Cauchy stress tensor, related to the small strain tensor ε(u s ) through a constitutive law σ(u s ) = C : ε(u s ). The elasticity tensor C is assumed to have usual ellipticity and uniform boundedness properties. In numerical experiments hereafter the Hooke law will be considered.

Coupling to the fluid domain: the hydro-elastic problem

The fluid is supposed to be irrotational, inviscid, homogeneous and incompressible, with density ρ f . We note u f the fluid displacement and p the pressure fluctuation around the rest state. Small harmonic solutions with angular frequency ω are introduced into the Euler equation. Following Morand and Ohayon [START_REF] Morand | Interactions fluides-structures[END_REF] a linearization procedure yields

-ρ f ω 2 u f + ∇p = 0 on Ω F , div u f = 0 on Ω F . (2)
The fluid motion being irrotational, there exists a fluid displacement potential ϕ such that u f = ∇ϕ and, according to Eq. (2), it is connected to the pressure through p = ρ f ω 2 ϕ. The mass conservation then implies that the potential obeys to the Laplace's equation within the fluid domain

-∆ϕ = 0 on Ω F . (3)
The fluid domain is bounded by the inflation inlet Γ l , where the pressure fluctuation is supposed to be zero, by the rigid support Γ m and by the interface Σ with the elastic solid. This implies for the fluid potential ϕ a Dirichlet condition on Γ l , a Neumann condition on Γ m , that accounts for sliding, and the transmission of the normal displacements and efforts on the interface Σ with the deformable container. The boundary condition on Γ l determines the reference value for the potential that would otherwise be defined up to an additive constant. As a result, the hydro-elastic vibrations problem can be formulated as: find ω ≥ 0, u s : Ω s → R d and ϕ :

Ω f → R such that ρ s ω 2 u s + div σ(u s ) = 0 on Ω s , (4a)
u s = 0 on Γ s , (4b)
σ(u s )n = 0 on Γ ext , (4c) σ(u s )n = ρ F ω 2 ϕ n on Σ. ( 4d
) -∆ϕ = 0 on Ω f , (4e)
ϕ = 0 on Γ l , (4f)
∇ ϕ • n = 0 on Γ m , (4g)
∇ ϕ • n = u s • n on Σ. ( 4h
)

Weak formulation and finite element method

In this section, following, e.g. [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF][START_REF] Morand | Interactions fluides-structures[END_REF], we introduce a weak formulation of the hydro-elastic problem described in Sec. 2 as well as a finite element discretization. Both are based on continuous and discrete added-mass operators that allow reduction to a simpler problem in which the only unknown is the solid displacement.

The set of infinitely derivable functions with compact support in a domain D is denoted C ∞ c (D). The Lebesgue space of square integrable functions on a domain D is denoted L 2 (D). Moreover we use the standard notation for Sobolev spaces (H s (D)) d , s ∈ R, d = 1, 2, 3, in one, two or three dimensions (see, e.g., [START_REF] Adams | Sobolev spaces[END_REF]). The Sobolev norm is denoted • s,D (dual norm if s < 0) and we keep the same notation whatever is the value of d.

Weak formulation of the hydro-elastic problem

We introduce first V s and V f , that are the spaces of kinematically admissible fields for the solid displacement and the fluid potential, respectively:

V s = v ∈ (H 1 (Ω s )) d : v = 0 on Γ s , V f = ψ ∈ H 1 (Ω f ) : ψ = 0 on Γ l .

Then Eq. (4a)-(4h) yield

Find ω ≥ 0, u s ∈ V s and ϕ ∈ V f such that

       K s (u s , v) -ω 2 M s (u s , v) = ω 2 C f s (ϕ, v) ∀v ∈ V s , M f (ϕ, ψ) = C f s (ψ, u s ) ∀ψ ∈ V f , (5)
with the solid stiffness and mass bilinear forms K s and M s defined, respectively, as

K s (u s , v) = Ω s σ(u s ) : ε(v) dΩ s and M s (u s , v) = Ω s ρ s u s • v dΩ s . (6)
The fluid mass and the coupling bilinear forms M f and C f s are defined, respectively, as

M f (ϕ, ψ) = Ω f ρ f ∇ ϕ • ∇ ψ dΩ f and C f s (ψ, v) = Σ ρ f ψ(v • n) dΣ. (7)
At this stage, we introduce the Neumann-to-Dirichlet operator for the Laplace equation on the fluid domain. For this purpose we define the trace space on the interface Σ, of functions in V f :

W f = {ψ| Σ : ψ ∈ V f }.
The space W f is endowed with the norm

ξ W f = inf ψ∈V f ,ψ| Σ =ξ ψ 1,Ω f ,
for ξ ∈ W f . The topological dual of W f will be noted W f (L 2 (Σ) being the pivot space), and the duality product between W f and W f will be noted • , • . The Neumann-to-Dirichlet operator is then defined as follows:

M : W f → W f η → ϕ η | Σ ,
where ϕ η ∈ V f is the unique potential field solution to

Ω f ρ f ∇ ϕ η • ∇ ψ dΩ f = ρ f η, ψ| Σ ∀ψ ∈ V f . (8)
It maps the prescribed Neumann condition η on Σ to the trace of the resulting potential ϕ η on this same interface. Note that it is sometimes referred to the Poincaré-Steklov operator [START_REF] Holmes | A spectral problem in hydroelasticity[END_REF][START_REF] Langer | Coupled finite and boundary element domain decomposition methods[END_REF] and is similar to the input impedance for unidimensionnal acoustic waveguides. Remark that u s • n ∈ L 2 (Σ) ⊂ W f , so from Eq. ( 5) and ( 8) there holds:

ϕ| Σ = M(u s • n),
for ϕ solution to Eq. ( 5). Therefore, for any v ∈ V s there holds:

C f s (ϕ, v) = Σ ρ f ϕ| Σ (v • n) dΣ = Σ ρ f M(u s • n)(v • n) dΣ.
This motivates the introduction of the added-mass bilinear form:

M a : V s × V s → R (u s , v) → Σ ρ f M(u s • n)(v • n) dΣ.
It describes how the fluid waves resulting from the interface motion influence back the vibration of the structure, The added-mass formulation of the spectral hydro-elastic problem now writes:

Find ω ≥ 0 and u s ∈ V s such that K s (u s , v) = ω 2 (M s (u s , v) + M a (u s , v)) ∀v ∈ V s . (9)
Note that at the continuous level formulations ( 5) and ( 9) are equivalent.

Finite element discretization

We discretize Problem (5) following the approach detailed in [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF][START_REF] Morand | Interactions fluides-structures[END_REF]. For this purpose, let V h s ⊂ V s (resp. V h f ⊂ V f ) be a conformal finite element discretization for the solid displacement u s (resp. for the fluid potential ϕ). In practice we choose Lagrange finite elements, that are continuous piecewise polynomials of degree k = 1, 2 for V h s and of degree 1 for V h f , built upon triangulations of each domain (h denotes the mesh size, see, e.g., [START_REF] Ciarlet | Basic error estimates for elliptic problems[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]). We denote by U s the column vector that contains all the degrees of freedom (hereafter abbreviated in dofs) associated to a discrete solid displacement u h s ∈ V h s , and similarly Φ the column vector that contains all dofs from a discrete fluid potential ϕ h ∈ V h f . We note respectively K s , M s , M f and C f s the matrices obtained from the bilinear forms K s , M s , M f and C f s after the standard finite element assembly procedure (see, e.g., [START_REF] Ern | Theory and practice of finite elements[END_REF]). The resulting discrete problem associated to Problem (5) then reads, in matricial form:

       K s U s = ω 2 M s U s + C f s Φ , M f Φ = C T f s U s . (10)
As noted in, e.g. [START_REF] Morand | Interactions fluides-structures[END_REF], the positive measure of Γ l implies that the fluid mass matrix M f is invertible. This allows to eliminate Φ in the first equation of Eq. ( 10) with help of the second equation, and we get:

K s U s = ω 2 (M s + M a ) U s with M a = C f s M -1 f C T f s . (11)
The added-mass matrix M a is a discrete counterpart of the added-mass operator M a . Note that no geometric conformity of the element nodes or edges at the fluid/solid interface Σ is required by this formulation, which simplifies slightly the meshing process. We recover a purely elastic problem from Eq. ( 11) if we simply omit the added-mass matrix (this corresponds to the case ρ f = 0 at the continuous level).

The numerical results are obtained with FreeFem++ finite element library [START_REF] Hecht | New development in FreeFem++[END_REF][START_REF] Hecht | FreeFem++ -Third Edition[END_REF] that provides a convenient open source environment for the specification of the weak formulation, its transformation into the discrete addedmass problem and numerical solving. Problem [START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF] is solved in FreeFem++ through an Implicit Restarted Lanczos Method implemented in the package ARPACK 1 . Tolerance for this iterative method has been fixed to 10 -10 . The functional interface of ARPACK is used: K s and (M s + M a ) are implemented as linear operators that apply to a vector (see [START_REF] Hecht | FreeFem++ -Third Edition[END_REF]Section 9.4]). This avoids the explicit computation of M -1 f to assemble the added-mass matrix M a . This is cheaper in terms of computational cost (memory and time), especially in that case where the inversion of the matrix M f implies the loss of sparsity of the added-mass matrix M a . However, as the coupling matrix C f s only involves dofs of U s that are located on the interface Σ, this density only corresponds to a reduced subset of dofs.

Regarding Dirichlet boundary conditions associated to the elastic solid, these are imposed on the stiffness matrix K s , but not on the mass matrix (M s + M a ). This enables to push to the highest frequencies the spurious modes arising from the treatment of Dirichlet boundary conditions with exact penalty (see [START_REF] Hecht | FreeFem++ -Third Edition[END_REF]Note 9.1]).

Characterization of the spectrum

The aim of this section is to characterize the spectrum of the hydro-elastic problem [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF], and to show it has a structure very similar to a pure elastic problem (i.e., uncoupled with a fluid): a countable increasing sequence of real positive eigenvalues of finite multiplicity. This is due to the simple model chosen to describe the fluid, that acts only as an added-mass. Note that for different fluid-solid coupled problems, with more general fluid models, such as incompressible Navier-Stokes, this characteristic may not be preserved anylonger (see, e.g., [START_REF] Fernández | Linear stability analysis in fluid-structure interaction with transpiration. I. Formulation and mathematical analysis[END_REF]). Mathematical characterization of the spectrum has already been carried out for very similar hydro-elastic problems, in, e.g., [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF][START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF][START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF][START_REF] Ryan | Eigenvalue and eigenfunction error estimates for finite element formulations of linear hydroelasticity[END_REF]. In the aforementioned works, the main difference in their model comes from the boundary condition on Γ l . In [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF][START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF] a free surface condition is chosen, whereas in [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF][START_REF] Ryan | Eigenvalue and eigenfunction error estimates for finite element formulations of linear hydroelasticity[END_REF] a completely enclosed fluid cavity is considered (Γ l = Γ m = ∅). As a result, we adapt the framework of Boujot [START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF] for the spectral characterization of our specific model. Before addressing the hydro-elastic problem, we recall the main result for the pure elastic problem.

Elastic problem

For the purely elastic problem the result below is well-known and can be found in, e.g., [START_REF] Babuška | Eigenvalue problems[END_REF].

Theorem 1. Let us consider the elastic problem derived from Eq. (5) by setting ρ f = 0. This spectral problem has a countable infinite sequence of real eigenvalues of finite multiplicity

0 < λ E 1 ≤ λ E 2 ≤ . . . → ∞, with the relationship ω E n = λ E n , n ≥ 1. Furthermore the corresponding eigenfunctions (u E n ) n≥1 satisfy the orthogo- nality properties K s (u E n , u E m ) = (ω E n ) 2 M s (u E n , u E m ) = (ω E n ) 2 δ nm , (12)
where δ nm denotes the Kronecker symbol.

Note that 0 is discarded as an eigenvalue because of the Dirichlet boundary condition (1b) that prevents rigid body motions.

Hydro-elastic problem

To characterize the spectrum of Problem (5), we need first a regularity assumption. Indeed to apply a compacity argument in the proof of Theorem 2, we need some extra regularity associated to Problem [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF]. Since in our situation, the Laplace equation is associated to mixed Dirichlet/Neumann boundary conditions, standard regularity results can not apply and the solution φ η is not expected to belong to H 2 (Ω f ) (see, e.g., [START_REF] Grisvard | Problèmes aux limites dans les polygones. Mode d'emploi[END_REF] and [START_REF] Ern | Theory and practice of finite elements[END_REF]Theorem 3.10,Remark 3.11]). Nevertheless we can make a weaker regularity assumption as follows:

Assumption 1. Let us suppose that there exists δ > 0 such that, given η ∈ L 2 (Σ), the solution φ η to Problem [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] belongs to H 1+δ (Ω f ), with the regularity estimate

φ η 1+δ,Ω f ≤ C η 0,Σ , where C > 0.
Assumption1 is valid for instance if Ω f is a simply connected domain in R 2 with polygonal boundary and interior angles less than 2π (no crack), and in that situation δ > 1 4 (see, e.g., [START_REF] Babuška | On the mixed finite element method with Lagrange multipliers[END_REF][START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF][START_REF] Grisvard | Problèmes aux limites dans les polygones. Mode d'emploi[END_REF]). Moreover this is true for a flat boundary (ω = π), as in Fig. 1, and in this case δ = 1 2 (see, e.g., [START_REF] Grisvard | Problèmes aux limites dans les polygones. Mode d'emploi[END_REF]). Moreover we need a compactness theorem in fractional-order Sobolev spaces, that is direct application of, e.g., [35, Lemma 10] (see also [START_REF] Di Nezza | Hitchhiker's guide to the fractional Sobolev spaces[END_REF]Theorem 7.1] for a more general result):

Lemma 1. Let Ω ⊂ R d , d ≥ 1, be a Lipschitz bounded open set, then, for every δ ∈ (0, 1), the injection of H 1+δ (Ω) into H 1 (Ω) is compact. Proof: Let us pick a bounded sequence (v n ) in H 1+δ (Ω). It means that (v n ) is bounded in H δ (Ω) and (∂ x i v n ) is bounded as well in H δ (Ω)
, where ∂ x i v n denotes the weak partial derivative of v n , respectively to direction x i , i = 1, . . . , d. We apply two times the Lemma 10 of [START_REF] Palatucci | Local and global minimizers for a variational energy involving a fractional norm[END_REF] (compactness of the injection H δ (Ω) ⊂ L 2 (Ω)) and we get a subsequence, still denoted

(v n ) such that v n → v in L 2 (Ω), ∂ x i v n → w i in L 2 (Ω), i = 1, . . . , d.
Furthermore, for any φ ∈ C ∞ c (Ω) there holds

Ω w i φ dΩ = - Ω v ∂ x i φ dΩ,
(this is obtained for instance using the writing w i -∂ x i v n + ∂ x i v n in the left hand side, and passing to the limit when n → +∞). Whence we obtain that v ∈ H 1 (Ω) and

w i = ∂ x i v. It results that v n converges to v in H 1 (Ω).
Then we state our main result:

Theorem 2. Let us consider that Assumption 1 is valid. Then Problem (5) has a countable infinite sequence of real eigenvalues of finite multiplicity

0 < λ H 1 ≤ λ H 2 ≤ . . . → ∞, with the relationship ω H n = λ H n , n ≥ 1. Moreover the corresponding eigenfunctions (u H n ) n≥1 satisfy the orthogonality properties K s (u H n , u H m ) = (ω H n ) 2 (M s (u H n , u H m ) + M a (u H n , u H m )) = (ω H n ) 2 δ nm , (13)
where δ nm denotes the Kronecker symbol.

Remark 1. We differentiate the eigenfunctions for the hydro-elastic problem using the notation (u H n ) n≥1 , whereas (u E n ) n≥1 corresponds to the eigenfunctions of the pure elastic problem. The relationship and correlation between these two families of eigenfunctions is investigated numerically in Sec. 6.3.

Proof:

The proof is carried out in two steps: 1) first, we reformulate Problem (5) into an equivalent problem, 2) then we check that the framework provided in [11, Sections 1.1-1.2] can be applied (we use slightly different notations).

1) We introduce the spaces

V = (H 1 (Ω s )) d × H 1 (Ω f ), H = (L 2 (Ω s )) d × H 1 (Ω f ),
endowed with the product norm as well as

V s f = y = (v, ϕ v ) ∈ V s × V f : M f (ϕ v , ψ) = C f s (ψ, v) ∀ψ ∈ V f ⊂ V. ( 14
)
Note that V ⊂ H with continuous injection. We introduce as well the bilinear form, for x = (u, ϕ) ∈ V and

y = (v, ψ) ∈ V: a(x, y) = K s (u, v).
Similarly we define, for x = (u, ϕ) ∈ H and y = (v, ψ) ∈ H:

b(x, y) = M s (u, v) + M f (ϕ, ψ).
Problem ( 5) can be recasted as:

Find ω ≥ 0 and x = (u s , ϕ u s ) ∈ V s f such that a(x, y) = ω 2 b(x, y) ∀y = (v, ψ) ∈ V s f . (15)
Note that ( 5) and ( 15) are equivalent (see [START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF]Section 2.2] for more details in the context of a similar problem).

2) Now let us remark that V s f is closed in V and we introduce H s f that is the closure of V s f in H for the norm of H. Note as well that a( • , • ) is a bilinear symmetric continuous form on V s f and b( • , • ) is a bilinear symmetric continuous form on H s f . For y = (v, ψ) ∈ H s f , there holds:

b(y, y) = ρ s v 2 0,Ω s + ρ f ∇ψ 2 0,Ω f ≥ ρ s v 2 0,Ω s + c ρ f ψ 2 1,Ω f , (16)
with c > 0. Here we applied the Poincaré-Friedrichs inequality (see, e.g., [START_REF] Ern | Theory and practice of finite elements[END_REF]Lemma B.66]) since any function ψ in V f satisfies a homogeneous Dirichlet boundary condition ψ = 0 on Γ l . It results that b( • , • ) is coercive on H s f . From Korn's inequality (see, e.g., [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF]) and the above property ( 16), there holds, for every λ > 0, and for y = (v, ψ) ∈ V s f :

(a + λb)(y, y) ≥ c( v 2 1,Ω s + ψ 2 1
,Ω f ), where c > 0 is a constant that depends on λ, ρ S , ρ F , on the ellipticity constant of the elasticity tensor C and on the Korn's and Poincaré-Friedrichs' constants. Thus (a + λb)( • , • ) is coercive on V s f for every λ > 0.

The last point to satisfy is the compact injection of V s f into H s f . For this purpose let us pick a bounded sequence

(y n ) = (v n , ϕ v n ) in V s f , n ≥ 1.
This implies that the sequence (v n ) is bounded in (H 1 (Ω s )) d : we can apply Rellich-Kondrachov Theorem (see, e.g., [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]Theorem 9.16]) and extract a subsequence, still denoted (v n ), that converges in (L 2 (Ω s )) d , with an associated subsequence in V f , still denoted (ϕ v n ). The elements ϕ v n ∈ V f are solution to:

Ω f ρ f ∇ ϕ v n • ∇ ψ dΩ f = ρ f v • n, ψ| Σ ∀ψ ∈ V f ,
that is Problem [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF] with boundary data η = v • n ∈ L 2 (Σ). We apply the Assumption 1 so that, for a given δ > 0, each ϕ v n is regular enough to belong to H 1+δ (Ω f ). In addition the regularity estimate combined to a trace inequality yields

ϕ v n 1+δ,Ω f ≤ C v n • n 0,Σ ≤ C v n 1,Ω s , which means that the sequence (ϕ v n ) is bounded in H 1+δ (Ω)
. Whence, applying Lemma 1 there is a subsequence, still noted (ϕ v n ) that converges into H 1 (Ω f ). To summarize we proved that there exists a subsequence (

y n ) = (v n , ϕ v n ) that converges in H s f .
We are ready now to apply the results of Ref. [START_REF] Boujot | Mathematical formulation of fluid-structure interaction problems[END_REF]Section 1.2]: there exists a sequence of numbers λ H m ≥ 0 growing to infinity as m → ∞, and a family of elements x m = (u m , ϕ m ) that are orthogonal in V s f for a( • , • ) and that form an orthonormal basis of H s f for b( • , • ). These are such that:

a(x m , y) = λ H m b(x m , y), ∀y ∈ V s f . Note that in fact λ H
1 > 0 due to the Dirichlet boundary condition (1b) that prevents rigid body motions. The pulsations ω H m are retrieved from the eigenvalues λ H m as ω H m = λ H m .

Parametric analysis

We aim at identifying the dimensionless parameters associated to Problem [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF]. Among other purposes, this analysis is achieved to highlight the size proportion effects due to the small thickness of the elastic solid in contrast with the fluid that occupies a larger domain. Such an analysis might help in the future to choose or design an appropriate reduced order structural model such as a plate or a shell.

As illustrated in Fig. 2, the first step consists in considering a solid membrane -of thickness e -interfaced with a squared-shape fluid domain -of size R × R (in 2D). At the bottom of the fluid domain there is a vent/opening -of width l -where pressure is imposed. Then, the following transformations are considered:

• in the solid domain: the reduced spatial coordinates x = x/R and y = y/e are introduced, mapping

Ω s = [0, R] × [0, e] to Ω s = [0, 1] × [0, 1]
. Moreover adimensional displacement fields u and v link to u and v via the scaling parameter U: u = Uu and v = Uv . Cartesian coordinates of u (resp. v ) are denoted by (u x , u y ) (resp. (v x , v y )).

• in the fluid domain:

Ω f = [0, R] × [-R, 0] maps to Ω f = [0, 1] × [-1
, 0] with x = x/R and y = y/R. The parameter Φ, homogeneous to a displacement potential, defines the adimensional fluid potential ϕ = ϕ/Φ and ψ = ψ/Φ. The bilinear operators defined in Sect. 3.2 can be reformulated within the new domains:

x y Γ ext Σ e -R R Ω f Ω s Γ l Γ m x y Γ ext Σ 1 -1 1 Ω s Ω f Γ l Γ m
• the solid stiffness operator K s = Ω s λ tr ε(u) tr ε(v) + 2µε(u) : ε(v) (with λ, µ the Lamé coefficients of the Hooke law) splits into three dimensionless operators:

K s (u, v) = U 2 ER e K ee + e R K eR + e 2 R 2 K RR with                              K ee (u, v) = Ω s 1 -ν (1 + ν)(1 -2ν) ∂ y u y ∂ y v y + 1 2(1 + ν) ∂ y u x ∂ y v x dΩ s , K eR (u, v) = Ω s ν (1 + ν)(1 -2ν) ∂ x u x ∂ y v y + ∂ y u y ∂ x v x + 1 2(1 + ν) ∂ y u x ∂ y v x + ∂ x u y ∂ x v y dΩ s , K RR (u, v) = Ω s 1 -ν (1 + ν)(1 -2ν) ∂ x u x ∂ x v x + 1 2(1 + ν) ∂ x u y ∂ x v y dΩ s .
with E the Young's modulus and ν the Poisson's ratio. The first term relates to the flexural deformation, the second to shear, and the third one to inplane membrane strain.

• the solid mass operator:

M s (u, v) = Ω s ρ s u • v dΩ s = ρ s eRU 2 M s (u , v ) with M s (u , v ) = Ω s u • v dΩ s .
• the fluid mass operator:

M f (ϕ, ψ) = Ω f ρ f ∇ ϕ • ∇ ψ dΩ f = ρ f Φ 2 M f (ϕ , ψ ) with M f (ϕ , ψ ) = Ω f ∇ ϕ • ∇ ψ dΩ f .
• the coupling operator:

C f s (u, ψ) = Σ ρ f ψ(u • n) dΣ = ρ f RUΦ C f s (ψ , u ) with C f s (ψ , u ) = Σ ψ (u • n)dΣ .
Gathering all these expressions, Problem (5) can be recasted as:

E ρ S e 2 K ee + e R K eR + e 2 R 2 K RR (u , v ) = ω 2 M s (u , v ) + ρ f ρ s Φ eU C f s (ϕ , v ) , ΦM f (ϕ , ψ ) = UR C f s (u , ψ ),
or, equivalently, using the added-mass operatorM a built from C f s and M f , that is the dimensionless counterpart of

M a : E ρ S e 2 K ee + e R K eR + e 2 R 2 K RR (u , v ) = ω 2 M s + ρ f ρ s R e M a (u , v ). (17)
In Eq. ( 17) four characteristic parameters can be identified:

1. the membrane aspect ratio ζ = R/e, 2. the ratio of fluid to solid densities η = ρ f /ρ s , 3. the Poisson's ratio ν that appears in the stiffness integrals K ee , K eR and K RR , 4. a scaling frequency Ω = E/(ρ s e 2 ), showing some similarity with the classic eigenfrequency √ K/M of a single d.o.f. mass-spring system (considering K = E/eS and M = ρ s eS for some surface area S of a membrane with thickness e).

The elastodynamics problem is recovered as the limit case ρ f → 0, i.e. η → 0. A last parameter can not be ignored: the measure of the surface Γ l is implicit in the definition of the kinematically admissible fluid fields, and can be quantized by the dimensionless vent width = /R. It is essential in the sense that it describes the relative measures of the fluid boundaries where Neumann and Dirichlet conditions apply, respectively. Besides, it plays a role similar to the measure of the free surface area for sloshing problems. In Sec. 6.2 we will present numerical results of the parametric study for a model of vocal fold replica, using these five characteristic parameters.

Example: a vocal fold replica

Experimental vocal fold replicas are currently used in speech production studies [START_REF] Cullen | Brass instruments: linear stability analysis and experiments with an artificial mouth[END_REF][START_REF] Gilbert | Artificial buzzing lips and brass instruments: Experimental results[END_REF][START_REF] Pelorson | An experimental replica of the vocal folds to study normal and pathological voice[END_REF][START_REF] Ruty | An in vitro setup to test the relevance and the accuracy of low-order vocal folds models[END_REF]. A typical design of these replica, as shown in Fig. 3, is made of a container -0.2 -0.5 mm thick latex -filled with a fluid -usually water. The residual stress is minimized through avoiding the stretching and the twisting of the latex membrane when it is mounted on its metallic support The mechanical frequency responses of the structure can be adjusted varying the internal pressure of the water filling the half-cylinder volume. The water volume inside the replica is connected to a wide-open water tank so that pressure fluctuations within the replica might not affect the elevation of the free surface in the tank. Pressure within the replica is determined thanks to the water height in the tank. This pressure control is not under the scope of the current paper, but explains the difference with other hydro-elastic formulations accounting for a free surface variable. Furthermore, the interaction between the replica and the glottal air flow is not investigated here.

On one hand, in-situ mechanical characterization of these replica and their comparison with numerical simulations are critical in order to validate and improve the physical and mathematical models. On the other hand, numerical simulations may reveal complex vibrational behaviours that are challenging to detect experimentally. The default configuration correspond to a 2D approximation of the vocal fold replica depicted in Fig. 3 and considers the numerical parameters values listed in Tab. 1. piecewise quadratic functions (P 2 elements, 5050 dofs) leads to errors on eigenfrequencies comparable to the ones obtained with a mesh of N = 32 layers for the piecewise linear finite element P 1 space (265518 dofs). For a given precision, choosing the P 2 finite element space thus reduces the memory requirement and the CPU time in the eigenproblem solving step. It is also known that P 1 elements suffer from Poisson locking near the incompressible limit, see e.g. Ref. [START_REF] Babuška | Locking effects in the finite element approximation of elasticity problems[END_REF]. We do know that plate or shell elements can overcome the increase of the number of solid dofs, but the current work is intended to be used as a reference for model reduction. Moreover, those models too are prone to shear and/or membrane locking that may overestimate the eigenfrequencies.

latex metal

No exact solution being known, the values are fitted on a power-law function to obtain an approximated solution and the order of convergence. We observe an almost second-order convergence when discretizing the displacement field using piecewise quadratic functions.

In the case of the hydro-elastic problem, the triangulations in the fluid and solid domains can be nonmatching accross the interface Σ. Within the convergence study, we chose a coarser mesh of Ω f with the size of the elements still varying with N. Table 3 illustrates the lowest hydroelastic eigenfrequencies using P 2 elements for the solid displacement u and P 1 elements for the fluid potential ϕ. The added-mass effect is evident, the fluid loading strongly lowers the eigenfrequencies. The same procedure is applied on the results to estimate the orders of convergence. The latter appears to be smaller than in the elastic problem and matches with the results reported by Bermudez et al. [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF].

Parametric study

The parametric study of the model is performed using the results of the analysis presented in Sec. 5. Fig. 5 shows the evolution of the first ten eigenfrequencies when independently varying the first four characteristic parameters previously identified in Sec. 5 (the influence of the width of the pressure inlet is object of a separate section, see Sec. 6.4).

Figs. 5(a-b) illustrates the influence of the solid stiffness on the eigenfrequencies for the elastic and the hydroelastic problems. As expected from Eq. ( 17), there is a linear dependency on the scaling frequency Ω and then a (1/2) power relation with the Young modulus E: ω ∝ √ E. Changing the container material while maintaining its density and its Poisson ratio would trivially scale its modal characteristics: systems with a stiffer material have higher natural frequencies, as one would intuitively expect. The fluid mass merely contributes with a uniform decrease of the frequencies, independently from the Young modulus. On Figs. 5(c-d) are shown the evolutions of the eigenfrequencies with the ratio of fluid to solid densities η = ρ f /ρ s . Obviously, the elastic problem does not consider the fluid so there is no variation and we recover the asymptotic behavior for very light fluids (η < 0.01). This is consistent with the theoretical formulation written in Sec. [START_REF] Barrientos | Analysis of a coupled BEM/FEM eigensolver for the hydroelastic vibrations problem[END_REF] Table 3: Computed eigenfrequencies (in Hz) of the first modes for the hydroelastic problem, when the number N of element layers accross the thickness of the membrane is varied. Piecewise quadratic continuous elements for u and piecewise linear continuous elements for ϕ.

the added-mass term becomes negligible if the fluid density is much weaker than the solid one. However, there is a heavy fluid regime (η > 1) where the eigenfrequencies decrease towards zero with (-1/2) power law relation with respect to η. In that case, the added-mass term becomes predominant in comparison to the solid mass term in Eq. ( 17) and ω is almost inversely proportional with √ η or, more precisely, with √ ρ f solely. Notice that the reference configuration lies in the heavy fluid domain, where one could almost neglect the solid mass term. Finally, in this figure, a curve does not conform to the common pattern and shows a heavy-fluid behavior even for smaller values of η. We will discuss this curve below, in Section 6.4. The effect of the aspect ratio ζ = R/e of the latex membrane is depicted in Figs. 5(e-f). It is possible to distinguish two regimes here too: first, for a thin membrane (ζ > 10), the eigenfrequencies are decreasing functions of the aspect ratio, with almost (-1) and (-3/2) power dependence for most of the solutions of the elastic and the hydro-elastic problems, respectively. This is not explicit from Eq. ( 17), but is consistent with the results obtained with the Reissner-Mindlin plate model (see, e.g., Ref. [START_REF] Durán | Approximation of the vibration modes of a plate by Reissner-Mindlin equations[END_REF]) where the frequencies of the bending modes of a thin plate (without fluid loading) scale with e. The (-3/2) dependency then comes from the heavy-fluid case, where the added-mass term dominates the solid mass term and has an extra e/R coefficient. It means that the eigenfrequencies have a higher sensitivity to the thickness in the presence of fluid than without it. On the contrary, for a thick membrane (ζ < 10) or even a "capillary" fluid vessel with thick wall (ζ < 1, not shown on the figures), the frequencies show an asymptotic independence with respect to the latex thickness. The reference configuration of the vocal fold replica lies in the thin membrane case. As in Fig. 5(f), this family of branches crosses with some others branches that exhibit different sensitivities with respect to the aspect ratio ζ. In the hydro-elastic case, the first of them has for example an almost (-1/2) power dependence, but does not appear in the elastic case, and may so be fluid-specific. On the contrary higher frequency solutions appear in both the elastic and the hydro-elastic problems. They may even be increasing functions of the aspect ratio.

Finally, Figs. 5(g-h) illustrate the evolution of the eigenfrequencies (normalized to the reference configuration)

along with the Poisson's ratio. It evidences that almost all the solutions scale in a common function of this coefficient, even if the resemblance is greater in the elastic case than in the hydro-elastic one. In the latter, there is at least one frequency ( f 4 ) with a variation that differs from the others.

Discussion on the modes shapes

Particular solutions branches were evidenced in the previous Figs. 5(d,f,h). It is important to investigate the associated mode shapes to understand their nature and to explain their different dependencies on the various dimensionless coefficients.

We first focus on the reference configuration, and compare the dofs vector sets (U E n ) n≥1 and (U H m ) m≥1 obtained in the elastic and the hydro-elastic problems, respectively. Table 4 provides values of the eigenfrequencies and the corresponding kinetic energies for the first hydro-elastic modes: the solid kinetic energy M s (v j , v j ), the fluid kinetic energy M a (v j , v j ) and the ratio of fluid to solid kinetic energies M a (v j , v j )/M s (v j , v j ).

(HE0) (HE1) (HE2) (HE3) (HE4) (HE5) (HE6) (HE7) (HE8) (HE9) Frequency (Hz) 5.7 A classic correlation descriptor used in the experimental modal analysis is the Modal Assurance Criterion (MAC) (see, e.g., Ref. [START_REF] Morales | Comments on the MAC and the NCO, and a linear modal correlation coefficient[END_REF]) that consists in a matrix defined as:

MAC n,m = U E n • U H m 2 U E n 2 U H m 2 ( 18
)
and that is robust to the convention of normalization used for the eigenvectors. Values close to 1 indicate that the dofs vectors U E n and U H m are almost colinear. On the contrary, values close to 0 implies that the vectors are not correlated. However, this criterion is not an orthogonality check in the sense of the displacement fields. This property is more conveniently obtained in the frame of the current numerical study, through the use of the solid mass bilinear operator or, equivalently, of the solid mass matrix (see, e.g., Ref. [START_REF] Morales | Comments on the MAC and the NCO, and a linear modal correlation coefficient[END_REF]). This is usually known as the normalized cross orthogonality (NCO):

NCO n,m = M s (u E n , u H m ) M s (u E n , u E n )M s (u H m , u H m ) = M s (u E n , u H m ) M s (u H m , u H m ) = (U E n ) T M s U H m (U H m ) T M s U H m ( 19
)
as the elastic modes are normalized with respect to the solid mass matrix M s . Fig. 6 compares the first ten dofs vector sets of the elastic (E) and hydro-elastic (HE) models. at first glance, there is an almost one-to-one correspondence of the solid modal displacements, at the exception of the hydro-elastic mode 4 (HE4) that will be discussed hereafter. The correspondence is not perfect: due to the fluid loading, the (HE) modes, namely mode (HE2) and higher, even if mainly correlated with one elastic mode, show as well small contributions of lower elastic modes. For instance mode (HE5) has strong correlation with mode (E4) and some small but not negligible similarity with modes (E0) and (E2). Symmetry considerations lead to the chess-like pattern visible in the lower part of the NCO matrix. The fluid loading thus couples the mode shapes that are orthogonal in the elastic case. Fig. 7 explicits the singularity of the mode (HE4). While the other hydro-elastic modes consider the vibrations of the fluid-loaded membrane, the mode (HE4) has kinetic energy that is essentially localized at the pressure inlet Γ d , i.e. with little participation of the membrane. The vibration of the membrane appears to be, in this mode, a mere side-effect of the main fluid motion: the ratio of the maximal fluid to solid displacements is about 15, while it is about 1 for the other modes, hereafter denoted as coupled modes. The ratio of fluid to solid kinetic energies is relevant too as it is bounded to 18 for all hydro-elastic modes but the mode (HE4) that exhibits a four times higher value (almost 83). The mode (HE4) thus ressembles the sloshing modes that appear in free surface problems accounting to pressure fluctuations related to the elevation of the surface (see, e.g., Refs. [START_REF] Morand | Interactions fluides-structures[END_REF][START_REF] Schotté | Various modelling levels to represent internal liquid behaviour in the vibration analysis of complex structures[END_REF]).

Another point of interest of Fig. 7 lies in the proportion of the fluid region involved in the vibration. While almost all the full fluid moves in the modes (HE0) to (HE2), the participation then seems to decrease for higher modes to a thin layer in the vicinity of the interface. This is also visible in Table 4. As the hydro-elastic solutions are normalized with respect to the total mass operator M s + M a (see Eq. ( 13)), there holds, for the denominator in Eq. ( 19):

M s (u H m , u H m ) = 1 -M a (u H m , u H m ) ≤ 1. (20)
For higher hydro-elastic modes, the fluid kinetic energy decreases towards zero due to the reduced dimension of the fluid vibration region. This increases the denominator in Eq.( 19), i.e. the solid kinetic energy, and thus lowers the values of the NCO coefficients as also visible in Fig. 6. Finally, another discrimination criterion is the order of convergence of the estimation of the singular mode eigenfrequency. As visible in Table 3, the mode (HE4) has a lower order than the modes (HE0) to (HE3). When considering gravity waves on the free surface, Ref. [START_REF] Bermúdez | A finite element solution of an added mass formulation for coupled fluid-solid vibrations[END_REF] evidences as well different behaviors for the singular modes and the coupled modes (called sloshing and hydroelastic modes, respectively).

Influence of the width of the pressure inlet

As exposed at the end of Sec. 5, the definition of the kinematically admissible fluid fields relies on a last parameter: the width of the pressure inlet that does not explicitly appear in the dimensionless problem. In order to emphasize its influence, the evolution of the first eigenmodes is computed for values of this parameter ranging from 0.1 mm to 8 mm.

From Fig. 8(a), the width seems to noticeably affect only the singular mode (black dotted curve). In fact, the representation of the modal shape ensures that the modes labelled (HE4) and (HE5) swap for /R ∼ 0.3 without any kind of coupling. The frequency of the singular mode increases along with the width of the inlet. This is consistent with the increase of the eigenfrequencies with the size of the free surface in sloshing problems [START_REF] Morand | Interactions fluides-structures[END_REF].

The singular mode can be discriminated from the first hydro-elastic eigenmodes by computing the ratio of the fluid and solid kinetic energies that takes higher values for the former than for the latter (see Fig. 8(b)). This ratio decreases for higher hydro-elastic coupled modes as the fluid domain involved in the vibration concentrates in the vicinity of the membrane for high-frequency vibrations. Another criterion for the discrimination bases on the kinetic energy flux at the pressure inlet defined as

Φ v j , v j = Γ l ρ f ∇ϕ v j • n • n 2 . ( 21
)
This compares with the potential energy of the free surface in sloshing problems. The energy flux is several order of magnitude higher for the singular mode than for the coupled modes, as visible in Fig. 8(c).

Conclusion

In this paper, a mathematical analysis of the added-mass formulation is proposed for the hydro-elastic vibration problem with prescribed pressure on a part of the fluid boundary. In addition a parametric analysis and a numerical study enhance the different possible configurations depending on the importance of the fluid loading on the membrane container. Furthermore the focus has been put on the distinction between coupled modes and a singular mode that has its kinetic energy localized in the vicinity of the fluid boundary where the Dirichlet condition applies.

This work is a first step in the study of the behavior of a vocal-fold replica as described at the beginning of Sec. 6. A first limitation concerns the 2D simulations, and future work will consist in the design and study of a 3D finite element model. In fact, the 2D case relies on the hypothesis of a plane strain, and could then realistically model the physical situation where the membrane is free at the extremities of the hemi-cylindrical structure. The 2D modal results correspond only to the first family of modes. In the 3D case, the latex should be constrained on these extremities with homogeneous Dirichlet boundary conditions to conform to the experimental situation. This should lead to a stiffer situation for the first modes, and the 3D eigenfrequencies should be then higher than the 2D ones.

Furthermore the current formulation is not able to explain the dependance of the eigenfrequencies on the static pressure prescribed at the inlet which is the primary experimental control of the replica. Thus the next steps concern the inflation of the solid container due to the finite static pressure and the vibration problem to be formulated on the inflated membrane, i.e. on a prestressed configuration. Preliminary results are given in Ref. [START_REF] Hermant | Hydro-elastic finite element model of a vocal fold replica[END_REF] and the present study still shows some relevance when considering a simple hyperelastic constitutive law such as the Saint-Venant-Kirchhoff law and assuming that the fluid is at rest in the inflated configuration.

Figure 2 :

 2 Figure 2: Adimensionalization of Problem (5).

Figure 3 :Figure 4 :

 34 Figure 3: Design of the vocal folds replica used at GIPSA-lab [37]. Left: metal support with water supply hole covered with the latex cylindric membrane. Right: sketch of the final assembly and of the water pressure control by the height of the free surface in the tank.

Figure 5 :

 5 Figure 5: Influence of the four characteristic quantities on the first ten eigenfrequencies ( f 0 to f 9 ) for the elastic (left column) and the hydro-elastic (right column) problems. The reference configuration is denoted by the vertical dotted lines.

Figure 6 :

 6 Figure 6: Representation of the NCO matrix for the first ten modes of the elastic and the hydro-elastic problems.

Figure 7 :

 7 Figure 7: First nine hydro-elastic modes in the reference configuration: solid and fluid displacements vector fields are represented with the same scale. Kinetic energy localization is also displayed (increasing values from yellow to red).

Figure 8 :

 8 Figure 8: Top: (a) evolution of the first eigenfrequencies for a varying width of the pressure vent and corresponding (b) ratios of fluid to solid kinetic energies and (c) energy flux at the pressure inlet. Bottom: Modes (HE4) and (HE5).

Table 2 :

 2 Computed eigenfrequencies (in Hz) of the first modes for the elastic problem, when the number N of element layers accross the thickness of the membrane is varied.

	, where

Table 4 :

 4 Frequency and energy indicators for the hydro-elastic solutions of the reference configuration.

		14.1	29.2	49.2	70.0	77.5	112.4 155.1 205.9 265.6
	Solid kinetic en-ergy M s (v j , v j ) Fluid kinetic en-ergy M a (v j , v j ) Ratio of fluid to gies solid kinetic ener-	0.053 0.067 0.086 0.103 0.012 0.124 0.143 0.163 0.183 0.203 0.947 0.933 0.914 0.897 0.988 0.876 0.857 0.837 0.817 0.797 17.9 13.9 10.7 8.7 82.8 7.1 6.0 5.1 4.5 3.9

http://www.caam.rice.edu/software/ARPACK/

Acknowledgments

This work was partially supported by grant VoFoCam (ANR-12-PDOC-0018) of the ANR french funding agency. The authors thank Xavier Pelorson for some helpful discussion and comments, and as well for his support.

Geometric parameters

Convergence study

We first investigate the convergence of the eigenfrequencies of the elastic problem when the mesh size is reduced. An example of the structured mesh used for the solid domain is shown in Fig. 4. The number N of layers accross the thickness of the membrane and the number of subdivisions of the boundaries Σ and Γ ext are varied accordingly to preserve an aspect ratio of the triangles close to 1.

Table 2 shows the first lowest frequencies f n for different successively refined meshes for the reference configuration, corresponding to parameters provided in Table 1. It is noticeable that a mesh with N = 2 layers accounting for


Tan Vu

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Luc Dugard
email: luc.dugard@gipsa-lab.grenoble-inp.fr

Peter Gaspar
email: gaspar@sztaki.mta.hu

Van Tan Vu

Optimal selection of weighting functions by genetic algorithms to design H∞ Anti

Keywords: H∞ control, Genetic Algorithms, Heavy vehicle, Active anti-roll bar control, Rollover, Roll stability

come   L'archive ouverte pluridisciplinaire

INTRODUCTION 1.1 Context

Rollover of heavy vehicle is an important road safety problem world-wide. Although rollovers are relatively rare events, they are usually deadly accidents when they occur. Moreover, the roll stability loss is the main cause of traffic accidents in which heavy vehicles are involved. In order to improve the roll stability, several schemes with possible active intervention into the vehicle dynamics were proposed. One of them employs active anti-roll bars, that is, a pair of hydraulic actuators which generates a stabilizing moment to counter balance the overturning moment [START_REF] Vu | Active anti-roll bar control using electronic servo-valve hydraulic damper on single unit heavy vehicle[END_REF].

On the other hand, the H∞ control design approach is an efficient tool for improving the performance of a closed-loop system in pre-defined frequency ranges. The key step of the H∞ control design is the selection of weighting functions. In many applications, the difficulty in choosing these functions still increases as performance specification is not accurately defined i.e., it is simply to achieve the best possible performance (optimal design) or to achieve an optimally joint improvement of more than one objective (multi-objectives design). So the weighting functions optimization to satisfy the desired performances is still an open problem. Recently, Genetic Algorithms were used in [START_REF] Alfaro-Cid | Optimisation of the weighting functions of an H∞ controller using genetic algorithms and structured genetic algorithms[END_REF] and [START_REF] Kitsios | H∞ controller design for a distillation column using genetic algorithms[END_REF], their formulation well suited for this type of problematic [START_REF] Do | Optimization of weighting function selection for H∞ control of semi-active suspensions[END_REF].

Related works

Some of the control methods applied to active anti-roll bar control on heavy vehicle are briefly recalled below: a-Optimal control: Sampson et al [START_REF] Sampson | An investigation of roll control system design for articulated heavy vehicles[END_REF], [START_REF] Sampson | Achievable roll stability of heavy road vehicles[END_REF], [START_REF] Sampson | Active roll control of single unit heavy road vehicles[END_REF] proposed a state feedback controller which was designed by finding an optimal controller based on a linear quadratic regulator (LQR) for single unit and articulated heavy vehicles.

The LQR was also applied to an integrated model including an electronic servovalve hydraulic damper model and a yaw-roll model of a single unit heavy vehicle. The input current of the electronic servo-valve is the input control signal [START_REF] Vu | Active anti-roll bar control using electronic servo-valve hydraulic damper on single unit heavy vehicle[END_REF].

b-Neural network control: A reinforcement learning algorithm using neural networks was used to improve the roll stability for a single unit heavy vehicle [START_REF] Boada | Active roll control using reinforcement learning for a single unit heavy vehicle[END_REF].

c-Robust control (LPV): Gaspar et al [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF], [START_REF] Gaspar | Reconfigurable control structure to prevent the rollover of heavy vehicles[END_REF], [START_REF] Gaspar | Prediction based combined control to prevent the rollover of heavy vehicles[END_REF] applied Linear Parameter Varying techniques to control active anti-roll bars combined with active brakes on a single unit heavy vehicle. The forward velocity was the varying parameter.

Paper contribution

Based on the H∞ active anti-roll bar control presented by the authors in [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF], this paper proposes the use of Genetic Algorithms to define automatically the weighting functions. Hence the following contributions are brought:

-The Genetic algorithms method is applied to define the weighting functions of the H∞ robust controller for active anti-roll bar system on the single unit heavy vehicle. Thanks to GAs, the conflicting objectives between the normalized load transfers and the generated torques are handled using a single high level parameter only, denoted α.

-The simulation results in frequency and time domains show the normalized load transfers at two axles when the tuning parameter value α moves from 0 to 1, compared with the case of passive anti-roll bar and the results of [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF]. In time domain, a cornering manoeuvre is used for the heavy vehicle. The forward velocity is considered up to 160 km/h to evaluate the roll stability and to determine the maximal velocity at which the normalized load transfers and the generated torques reach their limits.

The paper is organised as follows: Section 2 gives a brief introduction about multiobjective optimization using GAs. Section 3 presents the model of a single unit heavy vehicle. Section 4 develops the H∞ control synthesis to prevent rollover of heavy vehicles. Section 5 illustrates how to use the GAs to define the weighting functions of the H∞ active anti-roll bar. Section 6 presents some simulation results in frequency and time domains. Finally, some conclusions are drawn in section 7.

GENETIC ALGORITHMS AND MULTI-CRITERION OPTIMIZATION

Genetic algorithms

GAs are now widely used since the first study in [START_REF] Holland | Adaptation in natural and artificial systems, an introductory analysis with application to biology, control and artificial intelligence[END_REF], confirmed by a popular theory-oriented book [START_REF] Goldberg | Genetic Algorithms in Searching Optimisation and Machine Learning[END_REF] and an application-oriented book [START_REF] Davis | Handbook of genetic algorithms[END_REF]. The algorithms are based on the natural selection mechanism and have been proven to be very effective in optimization in many real applications such as finance and investment strategies, robotics, engineering design, telecommunications, etc. They are likely global optimization techniques (despite the high computational expense) using probabilistic, multi-points search, random combination (crossover, mutation) and information of previous iteration to evaluate and improve the population. A great advantage of GAs compared with other searching methods (for example gradient methods) is that they search regardless of the nature of the objective functions and constraints.

GAs initialize with a random population which evolves through genetic operations: selection, crossover and mutation. By using a selection process, the fittest individuals based on their fitness values are chosen; crossover and mutation are then applied to create the new population. The genetic operation on individuals continues until the optimization criterion is satisfied or a certain number of generations is reached.

Fitness function: The fitness of an individual is used to choose between "good" and "bad" individuals. An individual with a high fitness has a great chance to be selected.

Selection: This step is to sort and copy individuals by order of satisfaction of the fitness function. The higher the value of the fitness, associated to an individual, the greater the individual's chances to be selected to participate in the next generation. "Proportionate" [START_REF] Holland | Adaptation in natural and artificial systems, an introductory analysis with application to biology, control and artificial intelligence[END_REF] and "tournament" [START_REF] Miller | Genetic algorithms, tournament selection, and the effects of noise[END_REF] selections are the most popular methods.

Crossover: This main operation acting on the population of parents is an exchange of parts of chains between two selected individuals (parents) to form two new individuals (children). This exchange may be due either to a single or multiple points.

Mutation: Mutation operates on a single individual by randomly changing a part of it. In the binary coding case, it is done by reversing one or more bits in a chromosome.

Multi-criterion optimization

One well-known application of the GAs is to find the optimal solution for the multiobjective optimization problem involving multiple and conflicting objectives. This is a very popular problem in practice and can be described as follows:

1 2 ( ) ( ) min ( ) , . .2,
( )

obj obj x C n f x f x F x n f x ∈       =         M ( 1
)
where x is called the decision vector, C the set of possible decision vectors (or the searching space), and F(x) the objective vector.

The existence of an ideal solution x* that can simultaneously minimize all objective functions f1, f2, .., fnobj is in fact rarely feasible.

There are many formulations to solve the problem (1) like weighted min-max method, weighted global criterion method, goal programming methods... [START_REF] Do | Optimization of weighting function selection for H∞ control of semi-active suspensions[END_REF] and references therein. One of the most popular and simple approaches is the weighted sum method which converts the multi-objective problem into a single objective one. In this paper, one uses a particular case of the weighted sum method, where the multiobjective functions vector F is replaced by the convex combination of objectives:

1 1 min ( ), . , 1 obj obj n n i i i i i J f x s t x C α α = = = ∈ = ∑ ∑ (2)
The vector

SINGLE UNIT HEAVY VEHICLE MODEL

Fig 1 illustrates the combined yaw-roll dynamics of the vehicle modelled by a threebody system, where ms is the sprung mass, muf the unsprung mass at the front including the front wheels and axle, and mur the unsprung mass at the rear with the rear wheels and axle. The model variables are given in Table 1 and the parameters in [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF].

( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) s yf yr xz zz yf f yr r xx s xz s s f uf f uf ARf f r ur r ur ARr r yf uf uf uf uf uf tf uf f uf mv m h F F I I F l F l I m h I m gh m vh k b M U k b M U rF m v r h m gh k k β ψ φ φ ψ φ ψ φ β ψ φ φ φ φ φ φ φ φ β ψ φ φ φ φ + - = + - + = - + - = + + - - - - + + - - - - + + - = - + + - + -                ( ) ( )( ) ( ) ( ) f uf
b M U rF m v r h m gh k k b M U φ φ β ψ φ φ φ φ φ φ            + - + +   - = - + - - + - + - + +         (3.1) (3.2) (3.3) (3.4) (3.5)
The lateral tire forces Fyf and Fyr in the direction of velocity at the wheel ground contact points are modelled by a linear stiffness as:

yf f f yr r r F C F C µ α µ α =    =   (4)
where the tyre side slip angles are given as:

f f f r r l v l v ψ α β δ ψ α β   = -+ -    = -+     (5)
The moment of passive anti-roll bar impacts the unsprung and sprung masses at the front and rear axles as follows [START_REF] Vu | Active anti-roll bar control using electronic servo-valve hydraulic damper on single unit heavy vehicle[END_REF], [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF]:

A B A ARf AOf AOf uf A B A ARr AOr AOr ur t t t M k k c c t t t M k k c c φ φ φ φ  = -     = -   (6)
where kAOf , kAOr are respectively the torsional stiffness of the anti-roll bar at the front and rear axles, tA half the distance of the two suspensions, tB half the distance of the chassis and c the length of the anti-roll bars' arm.

Using the previous equation, the single unit heavy vehicle can be represented by the linear system in the state space form (7):

1 2

x Ax B w B u y Cx

  = + +  =    (7)
with the state vector:

uf ur x β ψ ϕ ϕ ϕ ϕ   =      
, the disturbance input:

f w δ   =   the control inputs: f r u U U   = 
 and the output vector:

uf ur y β ψ ϕ ϕ ϕ ϕ   =      
.

H∞ CONTROL SYNTHESIS OF ACTIVE ANTI-ROLL BAR ON HEAVY VEHICLES

Control objective, problem statement

The objective of the active anti-roll bar control system is to maximize the roll stability of the vehicle. Usually, an imminent rollover is detected when the calculated normalized load transfer reaches 1 (or -1), as explained hereafter. First, the lateral load transfer can be given by:

u u z w k F l ϕ ∆ = ( 8
)
where ku is the stiffness of tire, φu the roll angle of the unsprung mass and lw the half of vehicle's width. Then, the lateral load transfer can be normalized w.r.t. the total axle load Fz as follows:

z z F R F ∆ = (9)
The normalized load transfer R=±1 corresponds to the largest possible load transfer. In that case, the inner wheel in the bend lifts off.

While attempting to minimize the load transfer, it is also necessary to constrain the roll angles between the sprung and unsprung masses (φ-φu) so that they stay within the limits of the suspension travel (7-8deg), see [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF].

The performance characteristic which is of most interest when designing the active anti-roll bar, is then the normalized load transfer. The chosen control objective is to minimize the effect of the steering angle on the normalized load transfer R, in the H∞ framework. As explained later, the limitation of the torques Uf;r generated by the actuators is also crucial for practical implementation.

Background on H∞ control

The H∞ control problem is formulated according to the generalized control structure shown in Fig 2 [START_REF] Scherer | Linear matrix inequalities in control[END_REF], [START_REF] Skogestad | Multivariable Feedback Control[END_REF].

      =             (10)
  = = + -   (11)
The aim is to design a controller K that stabilizes the closed loop system and also reduces the signal transmission path from disturbances d to performance outputs z. This problem is then to find a controller K which minimizes γ such that ( , )

l F P K γ ∞ < (12)
By minimizing a suitably weighted version of ( , ) l F P K the control aim is achieved.

H∞ control synthesis for the active anti-roll bar of the single unit heavy vehicle model

In this section, the H∞ control design is presented for the active anti-roll bar system on a single unit heavy vehicle. Consider the closed-loop system given in Fig 3, which includes the feedback structure of the nominal model G, the controller K and the weighting functions Wij. In this diagram, Uf and Ur are the control inputs, y1 and y2 are the measured outputs, n1 and n2 are the measurement noises. δf is the steering angle considered as a disturbance signal, which is set by the driver. The variables z1, z2, z3, z4 and z5 represent the performance outputs.

A B B X X Z C D D W Y U C D D               =                    (13)
with the exogenous input (disturbance): The next section proposes an automatic and systematic way to obtain the optimal weighting functions used to solve the H∞ control design problem.

SELECTION OF WEIGHTING FUNCTIONS OF H∞ ANTI-ROLL BAR CONTROLLERS BY GENETIC ALGORITHMS

In industrial applications, multiple goals, often conflicting, have to be taken into account. Multi-criterion optimization (MCO) is then a powerful tool to find the best compromise solution balancing the conflicts, and is therefore of great importance in practice. In this section, MCO problem for the weighting function selection of the H∞ active anti-roll bar control on heavy vehicles is introduced and solved using the Genetic algorithms method. But first, the criterion for MCO problem must be defined.

Optimization objectives

The objective of the active anti-roll bar control system is to maximize the roll stability of heavy vehicles to prevent rollover in dangerous cases. However, such a performance objective must be balanced with the energy consumption of the anti-roll bar system due to the torques generation by the actuators. Therefore the objective function is selected as follows:

Normalized_load_transfer Torque

(1 )

f f f α α = + - (14
1 1 1 ( ) ( ) 2 1 1 ( ) ( ) 1 2 1 1 ( ) ( ) T T f r T T f r T T f r f R t dt R t dt T T U t dt U t dt T T f U t dt U t dt T T       = +                = +             ∫ ∫ ∫ ∫ ∫ ∫ (15)
where Rf;r are the normalized load transfers and Uf;r the torques generated at front and rear axles. Uf;rmax are defined when the optimal problem focusses only on the normalized load transfers (i.e. the torques are then not considered in the optimisation problem). In that case, 1 α = and Normalized_load_transfer f f = .

Multi-criterion optimization problem formulation

The weighting functions used in to the maximum expected value, corresponding to a 1 0 steering angle command. The weighting functions Wn1 and Wn2 are selected as:

1 2 0.01 n n W W = =
, which accounts for small sensor noise models in the control design. The noise weights are chosen as 0.01(m/s 2 ) for the lateral acceleration and 0.01( 0 /sec) for the derivative of the roll angle ϕ  [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF]. Note that other low pass filters could be selected if needed. The weighting functions Wzi represent the performance outputs (Wz1, Wz2, Wz3, Wz4 and Wz5). The purpose of the weighting functions is to keep small the control inputs, normalized load transfers and the lateral acceleration over a desired frequency range. The weighting functions chosen for performance outputs can be considered as penalty functions, that is, weights should be large in the frequency range where small signals are desired and small where larger performance outputs can be tolerated.

The weighting functions Wz1 and Wz2 corresponding to the front and rear control torques generated by active anti-roll bars are chosen as:

1 1 1 Z W Z = ; 2 2 1 Z W Z = (16)
The weighting functions Wz3 and Wz4 corresponding to the normalized load transfers at front and rear axles are selected as:

3 3 1 Z W Z = ; 4 4 1 Z W Z = (17)
The weighting function Wz5 is selected as:

52 53 5 51 54 55 Z Z s Z W Z Z s Z + = + (18)
Here, the weighting function Wz5 corresponds to a design that avoids the rollover situation with the bandwidth of the driver in the frequency range up to more than 4rad/s. This weighting function will directly minimize the lateral acceleration when it reaches the critical value, to avoid the rollover.

The parameters Zij are constant.

From equations ( 16) -( 18), the following variables will then be defined: Z1, Z2, Z3, Z4, Z51, Z52, Z53, Z54, Z55.

The MCO problem for the H∞ active anti-roll bar control can be defined as:

{ } _ _ 1 2 3 4 51 52 53 54 55 min ( ), ( ) : : [ , , , , , , , , ] |
T Normalized load transfer Torque p P T l u f p f p f f P p Z Z Z Z Z Z Z Z Z R p p p ∈   =   = = ∈ ≤ ≤ (19)
where f(p) is the vector of objectives, p the vector of weighting function parameters, p l and p u the lower and upper bounds of the weighting function selection.

The lower and upper bounds of the weighting function parameters are given in Table [START_REF] Boada | Active roll control using reinforcement learning for a single unit heavy vehicle[END_REF]. Besides the minimization of the objective function from equations ( 14) and (19), we also have to account for the limitations of the normalized load transfers, roll angle of suspensions as well as torques generated at each axle. These limitations are considered as the optimality conditions (binding conditions) shown in the Table [START_REF] Davis | Handbook of genetic algorithms[END_REF].

3 f R < 1 - 4 r R < 1 - 5 f U < 120000 Nm 6 r U < 120000 Nm

Genetic operation

The selection method used in this paper is the proportionate selection developed by Holland [START_REF] Holland | Adaptation in natural and artificial systems, an introductory analysis with application to biology, control and artificial intelligence[END_REF]. This method assigns a probability of selection to each individual, proportional to its relative fitness. Proportionate selection can be illustrated by a roulette wheel. The crossover happens with a probability of 0.9 and the mutation happens with a very small probability 0.095.

The proposed weighting function optimization procedure for H∞ active anti-roll bar control synthesis is as follows:

Step 1: Initialize with the weighting functions as in the previous paper [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF], the vector of weighting function selected as p=p0.

Step 2: Select lower bound, upper bound, scaling factor, offset and start point.

Step 3: Format optimal algorithm, select the vector of objectives with the variation of switch value from 0 to 1 and then solve the minimization problem.

Step 4: Select the individuals, apply crossover and mutation to generate a new generation: p=pnew.

Step 5: Evaluate the new generation by comparing with the binding conditions. If the criteria of interest are not satisfied, go to step 3 with p=pnew; else, stop and save the best individual: popt=pnew.

SIMULATION RESULTS

Optimization results

Thanks to the Genetic algorithms method, Table [START_REF] Do | Optimization of weighting function selection for H∞ control of semi-active suspensions[END_REF] synthesis the values of the variables Zi, Z5j in five cases for , the optimal problem focuses only on the torques generated.

Evaluation of optimization results in frequency domain

The limited bandwidth of the driver must be considered up to 4rad/s to identify any resonances in the response that may be excited by the driver [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF]. Therefore, it is necessary to consider the behavior of the heavy vehicle in a wider frequency range. In this section, the frequency response of heavy vehicle is shown in the nominal parameters case of the single unit heavy vehicle considered, characterized by the sprung mass ms= 12487kg, the forward velocity V at 70Km/h and the road adhesion coefficient µ= 1. , H∞ active anti-roll bar controllers reduce well the normalized load transfers compared to that of the case passive anti-roll bar and of the case of the previous paper [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF] which are chosen by trials and errors.

Evaluation of optimization results in time domain

In this section, the considered vehicle maneuver is cornering [START_REF] Gaspar | The design of a combined control structure to prevent the rollver of heavy vehicles[END_REF]. where the forward velocities are respectively 134, 131, 123, 104 Km/h. Note that in the previous paper [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF] we got 107 Km/h for the forward velocity. where the forward velocities are respectively 130, 125, 112, 97 Km/h. Note that in the previous paper [START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF], we got 104 Km/h for the forward velocity. In fact, the forward velocities of heavy vehicles are often considered up around 100 Km/h, so that choosing α between 1 and 0.65 is convenient.

CONCLUSIONS

In this paper, a weighting function optimization procedure using GAs for H∞ active anti-roll bar control on single unit heavy vehicle has been proposed. The conflicting objectives between the normalized load transfers and generated torques are handled using only one high level parameter, which is a great advantage to solve the multiobjective control problem. The simulation results in frequency and time domains have shown the efficiency of GAs in finding a suitable controller to satisfy the desired performance objectives.

Even if other structures for the weighting functions could be used, the ones used in this paper are shown to be simple enough while being efficient to solve the problem. For future works, the comparison with an LPV controller (scheduled by the forward vehicle velocity) will be of interest.

  of function J. By using various sets of α, one can generate several points in the Pareto set[START_REF] Do | Optimization of weighting function selection for H∞ control of semi-active suspensions[END_REF].

Fig. 1

 1 Fig. 1 Yaw-Roll model of single unit heavy vehicle [5].

Table 1 :φ

 1 Variables of yaw-roll model. Unsprung mass on the front axle Cf Tire cornering stiffness on the front axle mu;r Unsprung mass on the rear axle Cr Tire cornering stiffness on the rear axle m The total vehicle mass kf Suspension roll stiffness on the front axle v Forward velocity kr Suspension roll stiffness on the rear axle vwi Components of the forward velocity bf Suspension roll damping on the front axle h Height of CG of sprung mass from roll axis br Suspension roll damping on the rear axle hu;i Height of CG of unsprung mass from ground ktf Tire roll stiffness on the front axle r Height of roll axis from ground ktr Tire roll stiffness on the rear axle ay Lateral acceleration Ixx Roll moment of inertia of sprung mass β Side-slip angle at center of mass Ixz Yaw-roll product of inertial of sprung mass ψ Heading angle Izz Yaw moment of inertia of sprung mass ψ  Yaw rate lf Length of the front axle from the CG α Side slip angle lr Length of the rear axle from the CG φ Sprung mass roll angle lw Half of the vehicle width , u i Unsprung mass roll angle µ Road adhesion coefficient In the vehicle modelling, the differential equations of motion of the yaw-roll dynamics of the single unit vehicle, i.e. the lateral dynamics (3.1), the yaw moment (3.2), the roll moment of the sprung mass (3.3), the roll moment of the front (3.4) and the rear (3.5) unsprung masses, are formalized in the equations (3): 2 ( )

Fig. 2

 2 Fig. 2 Generalized control structure. with P partitioned as: 11 12 21 22 ( ) ( ) ( ) ( ) P s P s z d y P s P s u

Fig. 3 G

 3 Fig. 3 G-K control structure of H∞ active anti-roll bar control.

;

  where Uf /Ur are the torques at the front/rear axles. the performance output vector: and A, B1, B2, C1, D11, D12, C2, D21, D22 are matrices of appropriate dimensions.

  corresponding to the normalized load transfers and torques generated at two axles. They are defined as follows:

Fig

  problem focuses only on the normalized load transfers and when 0 α = , it means that Torque f f =

Fig. 4

 4 Fig. 4 shows the conflicting relation between the normalized load transfers and torques generated with some Pareto-optimal points, computed for the active anti-roll bar on heavy vehicle. They are generated for different values of α in the range of [0;1]. For α=0, Torquef

Figure. 4

 4 Figure. 4 The optimization results of some points in the Pareto frontier for the active anti-roll bar on heavy vehicle.

Fig 5 and 6

 6 Fig 5 and 6 show the normalized load transfers at the two axles Rf,r. They show that in the cases [1;0.85; 0.65] α =, H∞ active anti-roll bar controllers reduce well the normalized load transfers compared to that of the case passive anti-roll bar and of the case of the previous paper[START_REF] Vu | H∞ active anti-roll bar control to prevent rollover of heavy vehicles: a robustness analysis[END_REF] which are chosen by trials and errors.

Figure. 5

 5 Figure. 5 Frequency responses of the normalized load transfer at the front axle (Rf ) due to steering angle.

Figure. 6

 6 Figure. 6 Frequency responses of the normalized load transfer at the rear axle (Rr) due to steering angle.

Figure. 7

 7 Figure. 7 Steering angle δf during cornering maneuver.

Figure. 8

 8 Figure. 8 Time responses of the normalized load transfers (Rf,r) and torques generated at the two axles.

Figure. 9

 9 Figure. 9 Effect of the forward velocity on the normalized load transfer: front axle Rf.

Figure. 10

 10 Figure. 10 Effect of the forward velocity on the normalized load transfer: rear axle Rr.

Table 2 :

 2 Lower and upper bounds of the weighting functions.

		Wz1	Wz2	Wz3	Wz4			Wz5
		Z1	Z2	Z3	Z4	Z51	Z52	Z53	Z54	Z55
	Lower bound	50	50	0.1	0.1	0.5	1 3000	1	1	0.001

Table 3 :

 3 Binding conditions.

	No	Note	Maximum value	Unit
	1	uf φ φ -	< 7	deg
	2	ur φ φ -	< 7	deg

Table 4 :

 4 Optimization results for the weighting functions of H∞ active anti-roll bar.

	Controllers	Wz1 Z1	Wz2 Z2	Wz3 Wz4 Z3 Z4	Z51	Z52	Wz5 Z53	Z54	Z55
	SSSC2016 [18]	150000	200000	1	1	1	0.0005	50	100	0.01
	α =	1	258762.34 259996.91 0.91	0.64	0.97	0.46	392.54 801.38 0.22
	α =	0.85	273598.36 295104.17 0.45	0.26	1.17	0.80	334.15 968.50 0.16
	α =	0.65	112322.78 110837.28 0.72	0.75	0.63	0.54	139.23 97.46	0.02
	α =	0.5	166902.22 196036.53 1.09	0.92	0.97 0.0005 54.19 116.64 0.01
	α =	0	50	50	1.59	0.83	0.50	0.43	1	683.05 0.024