W Larbi

J.-F Deü

R Ohayon

A new finite element formulation for internal acoustic problems with dissipative walls

Keywords: dissipative acoustics, noise reduction, finite element method, symmetric formulation, wall normal fluid displacement field

This work concerns the variational formulation and the numerical computation of internal acoustic problems with absorbing walls. The originality of the proposed approach, compared to other existing methods, is the introduction of the normal fluid displacement field on the damped walls. This additional variable allows to transpose formulations in frequency domain to time domain when the fluid is described by a scalar field (pressure or fluid displacement potential). With this new scalar unknown, various absorbing wall models can be introduced in the variational formulation. Moreover, the associated finite element matrix system in symmetric form can be solved in frequency and time domain.

INTRODUCTION

In this paper, we are interested in the linear oscillations of an acoustic (i.e. inviscid, compressible, barotropic) fluid contained in a rigid cavity, with some or all of its walls covered by a thin layer of absorbing material able to dissipate part of the acoustic energy of the fluid.

Recently, considerable research has been devoted to the development and testing of noise reduction techniques by passive damping treatments. Many approaches have been proposed. For instance, porous materials have been widely used for noise and vibration control purposes. For example, a study of flexible porous damping materials is presented by Göransson [START_REF] Göransson | Acoustic finite element formulation of a flexible porous material-a correction for inertial effects[END_REF], where a finite element solution is proposed for the coupled acoustic-elastic wave propagation problem.

In a related approach, a wall impedance, characterizing the absorptive material, is used to take into account the effect of damping material. Let us mention the papers [START_REF] Kehr-Candille | Elastoacoustic damped vibrations-finite element and modal reduction methods[END_REF][START_REF] Bermúdez | Modelling and numerical solution of elastoacoustic vibrations with interface damping[END_REF][START_REF] Davidsson | Using the finite element method for control of sound fields in room acoustics[END_REF] where a wall impedance is introduced to describe the absorbing boundary conditions at dissipative walls. It can be noted that for a given material, the parameters of this impedance can be found experimentally [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF].

Different fields can be chosen to describe the behaviour of the fluid: pressure, fluid displacement, fluid potential displacement or fluid potential velocity (see, for example, [START_REF] Morand | Fluid-Structure Interaction[END_REF][START_REF] Ohayon | Structural Acoustics and Vibration[END_REF][START_REF] Ohayon | Reduced models for fluid-structure interaction problems[END_REF]).

In the particular case of dissipative acoustic problem, a pressure/displacement potential formulation in frequency domain of the fluid is given and numerically solved by finite element and modal reduction methods [START_REF] Kehr-Candille | Elastoacoustic damped vibrations-finite element and modal reduction methods[END_REF]. It should be noted that in this case, the frequency complex dependence does not permit to write the associated problem in time domain. In Reference [START_REF] Bermúdez | Modelling and numerical solution of elastoacoustic vibrations with interface damping[END_REF], a displacement formulation of dissipative acoustic problem is proposed. In this case, one has to pay attention to the discretization of the admissible class of irrotational motions. Special complex treatments are then required (Raviard-Thomas types of finite elements) and this vector formulation can lead to costly computation.

The aim of this paper is to develop a finite element formulation for damped acoustic cavity, which is valid in time and frequency domains using acoustic pressure as main variable. In order to achieve this objective, a new scalar variable is introduced to describe the damping interface namely the normal fluid displacement field. With this new scalar unknown, various interface damping models can be taken into account in the variational formulation. Moreover, the associated finite element matrix system is linear and can be solved in frequency and time domains. This new formulation leads to an unsymmetric matrix equation. A symmetric formulation can be obtained through the introduction of an intermediate unknown field, namely displacement potential of the fluid [START_REF] Morand | Fluid-Structure Interaction[END_REF]. Numerical examples are presented in order (i) to validate the new approach through comparison with modal results given in literature, and (ii) to evaluate the influence of the damping model on the dynamic response of the system. In this last case, a modal reduction approach is used in conjunction with a direct time integration method.

DISSIPATIVE ACOUSTIC BASIC LOCAL EQUATIONS

Consider an inviscid compressible barotropic fluid contained in a rigid cavity, with some of its walls covered by a thin layer of a viscoelastic material, in order to absorb part of the acoustic energy of the fluid. Let us denote by F the domain occupied by the fluid. The boundary of F is noted * F = A ∪ R ∪ P , with A the face covered by damping material and called 'absorbing boundary', R the rigid face and P the face where a given pressure p d is applied. The unit outer normal vector along * F is denoted by n. Additionally, the following notations for the quantities associated with the fluid: p for the pressure, u F for the fluid displacement field, F for the mass density and c for the speed of sound.

The local equations of the problem described in Figure 1 are

F * 2 u F *t 2 + ∇p = 0 in F (1a) p = -F c 2 div u F in F (1b) Ω F p d Γ R n Γ P Γ A Ω F Figure 1. Dissipative acoustic cavity. u F • n = 0 on R (1c) p = k I u F • n + d I *u F *t • n on A (1d) p = p d on P (1e)
Equation (1a) corresponds to the linearized Euler equations and Equation (1b) represents the constitutive equation of the barotropic fluid. The usual kinematic interface condition for rigid walls, Equation (1c), has been relaxed in the absorbing boundary to take into account the effect of the viscoelastic material represented by Equation (1d). On the right-hand side of this equation, the first term is proportional to the normal component of displacements and accounts for the elastic behaviour of the material, whereas the second one is proportional to the normal velocity and models the viscous damping. Finally, Equation (1e) represents the pressure force applied to the wall P of the cavity.

Starting from these equations, it is possible to obtain a displacement formulation by eliminating p in the previous system. The corresponding admissible space must verify curl u F = 0 in F in order to avoid, in the discretized case, the possible presence of spurious rotational eigenmodes. We refer to Reference [START_REF] Bermúdez | Finite element analysis of a quadratic eigenvalue problem arising in dissipative acoustics[END_REF] for a deep mathematical and numerical analysis of such formulation (throught the use of Raviart-Thomas finite elements).

DISSIPATIVE ACOUSTIC PROBLEM IN TERMS OF FLUID PRESSURE

In order to remove the problem of nonphysical modes and to arrive at a more compact system of equations, the acoustic problem is usually described either in terms of pressure or fluid potential displacement [START_REF] Morand | Fluid-Structure Interaction[END_REF]. However, in the case of dissipative acoustic cavity, it is not trivial to express the problem in time domain using only these variables. On the other hand, in frequency domain, and without considering the forcing terms (free motion), the problem can be written only in terms of pressure variable by elimination of u F in Equations (1) but leads to a nonlinear eigenvalue problem as described below. In this case, the pressure is expressed in the form p(x, t) = e t p(x) and the eigenvalue equations are given by p -

2 c 2 p = 0 in F (2a) ∇p • n = 0 on R (2b) ∇p • n = -F 2 k I + d I p on A (2c) p = 0 on P (2d)
Equation (2a) corresponds to Helmholtz equation in frequency domain (considering = -i ) deduced by the elimination of u F between Equations (1a) and (1b); Equation (2b) represents the rigid boundary condition on R ; and Equation (2c) represents the absorbing boundary condition on A written in terms of pressure field.

Thus, the variational formulation consists in finding and p belonging to the admissible space

C p = {p | p = 0 on p }, such that ∀ p ∈ C p : F ∇p • ∇ p dx + F 2 k I + d I A p p dx + 2 c 2 F p p d = 0 (3)
This formulation corresponds to that one proposed by Kehr-Candille and Ohayon [START_REF] Kehr-Candille | Elastoacoustic damped vibrations-finite element and modal reduction methods[END_REF], written here in terms of pressure instead of fluid displacement potential in Reference [START_REF] Kehr-Candille | Elastoacoustic damped vibrations-finite element and modal reduction methods[END_REF]. In the abovementioned paper, the equation is written in the form:

F ∇p • ∇ p dx -i F Z( ) A p p dx - 2 c 2 F p p d = 0 (4)
where is the angular frequency and Z( ) is the specific impedance of the absorbing wall, which is related to the damping parameters. For the Kelvin-Voigt model used in this work, and considering = -i , the specific impedance Z of the viscoelastic material is defined in terms of the coefficients k I and d I in the previous equation through the expression:

Z( ) = d I + i k I (5)
The real and imaginary parts of this specific impedance can be experimentally measured for any given frequency . These quantities are the specific resistance and reactance, respectively. It can be noted that in literature, where the time dependence of oscillating quantities is described by e i , the reactance is negative. For more details, the reader can refer to Reference [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF].

Equation [START_REF] Davidsson | Using the finite element method for control of sound fields in room acoustics[END_REF] shows that the dissipative acoustic problem written in terms of pressure leads to a nonlinear (nor even quadratic) eigenvalue problem and needs specific algorithms for its resolution. Frequency response functions can be obtained by sweeping frequency.

DISSIPATIVE ACOUSTIC FORMULATION IN TERMS OF FLUID PRESSURE AND FLUID NORMAL DISPLACEMENT AT THE ABSORBING BOUNDARY

It has been shown, in the previous section, that the dissipative acoustic problem can be expressed in terms of pressure unknown field in frequency domain and for free motion. In this section, we propose a new formulation, valid for the time and frequency domains, in terms of fluid pressure and fluid normal displacement at the absorbing boundary. In the following, this variable will be noted = u F • n. Thus, the dynamic equations are given in terms of this new variable and pressure p:

p - 1 c 2 * 2 p *t 2 = 0 in F (6a) ∇p • n = 0 on R (6b) ∇p • n = -F * 2 *t 2 on A (6c) p = k I + d I * *t on A (6d) p = p d on P (6e)
Equation (6a) corresponds to Helmholtz equation in time domain; Equation (6b) represents the rigid boundary condition on R ; Equation (6c) represents the absorbing boundary condition on A ; and Equation (6d) corresponds to the relation between fluid pressure and normal fluid displacement on A deduced from the behaviour of the viscolelastic material.

Multiplying Equation (6a) by a time-independent test function p ∈ C p , applying a Green's formula and dividing by F , we have

1 F F ∇p • ∇ p dx + 1 F c 2 F * 2 p *t 2 p dx + A * 2 *t 2 p d = 0 (7)
Similarly, let be the time-independent test function, associated to , belonging to the admissible space C . Multiplying Equation (6d) by ∈ C , we have

A p d -k I A d -d I A * *t d = 0 (8)
Thus, the variational (nonsymmetric) formulation of the acoustic cavity problem with absorbing wall consists, for given appropriate initial conditions, in finding p (with p = p d on P ) and such that Equations ( 7) and ( 8) are satisfied.

After discretizing by the finite element method, the bilinear forms of the variational formulation involved in Equations ( 7) and ( 8), we obtain the following matrix equation:

0 0 C T p K p Ḧ P + d I G 0 0 0 Ḣ Ṗ + k I G -C p 0 F p H P = 0 0 (9)
where H and P are the vectors of nodal values of and p, respectively (a subvector of P being given on P ), and the submatrices of Equation ( 9) are defined by

1 F F ∇p • ∇ p dx ⇒ P T F p P (10a) 1 F c 2 F * 2 p *t 2 p dx ⇒ P T K p P (10b) A * 2 *t 2 p d ⇒ P T C T p Ḧ (10c) A p d ⇒ H T C p P (10d) A d ⇒ H T G H (10e)
The damped vibration modes of the fluid are complex solutions of Equations ( 7) and ( 8), of the form p(x, t) = p(x)e t and (x, t) = (x)e t and setting p d = 0 on P . The corresponding matrix quadratic eigenvalue equation is deduced from Equation ( 9) with P = 0 on P (for sake of brevity, we have not changed the matrix notations):

k I G -C p 0 F p + d I G 0 0 0 + 2 0 0 C T p K p H P = 0 0 (11)
We can note that the fluid normal displacement at the absorbing boundary can be expressed in terms of the fluid pressure from the first line of Equation ( 11):

H = 1 k I + d I G -1 C p P (12)
Replacing H by its expression (12) in the second line of Equation ( 11), we obtain the following nonlinear eigenvalue problem in terms of P:

F p + 2 k I + d I C T p G -1 C p + 2 K p P = 0 (13)
Equation ( 13) can be directly obtained through discretizing the pressure formulation given in Equation ( 3). Moreover, it can be easily seen that if d I = 0 and k I tends to infinity, eigenvalues of acoustic problem in rigid cavity are obtained.

DISSIPATIVE ACOUSTIC SYMMETRIZED FORMULATION IN TERMS OF ( , p)

In the previous section, a nonsymmetric quadratic eigenvalue formulation has been obtained. The aim of this section is to establish a symmetric formulation of this dissipative acoustic problem in view of a direct treatment by finite elements. This symmetrization is obtained through introduction of an intermediate unknown field [START_REF] Morand | Fluid-Structure Interaction[END_REF]: a fluid displacement potential field , defined up to an additive constant, such that u F = ∇ . We shall see in the next section that this problem can be reduced, by means of an elimination procedure, to an eigenvalue problem in terms of ( , p).

The equations of the acoustic spectral problem in terms of ( , p, ) can be rewitten as follows:

F + p c 2 = 0 in F (14a) ∇ • n = 0 on R (14b) ∇ • n = on A (14c) = 0 on P (14d) p F c 2 = - 2 c 2 in F (14e) -2 F = k I + d I on A (14f)
The corresponding variational formulation, involving fluid normal displacement at the absorbing boundary, fluid potential displacement field, and pressure can be easily obtained as described below:

• Let be the test function, associated to , belonging to the admissible space C = { | = 0 on p }. Multiplying Equation (14a) by ∈ C and taking into account Equations (14c) and (14b), we obtain

-F F ∇ • ∇ x + F A d + 1 c 2 F p dx = 0 ∀ ∈ C ( 15
)
• Multiplying Equation (14e) by a test function p ∈ C p , we have

1 F c 2 F p p dx + 2 c 2 F p dx = 0 ∀ ∈ C p ( 16
)
• Multiplying Equation (14f) by a test function

∈ C , we have 2 F A d + d I A d + k I A d = 0 ∀ ∈ C (17)
Then, the variational formulation consists in finding and ( , p, ) ∈ (C , C p , C ) satisfying Equations ( 15)-( 17).

The discretization of the preceding formulation leads to the following symmetric matrix equation:

⎛ ⎜ ⎜ ⎝ k I G 0 0 0 K p 0 0 0 0 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P ⎞ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ d I G 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P ⎞ ⎟ ⎠ + 2 ⎛ ⎜ ⎜ ⎝ 0 0 A 0 0 B A T B T -F ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P ⎞ ⎟ ⎠= ⎛ ⎜ ⎝ 0 0 0 ⎞ ⎟ ⎠ (18)
where H, P and are the vectors of nodal values of , p and , respectively, and the submatrices of Equation ( 18) are defined by

F A d ⇒ H T A (19a) F F ∇ • ∇ dx ⇒ T F (19b) 1 F c 2 F * 2 p *t 2 p dx ⇒ P T K p P (19c) 1 c 2 F p d ⇒ P T B (19d) A d ⇒ H T G H (19e)
We can note that the fluid normal displacement at the absorbing boundary and the fluid pressure can be expressed in terms of fluid potential displacement from the first and the second line of Equation ( 18):

H = -2 k I + d I G -1 A (20a) P = -2 K -1 p B (20b)
Replacing H and P by their expressions, Equation (20), in the third line of Equation ( 18), we obtain the following nonlinear eigenvalue problem in terms of :

F + 2 k I + d I A T G -1 A + 2 B T K -1 p B = 0 (21)
This equation can be obtained directly by using a fluid potential displacement formulation. This formulation has been established by Kehr-Candille and Ohayon [START_REF] Kehr-Candille | Elastoacoustic damped vibrations-finite element and modal reduction methods[END_REF]. It is clear that this formulation presents the inconvenience that leads to a nonlinear eigenvalue problem which cannot be solved directly.

SYMMETRIC REDUCED MODEL FOR THE DISSIPATIVE ACOUSTIC EIGENVALUE PROBLEM

In this section, in order to obtain a symmetric reduced model for the dissipative acoustic problem, we propose to eliminate the potential fluid displacement in the symmetric formulation given in the last section (Equation ( 18)). First, the third line of Equation ( 18) is written as

A T H + B T P -F = 0 (22)
As the matrix F is not singular (due to = 0 on P ), we can express in terms of H and P from Equation ( 22):

= F -1 A T H + F -1 B T P (23)
Substituting by its expression (Equation ( 23)) in first and second lines of Equation (18), we obtain the following symmetric reduced system in terms of and p:

k I G 0 0 K p H P + d I G 0 0 0 H P + 2 AF -1 A T AF -1 B T BF -1 A T BF -1 B T H P = 0 0 (24)

Case in which P = ∅

In this case, Equation (14d) is removed from the local equations (14a)-(14f). The associated variational formulation is the same as Equations ( 15)-( 17) with admissible spaces without any constraints. In Equation (18), the matrix F is now singular and consequently cannot be eliminated. In order to reduce the problem, we use the procedure described by Morand and Ohayon for elasto-acoustic problem [START_REF] Morand | Fluid-Structure Interaction[END_REF] as detailed below.

Let us define a partitioning of by letting T = [ 1 T 2 ] where 1 denotes a particular component of (here, the first), which induces the following partitioning of the matrices A, B and F in Equation ( 18):

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ k I G 0 0 0 0 K p 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ H P 1 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ d I G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ H P 1 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ + 2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 A 1 A 2 0 0 B 1 B 2 A T 1 B T 1 -F 11 -F 12 A T 2 B T 2 -F 21 -F 22 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ H P 1 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (25)
F 22 being nonsingular, 2 can be expressed in terms of P, H and 1 using the fourth line of Equation ( 25):

2 = -F -1 22 F 21 1 + F -1 22 A T 2 H + F -1 22 B T 2 P (26)
Substituting this relation for the remaining, the third line of Equation (25) becomes:

(A T 1 -F 12 F -1 22 A T 2 )H + (B T 1 -F 12 F -1 22 B T 2 )P = (F 11 -F 12 F -1 22 F 21 ) 1 (27)
It can be shown that [START_REF] Morand | Fluid-Structure Interaction[END_REF]:

F 11 -F 12 F -1 22 F T 12 = 0 (28)
Thus, Equation (27) can be written as

a T H + b T P = 0 (29)
where matrices a and b are defined by

a T = A T 1 -F 12 F -1 22 A T 2 (30a) b T = B T 1 -F 12 F -1 22 B T 2 (30b)
Finally, replacing 2 by its expression (26) in the two first lines of Equation (25) and using Equation (29), the following matrix system is obtained:

⎛ ⎜ ⎜ ⎝ k I G 0 0 0 K p 0 0 0 0 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P 1 ⎞ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ d I G 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P 1 ⎞ ⎟ ⎠ + 2 ⎛ ⎜ ⎜ ⎝ A 2 F -1 22 A T 2 A 2 F -1 22 B T 2 a B 2 F -1 22 A T 2 B 2 F -1 22 B T 2 b a T b T 0 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎝ H P 1 ⎞ ⎟ ⎠= ⎛ ⎜ ⎝ 0 0 0 ⎞ ⎟ ⎠ (31)
It is important to note that previous reduction is not an approximation but a rigorous condensation procedure.

Reduced eigenvalue problem

The reduced eigenvalue problem (Equation ( 24) or (31) for the case P = ∅) corresponds in fact to an indirect (via the intermediate displacement potential unknown field) symmetric formulation in terms of and p. In a condensed form, these matrix equations can be written as

(K + D + 2 M)X = 0 (32)
In practice, this second-order eigenvalue problem can be solved by transforming it into a system of 2n first-order equations by complementing it with the equality MX -MX = 0, where n is the number of degrees of freedom (d.o.f.) in the system. The extended spectral equation takes the canonical form:

K 0 0 -M + D M M 0 X X = 0 (33)
The matrices M, D and K being real symmetric, complex eigenvalues come in conjugate pairs associated with conjugate pairs of mode shapes.

NUMERICAL EXAMPLES

In this section, numerical results are obtained via a specific Matlab code especially developed for this project.

Modal analysis of a square cavity with a damping wall

In order to validate the finite element implementation of the proposed damping formulation, we firstly consider the spectral problem of a square cavity of size (1 m × 1 m) completely filled with air and with an absorbing wall (Figure 2). In this example, initially proposed in Reference [START_REF] Bermúdez | Modelling and numerical solution of elastoacoustic vibrations with interface damping[END_REF], the following data are used: F = 1 kg/m 3 , c = 340 m/s, k I = 5 × 10 6 Pa/m and d I = 50 Pa s/m. The last two parameters are average impedance coefficients corresponding to a typical acoustic insulating fabric (a Johns Manville glass wool) in the frequency range (50-500 Hz) (see Figure 3). In Table I, the first five eigenfrequencies (in Hz) with uniform meshes (quadrangular element) and with increasing number of d.o.f. are presented.

The first and second columns present the eigenfrequencies of the rigid cavity computed from the pressure formulation and compared to the exact solution given in Reference [START_REF] Blevins | Formulas for Natural Frequency and Mode Shapes[END_REF] for a rectangular cavity of lengths A and B (here, A = B = 1 m):

f nm (Hz) = c 2 n A 2 + m B 2 , n,m= 0, 1, 2, . . . , (n + m = 0) (34)
This comparison enables us to validate the development of the fluid finite element. The three other columns correspond to the complex frequencies of the damping cavity computed from the ( , p)-symmetric reduced formulation (Equation (31)) and compared to exact solution (last column) given by Bermúdez and Rodríguez [START_REF] Bermúdez | Modelling and numerical solution of elastoacoustic vibrations with interface damping[END_REF]. A good agreement between exact and computed values can be observed even for the coarse mesh. In this example, the imaginary part of the frequencies derives from dashpot dissipation. Moreover, the difference between the real part (damped case) and the real value (undamped case) of the frequencies is due to the spring effect.

It can be noted that the first five mode shapes are different in the damped and undamped cases, except for the third one (Figures 456). This is due to the lack of symmetry appearing in the problem with one dissipative wall.

Modal analysis of a damping rectangular cavity

The second example concerns a 2D rectangular rigid cavity (A = 1 m, B = 0.75 m) filled with air ( f = 1 kg/m 3 , c = 340 m/s). As in the previous example, a thin layer of absorbing material covers the top wall of the cavity. Here, the parameters of the Kelvin-Voigt model are: Figures 789show the mode shapes and Table II gives the eigenvalues of the discretized problem given in Equation (31). We have computed them with different successively refined uniform meshes. In the last column of the table, the exact eigenvalues given by Bermúdez et al. [START_REF] Bermúdez | Finite element analysis of a quadratic eigenvalue problem arising in dissipative acoustics[END_REF] are included. As it can be seen, our solution converges to the exact one when refining the mesh. In this example, the imaginary part of frequencies is much more important than in the previous example. Here, the real and imaginary parts of the frequencies can have the same order of magnitude for some of the first modes. Table III shows the eigenvalues of the fluid for various boundary conditions: (i) rigid boundary conditions along all edges, (ii) null pressure on P and rigid boundary elsewhere, (iii) null pressure on P , absorbing wall on A , and rigid boundary on R . As it can be seen, the difference between frequencies in the first (1111 d.o.f.) and second (1100 d.o.f.) columns is due to the elimination of the fixed d.o.f. on the left edge where a pressure load is applied. In the third column, where additional normal fluid displacements are taken into account, the real parts of the complex frequencies of the dissipative pipe are close to the frequencies of the rigid pipe. Moreover, the imaginary parts, which are small compared to the real parts, are associated with the damping effect.

k I = 5 ×
Equation (18) can be transposed in time domain in a straithforward manner, leading to a standard time second-order matrix system. Let us partition the vector of nodal pressures into P (where the pressure is prescribed) and P (where the pressure is unknown). Similarly, we divide the fluid potential displacement into and , where is given from P using the following relation:

* 2 *t 2 = - p F in F (35)
deduced from the linearized Euler equations.

With this separation into known and unknown components, Equation (18) can be transposed in time domain as follows:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ k I G 0 0 0 0 0 K p K pp 0 0 0 K pp K p 0 0 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ H P P ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ d I G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Ḣ Ṗ Ṗ ˙ ˙ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 A A 0 0 0 B p B p 0 0 0 B p B p A T B T p B T p -F -F A T B T p B T p -F -F ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Ḧ P P ¨ ¨ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (36)
The third and the fifth lines of the previous matrix system are eliminated because they correspond to the known nodal values P and . Using the fourth line of Equation (36), ¨ can be expressed in terms of other variables in the following form:

¨ = F -1 (A T Ḧ + B T p P + B T p P -F ¨ ) (37)
Finally, replacing ¨ by its expression (Equation (37)) in the first and second lines of Equation (36) and writing all of the known variables on the right-hand side, we obtain the following symmetric matrix system:

k I G 0 0 K p H P + d I G 0 0 0 Ḣ Ṗ + A F -1 A T A F -1 B T p B p F -1 A T B p F -1 B T p Ḧ P = F F p ( 38
)
where F and F p are given by

F = (-A + A F -1 F ) ¨ -A F -1 B T p P (39a)
F p = (-B p + B p F -1 F ) ¨ -B p F -1 B T p P -K pp P (39b)
and with: ¨ =-P/ F (see Equation ( 35)).

It should be noted that the matrix system (38) can be directly integrated in time domain. However, it can be costly because the mass matrix is fully populated. Alternatively, it can used in its spectral form in a modal reduction procedure as described below.

Modal reduction.

The superposition method is applied here in order to reduce the finite element model proposed in Equation (38), which can be written in condensed form:

M Ẍ + D Ẋ + KX = F (40)
Therefore, Equation (40) has to be transformed into modal or generalized coordinates of much small degrees of freedom. The undamped modal data can be obtained by solving the real eigenvalue problem:

K{} + 2 M{} = 0 (41)
The modal transformation and truncation of matrix system (Equation ( 40)) is accomplished by using 'r' columns of the modal matrix [] ('r' represents the number of the first most significant vibration modes of the fluid). As a result, our differential system is reduced to 'r' d.o.f., which are solved using the Newmark numerical integration technique:

M R ẌR + D R ẊR + K R X R = F R (42)
where, the superscript R represents the reduced system. We then transferred the modal solutions back to the original coordinate system using the following transformation:

X = [] R X R (43)
7.3.2. Numerical results. In this last section, transient response of the damped pipe is computed using a direct time integration method applied (i) to the global unsymmetric ( , p) system (Equation ( 9)), and (ii) to the reduced symmetric system (Equation (42)). Two types of excitation are considered:

1. Pressure pulse: p d (t) = p 0 sin( t/T 1 ) for 0 t T 1 0 for t T 1 2. Harmonic excitation:

p d (t) = p 0 sin(2 t/T 1 )
where p 0 = 1000 Pa is the amplitude of the pressure load, T 1 = 20 s is the pressure pulse length (case 1) or the period of excitation (case 2). Figure 11 shows time and frequency responses of the dissipative acoustic pipe in a point M (2.5 m, 0.4 m) when a pressure pulse is applied to the left edge. As shown in Figure 11(a), the pressure levels of the undamped system have been significantly reduced due to the damping effect. Figure 11(b) presents the decrease in the peak levels that were obtained using the absorption wall at the right edge of the pipe. These figures clearly show that there is a substantial reduction in pressure amplitude due to the addition of the damping interface.

CONCLUSION

This paper presents a new formulation for analysis of internal acoustic problem with absorbing walls which is based on the fluid pressure field and the normal fluid displacement on the damped walls. Local equations and corresponding variational formulations leading to a reduced-order symmetric matrix system are presented. Numerical results are presented in order to illustrate the effectiveness of the proposed method. Extension to structural-acoustic problems are in progress [START_REF] Deü | Structural-acoustic vibration and transient problems with interface damping[END_REF] and will be presented in further papers.

Figure 2 .

 2 Figure 2. Fluid in square rigid cavity with an absorbing wall.

Figure 3 .

 3 Figure 3. Acoustic impedance for Johns Manville glass wool [2].

  10 4 Pa/m and d I = 200 Pa s/m.

Figure 4 .

 4 Figure 4. First five acoustic modes for the square rigid cavity.

Figure 5 .

 5 Figure 5. First five acoustic modes for the square cavity with an absorbing wall (real part).

Figure 6 .

 6 Figure 6. First five acoustic modes for the square cavity with an absorbing wall (imaginary part).

Figure 7 .

 7 Figure 7. First five acoustic modes for a rectangular rigid cavity.

Figure 8 .

 8 Figure 8. First five acoustic modes for the rectangular cavity with an absorbing wall (real part).

Figure 9 .

 9 Figure 9. First five acoustic modes for the rectangular cavity with an absorbing wall (imaginary part).

7. 3 .Figure 10 .

 310 Figure 10. Straight pipe with a pressure pulse applied at the left end.

Table I .

 I Frequencies (Hz) of a square rigid cavity with an absorbing wall.

	170.078 170.078 240.526 340.622 340.622	f 01 = 170 f 10 = 170 f 11 = 240.416 f 02 = 340 f 20 = 340	166.790-0.040i 168.737-0.023i 236.019-0.078i 337.365-0.168i 341.094-0.123i	166.233-0.039i 167.982-0.023i 235.073-0.079i 332.910-0.155i 335.327-0.141i	166.165-0.039i 167.919-0.023i 234.957-0.079i 332.370-0.154i 334.720-0.144i

Undamped Damped 961 d.o.f. Exact [10] 121 d.o.f. 961 d.o.f. Exact [3]

Table II .

 II Frequencies (rad/s) of a rectangular rigid cavity with an absorbing wall.

	475 d.o.f.	825 d.o.f.	1271 d.o.f.	Exact [9]
	267.34 + 323.69i	267.48 + 322.38i	267.54 + 321.78i	267.65 + 320.71i
	813.50 + 259.37i	813.41 + 259.30i	813.36 + 259.27i	813.29 + 259.21i
	1281.94 + 90.53i	1281.68 + 90.28i	1281.56 + 90.16i	1281.35 + 89.95i
	2187.21 + 298.66i	2184.56 + 298.03i	2183.33 + 297.74i	2181.15 + 297.21i
	2256.13 + 27.90i	2253.63 + 027.66i	2252.47 + 27.56i	2250.41 + 27.37i
	2415.15 + 239.82i	2412.54 + 238.45i	2411.34 + 237.82i	2409.21 + 236.70i
	3032.04 + 147.80i	3028.395 + 145.74i	3026.69 + 144.80i	3023.68 + 143.16i
	3302.15 + 13.24i	3293.36 + 12.99i	3289.29 + 12.88i	3282.07 + 12.69i
	3616.68 + 306.56i	3604.30 + 304.85i	3598.59 + 304.05i	3588.43 + 302.60i
	3765.36 + 281.50i	3753.29 + 278.85i	3747.71 + 277.61i	3737.81 + 275.41i

Table III .

 III Frequencies (Hz) of the pipe.

	Rigid cavity		Dissipative cavity
	(i)	(ii)	(iii)
	1111 d.o.f.	1100 d.o.f.	1111 d.o.f.
	34.00	17.00	16.62 -0.02i
	68.01	51.00	49.88 -0.17i
	102.04	85.02	83.23 -0.45i
	136.09	119.06	116.69 -0.82i
	170.17	153.13	150.28 -1.25i
	204.30	187.23	184.01 -1.70i
	238.48	221.38	217.87 -2.15i
	272.72	255.59	251.86 -2.59i
	307.02	289.86	285.96 -3.02i
	341.40	324.20	320.19 -3.42i

ACKNOWLEDGEMENTS

This research is partially carried out within the 6th Framework Programme of the European Commission (InMAR project) under Grant No. NMP2-CT-2003-501084.

Figure 12(a) illustrates the effect of the damping wall in the case of harmonic excitation. The addition of an absorbing boundary caused a small shift in some of the resonant frequencies of pipe except for the seventh, which corresponds to excitation frequency (200 Hz). This shift is due to the spring effect of the Kelvin-Voigt model. This figure also indicates that the damping significantly reduced the first resonant peaks, with the largest reductions achieved after the excitation peak. Figure 12(b) compares the noise levels for the damped pipe using a time integration method applied to the global unsymmetric ( , p) system or a modal reduction technique applied to the symmetric ( , p) formulation presented in Section 7.3.1 using only the 10 first fluid modes. The result improves the procedure of symmetrisation and the performance of modal reduction method.


W Larbi

J.-F Deu

R Ohayon

Vibration of axisymmetric composite piezoelectric shells coupled with internal fluid

Keywords: axisymmetric conical shells, piezoelectric laminates, fluid-structure interaction, free vibrations, finite element method

This paper presents the theoretical and finite element formulations of piezoelectric composite shells of revolution filled with compressible fluid. The originality of this work lies (i) in the development of a variational formulation for the fully coupled fluid/piezoelectric structure system, and (ii) in the finite element implementation of an inexpensive and accurate axisymmetric adaptive laminated conical shell element. Various modal results are presented in order to validate and illustrate the efficiency of the proposed fluid-structure finite element formulation.

INTRODUCTION

In the context of noise reduction techniques by active treatments, this work presents the theoretical formulation and the finite element implementation of the specific problem of piezoelectric composite shells of revolution filled with compressible fluid. The originality of this work lies (i) in the modal formulation of the fully coupled fluid/piezoelectric structure system, and (ii) in the development of an inexpensive and accurate adaptive axisymmetric laminated conical shell element.

For the resolution of structural-acoustic vibration problems, two classes of methods can be distinguished according to whether they are based on analytical or numerical discrete approaches. Analytical methods are classically restricted to simple geometries and boundary conditions. For example, the vibration of cylindrical thin shells filled with a compressible fluid has been studied using a wave propagation approach in [START_REF] Zhang | Coupled vibration analysis of fluid-filled cylindrical shells using the wave propagation approach[END_REF] for elastic structures, or using a state-space method in [START_REF] Ding | Free vibration of piezoelectric cylindrical shells filled with compressible fluid[END_REF] for piezoelectric structures. Three-dimensional exact solutions have recently been proposed for the free vibration analysis of arbitrary thick simply supported piezoelectric hollow cylinders filled with compressible fluid [START_REF] Chen | 3D free vibration analysis of a functionally graded piezoelectric hollow cylinder filled with compressible fluid[END_REF][START_REF] Deü | A state space method for free-vibration analysis of radially polarized laminated piezoelectric cylinder filled with fluid[END_REF]. Methods based on discretization of the structural-acoustic governing equations, such as the finite element method, are not restricted to specific boundary conditions but require a higher computational effort. For example, a mixed formulation based on a piezoelectric shell finite element for the structure and a three-dimensional boundary element for the fluid, is developed in [START_REF] Kaljevic | Steady-state response of acoustic cavities bounded by piezoelectric composite shell structures[END_REF] to calculate the coupled response of smart structural-acoustic systems. For problems involving axisymmetric geometries, semi-analytic approaches are frequently used because they are computationally more efficient and accurate than fully discrete methods approaches. These methods, which employed a discretization in Fourier series for the circumferential direction, have recently been developed for vibrations analysis of axisymmetric piezoelectric shells [START_REF] Hussein | Discrete layer analysis of axisymmetric vibrations of laminated piezoelectric cylinders[END_REF][START_REF] Benjeddou | Use of piezoelectric shear response in adaptive sandwich shells of revolutionpart 1: theoretical formulation[END_REF][START_REF] Correia | Active control of axisymmetric shells with piezoelectric layers: a mixed laminated theory with a high order displacement field[END_REF], as well as axisymmetric laminated composite cylinders filled with fluid [START_REF] Xi | Free vibration of a laminated composite circular cylindrical shell partially filled with fluid[END_REF][START_REF] Zhang | A comparative study of axisymmetric finite elements for the vibration of thin cylindrical shells conveying fluid[END_REF][START_REF] Kochupillai | Model reduction for parametric instability analysis in shells conveying fluid[END_REF].

In this paper, a semi-analytical approach is used for adaptive conical laminated shells filled with compressible fluid. The piezoelectric shell model is based on the Kirchhoff-Love theory and is supposed to be polarized in thickness direction. The associated finite element model combines an equivalent single-layer approach for the mechanical behaviour with a layerwise representation of the electric potential in the thickness direction. In this work, we consider that the fluid is inviscid, compressible and barotropic, gravity effect being neglected. Extension to vibration damping of more complex systems involving surface gravity effects in the inviscid incompressible case [START_REF] Schotté | Incompressible hydroelastic vibrations, finite element modelling of the elastogravity operator[END_REF] or internal gravity effects in the inviscid compressible case [START_REF] Andrianarison | Compressibility and gravity effects in internal fluid-structure vibrations: basic equations and appropriate variational formulations[END_REF] will be the purpose of further investigations. Considering the above assumptions, the fluctuating fluid forces are determined by means of the linearized Euler equations and reduced to the Helmholtz equation for the pressure. The dynamic formulation of the coupled system is derived from a variational principle involving shell displacements and rotation, electric potential and fluid pressure. Concerning the finite element discretization, a one-dimensional element with two nodes is used for the adaptive conical shell structure and a quadrangular element with four nodes for the fluid. Moreover, an appropriate interface element is developed to ensure the interaction between the fluid and structure. For the shell element, each node has (i) four mechanical degrees-of-freedom to describe axial, radial and circumferential displacements, as well as the rotation about the circumferential axis, and (ii) one electrical degree-of-freedom related to the sensor voltage per piezoelectric layer. For the fluid element, the pressure is the unique nodal unknown. It should be mentioned that the variables associated with these two elements are expanded in Fourier series in the circumferential direction. Finally, some numerical examples are presented in order to validate and demonstrate the effectiveness of the proposed formulation.

FINITE ELEMENT FORMULATION OF THE PIEZOELECTRIC STRUCTURE COUPLED WITH AN INTERNAL FLUID

Let us consider the linear vibrations of an elastic/piezoelectric structure completely filled with a homogeneous, inviscid and compressible fluid, neglecting gravity effects. We establish in this section the variational formulation of the coupled problem and the corresponding matrix equations.

Governing equations of the spectral-coupled problem

We consider a piezoelectric structure occupying the domain S at the equilibrium. The structure is clamped on a part u and traction free on the complementary part of its external boundary. The interior fluid domain is denoted by F and the fluid-structure interface by (see Figure 1). For the considered spectral problem, electric boundary conditions are defined by null electric potential on and null surface density of electric charge on the remaining part D . Thus, the total structure boundary is

* S = u ∪ ∪ = ∪ D with u ∩ ∩ = ∩ D = ∅.
The linearized deformation tensor is denoted by e and the corresponding stress tensor by r. Moreover, D denotes the electric displacement and E the electric field. S is the mass density of the structure and n S is the unit normal external to S .

Since the compressible fluid is assumed to be inviscid, instead of describing its motion by a fluid displacement vector field, which requires an appropriate discretization of the fluid irrotationality constraint [START_REF] Bermúdez | Finite element computation of the vibration modes of a fluid-solid system[END_REF], we use the pressure scalar field p. Let us consider c F as the constant speed of sound in the fluid, F the mass density of the fluid, and n the unit normal external to F .

The local equations describing the spectral fluid-piezoelectric structure-coupled problem are div r + 2 S u = 0 in S (1a)

r n S = 0 on (1b) u = 0 on u (1c) rn S = pn on (1d) div D = 0 in S (2a) D • n S = 0 on D (2b) = 0 on (2c) p + 2 c 2 F p = 0 i n F (3a) * p *n = 2 F u • n on (3b)
where is the angular frequency. Equation (1a) corresponds to the elastodynamic equation in the absence of body force; Equations (1b) and (1c) are the classical mechanical boundary conditions; Equation (1d) results from the action of pressure forces exerted by the fluid on the structure; Equation (2a) corresponds to the electric charge equation for a dielectric medium; Equations (2b) and (2c) are the previously described electric boundary conditions; Equation (3a) is the Helmholtz equation; and Equation (3b) is the contact condition for the fluid on .

The stress tensor r and electric displacement vector D are related to the linear strain tensor e and electric field vector E through the converse and direct linear piezoelectric constitutive equations

r = ce -e T E (4)
D = ee + dE (5)
where c is the elastic stiffness matrix, d is the dielectric permittivity matrix, and e is the piezoelectric coupling matrix. Moreover, we have the following gradient relations between the linearized strain tensor e and the displacement u, and between the electric field E and the electric potential :

e = 1 2 (∇u + ∇ T u) (6)
E = -∇ (7)
For a detailed derivation of those classical equations, we refer the reader, for example, to [START_REF] Tzou | Piezoelectric Shells: Distributed Sensing and Control of Continua[END_REF] for piezoelectric aspects and to [START_REF] Morand | Fluid-Structure Interaction[END_REF] for fluid-structure aspects.

Variational formulation in terms of (u, , p)

The local equations of Section 2.1 are expressed in terms of the chosen unknown fields of the piezoelectric structure/fluid boundary value problem, i.e. the structural mechanical displacement u, the electric potential in the structure , and the fluid pressure p.

In order to obtain the variational formulation associated with the local equations of the coupled fluid-piezoelectric structure system given in Equations ( 1)-(3), the test-function method is applied. We proceed in three steps, successively considering the equations relating to the structure (subject to fluid pressure actions), the electric charge equation for a dielectric medium, and the equations relating to the fluid (subject to a wall displacement). We recall that the chosen unknown fields are the structure displacement u, the electric potential and the fluid pressure p.

First, we introduce the space C * u of sufficiently regular functions u defined in S verifying u = 0 on u . Multiplying Equation (1a) by any function u ∈ C * u , then applying Green's formula, and finally taking Equations (1b) and (1d) into account, leads to S r(u, ) : e( u) dvpn

• u ds -2 S S u • u dv = 0 ∀ u ∈ C * u (8)
Secondly, we consider the space C * of sufficiently regular functions in S verifying = 0 on . Multiplying Equation (2a) by any function ∈ C * and integrating over S , we have

S div D(u, ) dv = 0 ∀ ∈ C * (9)
Finally, we consider the space C p of sufficiently regular functions p defined in F . Multiplying Equation (3a) by p ∈ C p , applying Green's formula, and taking Equation (3b) into account, we obtain

1 F F ∇p • ∇ p dv - 2 F c 2 F F p p dv -2 u • n p ds = 0 ∀ p ∈ C p ( 10
)
Using the piezoelectric constitutive Equations ( 4) and ( 5), and taking the electric boundary condition (2b) into account, Equations ( 8) and ( 9) can be rewritten in the following forms:

Finite element formulation

Let us introduce U, U and P corresponding to the vectors of nodal values of u, and p, respectively, and the submatrices corresponding to the various bilinear forms involved in Equations ( 10)-( 12) defined by

S ce(u) : e( u) dv ⇒ U T K u U (13a) S S u • u dv ⇒ U T M u U (13b) S e T E() : e( u) dv ⇒ U T C u U (13c) S ee(u) • E( ) dv ⇒ U T C T u U (13d) pn • u ds ⇒ U T C up P (13e) 1 F F ∇p • ∇ p dv ⇒ P T K p P (13f) 1 F c 2 F F p p dv ⇒ P T M p P (13g) S dE() • E( ) dv ⇒ U T K U (13h)
The variational equations ( 10)-( 12) for the fluid-piezoelectric structure-coupled problem can be written in discretized form

⎡ ⎢ ⎢ ⎣ K u -C u -C up C T u K 0 0 0 K p ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎣ U U P ⎤ ⎥ ⎦ = 2 ⎡ ⎢ ⎢ ⎣ M u 0 0 0 0 0 C T up 0 M p ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎣ U U P ⎤ ⎥ ⎦ (14)

Remarks

• The matrices involved in Equation ( 14) are unsymmetric. A symmetric formulation with mass coupling can be established by introducing the fluid displacement potential as additional variable [START_REF] Morand | Fluid-Structure Interaction[END_REF]. For the sake of brevity, this procedure is not detailed here and will be the subject of a future work. • Using the second line of Equation ( 14), the degrees-of-freedom associated with the electric potential can be expressed in terms of structure displacements as

U = -K -1 C T u U (15)
Thus, after substitution of U into Equation ( 14), we get the following spectral problem in terms of U and P:

K u + C u K -1 C T u -C up 0 K p U P = 2 M u 0 C T up M p U P (16)
It should be noted that the term C u K -1 C T u represents an 'added-stiffness' contribution due to the electromechanical coupling [START_REF] Galucio | A fractional derivative viscoelastic model for hybrid active-passive damping treatments in time domain-application to sandwich beams[END_REF].

• A standard unsymmetric system, associated to the fluid-elastic structure coupled problem, can be obtained from Equation ( 16) by choosing a zero elastic-electric coupling matrix C u (i.e. when piezoelectric coupling constants are zero)

K u -C up 0 K p U P = 2 M 0 C T up M p U P (17)

APPLICATION TO AXISYMMETRIC COMPOSITE CONICAL SHELLS FILLED WITH COMPRESSIBLE FLUID

In this section, the proposed finite element formulation is applied to the vibration analysis of axisymmetric conical elastic/piezoelectric shells coupled with an internal compressible fluid.

Elastic laminated conical shell

Consider a hollow truncated conical shell of length L, slope angle , having an arbitrary constant thickness h and made of N orthotropic layers (see Figures 2 and3). The particular geometry of the structure is used to reduce the dimensions of the problem through a semi-analytical procedure whereby variables are expanded as Fourier series in the circumferential direction. The orthogonality property of harmonic functions decouples the contributions of the components of the series and allows a separate analysis for each circumferential harmonic.

R(x) α L R 2 θ h e ,

Strain-displacement relations. Since a classical shell element based on Kirchhoff-Love

hypotheses is used in this work, the components of the displacement vector can be written as

u x = u(x, ) + z (x, ) (18a) u = v(x, ) + z (x, ) (18b)
u z = w(x, ) (18c)
u, v and w being the axial, circumferential and radial displacement components at the middle surface of the shell. and are the rotations about the and x axes, respectively, given by

= - *w *x (19)
= cos R v - 1 R *w * ( 20
)
where R denotes the local radius of the shell written by

R = x sin (21)
With these notations, the state of strain at any point x, (on the shell) and z (distance from the neutral surface), is equal to the sum of the middle surface strains (membrane) and the strains due to the changes of curvature (flexural). Denoting the middle surface extensional strains in the lengthwise and circumferential directions by x x and , the middle surface changes in curvature in the lengthwise and circumferential directions by x x and , and the middle surface in-plane shear strain and change in twist by x and x , the state of strain at any point (x, , z) within the shell can be written as

ε x x = x x + z x x (22a) ε = + z (22b) 2ε x = x + z x (22c)
where the factor 2 in the left-hand side of Equation (22c) comes from the definition of the strain tensor (Equation ( 6)) and the middle surface extensional strains are defined by

x x = *u *x (23a) = sin R u + 1 R *v * + cos R w (23b) x = 1 R *u * - sin R v + *v *x (23c)
and the middle surface curvatures are expressed by

x x = - * 2 w *x 2 (24a) = cos R 2 *v * - 1 R 2 * 2 w * 2 - sin R *w *x (24b) x = - 2 cos sin R 2 v + cos R *v *x + 2 sin R 2 *w * - 2 R * 2 w *x* (24c)
Equations ( 23) and ( 24) correspond to strain-displacement relations for thin conical shells in Love theory [START_REF] Soedel | Vibrations of Shells and Plates (3rd edn)[END_REF].

Fourier series expansion.

The geometric and material rotational symmetry of the structure allows a semi-analytical formulation for the shell finite element. Thus, the displacement field is expended in Fourier series in the circumferential direction according to the following expressions:

u = u 0 (x) + ∞ n=1 u n (x) cos n (25a) v = v 0 (x) + ∞ n=1 v n (x) sin n (25b) w = w 0 (x) + ∞ n=1 w n (x) cos n (25c)
where n identifies the harmonic component.

In fact, due to the axisymmetry, the modes for n 1 appear by pair, the two members of a pair being denoted as companion modes. Those companion modes have exactly the same frequency and their shapes have the same geometry, the only difference being in their angular position. Moreover, displacement components appear grouped in (u, *v/* , w) or in (*u/* , v, *w/* ) in Equations ( 23) and [START_REF] Zhang | Vibration analysis of cross-ply laminated composite cylindrical shells using the wave propagation approach[END_REF]. As a consequence, the displacement components of one companion mode are of the form (cos n , sin n , cos n ) or (sin n , cos n , sin n ). To obtain both companion modes, one would have to introduce both sine and cosine functions in each displacement component (u, v, w). However, for sake of brevity, only one of the companion modes is considered in Equations ( 25) and consequently in our simulations.

Using Equations ( 23)-( 25), middle surface extensional strains and curvatures can be written as

x x = 0 x x (x) + ∞ n=1 n x x (x) cos n , x x = 0 x x (x) + ∞ n=1 n x x (x) cos n = 0 (x) + ∞ n=1 n (x) cos n , = 0 (x) + ∞ n=1 n (x) cos n x = 0 x (x) + ∞ n=1 n x (x) sin n , x = 0 x (x) + ∞ n=1 n x (x) sin n
where the strain components can be rewritten, for any circumferential harmonic n, in the following matrix form:

e n = D m ⎡ ⎢ ⎣ u n v n w n ⎤ ⎥ ⎦ , j n = D b ⎡ ⎢ ⎣ u n v n w n ⎤ ⎥ ⎦ (26)
with

e n = [ n x x n n x ] T , j n = [ n x x n n
x ] T and the gradient operators

D m = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ * *x 0 0 S R n R C R - n R - S R + * *x 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , D b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 - * 2 *x 2 0 nC R 2 n 2 R 2 - S R * *x 0 -2 C S R 2 + C R * *x -2 nS R 2 + 2 n R * *x ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
where C = cos , S = sin and the subscripts m and b state for membrane and bending.

Degrees-of-freedom and shape functions.

The displacement components associated with each harmonic n are predicted by discretizing the shell with one-dimensional two nodes axisymmetric elements. Each node has four degrees-of-freedom to describe axial, circumferential and radial displacements, as well as the rotation about the circumferential axis. One element can thus be described in terms of the following nodal degrees-of-freedom vector:

U e n = [u n1 v n1 w n1 n1 | u n2 v n2 w n2 n2 ] T ( 27
)
where superscript e denotes an elementary quantity, and subscripts 1 and 2 correspond to nodes bounding the element. The longitudinal and circumferential displacements are assumed to vary linearly along the axial co-ordinate x, while the radial displacement is described by Hermite cubic polynomial in order to maintain C1 continuity due to the application of Kirchhoff-Love hypothesis. Thus, the elementary middle surface displacements and rotations of the shell are given in terms of the nodal degrees-of-freedom by

⎡ ⎢ ⎣ u e n v e n w e n ⎤ ⎥ ⎦ = N d (x)U e n , ⎡ ⎢ ⎣ e n e n 0 ⎤ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ - *w e n *x C R v e n + n R w e n 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = N r (x)U e n ( 28
)
where the interpolation matrices are defined by

N d = ⎡ ⎢ ⎣ N l 1 0 0 0 N l 2 0 0 0 0 N l 1 0 0 0 N l 2 0 0 0 0 N c 1 N c 2 0 0 N c 3 N c 4 ⎤ ⎥ ⎦
and

N r = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 - *N c 1 *x - *N c 2 *x 0 0 - *N c 3 *x - *N c 4 *x 0 C R N l 1 n R N c 1 n R N c 2 0 C R N l 2 n R N c 3 n R N c 4 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
with the following linear N l i (i = 1, 2) and cubic N c i (i = 1 . . . 4) shape functions given in the reference element

N l 1 = 1 2 (1 -), N l 2 = 1 2 (1 + ) N c 1 = 1 4 (2 -3 + 3 ), N c 2 = - L e 8 (1 --2 + 3 ) N c 3 = 1 4 (2 + 3 -3 ), N c 4 = - L e 8 (-1 -+ 2 + 3 )
in which is the elementary co-ordinate that varies from = -1 at node 1 to = 1 at node 2 and L e is the length of the element. Moreover, elementary extensional strains and curvatures can be expressed by

e e n = B m (x)U e n , j e n = B b (x)U e n ( 29
)
with the following discretized gradient operators:

B m = D m N d , B b = D b N d
Note that the third line of zeros in the N r expression is used in order to define the products N T d N r and N T r N d in the mass matrix defined below.

Elementary mass and stiffness matrices.

The interpolations of displacements and strains presented in the previous section are used to express the elementary mass and stiffness matrices of the composite laminated shell. The mass matrix is evaluated without neglecting rotational inertia. Using relation (13b) and combining Equations ( 18), ( 25) and (28), the elementary mass matrix is defined, for any harmonic n, by

M e u = (1 + 0n ) 1 -1 [I 0 N T d N d + I 1 (N T d N r + N T r N d ) + I 2 N T r N r ]R( )J S d
where J S is the Jacobian determinant which is simply given by J S = L e /2, is the Kronecker symbol and the zero, first and second mass moments of inertia are given by

I 0 = N k=1 k S (z k -z k-1 ), I 1 = 1 2 N k=1 k S (z 2 k -z 2 k-1 ), I 2 = 1 3 N k=1 k S (z 3 k -z 3 k-1 )
In the same way, using definition (13a) and combining Equations ( 22), ( 26) and (29), the elementary elastic stiffness matrix is given, for any harmonic n, by

K e u = (1 + 0n ) 1 -1 (B T m AB m + B T m BB b + B T b BB m + B T b DB b )R( )J S d
The matrices A, D and B are extensional, bending and extensional-bending coupling stiffness of the laminated composite (see for example [START_REF] Reddy | Mechanics of Laminated Composite Plates: Theory and Analysis[END_REF]). The components of these 3 × 3 symmetric matrices are given by

A i j = N k=1 Q k i j (z k -z k-1 ), B i j = 1 2 N k=1 Q k i j (z 2 k -z 2 k-1 ), D i j = 1 3 N k=1 Q k i j (z 3 k -z 3 k-1 )
where i, j = 1, 2, 6 and Q k i j represent the reduced material stiffness constants for each layer k in the global co-ordinate system deduced from the assumption of zero normal stress in the thickness direction.

Finite element modelling of the fluid domain

The aim of this section is to introduce the finite element discretization of the fluid domain considering the axisymmetry of the problem (see Figure 4). As for the structure displacements, fluid pressure can be expressed in terms of Fourier series by

p(r, , Z ) = p 0 (r, Z ) + ∞ n=1 p n (r, Z ) cos n ( 30
)
where p 0 and p n are, respectively, the axisymmetric and the nth harmonic component of the pressure.

The considered fluid domain is discretized into quadrilateral elements. Each element is bounded by four nodes and the nodal pressure is considered as the only unknown variable.

For any harmonic n, the elementary pressure p e n (r, Z ) can be expressed in terms of the vector of nodal pressures P e n p e n (r,

Z ) = N p (r, Z )P e n ( 31
)
where the interpolation matrix is defined by with the following linear shape functions, given in the reference element

N p = [N 1 N 2 N 3 N 4 ]
N 1 = 1 4 (1 -)(1 -), N 2 = 1 4 (1 + )(1 -), N 3 = 1 4 (1 + )(1 + ), N 4 = 1 4 (1 -)(1 + )
From definition (13g), and using Equations ( 30) and (31), the elementary mass matrix of the fluid can be expressed, for each harmonic n, by

M e p = (1 + 0n ) F c 2 F 1 -1 1 -1 N T p N p r ( , )J F d d ( 32
)
where J F is the Jacobian determinant of the transformation from the real to the reference element.

In order to compute the elementary stiffness matrix of the fluid, we have to evaluate the pressure gradient, for each harmonic n, using

∇p e n (r, Z ) = B p (r, Z )P e n ( 33
)
with the following discretized gradient operators:

B p (r, Z ) = * *r n r * *Z T N p (r, Z ) (34)
From Equations (13f), ( 30) and (33), and after transforming the discretized gradient operator in the reference co-ordinate system, the elementary stiffness matrix of the fluid can be expressed, for each harmonic n, by

K e p = (1 + 0n ) F 1 -1 1 -1 B T p B p r ( , )J F d d (35)

Elementary fluid-structure coupling matrix

Using the appropriate interface conditions at the common boundary between the fluid and the structure, the coupled fluid-structure matrix is derived. In first approach, a compatible mesh is considered at the interface so that the axisymmetric interface element can be easily defined. The interface element is bounded by two nodes, each having two degrees-of-freedom to describe, for any harmonic n, the normal structure displacement w e n and the fluid pressure p e n . The elementary normal displacement w e n is discretized by cubic shape functions and can be written, in terms of the nodal structure degrees-of-freedom, by

w e n = [0 0 N c 1 N c 2 0 0 N c 3 N c 4 ]U e n = N w U e n ( 36
)
The elementary pressure p e n is discretized by linear shape functions and can be written, in terms of fluid pressure in nodes located at the interface, by

p e n = [N l 1 N l 2 ] p n1 p n2 = N pi p n1 p n2 (37)
Thus, the elementary 8 × 2 coupling matrix can be written as

C e up = (1 + 0n ) 1 -1 N T w N pi R( )J S d (38)

Extension to piezoelectric conical shell

In this section, we present an extension of the laminated composite conical shell to the case of piezoelectric layers with radial polarization. The chosen multilayer model combines an equivalent single-layer assumption for the mechanical displacement and a layerwise representation of the transverse electric potential (the electric potential is assumed to vary linearly in the thickness of each piezoelectric layer). The advantages of this mixed laminate theory are linked to its effectiveness to model thin composite shells and to capture the through-thickness electric heterogeneity induced by the piezoelectric layers [START_REF] Saravanos | Mixed laminate theory and finite element for smart piezoelectric composite shell structures[END_REF][START_REF] Lammering | Multi-field variational formulations and related finite elements for piezoelectric shells[END_REF].

From the previous electric potential assumption and neglecting the in-plane components, the electric field is defined, for each piezoelectric layer k, by its transverse component

E k = -V k / h k (39)
where V k is the electrical potential difference between top and bottom surfaces (

V k = k -k-1
) and h k is the layer thickness (h k = z kz k-1 ). Instead of describing the electric state by the electric potentials k , we use in the finite element discretization the electric potential differences V k . Moreover, as for the displacement and pressure fields, the potential difference is developed in Fourier series by

V k (x, ) = V k 0 (x) + ∞ n=1 V k n (x) cos n (40)
Thus, the degrees-of-freedom of the axisymmetric shell element described in the previous section are augmented by the electric potential difference (voltage) of each layer. For a circumferential wave number n, the electric field of one multilayer piezoelectric shell element can be written in the following form:

E e n = -B U e n ( 41
)
with

E e n = [E 1 n . . . E k n . . . E N n ] T , U e n = [V 1 n . . . V k n . . . V N n ]
T , and

B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 h 1 0 • • • • • • 0 0 . . . 0 . . . . . . . . . 0 1 h k 0 . . . . . . . . . 0 . . . 0 0 • • • • • • 0 1 h N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ( 42
)
where N is the total number of layers. With these considerations, the elementary electromechanical coupling stiffness matrix, given by Equation (13c), can be written, for each circumferential wave number n, by

C e u = (1 + 0n ) 1 -1 (B T m FB + B T b GB )R( )J S d (43)
F and G being the membrane-electric and bending-electric coupling matrices, defined by

F = ⎡ ⎢ ⎢ ⎣ h 1 e 1 31 • • • h k e k 31 • • • h N e N 31 h 1 e 1 32 • • • h k e k 32 • • • h N e N 32 h 1 e 1 36 • • • h k e k 36 • • • h N e N 36 ⎤ ⎥ ⎥ ⎦ (44)
and

G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ z 2 1 -z 2 0 2 e 1 31 • • • z 2 k -z 2 k-1 2 e k 31 • • • z 2 N -z 2 N -1 2 e N 31 z 2 1 -z 2 0 2 e 1 32 • • • z 2 k -z 2 k-1 2 e k 32 • • • z 2 N -z 2 N -1 2 e N 32 z 2 1 -z 2 0 2 e 1 36 • • • z 2 k -z 2 k-1 2 e k 36 • • • z 2 N -z 2 N -1 2 e N 36 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ( 45
)
where e k 3 j ( j = 1, 2, 6) are the reduced piezoelectric constants for each layer k in the global co-ordinate system deduced from the assumption of zero normal stress in the thickness direction.

The previously defined matrices are derived from the global constitutive equations of the laminate, which connect, for any harmonic n, the resultant in-plane forces N n and bending moments M n , to the middle surface extensional strains e n and curvatures j n , and electric field

E n N n M n = A B B D e n j n - F G E n (46)
Finally, from Equation (13h), the electric matrix K e is given, for each harmonic n, by

K e = (1 + 0n ) 1 -1 B T HB R( )J S d ( 47
)
where H is defined by

H = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ h 1 d 1 33 0 • • • • • • 0 0 . . . 0 . . . . . . . . . 0 h k d k 33 0 . . . . . . . . . 0 . . . 0 0 • • • • • • 0 h N d N 33 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (48)
with d k 33 the dielectric permittivity in the thickness direction. The matrix H is derived by integrating the transverse electric displacement D k over the thickness of each layer such that

Q n = -F T e n -G T j n -HE n ( 49
)
where the electric displacement resultant vector Q n is defined by

Q n = z 1 z 0 D 1 dz • • • z k z k-1 D k dz • • • z N z N -1 D N dz T (50)
and where D k is given for each layer in the global co-ordinate system by

D k = e k 31 ε x x + e k 32 ε + 2e k 36 ε x + d k 33 E k (51)

NUMERICAL RESULTS

In this last section, numerical results, obtained with a Matlab program developed by the authors, are proposed in order to validate and analyse the previously described formulation. Section 4.1 presents the performances of the composite shell finite element for elastic or piezoelectric structures without internal fluid. Section 4.2 concerns elastic and piezoelectric cylindrical shells coupled with an internal compressible fluid. Table I presents a comparison of the non-dimensional frequency parameter = R 2 (1-2 )/E for various set of boundary conditions at both the base and top of the cone: free-free (F-F), simply supported-simply supported (SS-SS) v = w = 0, and clamped-clamped (C-C) u = v = w = = 0. Moreover, the results are compared with those given by Shu [START_REF] Shu | Free vibration analysis of composite laminated conical shells by generalized differential quadrature[END_REF] using a generalized differential quadrature technique, and those obtained by the finite element code Nastran using quadrangular shell elements. In this example, the circumferential wave number n ranges from 1 to 9. An excellent agreement can be observed for all modes. From this comparison, it can be concluded Table II. Mechanical and geometrical data of an elastic laminated cylinder [START_REF] Lam | Analysis of rotating laminated cylindrical shells by different thin shell theories[END_REF].

Material properties

Geometric data

E 1 = 19 GPa, E 1 /E 2 = 2.5 Inner layer thickness = h/3 G 12 = 4.1 GPa
Middle layer thickness = h/3 12 = 0.26

Outer layer thickness = h/3 = 1643 kg m -3 h/R = 0.002

Table III. Frequency parameter = R 2 /E 2 of a composite cylinder. that the present axisymmetric two-node shell element is numerically accurate and computationally efficient for thin isotropic conical shells.

L/R = 1 L/R =

Laminated composite cylinder.

This example gives a comparison of the frequency parameter = R √ /E 2 for a three-layer cross-ply [0, /2, 0] cylindrical shell with simply supported (SS-SS) boundary conditions. The mechanical and geometrical shell parameters are given in Table II. Moreover, 100 finite elements are used for the discretization of the cylinder.

The results presented in Table III show the excellent performance of our finite element model compared to the analytical solution proposed by Lam and Loy [START_REF] Lam | Analysis of rotating laminated cylindrical shells by different thin shell theories[END_REF] and the wave propagation approach used by Zhang [START_REF] Zhang | Vibration analysis of cross-ply laminated composite cylindrical shells using the wave propagation approach[END_REF]. It can be observed that for shells with small length-to-radius ratios, as in the case L/R = 1, the minimum frequency occurs for the larger n. Nevertheless, for L/R = 10, the natural frequencies decrease when n varies from 1 to 4, and increase for n larger than 4.

Laminated composite conical shell.

In order to validate our finite element model for composite non-cylindrical structures, the natural frequencies of a three-layer [0, /2, 0] conical shell are compared with those given by Nastran with quadrangular composite shell elements. The geometrical data, boundary conditions and spatial discretization are the same as those given for the isotropic conical shell example (see Section 4.1.1). Here, the thickness of each layer is h/3 and a graphite-epoxy material is used (see Table IV).

As in the previous case, a very good agreement between the frequencies can be observed in Table V, thus validating our composite conical shell element. It is also observed that the influence of boundary conditions on the variation of frequency with circumferential wave number n is significant. For the F-F conical shell, the frequencies increase as the circumferential wave number n increases. For the SS-SS and C-C conical shells, the frequencies decrease for circumferential wave number from n = 1 to 4 and increase from n = 5 to 9.

For illustration purpose, the mode shapes are presented in Figure 5. These three-dimensional deformations have been constructed from the one-dimensional element and using the Fourier series expansion of the displacement field in the circumferential direction.

Piezoelectric cylindrical shell.

In order to validate the formulation for piezoelectric shells, we present in this example a comparative study for a simply supported (SS-SS) piezoelectric cylindrical shell. The geometrical properties are L = 5 m, R = 1 m and h = 0.02 m. Moreover, the piezoelectric material is the PZT-5H (see Table IV). Table VI presents the frequencies calculated by our method and those given by an exact threedimensional solution proposed by Deü and Larbi [4]. This exact solution is based on a mixed state-space approach previously developed for the free vibration analysis of laminated piezoelectric plates actuated by transverse shear mechanisms [START_REF] Deü | Free-vibration analysis of laminated plates with embedded shear-mode piezoceramic layers[END_REF]. In this table, the modes are classified, for two values of circumferential harmonic (n), in order of increasing frequency. To specify the mode type, the corresponding axial wave number (m) obtained with the exact solution is also given. As it can be observed from this table, there is a good agreement between finite element (with 100 conical shell elements) and exact solutions for different electric boundary conditions corresponding to short-circuited ( = 0) or open-circuited (D = 0) piezoelectric materials. As expected, the natural frequencies are higher in the open-circuit case than in the closed-circuit one. It is noted that the error committed by the finite element approximation is lower than 3%. This comparison enables us to validate the piezoelectric aspect of the shell. It is important to note that this validation is restricted to the first modes, i.e. when a linear piezoelectric potential along the thickness of the piezoelectric layer is sufficient. Table VII gives the frequencies of (i) the fluid in rigid cavity, (ii) the empty shell, and (iii) the fluid-structure-coupled system. The structural frequencies in the coupled and uncoupled cases are compared to those obtained with a wave propagation method [START_REF] Zhang | Coupled vibration analysis of fluid-filled cylindrical shells using the wave propagation approach[END_REF]. This comparison enables us to check the validity of the fluid-structure proposed formulation. It can be observed that the coupled frequencies correspond to half of the empty shell frequencies. This is due to the added mass effect of the fluid which is, in this particular case of a liquid, quasi-incompressible for the first low-frequency modes. In order to illustrate the results in the coupled case, the mode shapes, in terms of fluid pressure and structure displacement, are presented in Figure 6.

Elastic case: comparison with an exact state-space solution.

In our second fluid-structure example, we consider a simply supported isotropic cylindrical shell containing a compressible fluid. The shell has the following geometrical and material properties: length L = 20 m, radius R = 1 m, thickness h = 0.002 m, Young's modulus E = 2.07788 GPa, Poisson ratio = 0.317756 and mass density S = 8166 kg m -3 . The considered fluid and the finite element discretization are the same as in the preceding example. In the case where the fluid is not taken into account, Table VIII shows a comparison, for n varying from 1 to 10, between our finite element results, an exact state-space solution based on the three-dimensional theory of elasticity [START_REF] Deü | A state space method for free-vibration analysis of radially polarized laminated piezoelectric cylinder filled with fluid[END_REF], and the results given by Loy et al. [START_REF] Loy | Vibration of functionally graded cylindrical shells[END_REF] using a Love shell theory and a Ritz method. A very good agreement between the three approaches can be observed. Table IX presents the natural frequencies for the fluid-structure-coupled problem. As can be seen in this table, the coupled frequencies obtained with our finite element approach agree very well with those obtained from the exact solution [START_REF] Larbi | Solution tridimensionnelle pour l'analyse des vibrations d'une structure cylindrique avec fluide interne[END_REF] for different values of circumferential (n = 1-4) and axial (m = 1-5) harmonics. It is also observed in Figure 7(a) that for all axial wave number m, the frequencies decrease as the circumferential wave number n increases except for m = 1. Moreover, Figure 7(b) shows, for all circumferential wave number n, that the frequencies increase as the axial wave number m increases. 4.2.3. Piezoelectric case: comparison with a state-space exact solution. In this last example, the piezoelectric cylinder studied in Section 4.1.4 is filled with a compressible fluid characterized by F = 1000 kg m -3 and c F = 1500 m s -1 . Table X presents the eigenfrequencies of the coupled system computed by our finite element method and those given by the exact state-space approach [START_REF] Deü | A state space method for free-vibration analysis of radially polarized laminated piezoelectric cylinder filled with fluid[END_REF]. The results, given for the two first circumferential harmonics, show a good agreement between the two methods. However, as previously mentioned, the difference between exact and finite element solutions can increase for higher modes, in particular due to the linearity assumption of the electric potential.

CONCLUSIONS

This paper concerns variational formulation and the finite element implementation of the fluid/ piezoelectric structure vibration problems. The coupled system consists of a piezoelectric structure (described by its field and its electric potential) containing a compressible fluid (described by its pressure field). With this description, the variational formulation of the coupled problem is established and the corresponding matrix equations are presented. This formulation is then applied to the vibration of axisymmetric conical elastic/piezoelectric composite shells coupled with an internal fluid. The piezoelectric shell is based on Kirchhoff-Love theory and is supposed to be polarized in thickness direction. The associated finite element formulation combines an equivalent single-layer approach for the mechanical behaviour with a layerwise representation of the electric potential in the thickness direction. The fluid is modelled by the Helmholtz equation in terms of pressure. The particular geometry of the structure is used to reduce the dimensions of the problem through a semi-analytical procedure whereby variables are expanded as Fourier series. Following these considerations, the shell is discretized by a onedimensional element with two nodes, and the fluid by a quadrilateral axisymmetric element with four nodes. Moreover, an appropriate interface element is developed to ensure the interaction between fluid and structure. Finally, numerical examples are presented validating the development of the composite axisymmetric shell element, the piezoelectric aspect and the fluid-structure finite element implementation. In this paper, we have used a non-symmetric fluid/piezoelectric structure formulation. The symmetrization can be carried out using procedures as described in [START_REF] Morand | Fluid-Structure Interaction[END_REF] and will be the purpose of a further publication.

Figure 1 .

 1 Figure 1. Fluid/piezoelectric structure-coupled system.

S

  ce(u) : e( u) dv -S e T E() : e( u) dvpn • u dsvariational formulation of the coupled fluid-piezoelectric structure spectral problem consists in finding ∈ R + and (u, , p) ∈ (C * u , C * , C p ) satisfying Equations (10)-(12).

Figure 2 .

 2 Figure

Figure 3 .

 3 Figure 3. Cross-sectional view of the laminated composite conical shell.

Figure 4 .

 4 Figure 4. Fluid domain: geometry and notations.

4. 1 .

 1 Free vibration of shells without internal fluid 4.1.1. Isotropic conical shell. Consider an isotropic single-layer conical shell of length L = 7.5 m, slope angle = /6, thickness h = 0.05 m, and radius at its large edge R 2 = 5 m. The isotropic material used in this first example is aluminium, with Young's modulus E = 69 000 MPa, Poisson ratio = 0.3, and mass density = 2700 kg m -3 . Moreover, the shell is discretized with 50 axisymmetric elements.

Figure 5 .

 5 Figure 5. The five first mode shapes of composite conical shell: (a) free-free; (b) simply supported-simply supported; and (c) clamped-clamped.

4. 2 .

 2 Free vibration of cylindrical shells filled with a compressible fluid 4.2.1. Elastic case: comparison with a wave propagation approach. In order to validate the proposed formulation in a fluid-structure case, we consider an elastic isotropic cylinder filled with compressible fluid. The shell is clamped at both ends and has the following properties: length L = 20 m, radius R = 1 m, thickness h = 0.01 m, Young's modulus E = 210 GPa, Poisson ratio = 0.3, and mass density S = 7850 kg m -3 . The considered fluid is a liquid with a mass density F = 1000 kg m -3 and a speed of sound c F = 1500 m s -1 . Moreover, 40 conical shell elements and 200 quadrangular axisymmetric fluid elements are used in the simulations.

Figure 6 .

 6 Figure 6. Fluid-structure-coupled modes: fluid pressure and shell total displacement.

Figure 7 .

 7 Figure 7. Variation of frequency with: (a) circumferential wave number n and (b) axial wave number m.

Table I .

 I Frequency parameter = R 2 (1 -2 )/E of an isotropic conical shell for n = 1-9.

		F-F			SS-SS			C-C
	n	Nastran	Present	Nastran	Present	Shu [22]	Nastran	Present	Shu [22]
	1	-	-	0.679	0.679	0.695	0.752	0.753	0.752
	2	0.014	0.015	0.358	0.358	0.359	0.466	0.466	0.466
	3	0.037	0.037	0.210	0.211	0.212	0.315	0.316	0.315
	4	0.065	0.066	0.190	0.190	0.181	0.260	0.260	0.261
	5	0.096	0.096	0.216	0.216	0.216	0.264	0.265	0.265
	6	0.132	0.132	0.252	0.253	0.254	0.293	0.294	0.293
	7	0.173	0.174	0.298	0.298	0.298	0.334	0.335	0.335
	8	0.221	0.222	0.351	0.352	0.353	0.385	0.387	0.387
	9	0.275	0.276	0.413	0.415	0.414	0.446	0.448	0.448

Table IV .

 IV Properties of graphite-epoxy and PZT-5H materials.

	Properties	Graphite-epoxy	PZT-5H
	c 11 (GPa)	183.443	126
	c 12 (GPa)	4.363	79.5
	c 13 (GPa)	4.363	84.1
	c 22 (GPa)	11.662	126
	c 23 (GPa)	3.918	84.1
	c 33 (GPa)	11.662	117
	c 44 (GPa)	2.87	23
	c 55 (GPa)	7.17	23
	c 66 (GPa) e 15 (C m -2 ) e 24 (C m -2 ) e 31 (C m -2 ) e 32 (C m -2 ) e 33 (C m -2 ) d 11 (10 -10 F m -1 ) d 22 (10 -10 F m -1 ) d 33 (10 -10 F m -1 ) (kg m -3 )	7.17 0 0 0 0 0 153 153 153 1590	23.3 1 7 1 7 -6.5 -6.5 23.3 150.3 150.3 130 7500

Table V .

 V Frequencies (Hz) of a laminated [0, /2, 0] conical shell.

		F-F		SS-SS		C-C
	n	Nastran	Present	Nastran	Present	Nastran	Present
	1	-	-	93.531	93.562	95.536	95.636
	2	1.485	1.538	53.452	53.456	59.639	59.786
	3	3.936	3.970	35.277	35.291	42.951	43.090
	4	7.175	7.199	28.422	28.451	35.612	35.721
	5	10.870	10.892	28.519	28.560	34.345	34.438
	6	14.852	14.880	31.610	31.667	36.758	36.880
	7	19.292	19.332	35.995	36.085	41.151	41.337
	8	24.325	24.386	41.428	41.570	46.883	47.156
	9	29.988	30.079	47.822	48.036	53.701	54.089
	10	36.289	36.421	55.094	55.405	61.453	61.985

Table VI .

 VI Frequencies (Hz) of a piezoelectric cylindrical shell.

				Short circuited			Open circuited
	n	m	Present	Exact	Error %	Present	Exact	Error %
	1	1	87.806	88.231	0.482	93.370	93.790	0.448
		2	216.127	216.260	0.062	228.162		0.028
		0	278.711	275.956	-0.998	278.711	275.956	-0.998
		3	312.564	311.647	-0.294	331.697	330.649	-0.317
		4	368.995	366.968	-0.553	383.703	392.134	2.150
	2	1	36.705	37.316	1.639	39.489	40.109	1.546
		2	113.540	114.210	0.587	121.979	122.677	0.569
		3	193.357	193.791	0.224	207.472	207.916	0.213
		4	259.388	259.346	-0.016	278.508	278.434	-0.027
		5	308.942	308.464	-0.155	332.121	331.551	-0.172

Table VII .

 VII Frequencies of a cylindrical shell coupled with a compressible fluid.

	m, n	Rigid cavity	Empty shell	Empty shell [1]	Coupled pb	Coupled pb [1]
	1,2	651.04	12.17	12.17	4.77	4.93
	1,3	822.10	19.62	19.61	8.53	8.94
	2,3	953.94	23.45	23.28	10.24	10.64
	2,2	701.41	27.72	28.06	10.93	11.48
	3,3	955.90	32.46	31.98	14.25	14.66
	1,4	934.42	36.44	36.47	17.21	18.26
	2,4	1192.97	37.42	37.37	17.72	18.73
	3,4	1194.38	40.09	39.78	19.08	19.96

Table VIII

 VIII

	1	13.587	13.561	13.548
	2	4.637	4.595	4.592
	3	4.288	4.257	4.263
	4	7.234	7.211	7.225
	5	11.546	11.519	11.542
	6	16.899	16.863	16.897
	7	23.245	23.197	23.244
	8	30.574	30.511	30.573
	9	38.882	38.802	38.881
	10	48.168	48.069	48.168

. Frequencies (Hz) of an empty cylindrical shell for n varying from 1 to 10.

n FE Exact [4]

Loy et al.

[START_REF] Loy | Vibration of functionally graded cylindrical shells[END_REF]

Table IX .

 IX Frequencies (Hz) of an elastic cylindrical shell filled with a compressible fluid: comparison with an exact state-space solution[START_REF] Deü | A state space method for free-vibration analysis of radially polarized laminated piezoelectric cylinder filled with fluid[END_REF].

		n = 1	n = 2	n = 3	n = 4
	m	FE	Exact	FE	Exact	FE	Exact	FE	Exact
	1	2.408	2.441	0.899	0.912	0.947	0.968	1.779	1.839
	2	9.073	9.179	3.409	3.453	2.040	2.068	2.142	2.193
	3	18.633	18.817	7.484	7.559	4.243	4.280	3.254	3.290
	4	29.547	29.849	12.902	12.987	7.357	7.391	5.144	5.155
	5	40.536	41.147	19.405	19.470	11.283	11.283	7.701	7.666

Table X .

 X Frequencies (Hz) of a piezoelectric cylindrical shell filled with a compressible fluid: comparison with an exact solution[START_REF] Deü | A state space method for free-vibration analysis of radially polarized laminated piezoelectric cylinder filled with fluid[END_REF].

				Short circuited			Open circuited
	n	m	Present	Exact	Error %	Present	Exact	Error %
	1	1	45.445	44.760	-1.531	49.491	48.271	-2.528
		2	108.920	107.390	-1.424	115.311	115.614	0.263
		3	154.565	154.131	-0.281	166.237	166.065	-0.103
		4	191.016	189.327	-0.892	209.635	204.211	-2.656
		5	220.018	217.190	-1.302	236.253	234.469	-0.761
	2	1	19.637	19.916	1.404	21.001	21.464	2.157
		2	62.497	62.597	0.159	66.928	67.623	1.027
		3	108.109	108.575	0.430	116.416	117.248	0.710
		4	148.610	149.022	0.276	161.251	160.912	-0.211
		5	182.892	182.817	-0.041	197.330	197.400	0.036

ACKNOWLEDGEMENTS

This research is partially carried out within the 6th Framework Programme of the European Commission (InMAR project) under Grant No. NMP2-CT-2003-501084.