Charbel Farhat

Edmond Kwan-Yu Chiu

David Amsallem

Jean-Sébastien Schotté

Roger Ohayon

Modeling of Fuel Sloshing and its Physical Effects on Flutter

The sloshing effects of an internal fluid on the flutter envelope of an aeroelastic system have received little attention in the open literature. This issue is nevertheless relevant for many aircraft, especially high-performance fighter jets carrying stores. This paper addresses some aspects of this problem as well as related modeling and analysis issues. These include the importance or insignificance of accounting for the hydroelastic effect when modeling an internal fluid and its container as well as accounting for that container when modeling the aerodynamics of the overall aeroelastic system. The paper also reports on the findings of four independent sets of flutter analyses performed for a wing-store test configuration and various fuel fill levels in the subsonic, transonic, and early supersonic regimes. Two of these sets of numerica l experiments relied on a computational-fluid-dynamics-based computational technology, and two of them on the doublet-lattice method or a supersonic lifting-surface theory, where applicable. The geometry of the chosen test configuration is that of the AGARD Wing 445.6 with a blunt-nosed store. The obtained computational results show that, at least for the considered wing-store configuration, ignoring the aforementioned hydroelastic effect tends to overestimate the added-mass effect and underestimate the critical pressure and flutter speed. They also reveal that, whereas the aerodynamics of the store may be neglected in the subsonic regime, they cannot be ignored in supersonic air streams. Finally, the performed computational study suggests that, in general, the critical pressure and flutter speed decrease with an increasing fuel fill level. T HE numerous effects of an internal fluid on the overall behavior of an aircraft or spacecraft have long been recognized by the aircraft design community [START_REF] Cazier | Structural Dynamic and Aeroelastic Considerations for Hypersonic Vehicles[END_REF][START_REF] Runyan | Flutter at Very High Speeds[END_REF]. In particular, criteria for designing launch vehicles to accomodate the sloshing effect of an internal propellant [START_REF] Braslow | Propellant Slosh Loads[END_REF] have been documented. Unfortunately, whether derived from experimental or computational data [START_REF] Braslow | Propellant Slosh Loads[END_REF][START_REF] Bauer | Theory of Liquid Sloshing in Compartmented Cylindrical Tanks Due to Bending Excitation[END_REF], these criteria did not in general incorporate the flexibility effect of the fluid container. For this reason, these criteria do not address the flutter issue, which is often finally resolved by a flutter flight test.

The cost of a flutter flight test campaign increases rapidly with the number of configurations requiring certification [START_REF] Dawson | F-16 Limit Cycle Oscillation Testing with Asymmetric Store[END_REF][START_REF] Harris | Accelerated Flutter Testing: A Possible Technique to Expedite Flutter Testing for New Aircraft[END_REF]. In the presence of stores containing an internal fluid, the fill level becomes an additional configuration parameter. Perhaps because of the already large number of other parameters to consider and the already substantial cost associated with flight testing, it seems that there is very little, if any, focus on the effects of the interaction between an internal fluid and its containing structure during flutter flight tests.

There have been, however, some experimental attempts to investigate the effects of an internal fluid on flutter. In particular, it was observed in [START_REF] Peterson | Nonlinear Fluid Slosh Coupled to the Dynamics of a Spacecraft[END_REF] that the dynamic response of a coupled fluid-tank system can exhibit modal resonances that are not predicted by any uncoupled analysis. It was also reported in [START_REF] Ferman | Fluid-Structure Interaction Dynamics in Aircraft Fuel Tanks[END_REF] that the presence of fuel between wing panels can alter their flutter characteristics. Little computational research has been performed, however, for predicting the effects of the presence of an internal fluid on flutter [START_REF] Amsallem | Stabilization of Projection-Based Reduced-Order Models[END_REF][START_REF] Firouz-Abadi | Effect of Fuel Sloshing in the External Tank on the Flutter of Subsonic Wings[END_REF].

Therefore, the main objective of the work reported here is to numerically predict the effects on flutter of the interaction between an internal fluid and its containing flexible structure. To this effect, two different representations of the internal fluid are considered. The first one aims at capturing accurately the hydroelastic vibrations of such a fluid and its store. The second one aims instead at reducing model preparation time and leads to considering other assumptions such as, for example, ignoring the aerodynamics of the store. In all cases, different computational technologies are considered for representing the external fluid. These are the versatile linearized computational fluid dynamics (CFD) approach, and the doublet-lattice method and a lifting-surface theory, which remain popular for compressible subsonic flows and supersonic ones, respectively, particularly in the design environment. Consequently, a secondary objective of this work is to report on the impact of these modeling issues on the prediction of the sloshing effects on flutter.

To this effect, the remainder of this paper is organized as follows.

In Sec. II, two different computational models are described for representing an internal fluid. In Sec. III, a higher-fidelity computational framework for aeroelasticity based on linearized CFD and lower-fidelity alternatives based on the linear aerodynamic theory are overviewed. In Sec. IV, the results of the flutter analysis of a wing-store geometry based on the AGARD Wing 445.6 and the various modeling approaches outlined previously are reported and discussed for various fuel fill levels. Finally, Sec. V concludes this paper.

Remark: In general, slosh refers to the movement of an internal fluid with a free surface in the presence of gravity [START_REF] Morand | Fluid Structure Interactions: Applied Numerical Methods[END_REF]. However, such a slosh occurs at frequencies that are too low to be of concern to flutter. Therefore, in this work, the effect of gravity on the movement of the internal fluid is ignored, but for the sake of simplicity, the movement of that fluid is still referred to as slosh.

II. Mathematical Models for an Internal Fluid

Here, two approaches of different fidelity levels are considered for accounting for the effects of an internal fluid in flutter prediction. The first approach is based on the theory of hydroelastic vibrations of an elastic structure containing an inviscid fluid. It results in an addedmass model that captures the hydroelastic effect [START_REF] Morand | Fluid Structure Interactions: Applied Numerical Methods[END_REF]. The second approach accounts only for the mass and inertia of an internal fluid and leads to a trivial added-mass model [START_REF] Schotté | Various Modelling Levels to Represent Internal Liquid Behaviour in the Vibratory Analysis of Complex Structures[END_REF]. It minimizes model preparation time. It also leads to considering the possibility of neglecting the aerodynamics of the container, which can be attractive for parametric store-configuration studies. Both of these models can, however, be easily integrated in existing aeroelastic fluttercomputation procedures.

A. Hydroelastic Vibrations Model

Here, the flexible store and the fluid it may contain are viewed as a coupled hydroelastic system. The internal fluid occupying a domain Ω F with a free surface Γ and a fluid-structure interface boundary Σ (see Fig. 1) is assumed to be inviscid, incompressible, homogeneous, weightless, and without surface tension. Hence, its dynamic behavior is governed by the linearized Euler equations

∇p -ρ F ∂ 2 u F ∂t 2 (1)
where p, ρ F , and u F denote the Eulerian pressure fluctuation, density, and displacement of the fluid, respectively, and t denotes time. Furthermore, assuming that it is characterized by small displacement amplitudes, the fluid motion can be assumed to be harmonic and is therefore represented by

u F u F a cos ωt (2a) p p a cos ωt (2b)
where the subscript a designates the amplitude, and ω denotes the circular frequency. From Eqs. ( 1), (2a), and (2b), it follows that

p ρ F ω 2 ϕ π (3)
where ϕ is the potential function satisfying ∇ϕ u F , and π is an arbitrary constant that is set here to π 0. After some manipulations, the incompressibility condition ∇ • v F x 0 ;t0, where v F denotes the fluid velocity and x 0 x u F x; t, leads to

Δϕ 0i nΩ F ( 4
)
The structural vibrations are also assumed to be harmonic with small amplitudes. Hence, they are represented by the dynamic equations of equilibrium σ S ij;j u S ρ S ω 2 u S i 0i nΩ S (5a)

σ S ij u S n S j F i on ∂Ω S \ Σ (5b)
where σ S ij denotes the stress tensor components (i 1; 2; 3, j 1; 2; 3) of the structure Ω S ; ρ S and u S denote the density and displacement of the structure, respectively; a comma denotes a partial spatial derivative; the subscripts i and j designate the ith and jth components of a quantity in three dimensions, respectively; n S denotes the normal to the surface of the structure ∂Ω S ; and F denotes the prescribed external force field acting on the structure. The transmission conditions coupling the fluid and structure components of the hydroelastic system at the fluid-structure interface Σ can be written as

σ S ij u S n Sj ρ F ω 2 ϕn F i ( 6
)
where n F denotes the normal to the fluid-structure interface Σ from the fluid side, and

∂ϕ ∂n F u S • n F (7)
The first of the previous two equations expresses the conservation of momentum at Σ; the second one is equivalent to a nonpenetration condition at this interface. Finally, on the free surface Γ of the fluid (see Fig. 1), pressure invariance yields the following boundary condition:

p ρ F gu F 3 0 (8)
(as ρ F g 0 due to the weightless assumption) where the subscript 3 designates the coordinate in the vertical direction. This yields ϕ 0 on the surface; see Eq. ( 3). In summary, the boundary-value problem (BVP) governing the vibrations of a system comprising a store and an internal fluid is given by

Δϕ 0i nΩ F ϕ 0o nΓ ∂ϕ ∂n F u S • n F on Σ σ S ij u S n Sj ρ F ω 2 ϕn F i on Σ σ S ij;j u S ρ S ω 2 u Si 0i nΩ S σ S ij u S n Sj F i on ∂Ω S \ Σ
and the pressure p of the internal fluid is obtained from Eq. (3). Note that, with the choice π 0, the previous BVP remains well-posed [START_REF] Morand | Fluid Structure Interactions: Applied Numerical Methods[END_REF] for ω 0.

The finite-element discretization of the previous BVP leads to the following system of matrix equations:

K -λM S -λT -T T M F u Φ 0 ( 9
)
where λ ω 2 , K, and M S denote the standard structural stiffness and mass matrices; M F is the standard internal fluid mass matrix; T is the matrix arising from the finite-element discretization of the transmission conditions [Eqs. ( 6) and ( 7)]; u and Φ denote the vectors of structural displacements and fluid potential degrees of freedom, respectively; the subscript specifies that the finite-element space in which ϕ is approximated satisfies the Dirichlet boundary condition ϕ 0 on Γ; and the superscript T designates the transpose operation. Because M F is nonsingular, the second row of the previous matrix equation gives

Φ M -1 F T T u ( 10
)
which transforms the system of Eq. ( 9) into the following simpler equation: Ku -λM S M A u 0

where

M A T M -1 F T T (12)
is the "added-mass " finite-element matrix. In the remainder of this paper, Eq. ( 11) is referred to as the consistent or hydroelastic addedmass model of the flutter problem in the presence of an internal fluid.

B. Frozen-Mass Model

A trivial approximation of the added-mass effect of an internal fluid is obtained by lumping at one or more points the mass and inertia of this fluid. Compared to the previous modeling approach, this simpler one reduces model preparation time because it eliminates the need to create a new mesh for the internal fluid for every fill level of interest. It also leads to considering the possibility of neglecting the aerodynamics of the store in the flutter analysis. The latter assumption enables the application of popular supersonic liftingsurface theories that nevertheless lack the capability of handling slender bodies to modeling the unsteady external aerodynamics of the given aeroelastic system; for example, see Sec. III.B.2. More importantly, this further assumption simplifies the flutter analysis of parametric store configurations and therefore the modeling and simulation of store clearance, and it significantly reduces the associated model preparation time.

For convenience, the previous approach is implemented in this work by treating the fluid as an almost-incompressible medium with the same total mass and density (and hence volume) as the original body of fluid. The shape and location of this medium are chosen to be those assumed by the internal fluid when the store is completely horizontal and at rest. The Poisson ratio of this medium is chosen to be close enough to the incompressible limit of 0.5, while still preventing finite-element locking during any phase of the numerical analysis. Its elastic modulus is chosen to be significantly lower than that of the containing structure.

III. Computational Frameworks for Flutter Analysis

Two well-proven computational frameworks of different level of fidelity are also considered for analyzing the effect of fuel sloshing on flutter. The first one is a higher-fidelity framework because it is based on the Euler equations for modeling the external flow and the linearized three-field arbitrary Lagrangian-Eulerian (ALE) [START_REF] Lesoinne | A Linearized Method for the Frequency Analysis of Three-Dimensional Fluid/ Structure Interaction Problems in All Flow Regimes[END_REF][START_REF] Lesoinne | A CFD Based Method for Solving Aeroelastic Eigenproblems in the Subsonic, Transonic, and Supersonic Regimes[END_REF][START_REF] Lieu | Reduced-Order Fluid/Structure Modeling of a Complete Aircraft Configuration[END_REF][START_REF] Farhat | Mixed Explicit/Implicit Time Integration of Coupled Aeroelastic Problems: Three-Field Formulation, Geometric Conservation and Distributed Solution[END_REF] formulation of fluid-structure interaction for modeling aeroelastic effects. The second computational framework is a lower-fidelity but more computationally economical one. In the subsonic domain, it is based on the standard doublet-lattice method (DLM), which uses the linearized aerodynamic potential theory [START_REF] Albano | A Doublet-Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows[END_REF][START_REF] Rodden | Refinement of the Nonplanar Aspects of the Subsonic Doublet-Lattice Lifting Surface Method[END_REF] to predict the aerodynamic loads on lifting surfaces and slender bodies. In the supersonic regime, it is based on a similar linear theory [START_REF] Liu | Further Studies of Harmonic Gradient Method for Supersonic Aeroelastic Applications[END_REF] that is, however, limited to planar lifting surfaces and therefore cannot handle stores. This theory and associated software are most convenient for studying here the effect of neglecting the aerodynamics of a store on an overall flutter analysis. The higher-fidelity computational framework is reviewed in more detail than the lower-fidelity one because it is not as popular yet.

In both computational frameworks, finite-element modal reduction is used for representing a structure.

A. Linearized Computational-Fluid-Dynamics-Based Computational Framework

This framework is considered here for four main reasons. 1) It is based on a perturbation theory, and flutter is essentially a perturbation problem.

2) Unlike linear methods, it incorporates the effects of shocks in the transonic regime and is valid in all of the subsonic, transonic, and supersonic regimes.

3) Unlike nonlinear methods, it is fast enough to appeal to production environments.

4) It is compatible with the added-mass model of the internal fluid presented in Sec. II.A because the hydroelastic model itself is also a linearized model based on a perturbation theory [START_REF] Morand | Fluid Structure Interactions: Applied Numerical Methods[END_REF].

Three-Field Formulation of Aeroelastic Problems

The three-field formulation of nonlinear computational aeroelasticity introduced more than a decade ago [START_REF] Lesoinne | Stability Analysis of Dynamic Meshes for Transient Aeroelastic Computations[END_REF] models a fluidstructure interaction problem by three coupled partial differential equations: one governing the fluid subsystem and written in an ALE coordinate system to facilitate the treatment of moving boundary surfaces, one governing the motion of the fluid mesh induced by the motion of boundary surfaces, and one governing the dynamic equilibrium of the structural subsystem driving the motion of these boundary surfaces. In the case of inviscid flows, these three equations can be written as

∂Jw ∂t ξ J∇ x • F w - ∂x F ∂t w 0 (13a) ρ S ∂ 2 u S ∂t 2 -divσ S ϵ S u S b (13b) ρ ∂ 2 x F ∂t 2 -div Ẽ∶ εx F 0 (13c)
where ξ and x F t denote the reference position of a fluid grid point and the time-dependent displacement of a fluid grid point, respectively; J detdx F ∕dξ is the Jacobian of the frame transformation ξ → x F ; F denotes the convective flux associated with the fluid state vector w; ϵ S is the strain tensor of the structure; b denotes the body forces acting on the structure; and ρ, Ẽ, and ε denote the ficticious density, elasticity properties, and strain tensor of the pseudostructural system adopted for modeling the moving fluid grid, respectively. All other variables have the same meaning as before.

The finite volume (fluid) and finite-element (structure) semidiscretizations of the previous equations lead to Ax F w ;t Fw; x F ; x F;t 0 (14a)

M S u ;tt f int u; u ;t -f ae u; w0 (14b) Kx F -Kc u 0 (14c)
where A denotes the diagonal matrix of cell volumes; w and x denote the fluid state and fluid grid position semidiscrete vectors, respectively; a subscript ;tdesignates a time derivative; F denotes the vector of numerical convective flux functions; M S is the finiteelement mass matrix of the structure; u is the semidiscrete displacement vector; f int and f ae denote the vectors of internal and external aerodynamic forces acting on the structure, respectively; K and Kc are pseudostiffness matrices; and Kc accounts for the effect of the structural motion on the fluid mesh at the interface between the structure and the external fluid [START_REF] Farhat | Load and Motion Transfer Algorithms for Fluid/Structure Interaction Problems with Non-Matching Discrete Interfaces: Momentum and Energy Conservation, Optimal Discretization and Application to Aeroelasticity[END_REF].

Linearization Around an Equilibrium Configuration

Consider an aeroelastic system in an equilibrium configuration (w o , w ;t o 0, x Fo , x F;t o 0, u o , u ;t o 0, u ;tt o 0) designated by the subscript o. To rapidly compute the response of this system to a perturbation of the form (δw, δw ;t , δx F , δx F;t , δu, δu ;t , δu ;tt ), which can be used to identify the aeroelastic parameters of this system and determine its stability or flutter, the problem unknowns are expanded as follows:

w w o δw; w ;t δw ;t x F x F o δx F x F;t δx F;t (15) u u o δuu ;t δu ;t u ;tt δu ;tt (16)
and Eqs. (14a-14c) are linearized as follows.

First, the fluid equations are linearized around the equilibrium configuration and rewritten in dimensionless form following the approach first described in [START_REF] Lesoinne | A Linearized Method for the Frequency Analysis of Three-Dimensional Fluid/ Structure Interaction Problems in All Flow Regimes[END_REF][START_REF] Lesoinne | A CFD Based Method for Solving Aeroelastic Eigenproblems in the Subsonic, Transonic, and Supersonic Regimes[END_REF]. This transforms Eq. (14a) into

A o δ _ w H o δ w E o C o δ _ x F G o δ x F 0 (17)
where

A o A x F o H o ∂ F ∂ w w o ; x F o E o ∂ A ∂ x F x F o w o C o ∂ F ∂ _ x F w o ; x F o G o ∂ F ∂ x w o ; x F o
Here, and throughout the remainder of this section, the bar notation indicates that a quantity is dimensionless, and a dot designates the derivative with respect to the dimensionless time:

τ V ∞ L o t ( 18
)
where L o denotes a reference length, and V ∞ is the freestream velocity.

The

M S δu ;tt D o δu ;t K 0 δu P o δ w (19)
where

K o ∂f int ∂u u o K 0 K o - ∂f ae ∂u w o ; u o D o ∂f int ∂u ;t u o P o ∂f ae ∂ w w o ; u o
To keep the notation as compact as possible, the subscript o and the prefix δ are dropped in the remainder of this paper. Hence, the same variables w, x F , and u are used next to denote the perturbations around the chosen equilibrium configuration of the fluid state, fluid mesh motion, and structural motion vectors, respectively. Then, the dimensionless fluid mesh position and velocity variables x and _

x are eliminated from the coupled system of linearized equations by introducing

KF K-1 Kc (20) so that x F KF u _ x F KF _ u ( 21
)
where u u L o ( 22
)
and u is measured with respect to the equilibrium configuration. The previous algebraic manipulations allow rewriting Eq. ( 17) as

A _ w H w E C KF _ u G KF u 0 (23)
Next, neglecting the effects of ∂f ae ∂u j o and D o , that is assuming that K 0 ≈ K o and no structural damping, Eq. ( 19) is projected onto a modal basis V of n m dry, natural, structural modes and therefore is transformed into

Iu m;tt Ω 2 u m P m w (24)
Here, I denotes the identity matrix, Ω 2 is a diagonal matrix containing the squares of n m natural circular frequencies of the structure, and the subscript m designates generalized modal quantities, e.g., u m denotes the generalized modal displacement coordinates and P m w the generalized aerodynamic forces satisfying

u L o u Vu m P m V T P (25)
In summary, Eqs. ( 23) and ( 24) constitute the governing linearized equations of motion of an aeroelastic system. They can be solved either in the time or frequency domain to determine the flutter characteristics of this system. Both computational approaches are briefly reviewed next.

Flutter Analysis in the Time Domain

To perform the flutter analysis of a given aeroelastic system, a pertinent initial perturbation can be defined, and a numerically stable and provably second-order time-accurate staggered solution procedure [START_REF] Lesoinne | A Higher-Order Subiteration Free Staggered Algorithm for Nonlinear Transient Aeroelastic Problems[END_REF][START_REF] Farhat | Two Efficient Staggered Procedures for the Serial and Parallel Solution of Three-Dimensional Nonlinear Transient Aeroelastic Problems[END_REF][START_REF] Farhat | Provably Second-Order Time-Accurate Loosely-Coupled Solution Algorithms for Transient Nonlinear Computational Aeroelasticity[END_REF] can be applied to the time integration of the coupled fluid [Eq. ( 23)] and structural [Eq. [START_REF] Farhat | Provably Second-Order Time-Accurate Loosely-Coupled Solution Algorithms for Transient Nonlinear Computational Aeroelasticity[END_REF]] semidiscrete equations of equilibrium. In this case, analyzing the predicted response of the aeroelastic system by a modal identification procedure, for example, an eigenrealization algorithm such as ERA [START_REF] Juang | An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction[END_REF], determines the flutter characteristics of the aeroelastic system.

Flutter Analysis in the Frequency Domain

Equations ( 23) and ( 24) can be transformed into the frequency domain and manipulated to define a generalized aerodynamic force matrix. Then, the flutter analysis of the aeroelastic system of interest can be performed in the frequency domain by substituting this matrix into the frequency domain counterpart of Eq. ( 24) and applying the pk method [START_REF] Frazer | The Flutter of Aeroplane Wings[END_REF][START_REF] Hassig | An Approximate True Damping Solution of the Flutter Equation by Determinant Iteration[END_REF] to the resulting dynamic equation as described next.

Let

k L o V ∞ ω ( 26
)
denote the dimensionless reduced frequency. From this definition and that of the dimensionless time [Eq. [START_REF] Rodden | Refinement of the Nonplanar Aspects of the Subsonic Doublet-Lattice Lifting Surface Method[END_REF]], it follows that

ωt kτ (27)
Let also

f ae m P m w (28)
denote the generalized aerodynamic forces appearing in Eq. ( 24). In the frequency domain, w, u m , and f ae m can be expressed as

w w a e ikτ u m u a m e ikτ f ae m f ae a m e ikτ ( 29
)
where i is the pure imaginary number satisfying i 2 -1; and w a , u a m , and f ae a m denote the amplitudes of w, u m , and f ae m , respectively. Substituting the first two of Eq. ( 29) into Eq. ( 23) and using the first of Eq. ( 25) leads to

w a k- 1 L o ik A H -1 ik E C G KF Vu a m ( 30
)
From Eq. ( 28), the second of Eq. ( 25), the third of Eq. ( 29), and Eq. ( 30), it follows that

f ae a m kV T P w a k- 1 L o V T Pik A H -1 × ik E C G KF Vu a m Q m u a m ( 31
)
where

Q m k- 1 L o V T Pik A H -1 ik E C G KF V (32)
is the generalized aerodynamic force complex matrix of size n m equal to the dimension of the structural modal basis. Now, from definition [START_REF] Rodden | Refinement of the Nonplanar Aspects of the Subsonic Doublet-Lattice Lifting Surface Method[END_REF], the second of Eq. ( 29), the definition of the reduced frequency [Eq. ( 26)], and Eq. ( 27), it follows that

u m ;t iωu m (33)
Next, decomposing Q m into its real part Q R m and imaginary part Q I m and using the previous identity leads to

P m w f ae m Q m u m Q R m u m iQ I m u m Q R m u m 1 ω Q I m u m ;t (34)
which transforms Eq. ( 24) into

Iu m;tt - 1 ω Q I m u m ;t Ω 2 -Q R m u m 0 (35)
In summary, starting from Eqs. ( 23) and [START_REF] Farhat | Provably Second-Order Time-Accurate Loosely-Coupled Solution Algorithms for Transient Nonlinear Computational Aeroelasticity[END_REF], which are based on the linearized CFD technology outlined in Sec. III.A.2, Eq. ( 35) is derived. Then, the flutter analysis of a given aeroelastic system can be performed in the frequency domain simply by applying the pk method [START_REF] Frazer | The Flutter of Aeroplane Wings[END_REF][START_REF] Hassig | An Approximate True Damping Solution of the Flutter Equation by Determinant Iteration[END_REF] to this equation.

Incorporation of the Internal Fluid Models

Equation ( 23) governs the external flow and therefore is valid independently from the presence or absence of an internal fluid in the structure. However, if the structure contains an internal fluid, Eq. ( 24) can be modified as follows. First, the added-mass matrix M A [Eq. ( 12)] or its counterpart obtained with mass and inertia lumping can be assembled into the structural mass matrix M S . Then, the first n m natural modes of the structure containing the internal fluid, that is the structure characterized by the finite-element matrices M S M A and K, are computed, and Eq. ( 19) is projected onto the modal basis defined by these modes to obtain a governing equation similar to Eq. ( 24).

Implementation in AERO

Both time domain and frequency domain aspects of the linearized CFD-based computational framework for flutter analysis described previously are implemented in the CFD-based aeroelastic software platform AERO [START_REF] Geuzaine | Aeroelastic Dynamic Analysis of a Full F-16 Configuration for Various Flight Conditions[END_REF][START_REF] Farhat | Application of a Three-Field Nonlinear Fluid-Structure Formulation to the Prediction of the Aeroelastic Parameters of an F-16 Fighter[END_REF]. This software has been validated for several wind-tunnel and flight-test configurations, including the AGARD Wing 445.6 and the F-16 Block 40 aircraft [START_REF] Geuzaine | Aeroelastic Dynamic Analysis of a Full F-16 Configuration for Various Flight Conditions[END_REF][START_REF] Farhat | Application of a Three-Field Nonlinear Fluid-Structure Formulation to the Prediction of the Aeroelastic Parameters of an F-16 Fighter[END_REF]. Its time-domain capability is used in this work to perform all reported flutter analyses.

Implementation in REELC

The frequency domain component of the linearized CFD-based computational framework for flutter analysis described previously is also implemented in the code REELC [START_REF] Mortchéléwicz | Aircraft Aeroelasticity Computed with Linearized RANS Equations[END_REF]. Hence, this code is also used in this work to perform and verify all reported flutter analyses.

B. Non-Computational-Fluid-Dynamics-Based Computational Framework

For both purposes of verification and compare-and-contrast, the doublet-lattice method (DLM) and ZONA51 lifting-surface theory are also considered in this work for predicting the sloshing effects of an internal fluid on the flutter envelope of an aeroelastic system. These two methods, which operate in the frequency domain, are complementary representatives of non-CFD-based computational frameworks. They are equipped here with the pk method for flutter analysis [START_REF] Frazer | The Flutter of Aeroplane Wings[END_REF][START_REF] Hassig | An Approximate True Damping Solution of the Flutter Equation by Determinant Iteration[END_REF]. The DLM is commonly used in the aircraft industry for load computations in the subsonic regime, and the ZONA51 lifting-surface theory is typically used for similar computations in supersonic air streams. Both methods can account for interference among multiple lifting surfaces and bodies. However, ZONA51 assumes planar lifting surfaces. Because both methods are well covered in the literature [START_REF] Liu | Further Studies of Harmonic Gradient Method for Supersonic Aeroelastic Applications[END_REF][START_REF] Rodden | Aeroelastic Analysis -User's Guide[END_REF][START_REF] Blair | A Compilation of the Mathematics Leading to the Doublet Lattice Method[END_REF], they are not overviewed next. Instead, only the ability or inability of these methods to handle aeroelastic configurations with stores and fuel tanks is briefly discussed.

Doublet-Lattice Method for Subsonic Flows

To represent the lifting characteristics of bodies like a fuselage, nacelle, or an external store, the slender-body theory has been associated with the DLM [START_REF] Giesing | Subsonic Steady and Oscillatory Aerodynamics for Multiple Interfering Wings and Bodies[END_REF] in the subsonic domain. The primary wing-body interference is approximated by a system of images of the DLM trailing vortices and doublets within a cylindrical interference body that circumscribes each slender body. The secondary wingbody interference that results from the DLM-bound vortices and doublets is accounted for by a line of doublets located on the longitudinal axis of each slender body. The boundary conditions of no flow through the lifting surfaces or through the body (on the average about the periphery) lead to the equations for the lifting pressure on the surfaces and longitudinal (and/or lateral) loading on the bodies in terms of the normal washes on the wing-body combination.

ZONA51 Lifting-Surface Theory for Supersonic Flows

Unlike the DLM, ZONA51 has no extension to slender bodies. Consequently, it cannot handle configurations with stores and fuel tanks. For this reason, it is used in this work primarily for the purpose of verification when assessing the validity of ignoring the store aerodynamics during the flutter analysis of a wing-store configuration with an internal fluid.

Implementation in NASTRAN and Incorporation of the Internal Fluid Models

The DLM and ZONA51 packages used in this work are those implemented in the SOL145 "Dynamic Flutter Analysis" of the commercial MSC NASTRAN software. For the purpose of this study, the hydroelastic and frozen-mass models described in Secs. II.A and II.B, respectively, were incorporated in both methods by using NASTRAN's programming language DMAP (from "Direct Matrix Abstraction Program") to superpose the added-mass matrix of the internal fluid to the structural mass matrix. Alternatively, the standard virtual mass approach available in NASTRAN could have been used for the same purpose. It was recently verified (see [START_REF] Schotté | Effect of internal liquids on the vibrations of aerospace structures[END_REF]) that both aforementioned approaches give the same results. For both DLM and ZONA51, MSC NASTRAN's implementation of the pk method is used for flutter analysis.

IV. Flutter Analysis of a Wing-Store Configuration at Different Fill Levels

To develop an insight into the effect of fuel sloshing on the flutter envelope of an aircraft, the flutter analysis of a hypothetical wing-store configuration filled with JP-8 is considered in this section. The numerical results obtained for different flight regimes and various store fill levels are reported and discussed.

A. Wing-Store Test System

The hypothetical aeroelastic system considered herein consists of the assembly of the 2.5 ft weakened model 3 of the AGARD Wing 445.6 [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF] and a store with a 3∕32-in.-thick aluminum (E 1.44 × 10 10 lb∕ft 2 , ν 0.33, ρ S 5.24 slugs∕ft 3 ) shell. The internal fluid has the same density as that of JP-8 fuel (ρ F 1.5579 slugs∕ft 3 ). The wing and the store are connected by two sets of linear multipoint constraints (see the Appendix for a detailed description of the geometry).

An undamped finite-element model of the wing with 800 triangular composite shell elements and 2646 DOF is constructed using the information given by Yates [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF]. It predicts natural mode shapes and frequencies that are similar to those derived experimentally. The hypothetical store is discretized into 1408 triangular shell elements and 4188 DOF. Figure 2 shows the finiteelement mesh of the geometry. Between parentheses: percent change with respect to the 0% store fill level. Between parentheses: percent change with respect to the 0% store fill level.

B. Simulated Ground Vibrations

Because in both computational frameworks adopted in this work the structure is represented by its natural modes, an eigenanalysis of the aeroelastic system described previously is first performed for various store fill levels. In general, the results obtained using AERO and NASTRAN are in good agreement, with relative discrepancies lower than 7% for all frequencies. For this reason, the analysis reported next refers only to the AERO computational results (see Tables 1 and2). However, the NASTRAN counterparts are also provided for reference (see Tables 3 and4).

Table 1 reports the first natural frequencies computed by AERO and listed using their physical description. The reader can observe that these two models of the internal fluid lead to a general trend of decreasing frequencies with increasing store fill levels. However, both models lead to results that also differ in an interesting manner. When the hydroelastic model of the internal fluid is used, the computed frequencies decrease monotonically with an increasing fill level. However, when the frozen-mass model is used, some frequencies, namely those associated with the third flexional and torsional modes, reach their minima at intermediate fill levels. Moreover, one can see that the store fill level affects some modes more significantly than others. As this fill level increases from 0 to 99%, the natural frequencies corresponding to the second bending mode and store-roll modes drop by more than 40%, whereas those corresponding to the other modes drop by as little as 3%.

Table 2 reports the relative discrepancies between those eigenfrequencies computed using the hydroelastic and frozen-mass models. Except for certain frequencies that are identically computed (within 0.5% discrepancy) using either model, the frozen-mass model tends to lead to lower frequencies than the hydroelastic model for the same fill level. The largest discrepancies are observed at intermediate fill levels (e.g., 30.9, 50, and 69.1%). In any case, the reader can observe that freezing the mass of the internal fluid instead of accounting for its small motion relative to the containing structure leads to the overestimation of the added-mass effect. Moreover, the hydroelastic model essentially predicts a monotonic decrease in the frequency for the first three torsional/bending modes, the coplanar mode, and the store-roll mode (which, as its name suggests, is predominantly characterized by the displacement of the store). This suggests that the effect of the internal fluid increases when the store fill level is increased. However, the frozen-mass model predicts that the presence of the internal fluid alters most significantly the frequency of the store-roll mode at the 50% fill level. It is furthermore observed that the largest differences between the frequencies obtained using the two internal fluid models occur for the store-roll mode. This is expected because the representation of the roll inertia of the fluid differs significantly from one model to the other.

C. Flutter Analysis

Next, the results of the flutter analysis performed on the wing-store test configuration using both computational frameworks described previously are reported and discussed. In both cases, the structure is represented by its first 10 natural modes. Flutter predictions are performed for various freestream Mach numbers spanning the subsonic, transonic, and early supersonic flow regimes, following the guidelines established by Yates [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF] for the stand-alone AGARD Wing 445.6 configuration (see Table 5).

AERO Results

First, the CFD-based results obtained using the higher-fidelity code AERO are presented. Table 6 reports, for both internal fluid models considered in this work and various freestream Mach numbers M ∞ , the variation with the store fill level of the critical Regardless of the internal fluid model used, the critical flutter speed, and therefore also the critical pressure, predicted by AERO are found to generally decrease with an increasing fill level, meaning that flutter onsets earlier at higher fill levels. Also, these critical quantities are found to change most significantly as the fill level varies from 30 to 70%. Exceptions are noticeable, however, at the two supersonic Mach numbers M ∞ 1.141 (see Fig. 3) and M ∞ 1.201, where AERO predicts that the aforementioned critical quantities significantly increase when the fill level is increased from less than 20 to about 31%, before they decrease when the fill level is further increased.

Overall, the reader can observe that the choice between the two internal fluid models considered in this work has little impact on the critical quantities predicted by AERO. In most cases, the difference between the critical velocity values, and therefore critical pressure values, computed using both models are within a few percent of the corresponding zero-fill reference values (see Fig. 3). Nevertheless, large relative differences are found for the case M ∞ 1.201, where the reference critical quantities are very small.

To characterize further the dependence on the flight conditions of the impact of the internal fluid model, the following quantity: Given that this model is an approximation of the hydroelastic model, the following comment can be reasonably added. From a design perspective, these findings imply that the frozen-mass model could be treacherous because, if relied on, it could lead to a design configuration that may experience unexpected flutter at certain flight conditions. Next, Fig. 5 reports the computed variations of the flutter speed index with the freestream Mach number for the half-full store and almost-full store cases, and both internal fluid models. Here, the flutter speed index is defined as

fV cr ∞ V cr ∞;frozen -V cr ∞;hydroelastic V cr ∞;hydroelastic
FSI V cr ∞ b s ω α μ p ( 36
)
where ω α is the first dry torsional mode of the structural system; b s is the semichord of the wing at the root; and μ m S ∕ρ ∞ V, where ρ ∞ is the freestream density of the test medium. In the absence of a store, V is the volume of a conical frustrum with lower and upper base diameters equal to the streamwise root and tip chords of the wing, respectively, and height equal to the wing semispan; and m S is the structural mass of the wing [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF]. In the context of this work, the volume of the smallest cylinder enclosing the store is added to V, and the structural mass of the empty store is added to m S . The reader can observe that, at certain flight conditions (e.g., M ∞ 0.96 and M ∞ 1.072), the flutter speed index is affected by the choice of the internal fluid model as much as, if not more than, taking into account the presence of an internal fluid in the flutter analysis.

REELC Results

Next, the computational results obtained for supersonic airstreams using the higher-fidelity code REELC are presented in Table 7. These results are in good agreement with those predicted by the code AERO (see Fig. 6).

NASTRAN Results

Finally, the flutter analysis results obtained using NASTRAN are summarized in Table 8. For subsonic airstreams, the comparison of the critical freestream velocities obtained using the DLM and AERO is performed in Fig. 7. The reader can observe that AERO and NASTRAN predict similar results that lead to the same conclusions.

In supersonic airstreams, the comparison of the results obtained using NASTRAN (ZONA51), AERO, and REELC is performed in Fig. 6. Here, the reader can observe that the lower-fidelity ZONA51 method cannot reproduce the results of the higher-fidelity CFDbased codes AERO and REELC. The discrepancies are due primarily to the absence of the aerodynamics of the tank from the ZONA51 computational model due itself to the inability of this method to deal with slender bodies. (This explanation was verified by repeating the flutter analysis in supersonic airstreams for the computational Between parentheses: percent change with respect to 0% store fill level. aeroelastic model excluding the aerodynamics of the tank using AERO and REELC, and reproducing the results obtained using ZONA51.) On the other hand, the reader can observe once again that AERO and REELC deliver relatively close predictions. This suggests that the aerodynamics of a store cannot be, in general, ignored when performing the aeroelastic analysis of a wing-store configuration.

Further Discussion of the Impact of the Chosen Model for the Internal Fluid

Figure 8 highlights the differences found between the critical freestream velocities predicted using a hydroelastic model for the internal fluid or a frozen-mass simplification and the various software packages considered in this work. These differences tend to be .9% fill 0.499 Hydroelastic 669.1 666.0 (-0.5) 652.3 (-2.5) 628.9 (-6.0) 606.6 (-9.3) 584.7 (-12.6) 0.499 Frozen mass 669.1 664.3 (-0.7) 647.6 (-3.2) 616.5 (-7.9) 583.8 (-12.7) 558.2 (-16.6) 0.678 Hydroelastic 919.8 916.0 (-0.4) 903.8 (-1.7) 881.3 (-4.2) 860.9 (-6.4) 838.6 (-8.8) 0.678 Frozen mass 919.8 913.9 (-0.6) 897.9 (-2.4) 871.3 (-5.3) 840.2 (-8.7) 809.9 (-11.9) 0.901 Hydroelastic 1156. smaller in the subsonic regime (below 6%) and larger in supersonic airstreams (above 10%). Interestingly, however, except for the case of the fill level of 99% at M ∞ 1.141, these differences have always the same negative sign, at least for the wing-store configuration considered in this work. This suggests that the frozen-mass assumption is most often a conservative one, in the sense that it leads to smaller critical freestream velocities than the hydroelastic model for the internal fluid.

V. Conclusions

Flutter computations were performed for a wing-store test configuration in subsonic, transonic, and early supersonic air streams, using both computational-fluid-dynamics-based and linear-theory-based computational frameworks. Various JP-8 fill levels of the store were considered. In each case, the fuel was represented using both a hydroelastic added-mass model, which accounts for the sloshing effect, and a frozen-mass model, which ignores it. Both independent sets of obtained computational results suggest that, in general, the critical pressure and flutter speed decrease with an increasing fill level. They also reveal that ignoring the hydroelastic vibration effect by representing the internal fluid as a simplified frozen mass leads to natural frequencies of the dry wing-store system that are more sensitive than otherwise to the store fill level. More importantly, they also suggest that, in general, but not necessarily at all speeds, ignoring the hydroelastic effect tends to overestimate the added-mass effect and underestimate the critical pressure and flutter speed. The obtained sets of computational results also reveal that, at least in the supersonic regime, excluding the aerodynamics of the store from the computational aeroelastic model, for example to accelerate parameteric studies, leads to erroneous numerical predictions.

  cell volumes of the fluid grid b = body forces acting on the structure b s = semichord of the wing at its root C = Jacobian matrix of the numerical flux function with respect to the fluid grid velocity c = arc center D = structural damping matrix d = distance det = determinant of a matrix div = divergence operator E = Young's modulus E = matrix resulting from the time derivative of the product of the fluid cells volumes and the discrete fluid state vector Ẽ = fictitious elasticity properties F = external force field acting on the structure F = numerical convective flux function vector F = fluid convective flux f ae = vector of external aerodynamic forces acting on the f int f m FSI G = vector of internal forces acting on the structure = generalized force vector = flutter speed index = Jacobian matrix of the numerical flux function with respect to the fluid grid position H = Jacobian matrix of the numerical flux function with respect to the discrete fluid state vector I = identity matrix i = pure imaginary number J = Jacobian of a frame transformation K = structural stiffness matrix K 0 = modified structural stiffness matrix K = pseudostiffness matrix Kc = matrix related to the pseudostiffness matrix K c = matrix accounting for the effect of the structural motion on the fluid grid at the fluid-structure interface KF = transmission matrix between the structural displacement vector and the fluid grid position vector k = dimensionless reduced frequency L = length L o = reference length M A = added-mass matrix M F = mass matrix of the internal fluid M S = structural mass matrix M ∞ = freestream Mach number m empty store = mass of an empty store m S = structural mass of the wing n F = normal to the fluid-structure interface n m = number of natural eigenmodes n S = normal to the surface of the structure P = Jacobian matrix of the vector of external aerodynamic forces with respect to the discrete fluid state vector P m = generalized Jacobian matrix of the vector of external aerodynamic forces with respect to the discrete fluid state vector p = pressure fluctuation of the internal fluid Q m = generalized aerodynamic force matrix R of the internal fluid u m = generalized modal displacement vector u S = structural displacement field V = matrix of structural eigenmodes V = volume of a conical frustrum V ∞ = freestream velocity v F = velocity field of the internal fluid w = fluid state vector w = discrete fluid state vector x = position vector field x = fluid grid position vector x F = displacement field of a fluid grid point Γ = fluid free surface δ = pertubation of a quantity ∂Ω S = surface of the structure ε = fictitious strain tensor ϵ S = strain tensor λ = square of the circular frequency μ = ratio of the structural mass of the empty store and the product of the freestream density and the frustum volume ν = Poisson's ratio ξ = reference position of a fluid grid point π = pressure constant ρ = fictitious density ρ aluminium = aluminium density ρ F = density of the internal fluid ρ S = structural density ρ ∞ = freestream density Σ = fluid-structure interface boundary σ S = structural stress tensor τ = dimensionless time Φ = vector of fluid potential degrees of freedom ϕ = potential function for the internal fluid ψ = angle ω = (circular) frequency Ω = diagonal matrix containing the circular frequencies of the structure Ω F = internal fluid domain Ω S = structural domain ω α = first dry torsional mode frequency of the structure ∇ quantity LE = store leading-edge quantity o = equilibrium configuration t = time derivative = quantity satisfying the free surface Dirichlet boundary conditions

Fig. 1

 1 Fig. 1 Cut-away view of a typical store showing the definition of the various domains and interfaces.

Fig. 2

 2 Fig. 2 Wing-store system: hypothetical assembly of the 2.5 ft weakened model 3 of the AGARD Wing 445.6 and a store.

is plotted in Fig. 4 .

 4 At all Mach numbers where the critical speed decreases monotonically with the increasing fill level, two observations are noteworthy. First, f is almost always negative, meaning that the frozen-mass representation of the internal fluid underpredicts the flutter quantities or, in other words, overestimates the added-mass effects of the internal fluid, in comparison to the hydroelastic representation of the internal fluid. Second, in most cases, the absolute value of f increases with the store fill level. Therefore, the overprediction of the added-mass effects is exacerbated by an increasing fill level. At Mach 1.201 and fill level 15%, f becomes positive. This means that, at these conditions, the cruder frozen-mass model overpredicts the size of the flight envelope.

Fig. 3 Fig. 4

 34 Fig. 3 Variations of the critical freestream velocity with the store fill level predicted by AERO.

Fig. 5

 5 Fig. 5 Variations of the flutter speed index with the freestream Mach number predicted by AERO.

Fig. 6

 6 Fig.6Comparison of the variations of the critical freestream velocity with the store fill level predicted by AERO, REELC, and NASTRAN (ZONA51) for supersonic airstreams.

  2 1154.0 (-0.2) 1141.4 (-1.3) 1120.5 (-3.1) 1105.3 (-4.4) 1089.7 (-5.7) 0.901 Frozen mass 1156.2 1151.0 (-0.4) 1135.9 (-1.8) 1114.0 (-3.6) 1091.3 (-5.6) 1071.4 (-7.3) 0.960 Hydroelastic 1275.6 1272.0 (-0.3) 1259.2 (-1.3) 1234.0 (-3.3) 1221.1 (-4.3) 1204.6 (-5.6) 0.960 Frozen mass 1275.6 1268.9 (-0.5) 1252.5 (-1.8) 1228.6 (-3.7) 1206.8 (-5.4) 1186.8 (-7.0) 1.072 Hydroelastic 1376.5 1323.4 (-3.9) 1389.3 (0.9) 1294.5 (-6.0) 1313.4 (-4.6) 1209.9 (-12.1) 1.072 Frozen mass 1376.5 954.3 (-30.7) 717.1 (-47.9) 1294.1 (-6.0) 1249.8 (-9.2) 1232.5 (-10.5) 1.141 Hydroelastic 1447.3 1189.5 (-17.8) 1376.3 (-4.9) 1337.2 (-7.6) 1303.7 (-9.9) 1269.4 (-12.3) 1.141 Frozen mass 1447.3 862.3 (-40.4) 642.6 (-55.6) 1318.1 (-8.9) 1226.9 (-15.2) 1206.2 (-16.7) 1.201 Hydroelastic 1392.4 1264.4 (-9.2) 1829.0 (31.3) 1752.6 (25.9) 1689.0 (21.3) 1631.2 (17.1) 1.201 Frozen mass 1392.4 984.4 (-29.3) 681.9 (-51.0) 1685.3 (21.0) 1597.1 (14.7) 1537.8 (10.4) 1.272 Hydroelastic 1664.4 1544.5 (-7.2) 2399.5 (44.1) 2210.0 (32.8) 2070.1 (24.4) 1988.1 (19.4) 1.272 frozEn mass 1664.4 1274.7 (-23.4) 818.7 (-50.8) 2100.2 (26.2) 1924.3 (15.6) 1810.8 (8.8) a Between parentheses: percent change with respect to 0% store fill level.

Fig. 7

 7 Fig.7Comparison of the variations of the critical freestream velocity with the store fill level predicted by AERO and NASTRAN (DLM) for subsonic airstreams.

Fig. 8

 8 Fig. 8 Differences between the critical freestream velocities predicted using a hydroelastic model of the internal fluid or a frozen-mass simplification and AERO, REELC, and NASTRAN.

  matrices H o , G o , and C o result from a first-order Taylor expansion of the vector of numerical flux functions around the equilibrium configuration ( w o , x F o ). The matrix H o has, in general, a rank equal to the number of fluid degrees of freedom (DOF). The matrix E o results from a similar expansion of the time derivative of the product of the cell volumes and the fluid state vector.

Next, the structural subsystem [Eq. (14b)] is linearized around the equilibrium state (w o , x Fo , u o ) to obtain

Table 1

 1 Variations with the store fill level of the natural frequencies computed by AERO

	Frequency, Hz a

a

Table 2

 2 Relative differences between the frequencies computed by AERO using the hydroelastic and frozen-mass models

		Discrepancies in the frequencies computed using
		two different internal fluid models (as a percentage
		of the frequency computed using the hydroelastic
				added-mass model)
	Mode	0% fill 30.9% fill 50.0% fill 69.1% fill 99.9% fill
	First flexion	-	-0.23	-0.44	-0.66	-0.72
	Second flexion	-	-0.30	-3.91	-4.09	-4.34
	First torsion	-	-3.61	-1.58	-2.47	-2.40
	Second torsion	-	0.05	0.09	0.11	0.15
	Third flexion	-	-1.24	2.15	-5.46	-0.05
	Third torsion	-	-18.57	2.09	0.06	0.30
	Store axirotation	-	-13.31	-32.98	-17.13	0.01
	Coplanar	-	-1.82	-3.74	-5.33	-3.37

Table 3

 3 Variations with the store fill level of the natural frequencies computed by NASTRAN

	Frequency, Hz a

(-14.02)

119.40 (-2.13) 119.10 (-2.38) 119.00 (-2.46) Store axirotation 162.00 120.20 (-25.80) 88.40 (-45.43) 91.53 (-43.50) 93.40 (-42.35) Coplanar 176.00 171.20 (-2.73) 165.40 (-6.02) 159.30 (-9.49) 154.50 (-12.22) a

Table 4

 4 Relative differences between the frequencies (ordered according to corresponding mode shapes) computed by NASTRAN using the hydroelastic and frozen-mass models

		Discrepancies in the frequencies computed using different
		internal fluid models (as a percentage of the frequency
		computed using the hydroelastic added-mass model)
	Mode	0% fill 30.9% fill 50.0% fill 69.1% fill 99.9% fill
	First flexion	-	0.27	-0.92	0.38	-0.94
	Second flexion	-	0.81	-5.85	2.57	-0.36
	First torsion	-	2.17	0.31	3.08	-0.78
	Second torsion	-	1.57	0.00	1.23	0.21
	Third flexion	-	1.25	2.18	-2.49	0.23
	Third torsion	-	-10.49	1.19	1.79	-0.67
	Store axirotation	-	-18.01	-32.00	-12.16	0.43
	Coplanar	-	-4.84	-4.94	-6.29	-0.71

Table 6

 6 Variations of the critical freestream velocity with the store fill level predicted by AERO Hydroelastic 3141.4 3026.4 (-3.66) 2847.3 (-9.36) 2589.5 (-17.57) 2408.9 (-23.32) 2285.5 (-27.24) 1.272 Frozen mass 3141.4 3029.1 (-3.57) 2735.1 (-12.93) 2436.7 (-22.45) 2207.7 (-29.72) 2069.5 (-34.12)

	V cr ∞ , ft∕s a

a Between parentheses: percent change with respect to 0% store fill level.

Table 5

 5 Freestream parameters chosen for the flutter analysis[START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF] cr ∞ of the freestream velocity, that is its values at the onset of flutter.

	0.499 0.4278 ρ ∞ , slugs∕ft 3 0.000830 0.000404 0.000193 0.000123 0.000107 0.000152 0.000125 0.000115 0.678 0.901 0.960 1.072 1.141 1.201 1.272 M ∞ 0.2082 0.0995 0.0634 0.0551 0.0783 0.0644 0.0593 ρ ∞ , kg∕m 3

Table 7

 7 Variations of the critical freestream velocity with the store fill level predicted by REELC in supersonic domain

	V cr ∞ , ft∕s a

a

Table 8

 8 Variations of the critical freestream velocity with the store fill level predicted by NASTRAN

	V cr ∞ , ft∕s a

Acknowledgments

The Stanford University authors acknowledge partial support by a research grant from the King Abdulaziz City for Science and Technology, partial support by The Boeing Company under Contract Sponsor Ref 45047, and partial support by the U.S. Air Force Office of Scientific Research under grant FA9550-10-1-0539. Any opinions, findings, conclusions, or recommendations expressed in this document are those of the authors and do not necessarily reflect the views of the sponsors.

Appendix A: Store Geometry

The test geometry considered in Sec. IV.A consists of the 2.5 ft weakened model 3 of the AGARD Wing 445.6 [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF] and an axisymmetric store with smooth front and back ends. Using a coordinate system (in feet) centered at the leading edge of the root of the wing (the same system as in Yates [START_REF] Yates | AGARD Standard Aeroelastic Configuration for Dynamic Response, Candidate Configuration 1[END_REF]), the store is positioned so that its front tip is in the plane x -0.1667 and its centerline is at the intersection of the planes y 0.75 and z -0.375. The wing and the store are connected by two sets of linear multipoint constraints. The first set connects the point on the wing with coordinates (0.7966, 0.7500, 0.0000) to the two points on the store with coordinates 0.7395; 0.6783; -0.2018 and 0.7395; 0.8218; -0.2018, respectively. The second set connects the point (1.947, 0.7500, 0.0000) on the wing to the two points 2.000; 0.6783; -0.2018 and 2.000; 0.8218; -0.2018 on the store (see Fig. 2).

The store is a cylinder capped at both ends by identical tapered sections, each of which begins with a spherical tip then varies smoothly toward the cylindrical section. The following figure shows a quarter of the store cross section through the centerline (the remainder of the store profile can be obtained by symmetry) and highlights the parameters that completely define the geometry.

Each cross section of a tapered section has two arcs: the leadingedge arc s LE , and the arc s 2 , which connects the leading-edge arc to the cylindrical section. These two arcs are centered at c LE and c 2 , respectively. The arc s 2 is tangent to s LE at one end and tangent to the line representing the cross section of the cylindrical wall at the other end.

The values of the parameters shown in Fig. A1 are as follows: R cyl 0.1875 ft, L cyl 1.8333 ft, L cap 0.3333 ft, R LE 0.04249 ft, and ψ 47 deg. From these values and the tangency conditions stated previously, it follows that R 2 0.4099 ft and d 2 0.2224 ft.

The local coordinates are aligned so that the x 1 axis runs along the centerline of the store, from the root to the tip of the tapered section. The variation of the radius of the tapered section with the coordinate x 1 is given by

The wall thickness of the store t wall is relatively small. Hence, the mass of the empty store can be computed as follows:


Detection of manhole covers in high-resolution aerial images of urban areas by combining two methods

The detection of small objects from aerial images is a difficult signal processing task. To localise small objects in an image, low-complexity geometry-based approaches can be used, but their efficiency is often low. Another option is to use appearance-based approaches that give better results but require a costly learning step. In this paper, we treat the specific case of manhole covers. Currently many manholes are not listed or are badly positioned on maps. We implement two conventional previously published methods to detect manhole covers in images. The first one searches for circular patterns in the image while the second uses machine learning to build a model of manhole covers. The results show non optimal performances for each method. The two approaches are combined to overcome this limit, thus increasing the overall performance by about forty percent.

Introduction

Urban growth is an ongoing trend and one of its direct consequences is the development of underground utility networks. Over the past century it was common practice for public service providers to install, operate and repair their networks separately [START_REF] Rogers | Condition assessment of the surface and buried infrastructure-A proposal for integration[END_REF], so it is now very difficult to find accurate records of utility network maps in cities in both industrialized and developing countries. This crucial problem will worsen as cities expand and their networks increase in size and complexity [START_REF] Jamil | Underground utility mapping and its challenges in Malaysia[END_REF] [START_REF] Metje | Mapping the underworld: State of the art review[END_REF]. Urban works will thus be more prone to delays with concomitant additional costs [START_REF] Yung | Evaluation of subsurface utility engineering for highway projects: benefit-cost analysis[END_REF] [START_REF] Jeong | Evaluation of an emerging market in suburface utility engineering[END_REF]. However, some of these networks have surface access traps which may be visible on airborne or satellite images. If correctly detected, these elements could serve as indicators of underground utility networks. Furthermore, they can be used as landmarks in photogrammetric applications [START_REF] Drewniok | Exterior orientation-an automatic approach based on fitting analytic landmark models[END_REF] or in geotechnical works such as subsidence calculations [START_REF] Liu | Exploration of subsidence estimation by persistent scatterer InSAR on time series of high resolution TerraSAR-X images[END_REF]. We 1 put forward a methodology to detect small urban objects, namely manhole covers and grates, on very high-resolution aerial and satellite images. Two methods are tested. The first is based on a geometrical circular filter whereas the second uses machine learning to retrieve some patterns. The results are compared and combined in order to benefit from the two approaches.

Methods

Circular detection filter

r 2 r r 1 R 1 R 2 S 1 3 S 1 1 S 1 5 S 2 5 S 1 6 S 2 6 S 1 7 S 2 7 S 1 8 S 2 8 S 2 1 S 2 2 S 2 3 S 2 4 S 1 4 S 1 2
Figure 1: Circular filter (from [START_REF] Niigaki | Circular Object Detection Based on Separability and Uniformity of Feature Distributions using Bhattacharyya Coefficient[END_REF]). Left: R 1 and R 2 are the two main regions; middle and right: definition of the subregions (the phase shift is π/4).

The geometrical approach is based on the method proposed in [START_REF] Niigaki | Circular Object Detection Based on Separability and Uniformity of Feature Distributions using Bhattacharyya Coefficient[END_REF] for the detection of circular patterns in a noisy and low contrasted image. The authors propose a filter that consists of two annular regions R 1 and R 2 of radius r 1 to r and r to r 2 , each of which is divided into eight sub-regions (see Figure 1). The filter is applied to a grey-scaled image, obtained with the luminance formula. Three indices are computed to detect a circular pattern on a sliding window, using normalized histograms of each region/subregion.

The first one estimates the similarity between two statistical distributions using the Bhattacharyya coefficient:

S(R 1 , R 2 ) = N x=0 p 1 (x)p 2 (x) (1)
with p 1 (respectively p 2 ) being the normalized histogram of R 1 (respectively R 2 ) and N the maximal intensity of the two histograms. The result of this index is 1.0 when the two histograms are identical and 0.0 when they are completely different. The optimal value for this index in case of a circular pattern is thus the lowest one.

The second index is computed to avoid detection of linear patterns. It is based on a comparison between intensity distributions of R 1 and the eight subregions of R 2 :

S 8 = max j∈1...8 S(R 1 , S j 2 ) ( 2
)
where S j 2 denotes the j th subregion of R 2 . This index is low when all oriented similarity scores between R 1 and R 2 subregions are small. The last index assesses the uniformity inside the two main regions:

U(R i ) = min j,j ∈{1..8} S(S j i , S j i ) (3)
where i stands for the region and j, j for the subregions.

The three indexes are merged in a global index for circular pattern detection:

ζ = (1 -max {S(R 1 , R 2 ), S 8 }) • U(R 1 ) • U(R 2 ) (4)
The greater the value of ζ, the higher the likelihood of a circular pattern.

Machine learning approach

The appearance approach uses a machine learning algorithm which builds a model from data. This method is quite efficient for urban object detection [START_REF] Chaumont | Automatic localization of tombs in aerial imagery: application to the digital archiving of cemetery heritage[END_REF]. The learning step involves in three stages [START_REF] Viola | Rapid Object Detection using a Boosted Cascade of Simple Features[END_REF]. 1) We extract small images of manhole covers and small random images from the training database. Each small image is resized to a constant size in order to be robust to all scale. 2) To transform the data to a feature vector, we extract multiple histograms from oriented gradients (HOG) [START_REF] Qiang | Fast Human Detection Using a Cascade of Histograms of Oriented Gradients[END_REF]. Each HOG vector measures the distribution of the gradient angles within the image.

3) The final step consists in using a linear SVM classifier [START_REF] Fan | LIBLINEAR: A Library for Large Linear Classification[END_REF] to create the required model. During the evaluation step, a sliding window evaluates all the positions and scales in the entire image.

Merging the approaches

There is a marked difference between the two approaches. The first one uses only pixel intensity while the second one uses a robust model that is built from the gradient. We experimentally observed that the false positives varied from one approach to the other. Therefore, we combined the results from the two methods to increase the detection performance. The final score is obtained by merging the scores given by each approach.

Let f (x) be a function that returns the value of ζ for a pixel x (Eq. 4) rescaled between 0 and 100, and g(x) is a function that returns the probability of a pixel x being a manhole cover center by the machine learning approach. The product function between f and g is called h (Eq. 5) and is used as a score value for the presence of manhole covers.

Note that all the objects may not be detected by both methods. For instance, function f returns a null score for rectangular shapes. We thus suggest other functions for scoring the presence of a manhole cover, s(x) for the sum (Eq. 6), n(x) and m(x) for the min and max (Eqs. 8 and 7). Scoring functions are evaluated in the experimental section 3.

h(x) = f (x).g(x) (5) s(x) =    2.f (x) g(x) = 0 2.f (x) f (x) = 0 f (x) + g(x) otherwise (6) m(x) = max(f (x), g(x)) (7)
Figure 2: Sample image with at 4 cm resolution on which we performed learning step.

n(x) = min(f (x), g(x)) (8)
The results are compared in terms of precision and recall: precision = TP TP+FP [START_REF] Yung | Evaluation of subsurface utility engineering for highway projects: benefit-cost analysis[END_REF] recall = TP TP+FN [START_REF] Jeong | Evaluation of an emerging market in suburface utility engineering[END_REF] with TP: true positive (number of correctly detected manhole covers); FN: false negative (number of omitted manhole covers); FP: false positive (number of objects confused as manhole covers).

Experimental results

Data sets

The methods were tested on a orthorectified aerial photograph (5,300 x 5,500 pixels at 4 cm resolution) of the town of Gigean in the South of France (see Figure 2). The scene contained 125 manholes, with an average size of 80 x 80cm, 91 of which were used to build the training set and 34 were used for the testing step. Before any localisation treatment, the image was firstly segmented and only pixels from the road and street networks were used. applied to ζ varies. All manhole covers could not be detected by this method because about 20% of them had a rectangular shape and were undetectable. We also noticed a significant number of false positives. In fact, there were many circular patterns on the road. Which were oil spots or potholes and they were enhanced as the image is smoothed.

Evaluation with and without machine learning

The ROC curve of the machine learning method, which has higher precision than the previous method, is plotted in Figure 3.b. Nevertheless, poor results were obtained; for a recall of 50%, only 20% of the objects were manhole covers. The poor results could stem from the small size of the training database.

Combined method

Figure 4 illustrates the combined method, with functions H, M, N , which was more efficient than the separate approaches presented in section 2. All the fusing methods except the one using the s(x) function gave better results than each method taken separately, as can be seen in Figure 5. The poorer results obtained with the summing function s(x) were due to the addition of all detections, including false ones. Hence, the precision was worsened according to Eq. 9.

With the remaining functions (m(x), n(x) and h(x)), the precision was increased by more than 25%. For instance, a recall of 40% corresponded to a precision of approximately 43% for the machine learning approach (M L curve) and reached 68% for the min function (N curve).

We noticed that the maximum function m(x) gave a better score: nearly twice the precision of the machine learning method. Actually, if the classifier probability value is high, then this response is more efficient than the dot or minimum return. This could have two explanations. First, nearly ideal circular shapes may be detected with high precision by the circular filter, but may have a lower probability with the machine learning method, which has no information on the shape. In contrast, the circular filter gives a lower score for rectangular manhole covers, contrary to the machine learning approach, which is more robust.

The main shortcomings of this kind of approach is that, when using two methods, the combined recall value is always lower than the lowest score of each method taken separately.

Conclusion

The objective of this work was to put forward a methodology to detect small urban objects on high resolution images in order to reconstruct buried utility networks. The first case study on a sewage network is highly encouraging, especially considering that we had a small training set, which was used for the machine learning method. The preliminary results show that nearly 40% of manhole covers were detected with a precision of 80%. This compares well with the common practice of field surveyors in southern France, who are frequently asked to locate only one third of manhole covers in order to cut down operational costs. Many options remain unexplored and will be investigated in the very near future. For instance, a square filter could be added to the geometrical detection procedure, additional SVMs could be tested as well as other combined methods. The high resolution of the aerial images for the localisation of small urban objects is a highly promising research field.

Figure 3 .

 3 Figure3.a illustrates the performance of the circular detection filter method using a ROC curve (precision as a function of the recall) where the threshold

Figure 3 :

 3 Figure 3: Precision vs recall. a) with circular pattern detection; b) with the machine learning method.

Figure 4 :

 4 Figure 4: Results of the combined method: summing (s), dot (h), maximum (m) and minimum (n) functions. Comparison with both approaches (ML-machine learning) and (CF-circular filter)

Figure 5 :

 5 Figure 5: The two figures represent the results obtained with the circular filter on the left and the machine learning method on the right. The red rectangles are the true negatives, the blue rectangles are the true positives and the green rectangles are the ground truth. The lines connect the two common results between two approaches. The testing image is at 4 cm resolution.