domain: spi info

text:
Oumar Cissoko
email: oumarbcissoko@gmail.com

Bastien Béchadergue

Iman Tavakkolnia

JWOC 2023

Keywords: Optical wireless communications, energy harvesting, photovoltaic module

In this article, we present the simulation model of a photovoltaic module (PVM) used in an optical wireless communication system for both data reception and energy harvesting. A light-emitting diode placed at the center of a 4×4×2.85 m room transmits a 10 kbps Manchester-encoded onoff keying data stream, that is received and decoded by a PVM moving across a reception plane 2 m from the source. We obtained a null bit error rate whatever the PVM position, and a maximum harvested power of 0.42 mW.

I. INTRODUCTION

Communications needs are growing at an extremely rapid pace, and so is the deployment of new wireless technologies to meet them, including radio frequency, future terahertz and optical wireless communications (OWC) [START_REF]Block diagram of the simulated system with the equivalent electrical model of the PVM[END_REF]. In a context of strong environmental concerns [START_REF] Andrae | On Global Electricity Usage of Communication Technology: Trends to 2030[END_REF], such technologies must be sustainable, especially in terms of energy efficiency and power consumption. One solution addressing the later issue is to collect the energy carried by the data signals received by communication systems and the ambient energy: this is energy harvesting [START_REF] Sezer | A comprehensive review on the state-of-the-art of piezoelectric energy harvesting[END_REF][START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF][START_REF] Rahul | Experimental Validation of Optical Wireless Receiver using Solar Panel with Bandwidth Enhancement Circuit[END_REF][START_REF] Tavakkolnia | Organic photovoltaics for simultaneous energy harvesting and high-speed MIMO optical wireless communications[END_REF].

In the optical field, such a solution can be implemented using photovoltaic modules (PVM), designed to convert optical energy into electrical energy, rather than conventional photodiodes, although the latter may provide higher data rate as they usually have larger bandwidth [START_REF]Block diagram of the simulated system with the equivalent electrical model of the PVM[END_REF]. For example, in [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF], the use of a PVM is demonstrated experimentally and a data rate of 11.84 Mbps is achieved simultaneously with a harvesting power of 30 mW when the distance between the emitter and the PVM is 0.75 m. However, these results were obtained at a given distance without considering other crucial metrics like the communication coverage.

In this paper, we simulate with MATLAB an end-to-end OWC system using a state-of-the-art PVM for both energy harvesting and optical signal reception, and we study its performance in terms of data rate, energy harvesting and communication coverage. To that extent, we first detail the system and PVM models (Section II) and then the simulation results obtained (Section III), before concluding the paper (Section IV) with final remarks and future works.

II. STUDIED SYSTEM ARCHITECTURE AND MODELLING

Figure 1 shows the architecture of the simulated system. The data transmitted in the form of light radiation from the light emitting diode (LED) propagates through the OWC channel before being collected by the PVM. This data is transmitted using on-off keying (OOK) modulation with Manchester coding (i.e. a bit '1' corresponds to the symbol '10' and a bit '0' to '01'). The PVM, on its side, is composed of a photovoltaic panel along with a receiver circuit that recovers the data and the electrical power they carry.

To simulate the communication performance of such a system, it is necessary to model its frequency response, which links the optical signal emitted ( 𝜙 , in [W]) to the load terminal voltage of the PVM (𝑣(𝑓), in [V]), and is given by:

 𝑣(𝑓) = 𝜙 𝐻(0)ℛ𝐻 (𝑓), 
where 𝐻(0) (without units) models the OWC channel gain, ℛ (in [A/W]) is the PVM responsivity and 𝐻 (𝑓) (in [V/A]) models the PVM frequency response.

Then the energy harvesting performance of the system, i.e. the electrical power 𝑃 collected on the terminals of the load 𝑅 of the harvesting branch, can be evaluated by calculating the DC current across 𝑅 while assuming the AC source of the circuit is off (as the AC part is strongly attenuated by the energy harvesting branch). This leads to:

 𝑃 = ( ) ℛ ( )
,  with 𝑟, 𝑅 , 𝑅 and 𝑅 as defined in Table I.

In practice, we can model the OWC channel using Lambert's emission law, which relates the optical power 𝜙 emitted by an LED to the optical power 𝜙 incident to a photoreceiver of area 𝑆, and is given by [START_REF]Block diagram of the simulated system with the equivalent electrical model of the PVM[END_REF]:

 𝐻(0) = = (
) cos (𝜙) cos(𝜓),  with 𝑚 the Lambertian order, defined as 𝑚 = -ln 2 / ln(cos 𝜃 / ), 𝜃 / the semi-angle at half power of the LED source, 𝜙 the radiation angle, 𝜓 the incidence angle and 𝑑 the distance between the LED and the PVM.

In parallel, we can model the PVM frequency response

𝐻 (𝑓) as 𝐻 (𝑓) = ( ) ( )
, where 𝑖 (𝑓) is the AC part of the photocurrent generated by the PVM from the incident optical power 𝜙 . Although the literature lists three main electrical models for PVMs -with one, two or three diodes -, our study is restricted to the use of the one-diode model as it remains accurate [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF] and requires less simulation time when estimating parameters or evaluating the PVM performance. In Fig. 1, the PVM is thus represented following the equivalent one-diode model proposed in [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF], which includes multiple passive dipoles modelling each specific physical behaviors, as detailed in Table I. According to this model:

 𝐻 = ,  with 𝑅 = / an equivalent impedance.
Finally, when converting 𝜙 into 𝑣(𝑓), the PVM also adds shot and thermal noise, which expressions are not detailed here but were simulated following the same model as in [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF].

III. SYSTEM PERFORMANCE

Table I lists the parameters used in our simulations, defined as in [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF]. With these parameters, the PVM frequency response has a flat -3 dB bandwidth from around 350 Hz to 110 kHz, which strongly limits the achievable data rate, unless bandwidth enhancement or equalization techniques are used [START_REF] Rahul | Experimental Validation of Optical Wireless Receiver using Solar Panel with Bandwidth Enhancement Circuit[END_REF]. Therefore, we set the OOK data rate to 10 kbps for the rest of the simulations, which consist in evaluating the spatial distribution of the bit error rate (BER) and energy harvested when the LED is at the centre of the ceiling of a 4×4×2.85 m room, and the PVM moves by steps of 20 cm along the 𝑥 and 𝑦 axis of a reception plane 2 m form the source. Figure 2 shows on the one hand an example of transmitted optical signal (in blue) and the resulting voltage signal output by the PVM (orange). We observe the output signal has a very low amplitude because of the OWC channel path loss, and that its shape is not square anymore because of the low-pass and high-pass effects of the PVM frequency response. However, the noise level is so low that the signal-to-noise ratio (SNR) remains high enough at any position of the PVM in the room (around 40 dB), so that the BER spatial distribution remains null (after transmitting 1000 bits) and the whole room is thus correctly covered. Such a high SNR also shows that more advanced modulations could be used to reach higher data rate. On the other hand, Fig. 3 shows the spatial distribution of the power harvested by the PVM, which is higher when the PVM is closer to the transmitter. This is due to the effect of the OWC channel gain (3), which depends on the distance between transmitter and receiver. The maximum value is eventually 0.42 mW at a distance of 2 m (at 0.75 m, it reaches 21.3 mW, which is line with the 30 mW measured in [START_REF] Wang | On the Design of a Solar-Panel Receiver for Optical Wireless Communications With Simultaneous Energy Harvesting[END_REF]), but it drops quickly as the PVM moves away from the optical axis of the emitter.

IV. CONCLUSIONS AND FUTURE WORKS

In this paper, we evaluated the communication coverage and energy harvesting performance of an OWC system using a PVM as receiver. We concluded that at 10 kbps, the communication is error-less across the whole room, but that the energy harvested drops quickly as the PVM moves away from the source. In future works, we will investigate the interest of beamsteering (i.e. real-time alignment between the source and PVM) to optimize the energy collected and validate these results through experimental measurements.

Fig. 2 .

 2 Fig. 2. Example of time evolution of the transmitted optical signal (blue) and voltage signal output by the PVM (orange).

Fig. 3 .

 3 Fig. 3. Spatial distribution of the electrical power harvested (in [mW]).

TABLE I

 I

	.	PARAMETERS VALUES USED.
	Description	Symbol	Value
	Serial resistor	𝑅	1 Ω
	Shunt resistor	𝑅	20 kΩ
	Small signal equivalent resistor	𝑟	720 Ω
	Branch resistance for energy conversion	𝑅	1 kΩ
	Resistor of the data reception branch	𝑅	50 Ω
	PVM capacitive effect	𝐶	34 nF
	Capacitance to block the DC component of the photo-generated current	𝐶	1 µF
	Serial inductance	𝐿	120 nH
	Inductance to attenuate the AC component of the photo-generated current	𝐿	10 mH
	PVM sensitive area	𝑆	432 cm²


Roger Ohayon
email: ohayon@cnam.fr

Reduced symmetric models for modal analysis of internal structural-acoustic and hydroelastic-sloshing systems

Reduced symmetric models for modal analysis of internal structural-acoustic and hydroelastic-sloshing systems

Introduction

This paper deals with appropriate computational methods for modal analysis of elastic structures containing an inviscid ¯uid (gas or liquid). These methods, based on Ritz±Galerkin projection using appropriate functional basis, allow us to construct reduced models expressed in terms of physical displacement vector ®eld u in the structure, and generalized coordinates vector r describing the behavior of the ¯uid. Those reduced models lead to symmetric generalized eigenvalue matrix system involving a reduced number of degrees of freedom for the ¯uid. More precisely, we construct symmetric matrix models of the ¯uid considered as a subsystem, by considering the response of the ¯uid to a prescribed normal displacement of the ¯uid±structure interface. Two distinct situations are analyzed, namely, linear vibrations of an elastic structure completely ®lled with a compressible gas or liquid and linear vibrations of an elastic structure containing an incompressible liquid with free surface eects due to gravity. The ®rst case is a structuralacoustic problem with modal interaction between structural modes in vacuo for structure containing a gas or incompressible hydroelastic modes for structure containing liquids, with acoustic modes in rigid motionless cavity. Wall impedance can also be easily introduced in order to take into account ¯uid±structure interface dissipation, for further forced response studies. The second case is a hydroelastic-sloshing problem with modal interaction between incompressible hydroelastic structural modes with incompressible liquid sloshing modes in rigid motionless cavity, involving an elastogravity operator elastogravity operator related to the wall normal displacement of the ¯uid±structure interface. For the construction of reduced models, the static behavior at zero frequency plays an important role. This is why we start from ``well-posed'' variational formulations of the problem, in the sense that zero-frequency behavior must be well retrieved in the equations. It should be noted that the so-called ``quasi-static correction'' term plays a fundamental role in the Ritz±Galerkin procedure (error truncation). The general methodology corresponds to dynamic 1 substructuring procedures adapted to ¯uid±structure modal analysis. For general presentations of computational methods using appropriate ®nite element and dynamic substructuring procedures applied to modal analysis of elastic structures containing inviscid ¯uids (sloshing, hydroelasticity and structuralacoustics), we refer the reader to [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF], and for the case of various damping effects, we refer the reader to [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF].

Structural acoustic problem

We consider here the harmonic vibrations of an elastic structure completely ®lled with a homogeneous, inviscid and compressible ¯uid considered as weightless (linear acoustic ¯uid). Among the applications, we ®nd, for instance, the important problem of acoustic comfort of ``passengers'' related to the vibroacoustic environment such as in launcher payloads, aircraft, automobiles. The particular case of a compressible liquid with a free surface, negecting gravity eects, is also considered (for applications involving compressibility eects in huge cryogenic tanks).

After an original derivation of the linearized equations of the ¯uid±structure coupled system, introducing an appropriate linear constraint in order to obtain a well-posed problem at zero frequency, we discuss the basic ideas leading to a reduced model of the ``¯uid subsystem''. For this purpose, acoustic modes in rigid motionless cavity are introduced as Ritz projection vector basis, including mandatory zerofrequency (static) solution of the coupled system. Let us recall that the internal ¯uid±structure system is characterized by a strong resonant behavior. Therefore, a direct ®nite element description can be costly. That is why, starting from one of the possible variational formulations of the problem, reduced symmetric matrix models are constructed.

Structural acoustic equations

Structure subjected to a ¯uid pressure loading

We consider an elastic structure occupying the domain X S at equilibrium. The internal ¯uid±structure interface is denoted as R. The angular frequency is denoted as x. The chosen unknown ®eld in the structure domain X S is the displacement ®eld u. The linearized deformation tensor is denoted as ij u and the corresponding stress tensor is denoted as r ij u. We denote by q S the constant mass density at equilibrium and by n the unit normal, external to the structure domain X S . Let du be the test function, associated to u, belonging to the admissible space C u .

The weak variational formulation describing the response of the structure X S to given harmonic forces of amplitude F d on the external structure boundary oX S nR, and to ¯uid pressure ®eld p acting on the internal ¯uid±structure interface R is written as follows. For all real x and Vdu P C u , ®nd u P C u such that ku

; duÀx 2 Z X S q S u du dx À Z R pn du dr Z oX S nR F d du
In Eqs. ( 4) and ( 5), k G u; du represents the classical structural prestress geometric symmetric bilinear form in which r 0 ij denotes the prestress tensor, and k P 0 u; du represents an additional prestress symmetric bilinear form due to rotation of normal n, in which P 0 denotes the initial pressure existing in the reference equilibrium con®guration. Finally, n 1 u represents the variation of normal n between the reference con-®guration and the actual con®guration.

Fluid subjected to a wall normal displacement

Since the ¯uid is inviscid, instead of describing the small motion of the ¯uid by a ¯uid displacement vector ®eld u F (see for instance [START_REF] Bermudez | Finite element computation of the vibration modes of a ¯uid±solid system[END_REF]), which requires an appropriate discretization of the ¯uid irrotationality constraint curl u F 0, we will use the pressure scalar ®eld p. It should be noted that the small movements corresponding to x T 0 are obviously irrotational, but in the static limit case, i.e. at zero frequency, we consider only ¯uids which exhibit a physical irrotational behavior.

Let us denote by c the (constant) sound speed in the ¯uid, and by q F , the (constant) mass density of the ¯uid at rest. We denote as X F the domain occupied by the ¯uid at rest (which is taken as the equilibrium state). The local equations describing the harmonic response of the ¯uid to a prescribed arbitrary normal displacement un of the ¯uid±structure interface R are such that

rp À q F x 2 u F 0j X F ; 5 p Àq F c 2 ru f j X F ; 6 u F n unj R ; 7 curl u F 0j X F : 8
Eq. ( 5) corresponds to the linearized Euler equation in the ¯uid. Eq. ( 6) corresponds to the constitutive equation of the barotropic ¯uid. Eq. ( 7) corresponds to the wall slipping condition. Eq. ( 8) corresponds to the irrotationality condition, only necessary in order to ensure that when x 3 0, u F tends to static irrotational motion, which corresponds to the hypothesis that for x 0, we only consider irrotational motions (for simply connected ¯uid domain).

Relation between static pressure p s and un

For x 0, Eqs. ( 6) and ( 7) lead to the constant pressure static ®eld which is related to the normal wall displacement by the relation (see [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF][START_REF] Ohayon | The eect of wall motion on the governing equations of structures containing ¯uids[END_REF], and for the dissipative case [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF])

p s À q F c 2 jX F j Z R undr; 9
in which jX F j denotes the measure of the volume occupied by domain X F .

Equations in terms of p and un

The elimination of u F between Eqs. ( 5)±(8) leads to

Dp x 2 c 2 p 0j X F ; 10 op on q F x 2 unj R ; 11 with the constraint 1 q F c 2 Z X F p dx Z R undr 0: 12
Eq. ( 10) is the classical Helmholtz equation expressed in terms of p. Eq. ( 11) corresponds to the kinematic condition de®ned by Eq. ( 7). The linear constraint de®ned by Eq. ( 12) corresponds to the global mass conservation which ensures that the boundary problem de®ned by Eqs. ( 10) and ( 11) is equivalent to the problem de®ned by Eqs. ( 5)± [START_REF] Ohayon | Symmetric formulations for modal analysis of internal ¯uid structure systems[END_REF]. In the absence of the condition de®ned by Eq. ( 12), we would obtain a boundary value problem in terms of p which is not valid for x 0 and which does not allow us to retrieve the value of p s given by Eq. ( 9). This condition is usually omitted in the literature.

Variational formulation in terms of u; p

Let dp be the test function, associated to p, belonging to the admissible space C p . The weak variational formulation corresponding to Eqs. (10)±( 12) is obtained by the usual test-function method using Green's formula. The weak variational formulation corresponding to the modal analysis of the structural-acoustic problem is then obtained by setting F d 0 in Eq. ( 1). The eigenvalue structural-acoustic variational formulation is then stated as follows. Find x 2 and u P C u and p P C p , such that for all du P C u and d p P C p ,we have

ku; duÀ Z R pn du dr x 2 Z X S q S u du dx; 13 1 q F Z X F rp rdp dx x 2 q F c 2 Z X F pdp dx x 2 Z R undp dr; 14 with the constraint 1 q F c 2 Z X F p dx Z R undr 0: 15
The variational formulation de®ned by Eqs. (13)±( 15) is original due to the presence of the constraint de®ned by Eq. (15) which regularizes the u; p formulation for x 0 (see also [START_REF] Ohayon | Symmetric formulations for modal analysis of internal ¯uid structure systems[END_REF]). In eect, usually, only Eqs. ( 13) and ( 14) are written, and as pointed out above, are not valid for x 0. In the case of a direct ®nite element discretization of Eqs. (13)±(15), we obtain a matrix eigenvalue system of the type AY x 2 BY,i n which A and B are not symmetric. That is why various symmetric formulations using for the ¯uid pressure ®eld p and displacement potential u, de®ned up to an additive constant and such that u F ru, have been derived. The resulting symmetric formulations are then obtained by elimination of p or u (see [START_REF] Felippa | Mixed variational formulation of ®nite element analysis of acoustoelastic/slosh ¯uid±structure interaction[END_REF][START_REF] Liu | Variational approach to ¯uid±structure interaction with sloshing[END_REF]), or for detailed complete analysis, [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF]. In the present case, we are not considering a direct ®nite element approach of the variational formulation de®ned by Eqs. (13)±(15).

Symmetric reduced model

Let us consider the following two basic problems. The ®rst one corresponds to the acoustic modes in rigid motionless cavity and is obtained by setting u 0 into Eqs. ( 14) and (15). The calculation of these acoustic modes is generally done by using a ®nite element procedure. If we introduce the admissible subspace

C Ã p of C p C Ã p p P C p ; Z X F p dx 0 ; 16
the variational formulation of acoustic modes is stated as follows: ®nd x 2 > 0 and p P C Ã p such that, Vdp P C Ã p , we have

1 q F Z X F rp rdp dx x 2 1 q F c 2 Z X F pdp dx; 17 with the constraint Z X F p dx 0: 18
It should be noted that, in practice, if the constraint condition (18) is ``omitted'', we only add a ®rst nonphysical zero-frequency constant pressure mode, the other modes corresponding to x T 0 remaining the same as those de®ned by Eqs. ( 17) and ( 18). This zero-frequency mode must not be retained in any Ritz± Galerkin projection analysis. In addition, we have the following orthogonality conditions

1 q F c 2 Z X F p a p b dx l a d ab ; 1 q F Z X F rp a rp b dx l a x 2 a d ab : 19
The second basic problem corresponds to the static response of the ¯uid to a prescribed wall normal displacemnt un. The solution, denoted as p s un, is given by Eq. ( 9). For any deformation un of the ¯uid± structure interface, p s un belongs to a subset of C p , denoted as C un

C un p s P C p ; p s À q F c 2 jX F j Z R undr : 20
In the variational formulation de®ned by Eqs. (13)±(15), p is searched under the form p p s un X Np a1 r a p a ; 21 in which N p denotes the number of retained acoustic modes. The decomposition (21) is unique. In addition, it should be noted that, since each eigenvector p a corresponding to x a T 0, veri®es the constraint de®ned by Eq. (18), then, using Eq. ( 9), we deduce that p and un satisfy the constraint de®ned by Eq. ( 15). The decomposition de®ned by Eq. ( 21) corresponds to a decomposition of the admissible class C p into the direct sum of the admissible classes de®ned, respectively, by Eqs. (20) and ( 16)

C p C un È C Ã p : 22
Following Eq. ( 21), the test function dp is then searched under the following form:

dp p s dun X Np a1 dr a p a : 23
The variational formulation in du de®ned by Eq. (13) becomes ku; duk s u; duÀ

X Np a1 r a Z R p a n du dr x 2 Z X S
q S u du dx; 24 in which ku; du is de®ned by Eq. (2) and k s u; du is such that

k s u; du q F c 2 jX F j Z R undr Z R dundr : 25
If we consider a ®nite element discretization of the structure, the corresponding discretized form of Eq. ( 24) can be written as

K K s U À X n a1
C a r a x 2 MU; 26 in which symmetric matrices K and K s correspond to ®nite element discretization of stiness symmetric bilinear forms de®ned by Eqs. (2)±( 4) and (25), respectively. In Eq. ( 26), M denotes the structural-symmetric mass matrix and rectangular coupling matrix C a corresponds to the discretization of the coupling ¯uid±structure contribution R r pdundr. The discretized form of Eq. ( 14) in dp can be rewritten as x 2 a l a r a x 2 l a r a x 2 C T a U: 27

From Eqs. ( 26) and ( 27), we obtain the following symmetric matrix reduced model

K tot 0 0 Diagl a U r x 2 M tot D D T Diagl a =x 2 a U r ; 28
in which r denotes the vector of N generalized coordinates r a , with 1 6 a 6 N p , and

K tot K K s ; 29 M tot M X Np a1 1 x 2 a l a C a C T a ; 30
D a X Np a1 1 x 2 a C a : 31
Further diagonalization of Eq. ( 28) implies a projection of U on the solutions of the following eigenvalue problem

K tot U b k b M tot U b : 32 Setting U X Nu b1 q b U b ; 33
in which q b are the generalized coordinates describing the structure. Using the orthogonality conditions associated with the solutions of Eq. (32), i.e.

U T b H M tot U b l s b d bb H and U T b H K tot U b l s b k b d bb H , Eq. (28) be- comes Diag k b 0 0 Diag x 2 a q r x 2 I Nu C ba C ba T I Np q r : 34
It should be noted that in vacuo structural modes can also be used (instead of those de®ned by Eq. ( 32)), but the resulting matrix system would not be diagonal with respect to U.

Wall impedance condition

Referring to [START_REF] Kehr-Candille | Elasto-acoustic damped vibrations. Finite element and modal reduction methods[END_REF] or for more detailed modeling, to [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF], wall impedance condition corresponds to a particular ¯uid±structure interface modeling. This interface is considered as a third medium with in®nitesimal thickness, without mass, and with the following constitutive equation p jxZxun À u F n; 35 in which Zx denotes a complex impedance. Eqs. ( 7) and ( 11) must be replaced by Eq. ( 35), using op=on q F x 2 u F n.

Case of a liquid with a free surface

Let us consider a liquid with a free surface at rest denoted as C. If we neglect gravity eects, the boundary condition on C is such that p 0j C : 36

In this case, constraint condition (12) or ( 15) is replaced by Eq. (36). Eq. ( 9) is replaced by p s 0. Admissible space de®ned by Eq. ( 16) becomes C Ã p fp P C p ; p 0g. In this case, the static problem de®ned in Section 3.3 leads to a zero-pressure ®eld. Let us remark that in this case, the ``structural'' modal basis may be constituted by the hydroelastic incompressible modes using the classical added mass operator.

The reduced modal matrix models has been extended to the dissipative case (dissipative internal ¯uid with nonhomogeneous local impedance wall condition) in [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF].

Incompressible hydroelastic-sloshing problem

We consider the harmonic vibrations of an elastic structure partially ®lled with an homogeneous, inviscid and incompressible liquid, taking into account gravity eects on the free surface C. After an original derivation of the linearized equations of the ¯uid±structure coupled problem, introducing an appropriate linear constraint in order to obtain a well-posed problem at zero frequency, we discuss the main ideas leading to a reduced model of the ``liquid subsystem''. For this purpose, sloshing modes in rigid motionless cavity are introduced as Ritz projection vector basis, including mandatory zero-frequency (static) solution of the coupled system. Let us recall that the internal ¯uid±structure system is characterized by a strong resonant behavior. Therefore, a direct ®nite element description can be costly. That is why, starting from one of the possible variational formulations of the problem, reduced symmetric matrix models of this system are constructed.

Hydroelastic-sloshing equations

Structure subjected to a ¯uid pressure loading

The notations are the same that those de®ned in Section 3.1. The weak variational formulation describing the response of the structure X S to given variation F d of the applied forces with respect to the equilibrium state on the external structure boundary oX S nR, and to ¯uid pressure ®eld p acting on the internal ¯uid±structure interface R is written as follows.

For all real x and Vdu P C u , ®nd u P C u such that

b ku; duÀx 2 Z X S q S u du dx À Z R pn du dr Z oX S nR F d du dr; 37 in which b k k k R : 38
In Eq. (38), ku; du is de®ned by Eq. ( 2), and k R is the elastogravity symmetric bilinear form such that (see [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF], chapter 6)

k R u; duÀ 1 2 q F g Z R zn 1 u du u z dun dr Z R
zn 1 duu du z un dr : 39

Fluid subjected to a wall normal displacement

We assume that the liquid is homogeneous, inviscid and incompressible. Free surface C is horizontal at equilibrium. We denote by z the external unit normal to C, and by g the gravity. The notations are similar to those of Section 3.2. The local equations describing the response of the ¯uid to a prescribed arbitrary normal displacement un of the ¯uid±structure interface R are such that

rp À q F x 2 u F 0j X F ; 40 ru f 0j X F ; 41 u F n unj R ; 42 p q F gu F nj C ; 43 curl u F 0j X F : 44
Eq. ( 41) corresponds to the incompressibility condition. Eq. ( 43) is the constitutive equation on the free surface C due to gravity eects.

Relation between static pressure p s and un

For x 0, Eqs. (41)±(43) lead to the constant pressure static ®eld which is related to the normal wall displacement by the relation (see [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF][START_REF] Ohayon | The eect of wall motion on the governing equations of structures containing ¯uids[END_REF])

p s À q F g jCj Z R
undr; 45 in which jCj denotes the measure of the area of free surface C.

Equations in terms of p and un

The elimination of u F between Eqs. ( 40 44) . This condition is usually omitted in the literature.

Variational formulation in terms of (u; p)

Let dp be the test function, associated to p, belonging to the admissible space C p . The weak variational formulation corresponding to Eqs. (46)±( 49) is obtained by the usual test-function method using Green's formula. Recalling Eq. (37), the variational formulation of the hydroelastic-sloshing problem is then stated as follows. Find u P C u and p P C p , such that for all du P C u and dp P C p ,w eh a v e b ku; duÀx 2

Z X S q S u du dx À Z R pn du dr Z oX S nR F d du dr; 50 1 q F Z X F rp rdp dx x 2 q F g Z C pdp dx x 2 Z R undp dr; 51 with the constraint 1 q F g Z C p dr Z R undr 0: 52
The variational formulation de®ned by Eqs. (50)±( 52) is original due to the presence of the constraint de®ned by Eq. ( 49) which regularizes the (u; p) formulation for x 0 [START_REF] Ohayon | Symmetric formulations for modal analysis of internal ¯uid structure systems[END_REF].

Symmetric reduced matrix model

Let us consider the following two basic problems. The ®rst one corresponds to the sloshing modes in rigid motionless cavity and is obtained by setting u 0 into Eqs. ( 47) and (49). The calculation of these acoustic modes is generally done by using a ®nite element procedure. It should be noted that, in practice, if the constraint condition (55) is ``omitted'', we only add a ®rst nonphysical zero-frequency constant pressure mode, the other modes corresponding to x T 0 remaining the same as those de®ned by Eqs. (54) and (55). This zero-frequency mode must not be retained in any Ritz± Galerkin projection analysis. In addition, we have the following orthogonality conditions

1 q F g Z C p a p b dr l a d ab ; 1 q F Z X F rp a rp b dx l a x 2 a d ab : 56
The second basic problem corresponds to the static response of the ¯uid to a prescribed wall normal displacement un. The solution, denoted as p s un, is given by Eq. (45). For any deformation un of the ¯uid± structure interface, p s un belongs to a subset of C p , denoted as C un

C un p s P C p ; p s À q F g jCj Z R undr : 57
In the variational formulation de®ned by Eqs. ( 50 If we consider a ®nite element discretization of the structure, the corresponding discretized form of Eq. ( 60) can be written as

b K K s U À X n a1 C a r a À x 2 MU F d ; 61
in which symmetric matrices b K and K s correspond to ®nite element discretization of stiness symmetric bilinear forms de®ned by Eqs. ( 38) and (60), respectively. The discretized form of Eq. (51) in dp can then be written as

x 2
a l a r a x 2 l a r a x 2 C T a U: 62

From Eqs. ( 61) and (62), we obtain a symmetric matrix reduced model whose expression is similar to the one given by expression (28). Similarly to Section 3.6, further diagonalization can be obtained by setting

U X Nu b1 q b U b ; 63
in which q b are the generalized coordinates describing the structure and U b are the eigenmodes of an eigenvalue problem similar to the one described by Eq. (32). We then obtain a similar matrix system than the one described by Eq. ( 34)

Diag k b 0 0 Diag x 2 a q r À x 2 I Nu C ba C ba T I Np q r F d 0 : 64
It should be noted that we can also use the incompressible hydroelastic modes, i.e. the modes of the coupled system constituted by the elastic structure containing an incompressible liquid, with p 0o nC (through an added mass operator). In this case, the resulting matrix system is not completely diagonal with respect to U variables [START_REF] -P. Morand | Fluid±Structure Interaction[END_REF].

Conclusion

Using a new approach, we have presented convenient and versatile tools for low frequency modal analysis of elastic structures containing linear inviscid homogeneous ¯uids for structural-acoustics and incompressible hydroelastic-sloshing problems. Those formulations, together with dynamic substructuring techniques [START_REF] Ohayon | Dynamic substructuring of damped structures using singular value decomposition[END_REF], lead to symmetric reduced matrix systems expressed in terms of generalized coordinates for the ¯uid (and if necessary for the structure).

  If we introduce the admissible subspace C Ã p of C p

			Z
	C Ã p p P C p ;	p dr	0	;	53
			C
	the variational formulation of acoustic modes is stated as follows: ®nd x 2 > 0andp P C Ã p such that, for all dp P C Ã p , we have
	1 q F	Z X F	rp rdp dx x 2 1 q F g	Z C	pdp dr;	54
	with the constraint
	Z

C

p dr 0: 55

  p denotes the number of retained sloshing modes. The decomposition (58) is unique. In addition, it should be noted that, since each eigenvector p a corresponding to x a T 0, veri®es the constraint de®ned by Eq. (55), then, using Eq. (45), we deduce that p and un satisfy the constraint de®ned by Eq. (52). The decomposition de®ned by Eq. (58) corresponds to a decomposition of the admissible class C p into the direct sum of the admissible classes de®ned, respectively, by Eqs. (56) and (57), C p C un È C Ã p . The variational formulation in u de®ned by Eq. (50) becomes S u du dx; 59 in which b ku; du is de®ned by Eq. (38) and k s u; du is such that

	b ku; duk s u; duÀ	X N p	r a	Z	p a n du dr x 2	Z
					a1		R			X S
	k s u; du	q F g jCj	Z R	undr	Z	R	dundr	:	60

)±(52), p is searched under the form p p s un X Np a1 r a p a ; 58 in which N q


C Hoareau

J.-F Deü

R Ohayon

Hydroelastic linearized vibrations taking into account prestressed effects due to internal liquid weight: numerical vs. experimental results

Keywords: Liquid weight, Pressure prestress effect, Geometrical nonlinearity, Hydroelastic vibrations, Finite element method

This paper concerns the numerical simulation of the hydroelastic vibrations of an elastic tank containing a free surface liquid, around a geometrically nonlinear equilibrium configuration that takes into account the prestress effect due to the liquid weight. The computational methodology is the following. Firstly, the hydrostatic pressure on the fluid-structure interface is considered as a non-uniform follower force which depends on the liquid height. Therefore, the prestressed equilibrium at a given free-surface level, called current configuration in the paper, is computed as a geometrically nonlinear problem, in which the height of the liquid is considered as an evolution parameter. We start from an empty structure, defined as the reference configuration, which is progressively filled to reach the level of the chosen configuration. Therefore the calculation of hydroelastic vibrations of a tank filled with liquid as a function of a emptying rate is straightforward. The methodology is validated with experimental results extracted from literature.

Introduction

This study deals with a hydroelastic problem, namely the linearized vibrations of elastic structures, filled with an inviscid incompressible fluid, taking into account the prestress effect due to the fluid weight at the fluid-structure interface and neglecting the solid weight. Note that the weight effect in structural vibrations has already been the subject of several analyses (see e.g. [START_REF] Chen | Dynamic effects of weights on vibration reduction by a nonlinear energy sink moving vertically[END_REF]). Let us recall that linearized hydroelastic vibration analysis is of great interest in the study of liquid rocket launchers control and stability. Generally, in hydroelastic vibrations, the elasticity of the tank is of prime importance together with the kinetic energy of the liquid considered as incompressible (the called added mass effect [START_REF] -P. Morand | Fluid structure interaction[END_REF][START_REF] Van Brummelen | Added Mass Effects of Compressible and Incompressible Flows in Fluid-Structure Interaction[END_REF]). This means that compressibility effects of the liquid (its potential energy) as well as the effects of sloshing of the free surface (potential energy due to gravity and/or surface tension) can be neglected [START_REF] Ohayon | Novel formulation for the effects of sloshing with capillarity on elastic structures in linear dynamics[END_REF]. Of course, the gravity effects are present to define the equilibrium static configuration around which the vibrations are studied. Furthermore, those tanks can be pressurized, and effects of the pressure above the liquid can be significant especially for liquid propelled launchers. They have been taken into account through classical prestress geometrical terms considering specific linearization approaches [START_REF] -P. Morand | Internal pressure effects on the vibration of partially filled elastic tanks[END_REF][START_REF] Schotté | Linearized formulation for fluid-structure interaction: Application to the linear dynamic response of a pressurized elastic structure containing a fluid with a free surface[END_REF]. For general nonlinear fluid-structure interaction in the time domain, see for instance [START_REF] Bazilevs | Computational Fluid-Structure Interaction: Methods and Applications[END_REF]. According to authors knowledge, the weight of the liquid on the structural vibrations has never been taken into account in coupled three-dimensional fluid-structure computations models [START_REF] Amabili | Nonlinear Vibrations and Stability of Shells and Plates[END_REF].

This constitutes the first originality of the paper. Nevertheless, it should be mentioned that semi-analytical approaches exist for the study of pre-deformed membrane in a semi-infinite fluid medium [START_REF] Di Trolio | Modeling added mass effects on the vibrations of air-backed, predeformed membranes[END_REF] or for rigid cylindrical tank with an elastic bottom [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part II: Linear axisymmetric vibration analysis[END_REF].

In addition, for liquid propelled launchers, the evolution of the frequencies as a function of the emptying rate is of prime importance. The period of vibration of the liquid structure system being very small compared to the time needed to emptying the tank during the launch, the vibrations of the system are classically performed for various liquid height which is considered as an evolution parameter. One major interest of the method presented here is to compute the vibrations considering nonlinear static prestressed state obtained through a procedure considering the weight of the liquid, depending on its height, from an empty situation up to a required filling level. It allows the computation of the hydroelastic vibrations corresponding to the appropriate liquid height. This constitutes the second originality of the approach.

Usually, the effects of geometric nonlinearities of prestressed type, due to the presence of the liquid (effects of gravity linked to the hydrostatic pressure which induces a non-uniform following force) are generally neglected in fluid-structure vibrations problems. Therefore, it is important to assess these effects on the linearized vibrations analysis around the prestressed state. Studies concerning the estimation of the prestress, from a nonlinear static evolution problem due internal or external pressure, can be found in [START_REF] Leissa | Vibration of Shells[END_REF][START_REF] Amabili | Review of studies on geometrically nonlinear vibrations and dynamics of circular cylindrical shells and panels, with and without fluid-structure interaction[END_REF][START_REF] Amabili | Displacement dependent pressure load for finite deflection of doubly-curved thick shells and plates[END_REF] and notably for hydrostatic pressure in [START_REF] Haßler | On the static interaction of fluid and gas loaded multi-chamber systems in large deformation finite element analysis[END_REF][START_REF] Hoareau | Nonlinear equilibrium of partially liquid-filled tanks: A finite element/level-set method to handle hydrostatic follower forces[END_REF][START_REF] Narayanan | Monolithic and partitioned approaches to determine static deformation of membrane structures due to ponding[END_REF]. Throughout all the paper, in static as in dynamic, it should be noted that the equation in the structure and on the fluid-structure interface are considered through an appropriate reciprocal mapping on a reference configuration, defined as the structure empty structure (in vacuo). This constitutes also an original aspect of the computational approach.

The computational methodology used in this work is based on the finite element method for both the structure and fluid domains (see e.g. [START_REF] -P. Morand | Fluid structure interaction[END_REF]). For FSI problems with external fluid, other numerical approaches can advantageously be used, such as the immersed methods [START_REF] Griffith | Immersed methods for fluid-structure interaction[END_REF]. To the best authors knowledge, very few experiments exist in literature to quantify the prestressed effect investigated in this contribution. The computational results presented in this work can still be compared with the experimental results given in [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part I: Experiment[END_REF] in order to validate the numerical developments.

The outline of the paper is as follows. Following this introduction, in Section 2, the nonlinear static liquidstructure computation as function of the liquid height is introduced. In Section 3, the linearized harmonic hydroelastic problem is developed. In Section 4, the computational methodology is presented. In Section 5, comparison and validation with experiments are discussed. Finally, in Section 6 conclusion is given followed by references.

Nonlinear static liquid-structure computation as function of the liquid height

For general studies on theoretical and computational nonlinear structural mechanics, see references [START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF][START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Wriggers | Nonlinear Finite Element Methods[END_REF][START_REF] Borst | Nonlinear Finite Element Analysis of Solids and Structures[END_REF][START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF].

Definition of the parameterized prestressed configuration

The structure is supposed to be nonlinear elastic isotropic and relevant of Saint-Venant Kirchhoff constitutive relations (i.e. large displacements and small strain assumption). We define as the prestressed configuration, the stresses in the structure containing the liquid given in the called current configuration

Ω ′ s such that ∂Ω s = Σ ′ i ∪ Σ ′ u ∪ Σ ′ e
where Σ ′ i is the fluid-structure interface, Σ ′ u is the surface with prescribed displacements and Σ ′ e is the surface of application of external load (see Fig. 1). Let (e x , e y , e z ) be the orthogonal unit vectors which define the Cartesian coordinates system. Concerning the liquid with horizontal free-surface under constant gravity force g = -ge z , the equations are relevant of the classical hydrostatic equation due to gravity. We do not consider here surface tension of the free surface. Thus, the prestressed state in the structure is only due to the weight of liquid seen as non-uniform follower forces [START_REF] Hibbit | Some follower forces and load stiffness[END_REF][START_REF] Elishakoff | Controversy associated with the so-called "follower forces": critical overview[END_REF] In practice, this problem will be solved progressively by filling the tank from an empty state until the prestressed configuration. Thus, the problem is treated as a parameterized problem depending on the fluid height denoted by l.

Definition of the reference configuration

It is important to note (see [START_REF] Hoareau | Nonlinear equilibrium of partially liquid-filled tanks: A finite element/level-set method to handle hydrostatic follower forces[END_REF]) that we will consider as the reference configuration, the empty structure, at equilibrium, at rest in a Galilean frame of reference. More precisely, this means that the current structural configuration, with the fluid structure interface, will be treated through an appropriate reciprocal mapping on this reference configuration as seen in Fig 2.

Computation of the prestressed state as a function of the liquid height

The structural displacement field is denoted by u s = u s e x + v s e y + w s e z in a Cartesian coordinate system.

Thus, the current position x ′ = x ′ e x + y ′ e y + z ′ e z is defined by:

x ′ = x + u s (1)
where x = xe x + ye y + ze z is the reference position of a material point in the structure on its reference 80 configuration Ω s bounded by ∂Ω s = Σ i ∪ Σ u ∪ Σ e , which are respectively, the interface, the surface with prescribed displacement and the distributed load surfaces mapped on the reference configuration. It should be noted that the nonlinear static problem is expressed only in terms of structural displacement. The static pressure, which depends on the structural displacement, is defined on the current fluid-structure interface Σ ′ i , considered on a structural material point such that:

p s (l; u s ) =    0 if z ′ ≥ l -ρ f g(z ′ -l) if z ′ < l (2)
where

z ′ = (x + u s ) • e z .
Let the transformation gradient and the Green-Lagrange strain operator, with respect to the reference position x be defined as:

F (u s ) = I + Gradu s (3)
E(u s ) = 1 2 (F T (u s )F (u s ) -I) (4)
The determinant of the transformation gradient is denoted by J as:

J(u s ) = det(F (u s )) (5)
For sake of brevity, the structural body forces are neglected. The given force depends upon the liquid height l as seen in Eq. ( 2). We will consider here the structure as clamped on Σ u . The structural mechanical equations are given on the reference configuration Ω s , considering F = F (u s ), E = E(u s ) and J = J(u s ), such that:

-Div(F S) = 0 in Ω s (6)

D : E = S in Ω s (7)
F Sn = JF -T p s n on Σ i (8)
u s = 0 in Σ u (9)
where Eq. ( 6) is the balance of linear momentum, Eq. ( 7) is the Saint-Venant Kirchhoff constitutive relation, which links the second Piola-Kirchhoff stress S, consequently S = S(u s ), where D is a constant fourth order 90 tensor of elasticity in which the notation ":" correspond to the trace operator, Eq. ( 8) is the equilibrium of follower forces surface traction and Eq. ( 9) correspond to a fixed displacement boundary condition. We recall that the pressure p s depends of the unknown function u s and of the given data only (i.e. the height l).

Let C u be the admissible space of regular displacement field:

C u = {v "regular enougth", such that v = 0 in Σ u } (10)
Let introduce a virtual displacement field δu in C u . The static nonlinear problem consists in finding the displacement field u s in C u such that:

δW int (u s , δu) -δW ext (l; u s , δu) = 0, ∀ δu ∈ C u (11)
where δW int and δW ext are respectively the internal and external virtual works defined as:

δW int (u s , δu) = Ωs δE(u s , δu) : S(u s )dΩ (12) δW ext (u s , δu) = Σi -p s (l; u s ) J(u s )F -T (u s )n • δu dΣ (13)
where the term δE is the virtual Green-Lagrange matrix defined as follows:

δE(u s , δu) = 1 2 δF T (δu)F (u s ) + F T (u s )δF (δu) (14)
with

δF (δu) = Gradδu ( 15
)
The resolution of the problem Eq. ( 11) is based on the finite element approach. It should be noted that, for computational and software facilities, the computation of the static response of the structure, for each level of liquid, is carried out on the reference configuration defined in Section 2.2. A finite element discretization of the 3D domain Ω h s and its boundaries

∂Ω h s = Σ h i ∪ Σ h u ∪ Σ h
e is now considered (see Fig. 1 for the continuous representation on the current configuration). The discretized displacement field and its virtual form are denoted u h s and δu h such that:

u h s = N u u s and δu h = N u δu (16)
where N u is a known shape function matrix, u s is a unknown nodal vector and δu is the associated virtual nodal displacement. Considering Eq. ( 11), the discretized problem consists in finding u such that:

f int (u s ) -f ext (l; u s ) = 0 (17)
where f int and f ext are respectively the internal and external finite element forces defined as:

δW int (u h s , δu h ) = δu T f int (u s ) ( 18
)
δW ext (u h s , δu h ) = δu T f ext (u s ) (19)
The resolution of this nonlinear problem can be done considering a Newton-Raphson type algorithm. The geometrical, the material and the follower forces tangent stiffness matrices, respectively denoted by K g , K m and K f defined as:

∂f int ∂u s = K g + K m and ∂f ext ∂u s = K f ( 20
)
Those operators are needed in the Newton-Raphson loop to obtain a quadratic convergence rate of the nonlinear solution. We will see in the following that K m and K g are also useful to take into account the prestressed effect on the linearized hydroelastic vibrations of the fluid-structure system. If the structure is completely free (i.e. no prescribed displacement), the tangent stiffness matrices should be adapted (e.g. through a pseudo-inverse analysis) to the singularity involved with respect to the rigid body motions.

Linearized harmonic hydroelastic problem

Liquid-structure assumption of the linearized hydroelastic problem

Let assume linearized vibrations of the fluid-structure system around a known equilibrium state (defined as the prestressed state which depends on the fluid height). The circular frequency is denoted by ω. Regarding the liquid, it is assumed to be inviscid, homogeneous and incompressible. Thus, those small motions are irrotational. Of course, for large amplitude nonlinear vibrations, this irrotational assumption would be not relevant but this case is out of the scope of this article. Nevertheless for moderate nonlinear vibrations, not considered here, the hypotheses of irrotationality could be maintained. Consequently, the equations of the fluid can be expressed through a scalar field namely the pressure field or the displacement potential field.

We suppose that during those linearized vibrations, the potential energy of the free-surface, due to gravity and to surface tension, are negligeable with respect to the kinetic energy of the fluid and to the potential and kinetic energies of the structure. The consequent is that, for the fluid dynamic equation, the free-surface condition is expressed as a zero pressure condition (or zero potential of displacement condition). Let us recall that the gravity effect, on the fluid structure interface, has been of course taken into account for the static equilibrium and in the parameterized evolution of the static prestressed computations.

Linearized equations in the structure

The structural part of the system is supposed to be excited by a harmonic distributed load denoted by f . The response of the prestressed system is defined by knowing the dynamic structural displacement u, knowing the static displacement u s which is the results of the previous static nonlinear problems. Thus, the total displacement denoted by u tot is defined as:

u tot = u s + u (21)
The structure perturbation u is supposed to be small compared to the characteristic length of the structure t s (e.g. the tank thickness) such that:

|| u || L2 t s << 1 (22)
The current harmonic configuration is thus supposed to be the same configuration of the prestressed one (i.e.

called the current configuration from the static problem). We take a Taylor series expansion of the virtual works, in which we retain up to the first order term in u written as follows:

δW acc (u tot , δu) = 0 -ω 2 m(u, δu) (23)
δW int (u tot , δu) = δW int (u s , δu) + k t (u s ; u, δu) (24)
δW ext (h; u tot , δu) = δW ext (l; u s , δu) static -c(u s ; p, δu) + f (δu) dynamic ( 25
)
where the term p corresponds to the dynamic pressure perturbations of the fluid at the fluid-structure interface such that:

p tot = p s + p (26)
The terms involved in the linearized virtual acceleration and external works can be expressed on the reference configuration knowing the static displacement as follows:

m(u, δu) = Ωs ρ s u • δu dΩ (27) f (δu) = Σe f • δu dΣ (28) c(u s ; p, δu) = Σi pJ(u s )F -T (u s )n • δu dΣ ( 29
)
where m is the virtual work induced by the inertia of the structure, f is the external fluctuation load and c is the coupling operator between the fluctuation of pressure and the virtual displacement. The virtual potential work from the tangent stiffness is defined as the sum of the geometrical and the material stiffness contributions. Then, the expression of the geometrical and material virtual stiffness works k m and k g are given as:

k m (u s ; u, δu) = Ωs δE(u s , δu) : D : [ε(u) + γ(u s , u)] dΩ (30) k g (u s ; u, δu) = Ωs γ(u, δu) : D : S(u s ) dΩ (31)
k t = k m + k g (32)
where the operators ε and γ, such that γ(u, v) = γ(v, u), depend on the gradient defined for given displacement fields u and v as follows:

ε(u) = 1 2 [Grad T u + Grad u] (33) γ(u, v) = 1 2 [Grad T uGrad v + Grad T v Grad u] (34)
Finally, the linearized harmonic problem consists in finding u d ∈ C u , knowing u s such that:

k t (u s ; u, δu) -c(u s ; p, δu) -ω 2 m(u, δu) = f (δu), ∀ δu ∈ C u (35)
It can be seen that the operators k m and k g are symmetric in u and δu. For the same reasons indicated in the previous section, the variational formulation of the local equations of the structure is posed on the reference configuration (see Section 2.2). That is why the coupled term representing the virtual work done by the dynamic pressure fluctuation p on the virtual normal displacement is written here not on the current configuration e.g. the structure filled with a liquid with free surface at height l, but on the reference 130 configuration.

Linearized liquid equations

We are interested in the linearized vibrations of the liquid part. Since the liquid is irrotational, the displacement u f is derived from a potential displacement ϕ with respect to the liquid configuration at rest Ω ′ f (the static liquid configuration) as:

135 u f = gradϕ (36)
The relations in the fluid domain involving the fluctuation of dynamic pressure p and the fluid motion such that are given as:

div u f = 0 in Ω ′ f (37) gradp -ω 2 ρ f u f = 0 in Ω ′ f (38) u f • n ′ = u • n ′ in Σ ′ f ( 39
)
where the Eq. ( 37) is due to the incompressibility condition of the fluid motion, Eq.

         ∆ϕ = 0 in Ω ′ f gradϕ • n ′ = u • n ′ in Σ ′ i ϕ = 0 in Γ ′ (40a) (40b) (40c)
Eq. ( 40) involves only the geometry and shows that ϕ is a linear function of u f and that consequently, "the kinetic energy" of the incompressible liquid:

E kin = 1 2 Ω ′ f ρ f | ∇ϕ | 2 dΩ e.g. 1 2 Ω ′ f ρ f | u f | 2 dΩ , (41)
is a quadratic form of u f such that the discretized form is the classical added mass matrix operator [START_REF] -P. Morand | Fluid structure interaction[END_REF][START_REF] Van Brummelen | Added Mass Effects of Compressible and Incompressible Flows in Fluid-Structure Interaction[END_REF].

Let us take the variational formulation of Eqs. (40) in the current configuration of the fluid:

Ω ′ f gradϕ • gradδϕ dΩ ′ + Σ ′ i u • n ′ δϕ dΣ ′ = 0 ∀δϕ ∈ C ϕ ( 42
)
where C ϕ the space of regular functions define on Ω ′ f verifying Eq. (40c). It should be noted that the term expressed in Σ ′ i , will be written, through the reciprocal mapping, on the reference configuration, defined in Section 2.2, such that:

Σ ′ i δϕn ′ • u dΣ ′ = Σi δϕ J(u s )F -T (u s )n • u dΣ (43)
For sake of brevity, we will use the same notations for the unknown fields u and ϕ defined either on Ω ′ f and Σ ′ i or Ω f and Σ i . Finally, the problem expressed in term of the variable ϕ (instead of p) and u, can be written as follows:

h(δϕ, ϕ) + c(u s ; δϕ, u) = 0, ∀δϕ ∈ C ϕ ( 44
)
where the operator h correspond to the virtual work of the fluid kinetic energy and c corresponds to the liquid -structure coupling operator:

h(ϕ, δϕ) = Ω ′ f gradϕ • gradδϕ dΩ ′ (45) c(u s ; δϕ, u) = Σi δϕJ(u s )F -T (u s )n • u dΣ (46)
It should be noted that the discretization of the two terms of the variational formulation will be carried through two different aspects. The first contribution which involves only the displacement potential of the liquid will be done on this prestressed static configuration, around which we study the dynamic fluctuation.

The second contribution, i.e. the coupled term, will be discretized on the reference configuration through a reciprocal mapping method. This is not an issue due to the fact that the construction of the operator is based on the known static displacement at the interface. The prestressed static configuration of the fluid can be done: (i) through a computation of a pre-deformed fluid mesh or (ii) considering a re-meshing of the fluid domain. The coupling term here gives a relation between the structural acceleration at the interface and the virtual potential of displacement.

Finite element discretization of the coupled problem

The potential of fluid displacement and the fluctuation of pressure are linked considering Eq. (36) and Eq. ( 38) such that:

p = ω 2 ρ f ϕ (47)
Consequently, the linearized structural virtual equation is written as:

k t (u s ; u, δu) -ω 2 ρ f c(u s ; ϕ, δu) -ω 2 m(u, δu) = f (δu), ∀ δu ∈ C u (48)
As done in Section 2, let introduce the finite element discretization of the fields u h and ϕ h and their respective virtual fields. The discretized structural equation of Eq. ( 48) is given as:

K t u -ω 2 Mu -ω 2 ρ f Cϕ = f (49)
Concerning the discretization of Eq. ( 44), the volumic term h(ϕ, δϕ) is discretized using the finite element discretization of Ω ′h f . On the other hand, the second term of Eq. (44) which is the liquid-structure coupling term is discretized on the reference configuration (let us recall that the chosen configuration has been previously defined in Fig. 2 with no prestressed effect) using Eq. (46) :

Hϕ + C T u = 0 (50)
All the discretized operator are given as follows:

k t (u h s ; δu h , u h ) ⇒ δu T K t (u s )u (51) m(δu h , u h ) ⇒ δu T Mu (52) c(u h s ; δu h , ϕ h ) ⇒ δu T C(u s )u (53) f (δu h ) ⇒ δu T f (54) h(δϕ h , ϕ h ) ⇒ δϕ T Hϕ (55) c(u h s ; δϕ h , u h ) ⇒ δϕ T C(u s )u (56)
Taking into account the free-surface prescribed mandatory equation (40c), and for sake of simplification, using the same notations for C and H, the coupled matrix system, resulting from Eq. ( 49) and Eq. ( 50)

writes:     K t O O O   -ω 2   M ρ f C ρ f C T -ρ f H       u ϕ   =   f 0   ( 57
)
or equivalently

K t u -ω 2 (M + ρ f M a )u = f with M a = CH -1 C T (58)
Solving Eq. ( 58), we deduced directly ϕ and the discretized pressure p as follows:

ϕ = -H -1 C T u and p = ω 2 ρ f ϕ (59)
in which K t is symmetric and H, initially of kernel of dimension 1, is invertible using Eq. (40c), therefore the added mass matrix M a is a full, symmetric positive definite matrix, involving only structural fluid-structure interface degrees of freedom. In a practical situation, the matrix H -1 is never computed directly.

Computation methodology

The aim is to compute the linearized hydroelastic vibrations of an elastic tank prestressed by a liquid weight for various given liquid free-surface height l i ∈ [l min , l max ]. At first, for each values l i , static nonlinear solutions u (i) s = u s (l i ) and the associated tangent stiffness matrix K (i) t = K t (l i ) are then computed step by step from l min to l max (see Section 2), which is an originality of the approach. We recall here that the discretized matrices are computed on the previously defined reference configuration in Section 2.2. Those tangent stiffness operators K (i) t are stored and will be used for the calculation of the prestressed hydroelastic vibrations for various levels of liquid presented in Eq. ( 58). As mentioned in the previous section, the added mass matrix M a will never be constructed directly. Eq. ( 58) shows that the coupled eigenvalue hydroelastic problem is expressed only in terms of structural field unknown u (the potential of displacement and the pressure are then calculated on a straightforward manner by Eq. ( 59)), for a given level. The solution of Eq.

(58) can be searched using an appropriate basis of the structural admissible structural space C u . Therefore, we can choose as a possible basis of C u the eigenvectors of the prestressed structure without taking into account liquid inertia:

[K (i) t -ω 2 M]u = 0 (60)
The first N eigenvectors denoted by u

B (i) u = u (i) α α=1...m (61)
A remark has to be addressed at this point concerning the positivity of the tangent stiffness operator. For the practical application, we suppose that the action of the liquid weight in the structure tends to stiffen the system. Thus, we do not consider here the dynamic buckling, which can be the case for immersed structures.

In the following, the vector of the basis B (i)

u will called the prestressed dry modes. Of course question of truncation and related validations will be addressed later in the numerical application. Let us stress here that the computation of the added mass matrix no longer necessitates the direct inverse of a matrix. We proceed by considering known matrices H (i) = H(l i ) and C (i) = C(l i ) associated with the various liquid height. For each dry modes u (i) α , the computation of the liquid response in terms of displacement potential field can be obtained by solving the linear system of equation:

H (i) ϕ (i) α = -C (i)T u (i) α (62)
An approximation of u(l i ) = u (i) denoted by u

(i)
approx is constructed as a linear combination of the truncated basis u (i) α :

u (i) approx ≃ m α=1 q (i) α u (i) α = B (i) u q (i) with B (i) u =      . . . . . . u (i) 1 . . . u (i) m . . . . . .      and q =      q (i) 1 . . . q (i) m      (63)
Then, we consider an approximation of the displacement potential ϕ(l i ) = ϕ (i) denoted by ϕ

(i)
approx as a linear combination of the computed ϕ α (l i ) = ϕ (i) α as:

ϕ (i) approx = m α q (i) α ϕ (i) α = B (i) ϕ q (i) with B (i) ϕ =      . . . . . . ϕ (i) 1 . . . ϕ (i) m . . . . . .      (64)
Considering the left hand side of Eq. ( 49), and a left multiplication by B (i)T u , the reduced hydroelastic problem projected of the prestressed dry modes is written as:

B (i)T u K (i) t B (i) u K (i) r q (i) -ρ f ω 2 B (i)T u C (i) B (i) ϕ M (i) ar q (i) -ω 2 B (i)T u M (i) B (i) u M r = I q (i) = 0 (65)
where the matrices K (i) r and M (i) ar are respectively the reduced tangent stiffness and added mass matrices which are (N × N ) matrices defined as:

K (i) r =      . . . ω 2(i) α . . .      and M (i) ar =      . . . . . . m (i) αβ . . . . . .      (66)
where the reduced stiffness and mass matrices are diagonal matrices and the terms of the reduced added mass matrix are defined as:

m (i) αβ = u T(i) α C (i) ϕ (i) β = -u T(i) α C (i) H -1(i) C T(i) u (i) β (67)
Of course, the resulting reduced added mass matrix is symmetric positive definite.

The eigenvalue analysis of the reduced hydroelastic problem writes:

[K (i) r -ω 2 (M (i) ar + I)]q = 0 (68)
and gives the hydroelastic eigenmodes of the liquid structure problem.

Comparison and validation with experiments

Hydroelastic vibrations of a clamped circular plexiglass plate under a column of water

The computation methodology is applied to experiments carried out by M. Chiba [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part II: Linear axisymmetric vibration analysis[END_REF]. It consists in evaluating the linearized hydroelastic vibrations of a thin circulate plate made of plexiglass under a column of water. The great interest of this example is that it highlights the prestressed effects due to the liquid weight on the linearized hydroelastic vibrations. The procedure is illustrated in Fig. 3. It consists in (i) computing the nonlinear static solutions of the elastic structure subjected to the hydrostatic pressure forces and (ii)

establishing the linearized vibrations of the coupled system. The liquid is in contact with the structure upon the upper surface of the clamped plate, while the lower surface is free of charge. The liquid height is varying from l min = 0 mm to l max = 250 mm. The solid material parameters needed for the Saint-Venant Kirchhoff constitutive relations are the Young modulus E = 6.9 × 10 9 Pa, the Poisson ratio ν = 0.38 and the solid mass density ρ s = 1.4 × 10 3 kg.m -3 . The liquid mass density is ρ f = 1.0 × 10 3 kg.m -3 .

Preliminary nonlinear static evolution solution

The Fig. 4 illustrates the mesh parameterization of the elastic plate to compute the nonlinear static problem. The computation has been done with quadratic hexahedron elements, which has been chosen in order to obtain an accurate solution even if the aspect ratio is rather small (e.g. the ratio between the smaller and the larger edges of an element).

In Fig. 5, the quantity of interest considered is the bottom displacement q at the center of the plate in function of the liquid height. At l = 250 mm, the absolute value of the displacement is almost 11 times higher than the thickness of the structure. Considering this order of magnitude the static solution can be considered in a nonlinear regime. An "overkill" solution denoted by q ref , for a liquid height l = 250 mm, has been computed considering a large number of degrees of freedom. Then, simulations have been performed to compare the value of q to the reference as function of the number of elements showing from which the static 225 solution can be considered as independent of the mesh. The table 1 contains all mesh parameters described in Fig. 4. In [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part I: Experiment[END_REF], the deflection value is not measured nor given explicitly given. At this point of the paper, it is not possible to known if the results from the computation are in good agreement with the experiment.

(n θ ,n r ) (2,
Those considerations will be presented later by considering the natural coupled frequencies function of the liquid height. In order to conclude about the static computation, at each liquid height, the nodal static displacement solutions u

s are known. Thus, the geometrical stiffness matrices K (i) g , the material stiffness matrices K (i) m and the coupling matrices C (i) needed for the hydroelastic problem can be constructed and stored at those liquid heights. In the sing, we show how to compute the operators H (i) from a liquid mesh coincident with the solid mesh at the liquid-structure interface.

Prestressed hydroelastic vibrations

The computation of the hydroelastic natural frequencies needs to mesh the liquid at equilibrium for each liquid height. Here, an associated fluid mesh coincident with the current configuration of the structure is generated as follows. At first, the Fig. 6 illustrates the liquid mesh parameterization that we use to generate a liquid mesh with an appropriate liquid height. Then, the given liquid mesh is deformed using a vector Laplacian equation with a prescribed boundary condition at the liquid structure interface. The risk of distortion is really low because the liquid mesh deformation near the liquid-structure interface is negligible.

Considering the liquid mesh domain, the operators H (i) are generated for each liquid height l i . The next step consists in expressing and solving the hydroelastic problem by using the computation methodology presented in Section 4. Considering the geometry of the circular plate, the modal shapes can be characterized by given a couple number (n, m) and a symbol presented in Fig. 7. The number n corresponds to the number of radial waves and m to the angular waves. Those characterizations can still be used for the modal shape of the structural part of hydroelastic modes. In Fig. 8, the hydroelastic natural frequencies are given as function of the liquid height in a frequency range between 0 Hz and 80 Hz. Qualitatively, the graph shows a decrease and then an increase of the natural hydroelastic frequencies for all the modes, except for the first one (0, 1) which have an opposite behavior. The curves are in agreement with the experimental observations presented in reference [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part I: Experiment[END_REF]. We can see that the natural frequencies of the prestressed dry modes tend to increase the natural frequencies in function of the liquid height. Conversely, for an elevation of the liquid height, it is well known that the added mass effect tends to lower the natural hydroelastic frequencies. Therefore, in this example, we have a balance between the prestressed state due to the weight of the liquid and the added mass effect. At first, for an almost empty structure, a decrease of the natural frequencies can be observed, except for the first mode. At a given height, the natural frequencies reach a minimum value before rising again due to the prestressed effect. In Fig. 9, sectional views of the non-dimensional displacement potential fields are represented respectively for the first and second modes. For the second mode, the spatial fluctuation with the behaviour of the evolution of hydroelastic natural frequencies.

The Fig. 10 shows a quantitative comparison between the experiments and the computed natural frequencies which are in very good agreements. When the liquid height l is near zero, differences can be observed.

Many uncontrolled parameters can be the sources of those differences : (i) the prestressed effect due to the plate weight without the liquid could be significant at this stage, (ii) the presence of uncertainties of material parameters used in the simulation or (iii) the sloshing effect at the free-surface. Nevertheless, once the liquid height is elevated, a very good agreement between the simulations and the experiments are observed for the first and the third modes. Concerning the second mode, the behavior is overall the same even if a small discrepancy is observed. Numerical results for three particular hydroelastic frequencies are given for various liquid heights in Table 2. A sharp evolution of frequency between an empty state (l = 0 mm) and a nearly empty state (l = 1 mm) can be observed for the first mode. This behavior is due to the prestress effect as described in [START_REF] Chiba | Nonlinear hydroelastic vibration of a cylindrical tank with an elastic bottom, containing liquid. part II: Linear axisymmetric vibration analysis[END_REF].

Conclusion

In this paper, we have proposed a numerical strategy for the evaluation of hydroelastic vibrations of an elastic tank containing a free surface liquid, around nonlinear prestressed equilibrium configurations due to the liquid weight. The main originality of the approach lies in the resolution of the problem by considering the liquid height as an evolution parameter. The methodology has been validated by comparison of the numerical solutions with experimental results from the literature. The computation uses the prestressed dry modes on a reference configuration. The numerical methodology presented in this paper is valid for compatible fluid-structure interface meshes. The extension of the methodology for incompatible meshes will be the purpose of future investigations. This would allow more complex parametric situations (e.g. influence of ring stiffeners geometry for example). Deeper convergence analyses are under investigations through appropriate hydroelastic projection-based reduced order models for arbitrary complex 3D shapes.

as shown in Fig 1 .

 1 The liquid pressure on the interface only depends on the wetted interface position for a given configuration.

Figure 1 :

 1 Figure 1: Various filling liquid levels in an elastic tank considering an evolution of the liquid height l.

  75

Figure 2 :

 2 Figure 2: Reciprocal mapping of the liquid-structure interface from the current configuration on the reference configuration.

  (38) correspond to the linear Euler equations and Eq. (39) corresponds to the normal displacement continuity between the normal fluid displacement field u f • n ′ and the normal solid displacement field u • n ′ on the liquid-structure interface denoted by Σ ′ i , in which n ′ denotes the external unit normal with respect to the solid. Considering the previous equations and the boundary conditions, expressed in term of potential displacement field ϕ as 140 follows:

α

  = u α (l i ), normalized by the mass (i.e. u associated with the eigenvalues denoted by ω 2(i) α , generate a basis denoted by B (i) u defined as:

Figure 3 :

 3 Figure 3: (a) Computation of the prestressed state of the circular plate due to the liquid weight for various liquid height l i and (b) liquid and structure domains of the linearized hydroelastic vibrations of the prestressed plate.

Figure 4 :

 4 Figure 4: (a) Parameterized geometry with the thickness t = 0.35 mm, the internal radius R int = 144 mm; (b) Parameterized mesh by the numbers n θ and nr, the number of elements in the thickness is fixed as two elements.

Figure 5 :

 5 Figure 5: (Left) Evolution of the deflection bottom function of the liquid height with a thickness ts = 0.35 mm and (Right) convergence analysis of the quantity of interest at l = 250 mm in function of the number of structural degrees of freedom.

Figure 6 :

 6 Figure 6: (a) Parameterized geometry of the fluid and the solid of thickness ts, the internal radius R int and the fluid height l; (b) Parameterized mesh with the numbers n θ and nr.

  255

Figure 7 :

 7 Figure 7: Characterization of the six first modal deformed shapes of the elastic bottom.

Figure 8 :

 8 Figure 8: Evolution of the computed hydroelastic natural frequencies as function of the liquid height.

Figure 9 :

 9 Figure 9: Sectional views of the non-dimensional displacement potential responses for two first modes at various liquid heights.

Figure 10 :

 10 Figure 10: Hydroelastic frequencies as function of the fluid height compared to experimental values extracted from [18].

Table 1 :

 1 Results and mesh parametrization.

		2)	(5, 5)	(10, 10)	(20, 20)	(50, 50)
	n elem	40	250	1000	4000	25000
	n dof	771	4380	16995	66965	264009
	q(l = 0.250 mm)	-3.12 mm	-3.74 mm	-3.84 mm	-3.86 mm	q ref = -3.87 mm
	100× | q-q ref | /q ref	19.3 %	3.4 %	0.77 %	0.25 %	-

Table 2 :

 2 Evolution of hydroelastic frequencies (n, m) with prestress effect from liquid weight for various liquid heights.

	, 1)	17.5 Hz	27.6 Hz	20.7 Hz	19.4 Hz	18.5 Hz	17.8 Hz	16.6 Hz
	(1, 1)	36.6 Hz	34.4 Hz	24.3 Hz	24.0 Hz	25.6 Hz	28.0 Hz	32.6 Hz
	(2, 1)	60.2 Hz	46.0 Hz	31.7 Hz	33.1 Hz	38.3 Hz	43.3 Hz	51.2 Hz

of potential remains localized near the liquid-structure wall. On contrary, the spatial fluctuation of fluid response for the first mode seems to affect the overall fluid domain. Those observations are in agreement


Integration of Integro-differential Equations and Deep Learning

François Lemaire, Louis Roussel CFHP team, CRIStAL Laboratory, Université de Lille, France

Integral equation modelling

Dynamical system ẋ1 (t) = x 2 ẋ2 (t) = θx 1 (t) y(t) = x 2 (t) (known) y(t) -y 0 = θ t 0 τ 2 0 y(τ 1 )dτ 1 dτ 2 + θx 1 (0)t ÿ(t) = θy(t) Integral Elimination Differential Elimination [START_REF] Boulier | The DifferentialAlgebra project[END_REF] [2, Integrate] Algorithm Many models (such as the SIR epidemiology model) consist of non-linear differential equations, involving unknown numerical parameters. Recent work investigates the treatment of integro-differential equations i.e. equations involving both derivatives and integration operators. The introduction of integral equations helps to increase the expressiveness of the models, improve the estimation of parameters and reduce the size of the intermediate equations.

Integration of equations with deep learning [3]

Through learning techniques used for the translation of texts, Lample and Charton [START_REF] Lample | Deep Learning for Symbolic Mathematics[END_REF] succeed in calculating primitives of mathematical expressions, and solutions of differential equations.

The model used is called

Transformer and is widely used in natural language processing [START_REF] Vaswani | Attention Is All You Need[END_REF]:

• It takes a list of words as input and returns the list of translated words as output,

Adaptation of [3] to integro-differential equations

In order to integrate integro-differential equations, we rewrote from scratch the implementation of [START_REF] Lample | Deep Learning for Symbolic Mathematics[END_REF] (based on PyTorch) to fit our needs.

Equation f ′ Solution f 1 x - tx ( x) 2 t x (-a 0 x -1) xy + ÿ(t + a 0 x) + a 0 xy x + xyt -a 0 ẏx - ẏ (t + a 0 x) 2 ẏ + xy t + a 0 x

Retrieve

A B from A ′ B -AB ′ (and not A ′ B-AB ′ B 2

)

The aim of this experiment is to integrate equations that require the introduction of an integrating factor. Datasets are composed of pairs (A ′ B -AB ′ , A/B) where A and B are random irreducible polynomials. Here, we used 420 000 pairs for the training and 10 000 for the evaluation of the accuracy.

Conclusion

We achieve a quite good accuracy on the integration of a specific set of integrodifferential equations. Our next goal will be focused on the integral elimination of differential systems with deep learning.

Bibliography

Figure 1 :

 1 Figure 1: Modelling, elimination and [2, Integrate] algorithm.

Figure 2 :

 2 Figure 2: A signal, its derivative and its integral (with noise in red).

Table 1 : Accuracy of the model versus the accuracy of Maple.

 1 GenerationOur generation process creates unary-binary trees, where nodes are +, -, * , ÷ plus the ∂ and operators, and leaves are x(t), y(t), a 0 and t.

	f : a 0 ẋ(t) + y(t)

Table 2 : Equations that the model can integrate, but not Maple.

 2

Table 3 : Accuracy of the model versus the accuracy of Maple.

 3 -ẏ2 x -x 2 ÿ + 2 ẋ y + ẏ2 x y + y 2 ẏx + ẋ ẏy -y 2 ÿ y -ÿxy -ẋ ẏy y -ÿxy y

	Equation f ′

Table 4 : Equations that the model can integrate, but not Maple.

 4