Investigation of Sloshing Effects on Flexible Aircraft Stability and Response

The present paper provides an investigation of the effects of linear slosh dynamics on aeroelastic stability and response of flying wing configuration. The proposal of this work is to use reduced order model based on the theory of the equivalent mechanical models for the description of the sloshing dynamics. This model is then introduced into an integrated modeling that accounts for both rigid and elastic behavior of flexible aircraft. The formulation also provides for fully unsteady aerodynamics modeled in the frequency domain via doublet lattice method and recast in time-domain state-space form by means of a rational function approximation. The case study consists of the so-called body freedom flutter research model equipped with a single tank, partially filled with water, located underneath the center of mass of the aircraft. The results spotlight that neglecting such sloshing effects considering the liquid as a frozen mass may overshadow possible instabilities, especially for mainly rigid-body dynamics.


Introduction
Large passenger aircrafts are likely to have wings that are highly flexible structures, carrying an amount of fuel comparable in weight to that of their structural components, which can deform significantly when atmospheric turbulence or gust is encountered. Then, the dynamic loads that act on the structure due to the oscillatory motion of the heavy liquid mass may represent a critical issue seriously influencing aircraft performances or even its safety.
Sloshing means any motion of the free liquid surface caused by any disturbance to partially filled liquid containers. The basic problem of sloshing involves the estimation of pressure distribution, forces, and moments applied by the liquid to the tank.
Although sloshing can assume nonlinear behaviors experiencing possible instabilities due to tank high acceleration values [1], the main problem coped in this work is the estimation of the effects of linear sloshing onto aircraft aeroelastic stability and response to little perturbations. There are several analytic methods that account for sloshing for simplified tank geometries (as parallelepiped or cylinder): according to Refs. [2,3], a realistic representation of the liquid dynamics inside closed containers can be approximated by an equivalent mechanical model (EMM), whose parameters can be suitably related with the physical quantities obtained from the linearized potential flow theory [4].
Among the mentioned literature [3,4] which was stimulated by specific needs of aerospace community, several contributions on sloshing in aeronautic applications can be found. In particular, the effects of sloshing on aircraft aeroelastic flutter stability was considered in Refs. [5,6] where liquid dynamics were modeled by means of EMMs in Ref. [5] and frozen fluid approach in Ref. [6]. A detailed formulation of fluid-structure interaction is also present in Ref. [7] that aim at integrating flight dynamics, aeroelasticity and sloshing dynamics of flexible launch vehicles. Specifically, a reduced order model of sloshing forces were obtained by means of boundary element method. The fuel sloshing was also investigated in Refs. [8,9] to study its effect on aircraft component loads.
The final perspective of the present work is the integration of a reduced order model (ROM) based on the EMMs into the flexible aircraft flight physics, where flight dynamics, 1 3 aeroelasticity and controls are involved. We will rely on the approach mainly described in Ref. [11], where the fully coupled equation of motion have been derived by assuming practical mean axes constraints. The equations of motion were linearized around aeroelastic trim conditions and recast in time-domain state-space form by approximating doublet lattice method fully unsteady aerodynamics via rational polynomial functions, thus obtaining a model that includes rigid-body, elastic, and aerodynamic state variables, as well as the variables related to the aforementioned fluid slosh dynamics. The only maneuver considered in the analysis is the stationary level flight condition, although the present formulation can be used for any maneuvering condition.
Employing as reference aircraft the body freedom flutter (BFF) with a parallelepiped tank (partially filled with water) placed underneath its center of mass, the effects of the sloshing dynamics on the stability and response have been investigated comparing the case in which the liquid is assumed as a frozen mass and the case in which the liquid is allowed to slosh. Finally, a simple proportional control law is assumed for system stability augmentation.
The paper is organized as follows: the equivalent mechanical model along with its integration into the integrated framework of flight dynamics and aeroelasticity is introduced in Sect. 2; the aeroelastic stability and response of the present test case, namely the considered flying wing research model, are illustrated in Sect. 3 along with the implementation of the control law; a section of concluding remarks ends the paper.

Theoretical Background
In this section, the reader is first addressed to analytical models in Sect. 2.1, generally used in engineering applications to account for force and moments generated by small perturbations of simplified tank geometries. After that, the reduced order model obtained for the description of the sloshing is integrated into the aeroservoelastic model of a flexible aircraft in Sect. 2.2.

Lateral Sloshing Analytical Models
The analytical representation of linear sloshing for potential flows is here presented. Linear fluid slosh dynamics for almost rigid parallelepiped tank is featured by the excitation of lateral forces and moments when perturbing the system by means of a lateral acceleration and rotations. However, for such simplified geometries, vertical forces are not taken into considerations, since they arise only in the presence of more relevant nonlinear phenomena.
Starting from the basic assumptions generally adopted for such analytical models (mainly related to the fluid properties and tank geometry), the procedure here presented considers the formulation in Refs. [2][3][4] as a starting point with reference to a rectangular tank geometry. Indeed, the analytical models provide a Laplace domain expression of lateral force and moment about the center of mass of the liquid denoted as G arising by lateral motion and rotation of a 2D tank as a function of the tank size and filling level. Considering a rectangular tank as in Fig. 1, the force and moment can be expressed as: where ~ represents the Laplace transformed variables, s is the Laplace variable, n are the natural frequencies of the sloshing dynamics, is the density of the fluid, g is the gravity acceleration, a and h f represent, respectively, the tank-edge length in y-direction and the height of the liquid mass and b is tank-edge length in x-direction. Equation 1 can eventually be recast as where is a linear sloshing operator derived by analytical models, hereafter referred as generalized-sloshing-forces matrix (GSF). Considering a limited number of sloshing dynamics N s and according to the equivalent mechanical models (EMMs), the above equation can be rewritten as where (1) and s is a diagonal matrix, having as i-th diagonal element, a term associated with the modal damping of sloshing phenomenon (see Ref. [10]). I f is the effective sloshing fluid inertia moment (see Ref. [3]) and This allows us to express the Eq. 2 as follows by defining as ̃ the sloshing modal state space vector (in the Laplace domain) that collects the sloshing states, each of them representing a different dynamics of the fluid contained into the tank. For a wide range of applications in the framework of linear sloshing of simplified tanks, these Laplace-domain models allow to obtain sloshing forces and moments by having as time-domain counterpart a set of linear ordinary differential equations.

Integration of Sloshing ROM into a Flexible Aircraft Model
The present work aims at including sloshing dynamics into the integrated modeling of flight dynamics and aeroelasticity. Therefore, the equations of motion of the flexible aircraft require to be updated by adding further equations and further variables. The present modeling of coupled aeroelasticity and flight dynamics is based on the formulation in Refs. [11,12], in which the rigid body degrees of freedom are associated with a set of pratical mean axes (PMAs) and the linearized structural dynamics is described as a combination of unconstrained aircraft mode shapes in the PMAs coordinate system that simplifies the inertial coupling description. The PMAs are featured by having the origin coinciding with the instantaneous center of mass of the aircraft, and the orientation of the principal axes remains constant in the deformed configuration.
The equations of motion of the unrestrained flexible aircraft are provided below by including the conservation of momentum and angular momentum along with the structural dynamics equations: where and are the centre of mass velocity and the angular momentum, respectively, m is the total mass, and are the total force and moment, respectively, whereas q n denotes the n-th modal coordinate (which, together with the n-th elastic mode shape of the unconstrained aircraft n , allows to write the elastic displacement of a generic point as = ∑ ∞ n=1 q n n ) and m n , k n , and f n are the n-th modal mass, stiffness and external force. The equations of motion in Eq. 7 are then recast with respect to a body frame of reference and linearized around a level flight aeroelastic trim condition: where is the inertia tensor, and is the angular velocity. It is worth noting that the inertial coupling terms in Ref. [11] are here neglected. Moreover, the dynamics of the control surfaces along with their effects on rigid-body and structural dynamics enriches the present modeling [13]. In the present formulation, the variables associated with a second-order dynamics are grouped in the following vector: where are, respectively, the perturbation vectors of the center of mass coordinates in inertial frame of reference, of the Euler angles and of the control surface rotations, with Δ a , Δ r and Δ e associated, respectively, to the aileron, the rudder and the elevator. Moreover, the perturbation vector of the modal coordinates is given by where N e is the number of elastic modes. Nevertheless, the equations of motion of the aircraft are expressed in a noninertial frame of reference. Indeed, the following vector is defined: where Δ = {Δu, Δv, Δw} and Δ = {Δp, Δq, Δr} are the translational and angular velocities in the PMAs coordinate system, respectively. The linearized relation between Δ in Eq. (12) and Δ is expressed as with which allows to highlight the link between the variables expressed in the PMAs and those defined in the inertial reference system (being * 1 a square matrix with dimension ( 9 + N e )).
As far as it concerns the small-disturbance aerodynamics, in this work it is modeled via the doublet lattice method (DLM) available in MSC Nastran for steady and unsteady linear aeroelastic analyses [14]. In particular, the MSC Nastran flutter solver provides the so-called generalized aerodynamic force (GAF) matrix . This provides a frequencydomain fully unsteady representation of small-disturbance aerodynamics due to perturbations of rigid body and elastic degrees of freedom. Specifically, the aerodynamics is provided as a function of the reduced frequency that coincides with Strouhal number and is expressed as Although is generally provided by assuming modes defined in the inertial frame of reference, it can be easily transformed as done in Refs. [11,12] to be applied when body-frame coordinates are used instead. The GAF matrix is then approximated via following rational-polynomial interpolation in the reduced frequency domain Then, analytical expansion allows to express unsteady aerodynamic forces in time domain as it follows: where q D = 1 2 ∞ U 2 ∞ is the reference dynamic pressure and is the vector of N a additional aerodynamic variables representing the aerodynamic wake dynamics.
The aeroservoelastic model in level flight can be expressed in the Laplace domain as follows: where are the mass, damping and stiffness matrices, respectively, which are such as to include in their definition the matrices ̂ 2 ,̂ 1 and ̂ 0 deriving from the finite states interpolation of the GAF matrix.
This interpolation leads to the definition of a representative vector of the wake dynamics, indicated with Δ , which can be expressed in the Laplace domain as where p ∶= sb∕U ∞ is the non-dimensional Laplace variable (b is the reference half-chord). The mass e , damping e , and stiffness e matrices of the integrated system are the same as in Refs. [11,12]. Moreover, the description of the aircraft motion in the PMA non-inertial reference requires accounting for the projection of the weight force on the aircraft body reference. Under the assumption of small perturbation with respect to the trimmed configuration, such a contribution was modeled as an additional stiffness term indicated as . Finally, ̃ in Eq. (16) is the vector collecting the generalized modal forces applied at the control surfaces. To integrate sloshing into the presented model, it is worth noting that the sloshing forces in Eq. (6) are expressed with respect to the center of mass of the undisturbed fluid that depends on the filling level of the tank. Generally, for finite element applications, it is useful to express forces and moments with respect to a fixed point, which in this case can be associated with the geometric center of the tank. This needs to introduce a transformation matrix, indicated as , that is given by Subsequently, other transformation matrix indicated as t is introduced to express the state variables of the tank (that now refer to the geometric center of the tank) with respect to the perturbation variables of the complete aircraft. t is defined as where ( ) is the rotation matrix , related to the tank, which allows to pass from a "tank" frame of reference to the FE (17) model frame of reference. The angle measures the deviation between these two systems. t instead, represents a 6 × (9 + N e ) matrix whose columns are the modal shape vectors associated with the rigid and structural modes (assuming empty tank) and referred to the node where the tank is located. Therefore, it is now possible to express the system in Eq. (16), including sloshing liquid inside the tank, as follows where ̃ is the generalized sloshing load, and its expression is given by The sloshing modal coordinate Δ̃ has been assumed equal to having considered the following transformation: The vectors Δ and Δ can be subsequently augmented by taking into account the modal coordinates associated with sloshing dynamics The equations of motion of the global system, defined in the Laplace domain and in a second-order form, are given by the following expression (21) (s + )Δ̃+ Δ̃= q D̂ Δ̃ +̃ +̃ , (26) s̆ +̆ ̃s +̆ Δ̃s = q D̂ Δ̃ +̃ , where the mass, stiffness and damping matrix associated with the final updating system are defined as follows:  It is worth noticing that the introduction of sloshing terms recalls the dynamic condensation (i.e., Craig-Bambton reduction, Ref. [15]), where the finite states, Δ̃ , define the sloshing substructure dynamics. In this framework, the master DoFs have been projected onto the vibration modes of structure with the empty tank. Moreover, the matrix ̄ arises because of volume forces that depend on the mass displacement. In the definition of the matrix ̆ , the term −̄ * 1 comes from the kinematic relationships expressed in Eq. (13) that are used to complete the integrated model. The state space vector ∈ ℜ (2(9+N e +N s )+N a ) can be expressed as follows Finally the global integrated and linearized system can be recast as = m a , m r , m e is the input vector that lists the moments applied at the different control surface hinges (see Ref. [13]). The state-space matrix and the input matrix H are where * ,̆ and ̆ have been created to fit the dimension of the state-space vector (including the sloshing dynamics), considering respectively as first block, matrices * 1 ,̂ and ̂ , and all the remaining elements equal to zero and ̂ is a selection matrix having dimension [(9 + N e + N s ) × 3].

Sloshing Integrated Aeroelastic Analysis
In the present work, the finite element model presented in Ref. [13] has been employed to obtain an aeroservoelastic model to be integrated with the sloshing reduced order model. The model mentioned above is the body-freedomflutter (BFF) research model which is an unmanned flyingwing research aircraft, developed for studying the effects of mutual coupling between the rigid-body motion and the structural vibration of the aircraft. The first six elastic modes, with their respective frequencies, are shown in Fig. 3. To study the effects of the sloshing dynamics on the stability and the response of the complete aircraft, it is possible to consider different configurations that foresee the presence of one or more tanks. The configuration considered in the present work is equipped with a rectangular tank (containing water) underneath the center of mass of the aircraft. This configuration is represented in Fig. 2.
The geometrical parameters of the tank are reported in Table 1, whereas the two natural frequencies of sloshing are Fig. 4 Computational framework for integrated stability and response  Table 2. In the present analysis, only one sloshing mode has been chosen for both the dynamics along x and y. However, the performed sensitivity analysis did not show any particular difference in terms of stability and response. From the point of view of the FEM solver, the presence of the tank was managed by creating a node located in the geometrical center of the tank, linked to the body of the aircraft by means of the rigid body element (RBE). A concentrated mass (assigned in Nastran by CONM2 element) is considered to represent the structural mass of the tank. A numerical code has been developed to perform the analysis concerning the stability and the response of the integrated system. Its flowchart is shown in Fig. 4.
The inputs to the developed code are the aircraft FEM and DLM models along with a set of parameters related to the tank geometry and the stowed fluid. Based on these inputs, the code collects results from MSC Nastran linear structural and aeroelastic analyses and from the equivalent mechanical model of the fluid to build a database of modal steady/unsteady aerodynamics and sloshing characteristics, that, in turn, are used to compute the state and input matrices in Eq. 30. Trim analysis plays a key role in all this process. Its determination involves calculating the forces and the trim variables in level flight. Also slosh dynamics plays an important role in the computation of trim: it contributes through the presence of the matrices ̄ and ̄ defined in Eq. 22. Finally, once the aircraft database is built, the code is able to perform both stability and response analyses of the integrated system of flight dynamics and aeroservoelasticity coupled with slosh dynamics.

Stability Analysis
In this section, the results obtained from the stability analysis of the integrated model are presented highlighting the differences between the case in which the fluid contained in the tank is frozen and the case in which the fluid is free to slosh side to side of the tank. Frozen fluid mass means that no oscillatory behavior is present. So in this specific case, the fluid contribution corresponds simply to a ballast. The traditional approaches employed to describe fuel in the aeronautical tanks consider the fuel as a frozen mass. It is therefore evident that it is important to highlight the variations introduced in the study of stability and response by the introduction of such sloshing effects. The stability analysis consists in determining the eigenvalues of the state matrix expressed in Eq. 30, varying the free stream velocity U ∞ (from 15 up to 30 m/s). In this framework, the Mach number is kept fixed at 0.2. The results obtained in the aforementioned cases are shown in Figs. 5 and 6, showing the stability scenario in frozen fluid configuration (Fig. 5) and sloshing fluid case (Fig. 6).
As it can be noted from Fig. 5b spiral, phugoid and dutch roll modes are stable, whereas Fig. 5a shows flutter occurring with the pole originated from short period coupling with the bending mode dynamics. By defrosting the fluid mass, some differences arise. In Fig. 6b it is possible to identify two poles associated with the sloshing modes. Slosh dynamics clearly influences the lateral directional stability of the aircraft, since the dutch roll mode becomes unstable. Furthermore, the system gains stability by increasing the free stream velocity, The considered flutter mode involves mainly the sloshing mode and this triggers phugoid to be less stable than the previous case. So the integration of the sloshing in the aircraft model, by means of a tank placed underneath the aircraft center of mass, has a relevant influence on the lateral dynamics of the aircraft.

Response Analysis
The previous section provided that the most critical problem for the present case study consists in the lateral directional stability. Thus, the present section presents the response of the aircraft to a command to the ailerons by means of a specified hinge moment. This kind of command has been where the command period is T ctr = 1 m/s. The maximum amplitude of the impulse is assumed to be equal to 0.5 Kg m. The reference velocity is fixed at 15 m/s. For this value of the aircraft speed, the system is stable only in the case of frozen fluid. The results are shown in the Fig. 7.
Again, It is evident the influence of the sloshing dynamics on the response of the considered aircraft. After the application of the aileron hinge moment, the aircraft starts moving laterally (see Fig. 7a, d), increasing its velocity in direction y (see Fig. 7b, e). This causes the fluid to slosh inside the tank. Specifically, the fluid is excited by the translational motion in y-direction and the rotation about the x-axis of the tank, so causing a positive feedback generating of a moment which amplifies the lateral motion. From Fig. 7e-g, the coupling between the lateral-directional dynamics and the sloshing mode (in y-direction) is evident. Indeed, velocity in direction y, v G , roll angle , sloshing mode r 1 and, slightly less, the yaw angle present an evident oscillatory behavior, with a divergent amplitude. It should be noted that this representation does not take into account nonlinear effects that arise with the free surface breakage on the tank walls, which may increase the damping as well as limiting the response of the aircraft to a limit cycle oscillation. An extension of the present approach may consider the modeling of nonlinear damping associated with higher free surface oscillations.

Implementation of a Proportional Feedback Control Law Using Measured Outputs with Sensor
The goal of this section is to provide a feedback proportional control to the aileron to control the instabilities arisen in the previous section. The control strategy takes advantage of the measurements obtained by a virtual inertial measuring unit (IMU) placed on an FE node close to the tank. Depending on the node of positioning of the IMU, different outputs in terms of accelerations and rotations can be obtained, and in turn, different control strategy must be considered. This arises because different nodes of the FE model are characterized by having different modal participations. The N out outputs measured by the sensors can be expressed as follows: where a , v and d are the selection matrices able to express the measurements by means of vector Δ s and its derivatives, and ̂ is the measurement noise (which will be considered = 0.5 a sin 2 null in the present analysis, assuming ideal measurements). The proportional feedback control action implemented is defined as =̂ y , where y is a matrix containing the proportionality coefficients having dimension 3 × N out . So, by including this control action into the state space form of Eq. (26), a new closed loop system having the following updated mass, damping and stiffness matrices is obtained: In this simplified control modeling, y matrix is basically a zero matrix with non-null elements in [1,2] and [3,5], and respectively equal to 0.05 and 1. To perform the stability and response analysis, the mass, damping and stiffness matrices in Eq. (30) have to be replaced by the ones defined in Eq. (32). An updated state matrix is thus obtained, which provides the stability scenario in Fig. 8 represented by varying the free stream velocity from 15 up to 30 m/s. It can be easily noted that in this case, the control action guarantees that, for the considered speed range, all the rigid modes as well as the sloshing modes are stable. Furthermore, after the application of such control, flutter occurs due to the crossing of the imaginary axis by the poles originated by the I bending mode.

Concluding Remarks
In this work, an investigation on the sloshing-tank effects on integrated aircraft modeling for aeroelasticity and flight dynamics has been carried out. Equivalent-mechanical models (EMMs) theory has been used to perform its integration into the aeroservoelastic model of the flexible aircraft, aiming at representing the complete aircraft system that takes into consideration also the sloshing dynamics. The reference case study consisted of the body-freedom flutter (BFF) research aircraft, to which a parallelepiped tank, partially filled with water, was positioned underneath its center of mass. Neglecting the inertial coupling between rigid body and elastic body dynamics and considering the leveled flight as reference maneuver, stability and response analyses were performed to show the influence of the sloshing dynamics on the global system.
The results obtained were compared to the ones related to the case in which the tank fluid is considered as frozen, and therefore a simple ballast. The stability scenario and the control time response highlighted the coupling between the lateral-directional dynamics of the aircraft and the sloshing of the fluid inside the tank, causing a loss of stability of the dutch roll mode. Taking advantage of the use of a IMU that the FE formulation allows to place in any grid point of the structure FE model, it has been possible to implement The natural extension of this work is the possibility of exploiting data derived from the implementation of high fidelity numerical codes to describe more general sloshing phenomena, taking into account any shape of tank, and possibly also tank deformations. Indeed, they may provide also a consistent database to build a CFD-based reduced order model for sloshing, whose mathematical structure may be identical to that used for EMMs. The final challenge will be to accurately model the phenomenon of sloshing in all its aspects. For this aim, it will be necessary to perform fully coupled simulation of different physics, as that concerning violent fluid sloshing. By exploiting different ROM techniques, as for example the response surface methodology (RSM), a simplified descriptive model of the nonlinear vertical sloshing can be obtained, allowing its subsequent inclusion in the current integrated formulation.