A multi-physics computational model of fuel sloshing effects on aeroelastic behaviour

A multi-physics computational method is presented to model the effect of internally and externally-carried fuel on aeroelastic behaviour of a pitch–plunge aerofoil model through the transonic regime. The model comprises three strongly coupled solvers: a compressible finite-volume Euler code for the external flow, a two-degree of freedom spring model and a smoothed particle hydrodynamics solver for the fuel. The smoothed particle hydrodynamics technique was selected as this brings the benefit that nonlinear behaviour such as wave breaking and tank wall impacts may be included. Coupling is accomplished using an iterative method with subcycling of the fuel solver to resolve the differing timestep requirements. Results from the fuel-structural system are validated experimentally, and internally and externally-carried fuel is considered using time marching analysis. Results show that the influence of the fuel, ignoring the added mass effect, is to raise the flutter boundary at transonic speeds, but that this effect is less pronounced at lower Mach numbers. The stability boundary crossing is also found to be less abrupt when the effect of fuel is included and limit cycles often appear. An external fuel tank is seen to exhibit a lower stability boundary, while the response shows a beating effect symptomatic of two similar frequency components, potentially due to interaction between vertical and horizontal fuel motion. & 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Fuel slosh has been shown to be important in the dynamics of a number of vehicle types, including spacecraft (Abramson, 1966;Slabinski, 1978;Vreeburg, 2005;Saturn Flight Evaluation Working Group, 1961), tanker vehicles (Sankar et al., 1992) and ships (Kim et al., 2007). At takeoff, fuel may comprise around 40% of the weight of a large aircraft, while military aircraft may operate with a range of fuel distributions at a wide variety of spanwise hardpoints. Consideration of the influence of fuel motion on loads is therefore critical, and alongside this the influence on aeroelastic motion should be understood. Although the dynamic effect of fuel motion is often assumed to be equivalent to damping, this has only recently been explored numerically. The objective of this work is to develop a multi-physics model to allow determination of the influence of fuel motion for a series of two-dimensional test cases, using more widely applicable modelling assumptions than applied so far, and to a range of larger amplitude motions.
The influence of slosh on rigid body dynamics and control systems has been extensively studied (Bauer, 1963), partly owing to the close correspondence of the sloshing frequencies to those of the controlled body, both with and without gravity (Peterson et al., 1989). For these problems it is important to solve the complete coupled model as dynamic-slosh Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/jfs details of the mathematics, but the process is to take a function that approximates the Dirac delta function and convolve it with the function to be represented as in Eq. (1). The angle brackets denote the integral representation, and the objective is to approximate ∇p=ρ in order to be able to integrate forward (in inviscid flow) Du=Dt ¼ À∇p=ρ. The process starts with an integral representation Here f ðxÞ is a function to be represented by integration over a volume Ω, Wðx À x 0 ; hÞ is a smoothing function that depends on the difference between the location of interest and the points over which the integral is evaluated. h is a characteristic length of W and is known as the smoothing length. As Price (2012) shows Eq. (1) is second order accurate providing W is symmetric and normalised to give R x W dx ¼ 1, which may be seen from expanding (assuming gradients of f are known) Symmetry implies that the second term in integral, being antisymmetric, must vanish to leave the second order term as the lead contribution. This integral representation is then approximated by summation over N nearby particles, shown in Eq.
(3). This summation is possible because W is usually chosen to be a compact function which is non-zero only for a small fraction of the total domain. The radius over which W is non-zero is known as the smoothing length or support radius, thus N is the number of particles within the support radius In replacing the integral with a discrete summation, additional errors are incurred that depend on how closely the discrete system satisfies the normalisation and symmetry requirements (Price, 2012) Ideally it would be true that P ΔV j Wðx À x j Þ ¼ 1 and P ΔV j ðx j À xÞWðx À x j Þ ¼ 0, but this depends on the particle distribution. In practice, this means that the order of accuracy of the method is dependent on the particle spacing within the smoothing length. Second order accuracy is only retained if the particles are evenly distributed within the smoothing length, in much the same way as central difference finite volume schemes only retain their second order accuracy on smoothly stretching meshes. If each particle is assigned a mass and a density then Eq. (3) can be written as Derivatives can be approximated in the same way However, while the derivative of f may not be known the derivative of the smoothing function is, so the following identity may be applied: Hence And by applying the divergence theorem the first volume integral becomes a surface integral over the problem domain S, Usually the support domain is entirely within the problem domain such that the surface integral is zero (as W¼0 at the edge of the support) and the representation can be expressed entirely in terms of the volume integral. This is not true in the case where the support domain overlaps the boundary (either a free surface or a fixed boundary).
The particle approximation can then be performed to give Note that in Eq. (10) the gradient of the kernel on the right hand side is evaluated at the jth particle and as most smoothing functions are symmetric (this is required for second order accuracy) Eq. (10) can equivalently be written as Eq. (11), where the subscript i denotes evaluating the gradient at the ith particle. The sign changes as for a symmetric kernel the gradient will be antisymmetric. With this information it is possible to construct the governing equations in the SPH formulation to give Using gives where on the right hand side, f ¼ p=ρ and f ¼ ρ can be identified and used to construct the approximation from Eq. (11).
Alternative identities can be employed giving different formulations to Eq. (13) (Price, 2012). At present the SPH formulation used is Eq. (14) is the quintic kernel suggested by Wendland (1995). The equation of state, Eq. (17), is commonly used in SPH to remove the need to solve a pressure Poisson equation, and an equation of this form was first suggested by Cole (1948) for water, before Monaghan (1994) then suggested setting B such that the speed of sound in the fluid is an order of magnitude greater than the fastest particle velocity. γ here is set to 7 following Cole's (1948) work. This is known as the weakly compressible method, and is able to predict free surface motion (Monaghan, 1994;Lee et al., 2008). The second term in Eq. (15) is a viscosity term proposed by Morris et al. (1997) which helps prevent non-physical particle interpenetration. The question of how to deal with boundary conditions in SPH is complex and a number of schemes exist. Monaghan (1988) suggests three ways of enforcing a rigid boundary: 1. Forces of length scale equal to the smoothing length. 2. Reflecting particles that cross the boundary without energy loss. 3. Fixing particles on the boundary to provide a repulsive force.
In this paper the third method has been used with the modification that the solid boundary particles are fixed only relative to the others in that rigid body. Force and moment on a rigid body are therefore calculated by summing the contribution from each of the solid boundary particles.

Time integration
The Newmark-beta method (Huebner et al., 1995) with γ ¼ 1=2, β ¼ 1=4 (these constants give the constant average acceleration integration and a scheme that is unconditionally stable) is used to time integrate the SPH code by applying it to Eqs. (15) and (16) ðwhere € s ¼ Du=DtÞ: The Newmark-beta method is implicit, so the particle acceleration and density rate of change from Eqs. (15) and (16) are substituted into Eqs. (18) and (19) to perform the time integration, and fixed point iteration carried out to converge to a solution for time level n þ1.

Aeroelastic model
The objective of this work is to explore the chordwise/pitching influence of fuel on the aeroelastic system. Although for a real aircraft wing spanwise fuel motion is important, for clarity and to enable comparison to established results, this complexity is not incorporated for this analysis. A two degree of freedom model still models the pitch-plunge flutter behaviour, and can incorporate the consequences of pitch and plunge motions on the fuel. Isogai's (1980) model, previously used for calculating flutter boundaries via linear stability and CFD methods, is adopted where the governing coupled differential equations are given by Eq. (20) and parameters illustrated in Fig. 1: Here m is the mass of the wing (kg/m), S the static imbalance (kg), I α the moment of inertia (kg m) about the elastic axis, K y the plunge spring stiffness (N/m/m), K α the pitching spring stiffness (N) and L and M are the vertical force and pitching moment about the elastic axis respectively. L is positive upwards and M is positive nose up. This system is usually expressed in a non-dimensional form given by (Alonso and Jameson, 1994) Here dots rather than dashes denote the derivative with respect to non-dimensional time τ ¼ ω α t. The non-dimensional parameters are as defined by Isogai (1979) and are The mass ratio μ ¼ m=πρb 2 , which expresses the ratio of section mass to the mass of air. Non-dimensional static imbalance x α ¼ S=mb, indicating the level of pitch-plunge inertial coupling.   The non-dimensional radius of gyration r α ¼ which gives the ratio of aerodynamic to stiffness forces.
Reduced frequency k c ¼ ω α b=U 1 , which implies the level of aerodynamic unsteadiness.
The Newmark-beta scheme, also used for the SPH integration, is used for time integration of the structure. It is useful to rearrange Eq. (21) for the accelerations in terms of the positions at time t þΔt and then substitute this into Eq. (19) to solve for the positions; the accelerations and velocities can then be simply found.
It has been shown (Piperno and Farhat, 2001) in aeroelastic calculations where the external flow and structure are simulated that energy conservation has a large impact on the accuracy of aeroelastic stability predictions as any energy loss acts as a damper and a gain acts as negative damping. It may also be shown that given an energy conserving time integrator (such as the Newmark scheme), the global energy of a system consisting of a rigid body represented by boundary SPH particles and a fluid is conserved if the force and moment on the rigid body are found by summation over the forces on the boundary particles. Hence, total force and moment of the fuel are found here through a summation of forces acting on the boundary particles.
The external flow is governed by an inviscid arbitrary Lagrangian-Eulerian model (Hirt et al., 1974) and solved using a cell-centred, finite volume (Jameson et al., 1981), solver of the Jameson-Schmidt-Turkel type, as shown in Fig. 2. Time integration is implicit in real time using a 2nd order backward difference, but explicit in pseudo-time via a Runge-Kutta integration, and mesh motion is accomplished using a radial basis function method (Rendall and Allen, 2008). In this work a geometry preserving cut-cell mesh is used, with 12 refinement levels from farfield to the surface.

Inclusion of sloshing and coupling of models
The equation of motion may be modified to include sloshing and again non-dimensionalised using τ ¼ ω α t, η ¼ y=b, where F is the vertical force from the fuel and G is the moment (note the deliberate analogy to lift and moment coefficients by replacing dynamic pressure with 1 2 ρ f gh2b, the hydrostatic head for semi-chord depth), to give where the speed index S i ¼ V 1 =ω α b ffiffiffi ffi μ p , g is the acceleration due to gravity, b is the semi-chord of the wing, the slosh index  sloshing reduced frequency k sl ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω 2 α b=g p is implied that may be related to its aerodynamic equivalent by k sl ¼ k c V 1 = ffiffiffiffiffi ffi bg p . gb=c 2 is a ratio of characteristic speeds, as c is the speed of sound in the external flow and gb is proportional to the characteristic wave speed of shallow water waves.

Time marching
For time accurate aeroelastic simulation it is preferable to synchronise the fluid and structural equations every timestep. It is possible to solve both the fluid flow and the structural problem together using a monolithic scheme (Bendiksen, 1991;Blom, 1998), but this is impractical for many problems, lacks flexibility and limits application of the wealth of experience that exists in the computational fluid dynamics (CFD) and FE fields. It is more usual to apply the partitioned method (Piperno et al., 1995) and solve each of the fluid and structural domains separately using well established techniques for the field in question and then couple the solutions by the exchange of forces and displacements. In this work the partitioned scheme is used owing to these advantages and also the difficulty of envisaging a monolithic scheme where part of the solution is calculated using a meshless method and part using a mesh based method.
For strong coupling, sub-iteration is often used to achieve convergence over a timestep (Alonso and Jameson, 1994). In comparison, weak coupling alternates between structural and fluid calculations without regard for synchronisation (Farhat Fig. 11. System in motion, Tα is the torsional period, 67% fill. and Lesoinne, 2000). A strongly coupled partitioned approach requires the exchange of information between solvers and this is an important and well studied problem with a number of methods for exchanging information across the different mesh geometries in each domain (de Boer et al., 2007;Farhat et al., 2003;Rendall and Allen, 2008). However, in this work the structural model is such that integration of forces around the boundary is sufficient for the transfer of information, but owing to the different timesteps the most important consideration is how force and displacement information is exchanged in time.
Three separate codes must be coupled together; the Euler code for the external aerodynamics, the structural code and the SPH for the internal fuel. The most natural way to do this is through the structural solver; forces on the aerofoil due to the aerodynamics and the internal fuel can be applied as source terms. Once Eq. (22) has been solved the pitch and plunge are communicated to the SPH and Euler codes which then adjust their boundary geometry.
However, due to the wall boundary condition (in the sense that a timestep integration must not move a particle through a wall) and the stiffness of the weakly compressible pressure/density relation the SPH code requires in general much smaller timesteps than the Euler or structural models. Typically the SPH timestep is in the region of ten or twenty times smaller than the Euler step, so to accommodate this mismatch a number of SPH timesteps are taken for each Euler step. During these substeps the fuel tank is forced in such a way that the position at one of these real substeps is calculated by a linear interpolation between the position at the Euler/structural timesteps, as shown in Fig. 3. The force coefficients are extracted from the SPH code only at the final substep.     The three solvers are therefore strongly coupled through iteration following the process outlined in Fig. 4. The SPH and Euler solvers are both advanced to the next timestep where the force coefficients on the structure are calculated. These force coefficients are then used to calculate new pitch and plunge positions at the end of the timestep, which are then communicated to the Euler and SPH solvers to iterate towards an improved solution. This process cycles until a convergence criterion is met after which the timestep advances; a solution is considered to be converged when the changes over each iteration in position, velocity and density are three orders of magnitude smaller than the initial change.

Sub-iteration of SPH tool
In order to validate the concept of running each model at a different timestep a comparison was made for a single time history, as shown in Fig. 5. The baseline was obtained by running the CFD code at every 16th SPH step, which means the CFD code is running at a smaller timestep than required for accuracy. Indeed, for this case the overhead in simulation time for running the CFD code at every 16th SPH timestep instead of every 64th SPH step was 312%. Little difference is seen as a result of the larger Euler step, and so this technique allows the aerodynamic timestep to be set independently from the SPH timestep. Final simulations used 32 SPH steps per Euler step.

CFD-structure
Previous studies have been conducted using the two degree of freedom model used in this paper. Alonso and Jameson (1994) calculate a flutter boundary for the NACA-64A010 aerofoil using the structural parameters suggested by Isogai (1979) to model a swept back wing. Using the notation set out in Section 3 these are μ ¼60, x α ¼ 1:8, r α ¼ 3:48 and ω h =ω α ¼ 1.
To find the flutter boundary the aerofoil was given an initial impulsive angular velocity and the resulting time histories were examined to find the flutter index at which the perturbation began to grow. The flutter boundary produced using the time marching method presented in this work, and also using a p-k analysis (Hassig, 1971), is compared to Alonso and Jameson's (1994) in Fig. 6. The good agreement between the codes suggest that the flow-structure coupling scheme used here is accurate, while the slightly lower boundary prediction of the p-k analysis is consistent with a component of nonlinear damping from shock motion not being completely represented in the harmonic coefficients used in the p-k method.

Numerical test
In order to provide a comparison between the fuelled and unfuelled systems that is independent of inertial effects, the total mass, moment of inertia and centre of mass of the wing and fuel at equilibrium are kept constant whatever the fuel fill level. This means for a full fuel tank the dynamics of the system will be very similar to the system with no fuel, and that if the fuel were to be 'frozen', the systems would be identical. Fig. 7 shows the comparison for the same initial perturbation from equilibrium. The system with a full, free-flowing fuel tank has a different initial transient and smaller amplitude compared to the original model. This is expected as the perturbation of the structure sets up flow in the fuel and damps the system reducing the amplitude.

Experimental comparison
To validate the SPH-structure coupled model, an experimental rig was constructed consisting of a hinged box. As shown in Fig. 8 the frame was built of aluminium sections, between which perspex panels were mounted and sealed to contain the water. The empty model mass and centre of mass location were measured and the moment of inertia found through comparison of the lumped mass frequency to the actual frequency, while damping in the bearing was characterised for the axle loads of interest here.
Tests were performed with the tank a fill level of 1 8 and time histories measured for a range of starting amplitudes (51,101,151,201,251,301) using an angular encoder. In this work, comparison is only made to the 201 starting amplitude case. Fig. 9 illustrates experimental and numerical angular time histories. The SPH-structure model replicates the measured behaviour well, while the variations further on in the history at smaller amplitudes are likely to be linked to variability in the axle damping. Free surface shapes are compared in Fig. 8 at a series of different times during the motion, and the agreement is commensurate with that seen for rigid dam break simulations (Gomez-Gesteiraa et al., 2012).

Results
For the following results gb=c 2 ¼ 5:435 Â 10 À 4 and ρ f =ρ g ¼ 833:3. These values were chosen based upon reasonable values for a full scale aircraft, and the geometries of the internal and external cases were selected to be representative of aircraft configurations.

Internal tank
Tests were conducted at fill levels of 40% and 67% for an internal fuel tank in a NACA-64A010 aerofoil with structural parameters suggested by Isogai (1979) for a swept back wing. Fig. 10 shows the initial configuration of the test for 67% fill by depth. The fuel is held in the internal tank and the free stream Mach number is 0.85. Fig. 11 shows the system in motion at times given in terms of the pitching natural period.  Figs. 12 and 13 share a similar trend. Addition of the fuel has the effect of raising the flutter boundary, especially for transonic Mach numbers, confirming an extra damping effect from the fuel motion, most likely as a consequence of k sl ¼ kV 1 = ffiffiffiffiffi ffi bg p . An increase in the Mach number therefore increases the sloshing non-dimensional frequency relative to that of the aerodynamic system, which introduces higher damping from higher sloshing frequencies. An exception to this is shown in Fig. 13 for the 67% fill case, where a single point does exist with a lower flutter speed index. The conclusion is therefore consistent with Farhat et al. (2013) where it was determined that ignoring the fuel motion tends to underestimate flutter speeds, but exceptions, as noted, may exist in specific cases. In the case of the structure alone the system tends to go straight from stable to unstable as the flutter speed index is raised, however with fuel present Figs. 12 and 13 show that there is a large region where the system enters a limit cycle, giving a softer stability crossing. This may be due to the fuel damping effect increasing as the motions become more violent. Of the two fill levels, Fig. 14 shows that the 67% level is slightly less stable, possibly as a consequence of a larger component of the mass being able to move freely as fuel.

Results in
The variation of the flutter boundary with different values of the fuel index defined in Section 4 was studied for the internal case. Results are shown in Fig. 15 for a fixed external Mach number of 0.825. The graph shows that for this Mach number the critical fuel index remains nearly constant over a range of flutter speed index values. For both high and low flutter speed indices the fuel index at flutter is greatly and abruptly reduced and past the point shown in Fig. 15 the system is either unconditionally stable with variation in fuel index or unconditionally unstable and as such, no flutter boundary exists.

External tank
A 30% full external fuel tank of the under wing type was also considered. Again the NACA-64A010 aerofoil was used with the same structural parameters. The aerodynamic effects of the fuel tank were not considered; it only interacted with the wing through the forces due to the fuel inside. The tank was 30% full. Fig. 16 shows the initial set up for a free stream Mach number of 0.8 and Fig. 17 shows snapshots in time of the system in motion.
The time marching method for establishing the flutter boundary was used and Fig. 18 shows the flutter boundary obtained. Fig. 18 suggests a damping effect which is on the whole is less than for the internal tank, while behaviour of the system as it becomes unstable is also different. In the internal tank case as the flutter speed index increases the system reaches a limit cycle, which it retains as the speed index increases further, until finally becoming fully unstable and experiencing growth. When the external tank is used the system is stable until the flutter boundary, at which point the beating type phenomenon shown in Fig. 19 is seen together with a larger region of limit cycles. This beating indicates that there are two frequency components that are close together, which both influence the motion. A probable cause is the interaction between vertical fluid motion and surface wave motion; when the system is beating the minima of the amplitude envelope coincide with wave impacts on the tank top. Fig. 20 shows the _ α vs. α phase plots for the two internal tank cases at varying speed indices. The initial limit cycle is clear in Fig. 20(a) and (b), and as the speed index grows these LCOs expand and assume a more complicated structure. This takes place because the increase in speed index raises the aerodynamic forcing in comparison to the spring stiffness, and excites more violent motions of the contained fuel. A key change appears when the accelerations of the motion approach gravity, and the fuel moves away from the walls. The subsequent slamming impacts between the fuel and walls lead to nonlinear behaviour with large, abrupt changes in the forces, and thus the turning points in the phase plots (Figs. 20(e) and (f)). The time histories for these plots shown in Fig. 21 reveal that the turning points are linked to rapid changes in the motion near the peak amplitudes, which is when the the acceleration crosses through gravity and slamming is likely. Figs. 22 and 23 confirm the link to slamming behaviour, as time snapshots just before and just after show clear preslamming and post-slamming configurations of the fuel.

Further analysis of motion characteristics
To help assess the nature of these motions the largest Lyapunov exponents were calculated using the technique of Rosenstein et al. (1993). This approach identifies a log-linear region in the divergence of nearby trajectories plotted against time. For some of the LCOs no log-linear region was observable, indicating non-chaotic behaviour, but for the higher speed indices positive exponents were found, indicating a chaotic system. As shown in Table 1, these ranged from a low of 3:62 Â 10 À 3 up to 7:56 Â 10 À 3 . The presence of a chaotic response would be consistent with the more intricate phase plots, and can be attributed to the slamming behaviour of the fuel.
The amplitude of the LCOs was also explored across the speed index range for the 40% and 67% internal fill cases at M ¼0.825, as shown in Fig. 24. Once the flutter boundary is crossed, LCO amplitude increases as the speed index is raised further. Since the LCO behaviour does not always yield a response of constant amplitude (see phase plots of Fig. 20), clear identification of a single amplitude is not possible and this accounts for the element of noise in Fig. 24. The 67% fill case shows a lower flutter point, but the amplitude growth with speed index is very similar for both fill levels, except at the higher speed indices which begin to approach the point of unconstrained growth.

Conclusions
An SPH code, a structural model and an external aerodynamics solver have been developed and coupled in a timeaccurate fashion. This multi-physics technique allows the non-linear behaviour of the fuel, such as wave breaking and tank top impacts during aeroelastic motion, to be modelled, along with the effects of that fuel motion on aeroelastic behaviour. A flutter boundary has been calculated for different fuel fill levels and it is seen that the addition of fuel usually raises the flutter boundary even when the effect of the fuel added mass is neglected (see also Farhat et al., 2013). Additionally, a region where the system enters a limit cycle is observed. A contrast is also seen between an internal and external tank; the external tank provides less damping than the internal tank. A beating type time history is seen in the limit cycle region of the  3.62 Â 10 À 3 1.10 NC 4.77 Â 10 À 3 7.56 Â 10 À 3 external tank response and it is suggested that this is due to an interaction between the vertical and horizontal fluid motion. It can be seen that there is a critical fuel index value that leads to flutter instabilities. The time marching technique used here is robust and can easily be applied to many geometries of tank including ones with complex baffle configurations. Practical fuel tanks will not have such a simple geometry and allowing the fuel to move in an unrestrained manner may cause high impact loads or stability problems. Future work will therefore consider evaluating a range of tank and baffle geometries for internal and external arrangements.