Sloshing reduced-order model trained with Smoothed Particle Hydrodynamics simulations

The main goal of this paper is to provide a Reduced Order Model (ROM) able to predict the liquid induced dissipation of the violent and vertical sloshing problem for a wide range of liquid viscosities, surface tensions and tank filling levels. For that purpose, the Delta Smoothed Particle Hydrodynamics (δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}-SPH) formulation is used to build a database of simulation cases where the physical parameters of the liquid are varied. For each simulation case, a bouncing ball-based equivalent mechanical model is identified to emulate sloshing dynamics. Then, an interpolating hypersurface-based ROM is defined to establish a mapping between the considered physical parameters of the liquid and the identified ball models. The resulting hypersurface effectively estimates the bouncing ball design parameters while considering various types of liquids, producing results consistent with SPH test simulations. Additionally, it is observed that the estimated bouncing ball model not only matches the liquid induced dissipation but also follows the liquid center of mass and presents the same sloshing force and phase-shift trends when varying the tank filling level. These findings provide compelling evidence that the identified ROM is a practical tool for accurately predicting critical aspects of the vertical sloshing problem while requiring minimal computational resources.


Abstract
The main goal of this paper is to provide a Reduced Order Model (ROM) able to predict the liquid induced dissipation of the violent and vertical sloshing problem for a wide range of liquid viscosities, surface tensions and tank filling levels.For that purpose, the Delta Smoothed Particle Hydrodynamics (δ-SPH) formulation is used to build a database of simulation cases where the physical parameters of the liquid are varied.For each simulation case, a bouncing ball-based equivalent mechanical model is identified to emulate sloshing dynamics.Then, an interpolating hypersurface-based ROM is defined to establish a mapping between the considered physical parameters of the liquid and the identified ball models.The resulting hypersurface effectively estimates the bouncing ball design parameters while considering various types of liquids, producing results consistent with SPH test simulations.Additionally, it is observed that the estimated bouncing ball model not only matches the liquid induced dissipation but also follows the liquid center of mass and presents the same sloshing force and phase-shift trends when varying the tank filling level.These findings provide compelling evidence that the

Introduction
Aircraft wings can deform up to 20% of the wing span when they encounter an atmospheric gust or turbulence.The fuel tanks are located inside of the wings and their liquid weight is comparable to the one of the structural components.The standard engineering practices for wing design do not consider the effect of fuel slosh within the tanks for the determination of the aircraft design loads.The two main tests to identify and verify the dynamic characteristics of the aircraft structure (Ground Vibration and Flight Flutter tests) only include small amplitude excitations which do not significantly excite the main sloshing modes.
Depending on the excitation amplitude sloshing in confined tanks can be a highly nonlinear phenomenon as some early theoretical studies indicate, see for example the classical Refs.[16][17][18].Within the framework of the SLOWD project, the initial experimental test involved examining the free response of a cantilever beam with partially water-filled tanks mounted on it.This investigation successfully demonstrated the presence of additional damping induced by sloshing (see Refs. [19][20][21]).Following this result, Single-Degree-Of-Freedom (SDOF) experiments were developed as part of the SLOWD project to isolate the vertical component of tank motion in order to the study the dissipative behavior of sloshing.Both free decay response (see Refs. [3,4,22]) and forced vertical harmonic excitation (see Refs. [5,23]) test setups were developed for this purpose.
Computational methods are another way of studying violent sloshing problems.There have been recent numerical studies making use of mesh-based techniques, such as finite volumes (see Ref. [24]), the particle-Finite-Element-Method (p-FEM) (see Ref. [25]) or the Consistent Particle Method (CPM) (see Ref. [26]).However, in the present paper, the Lagrangian Smoothed Particle Hydrodynamics (SPH) method is used.Taking this approach and neglecting the influence of the gas phase in the problem, the liquid phase is discretized in a set of particles that interact with one another.The SPH method has been broadly used in engineering and research applications, such as free surface flows (see Ref. [27]) where the method presents and advantage simulating complex geometries (see Refs. [28,29]), or fluid-structure interaction problems (see Refs. [30,31]) and violent sloshing problems (see Ref. [32]).The Weakly-Compressible SPH formulation has recently been applied to the vertical and violent sloshing problem in Refs.[10,11].In these studies, the tank-liquid system is first excited with a sinusoidal imposed motion and the global energy balance of the problem is examined.Then, a validation against the experimental measurements performed in Ref. [22] is carried out coupling the fluid solver with a one degree of freedom system that represents the mechanical behaviour of the experimental tests.Further validations of the experiments carried out in Ref. [22] can be found in Refs.[6,33].Additionally, a thorough comparison between 2D and 3D SPH simulations was carried out and validated against imposed sinusoidal experiments in Ref. [34].The SPH-based method can be con-sidered as an efficient and accurate numerical method when compared to mesh-based CFD solvers to simulate confined and vertical violent sloshing flows.In the present work, the classical SPH formulation has been modified by adding a diffusion term in the continuity equation which helps regularizing the pressure field of the simulations.This formulation is called the Delta Smoothed Particle Hydrodynamics (δ-SPH) formulation as presented in Ref. [35].Further improvements to the delta-SPH formulation have been introduced in Ref. [36] by implementing a Background Mesh scheme.
As part of the SLOWD project, this work involves the use of acceleration data from fluid-structureinteraction (FSI) simulations, performed by employing the δ-SPH formulation (see Ref [35]), describing the SDOF sloshing experiment carried out in the laboratories of the Universidad Politécnica de Madrid (UPM) (see Ref. [22]).The experimental object consists of a single degree of freedom mass-spring system coupled with the slosh dynamics within a Froude scaled tank.The simulated numerical data is utilized to train a reduced-order model capable of accurately replicating the dissipative characteristics observed in vertical sloshing, while accounting for variations in the physical properties of the liquid and the filling level.To ensure the identification of an accurate and comprehensive model for diverse liquid types, an extensive series of numerical simulations is conducted.These simulations involve the systematic variation of parameters associated with viscosity, surface tension, and filling level.The classical ROMs designed to tackle sloshing problems are intrinsically based on potential fluid theory or small lateral perturbations (see Refs. [37][38][39]).In Refs.[18,38], for example, several analytical methods were proposed to study slosh dynamics in simplified tank geometries.Specifically, a realistic representation is provided by the Equivalent Mechanical Models (EMMs), whose parameters can be suitably related with the physical quantities obtained from the linearized potential flow theory [37].However, these models are unable to describe complex nonlinear phenomena such as vertical sloshing, being limited to cases of sloshing caused by small perturbations imposed on the tank.In other words, they fail to reproduce the sloshing forces generated during impacts between liquid and tank.Thanks to the research carried out in the SLOWD project, several reduced-order models have been developed that are able to compensate for these shortcomings.Among them is the bouncing ball EMM proposed in Ref. [4].It is designed by assuming that the fluid behaves like a single solid mass that follows a curved path in free flight, and then bounces elastically off the walls of the tank.When the fluid reaches the top or bottom of the tank, it experiences a smooth visco-elastic force that acts like a spring and a damper.The parameters of this model required tuning via experimental data or full order CFD simulations.Reference [12] introduced another bouncing ball model, trained on the basis of the free response data measured in the one-degreeof-freedom experiment presented in Ref. [22].These models accurately reproduce the dissipative behaviour of sloshing, but have limited capabilities in generalising the estimation of sloshing forces as operating parameters vary (such as the frequency of tank motion).For this reason, a new surrogate model for sloshing was developed in Ref. [40], based on the identification of response surfaces capable of providing estimates of sloshing forces for a wide range of kinematic excitation conditions.Its identification process used direct interpolation of known experimental response.References [15,41] developed instead a data-based ROM consisting of a neural network trained directly with experimental data.Specifically, the model was identified by exploiting the experimental data obtained from the setup introduced in Ref. [5] to investigate the dissipative behaviour of the fluid sloshing inside a tank set in vertical motion.Despite the high performance of these models, their capabilities are limited to one type of liquid and one filling level.As already anticipated, in this work the aim is to overcome this limitation by going on to develop a ROM that also takes into account different liquid types and multiple fill cases.The bouncing ball introduced in Ref. [12] is used for this purpose.Specifically, for each SPH-based simulation case, a bouncing ball model is identified to emulate sloshing dynamics.Then, an interpolating hypersurface-based ROM is defined to establish a mapping between the considered physical parameters of the liquid (varied in the numerical simulations) and the identified ball models.The resulting hypersurface effectively estimates the bouncing ball design parameters while considering various types of liquids, producing results consistent with SPH test simulations.
The paper is organised as it follows.The δ-SPH formulation, together with the SDOF experiment description and numerical model validation are presented in Sect. 2. Section 3 presents the process of identifying the reduced-order sloshing model, including a description of the bouncing ball model and the simulation environment equivalent to the fluid-structure interaction problem under consideration.The article ends with a section presenting the results of tests performed with the identified ROM (Sect.4) and a section with concluding remarks.

δ-SPH formulation
In this section, the weakly compressible δ-SPH formulation that will be used to give the ROM a predictive capability is detailed.This formulation follows the classical SPH formulation with continuity, momentum and particle trajectory equations and adds a density diffusion term in the continuity equation that regularizes the pressure field.Additionally, the minor compressibility of the formulation is modelled through a barotropic equation of state.

SPH model
The flow evolution for a weakly-compressible barotropic fluid is governed by the Navier-Stokes equations and can be expressed in Lagrangian form as: where D Dt represents the Lagrangian derivative, p is the pressure field, g is a generic external volumetric force field, Ω st represents the surface tension force and V the viscous stress-tensor.The flow velocity, u, is defined as the material derivative of a fluid particle position, r.
In order to close the system of equations, the fluid is assumed to be weakly-compressible.This compressibility is modeled with a linear Equation of State (EOS) that relates the fluid density with the pressure variations such that: where ρ 0 is the reference liquid density and c s is the numerical speed of sound that has been chosen to satisfy c s ≥ 10 U ref where U ref is the maximum tank velocity of the simulation.In this way, it is ensured that the density variations of the fluid are below 1 % of ρ 0 with a Mach number below or equal to 0.1.For a Newtonian fluid, the divergence of the viscous stresstensor can be expressed as: where D is the strain rate tensor, D = (∇u + ∇u T )/2, I the Identity tensor and λ and μ the bulk and dynamic viscosity coefficients respectively.In Eq. 3 it is assumed that for a weakly-compressible fluid tr(D) = ∇•u ≈ 0. The δ−SPH formulation detailed in Ref. [35] is used in this work.The set of Eqs. 1 turns in SPH into: where subscript j represents a general neighbor of particle i.The volume of the neighbor particle is denoted as V j , defined for a fluid particle as V j = m j /ρ j .The kernel function is W i j = W (r j − r i , h) and h represents the smoothing length.Therefore, ∇W i j represents the spatial derivative of the kernel function with respect to the i − th particle.In the work presented in this paper, a Wendland C2 kernel with h/dr = 2 is used, being dr the distance between particles.The constant μ sph is the maximum between the physical dynamic viscosity μ and the artificial one as stated in Ref. [42].It can be expressed as: where K = 8 for 2D simulations and α is the artificial viscosity constant that is set to 0.01 for all simulations except for the ones where the variation of the dynamic viscosity of the liquid is analysed.In these cases the artificial viscosity term is set to zero (α = 0 and μ sph = μ) to retain the physical dynamic viscosity of the liquid.Furthermore, D i j and Π i j are the density and momentum diffusivity terms added to the continuity and momentum equations respectively which constitutes the standard δ-SPH model (see Ref. [35]).They are expressed as: where ∇ρ L i is the renormalized density gradient, defined by Ref. [43] and j * refers exclusively to fluid neighbor particles.The reader is encouraged to see Ref. [6], where the convergence property in terms of particle resolution of the density diffusivity term is evaluated for different liquids with surface tension effects.
The turbulence of the problem could be tackled using a LES subgrid model such as the one presented by Ref. [44] and later on applied in Refs.[10,11], however, it is not the scope of this work to deal with the turbulence of the vertical sloshing problem and it is assumed that the extra viscosity added by the Π i j term is broadly representative of the turbulent viscosity.Additionally, an analysis of the liquid viscosity in the sloshing induced damping is performed in Ref. [6].
In order to avoid numerical cavitation, also known as tensile instability a Particle Shifting Technique (PST) is used represented by the δu i term in the particle trajectory evolution equation.Following the PST formulation presented in Ref. [45], δu i is calculated as: Eq. 7 is corrected near the free surface according to the formulation presented in Ref. [46].An additional correction is applied to the PST term close to the free surface to avoid shifting velocities normal to the free surface.For this purpose, the SPH solver uses a free surface tracking algorithm following the work done by Ref. [6].
In the momentum equation of the set of Eqs. 4, a surface tension field Ω st i is added to model the free surface forces that have an influence in the free surface shape, specially at initial stages, when vertical Rayleigh-Taylor instabilities are triggered (see Ref. [47]).For this purpose, the formulation detailed in Ref. [48] is implemented in this work, and it reads: where n k is the normal at the free surface for a general particle k, calculated as: where λ represents the minimum eigenvalue of the renormalization matrix defined by: For the time integration of the set of Eqs. 4 a purely explicit formulation with a Leapfrog based predictorcorrector scheme is used, with the main advantage that only one time derivative is needed for each time step (see Ref. [49]).
For the solid walls of the tank, the numerical boundary integrals approach, as described by Ref. [50] is followed.In the boundary integrals formulation, the differential operators from Eqs. 1 are redefined near a solid boundary as: where f i is a generic field function, n j represents the normal of the surface pointing outwards, s j is the area of the surface element and γ i the Shepard renormalization factor, that is calculated geometrically through a semi-analytical approach as explained in Ref. [51].It should be noted that close to a boundary the neighbor particles j are divided into two groups: fluid particles j * and boundary particles j.In this work, the expression to compute all gradients will follow this subdivision procedure following the methodology presented in Refs.[52,53].

SDOF model
The numerical δ-SPH formulation described in the previous section is applied to the sloshing problem of the vertically accelerated SDOF experiment depicted in Fig. 1.The experimental setup consist of a Froudescaled tank filled with liquid up to a filling height h l which is attached to a mechanical guide and a set of springs that constraints the tank motion to be purely vertical.The tank features a base area of 60 cm 2 and a height of 6 cm.It is contained in a C-shaped wooden structure that is attached to a set of six springs (half on the upper side and half on the lower side).Each individual spring presents an stiffness of k i = 708.67N/m which corresponds to a total system stiffness of K = 4252.03N/m.The upper set of springs is attached to a metallic plate that is joined to the embedded load-cell.The lower set of springs are attached to the floor.The system is displaced a certain initial distance z T = z T 0 from its equilibrium position until it is fixed by the action of a pair of permanent solenoids.Then, when the current is turned on the solenoids release the system starting the oscillatory motion and the beginning of the experiment.Three quantities are measured from each test: the tank acceleration, the tank position and the force at the load-cell.Using these magnitudes, the vertical liquid motion can be derived following the methodology explained in Ref. [22].
The SDOF system can be modeled using the following equation: where z T is the vertical tank displacement.The structural mass is denoted by m s , B 0 = 0.38 N and B 1 = 1.73 kg/s account for the structural damping coefficients of the rig, K = 4252.03N/m is the total spring stiffness and F slosh the vertical liquid load which its calculation will be detailed later.The liquid weight m l g is subtracted from F slosh so only the dynamic part of the force is included in the system reported in Eq. 12.The liquid and structural masses of the system are balanced depending on the filling level and the type of liquid used in order to have a total mass m = m s + m l = 2.58 kg so that the natural frequency of the total system is fixed to 6.46 Hz.The sloshing force acting on the SDOF system is calculated from the SPH solver according to 123 Fig. where F p slosh and F V slosh are the corresponding pressure and viscous components of the vertical sloshing load.According to the boundary integrals formulation, these forces can be expressed in the SPH model as: where the second term of the viscous force F V slosh is accounted only in case no slip boundary conditions are enforced.A similar procedure to compute the forces is followed by Ref. [54].

Validation against experimental decay tests
The δ-SPH solver used is AQUAgpusph, an opensource code developed by CEHINAV (see Ref. [50]), which is validated against for the violent and vertical sloshing case that was experimentally assessed in Ref. [22].For the validation, 2D simulations have been carried out for 10 g, 50% filling level water and oil experiments.Figures 2 and 3 show the comparison between the SPH simulation and experimental data for the tank acceleration and vertical sloshing force in both water and oil tests.
The tank motion seems to be well captured by the SPH simulations for both water and oil tests in terms of tank frequency and envelope of the decay.Regarding the vertical sloshing force, the SPH simulations tend to present pronounced spikes close to the force maxima due to the noise coming from the integrated pressure field.Overall, the match between simulations and experiments is good.Further details on the validation of the AQUAgpusph code with respect to the vertical sloshing problem can be found in Refs.[6,33].
Figure 4 displays experimental and SPH snapshots of some of the characteristic instants of the flow kinematics.The experiment starts with a flat surface and when it is accelerated upward a ripple is formed in the free surface (t = 0.075 s) which then develops into Rayleigh-Taylor instabilities that impact the upper tank wall (t = 0.1 s).After that, the flow becomes highly turbulent and violent with bubble formation and many liquid-to-liquid and liquid-to-wall impacts (t = 1.43 s).At this stage, the flow kinematics is not captured to a great level of detail by the δ-SPH simulation since 3D and multiphase effects are involved.However, the relative motion of the bulk liquid mass position with respect to the tank motion is well captured.Finally, the liquid restores its initial state displaying a standing wave like regime (t = 3.31 s).The reader is encouraged to watch and download the videos of the experiments and SPH simulations in http://canal.etsin.upm.es/files/SLOWD/Fig. 2 Comparison between δ-SPH simulations and experimental data of the acceleration of the tank over time for water and oil decay tests at 50% filling level and 10 g of maximum acceleration Fig. 3 Comparison between δ-SPH simulations and experimental data of vertical sloshing force over time for water and oil decay tests at 50% filling level and 10 g of maximum acceleration videos/.At the bottom side of Fig. 4 the pressure fields at time instants t = 0.075 s, 0.1 s, 1.43 s and 3.31 s are shown.It can be seen that, despite the addition of the diffusion term in the continuity equation, the pressure field displays spurious acoustic noise which affects the accuracy of the integrated sloshing force.This source of error may potentially come from the explicit timestepping.To overcome this issue two approaches could be adopted: implement an implicit time integration scheme as done in Ref. [55] (in which case the use of artificial diffusion terms could also be avoided) or use Velocity-divergence Error Mitigating (VEM) and Volume Conservation Shifting (VCS) schemes that ensure a divergence-free velocity field (see Ref. [56]).

Sloshing reduced order model identification
This section presents a comprehensive procedure for identifying a reduced-order model that accurately replicates the dissipative behavior of sloshing in a fluidstructure interaction problem.The section begins by providing a detailed description of the bouncing ball model, with its constitutive parameters and underlying principles.Next, the simulation environment developed to reliably emulate the fluid-structure interaction problem is introduced.The bouncing ball model is integrated into this environment to reproduce the nonlinear dynamics of vertical sloshing.Finally, the section illustrates the ROM identification logic, which uses an interpolating hypersurface to establish a mapping starting Fig. 4 Experimental snapshots (top) of a vertical sloshing experiment with water at a 50% filling level and 10 g of maximum acceleration at t = 0.075 s, 0.1 s, 1.43 s and 3.31 s.Free surface evolution (middle) simulated with the δ-SPH method at the same instants.Nondimensional pressure field (bottom) of the δ-SPH simulation at the same instants from physical fluid parameters and fill level, resulting in the estimation of a bouncing ball model that accurately simulates the desired sloshing behavior for a generic liquid.

Bouncing ball model
The bouncing ball model used in this work was introduced in Ref. [12].It is an equivalent mechanical model able to emulate the mechanics of the impacts that can occur between the sloshing fluid and the walls of the tank, when the latter is excited violently in a vertical direction.In particular, the sloshing forces are replaced by the forces exchanged between the tank and a ball bouncing inside the rigid tank.Figure 5 shows the design scheme of the bouncing ball.
A part of the total liquid mass m l = m b + m f is associated to the ball m b , whereas another portion is considered frozen m f , namely attached to the wall.The equation of motion of the bouncing ball can be thus summarized as follows where F b is the force exchanged between the tank and the ball, z T (t) represents the motion of the tank whereas z b is the absolute vertical motion of the bouncing ball.The frozen mass m f is expressed as m f = (1 − β) m l and β is a nondimensional parameter that can take values between zero and 1 and regulates how much mass of the liquid is considered frozen.A variable s(t) that represents the relative motion of the ball with respect to the tank is defined as: e l s e w h e r e (16) where r 0 is the radius of the rigid ball and h is the height of the tank.When the condition, z b (t) < z T (t) + r 0 occurs, the ball is in the impact region with the floor of the tank.Whereas, the condition z b (t) > z T (t)+h −r 0 is related to impact with the ceiling.The new variable s(t) is null when the ball finds itself in the suspension or floating phase, where there is no impact.Moreover it is introduced the relative velocity of the ball with respect to the tank ν = żb − żT .It follows that we can define the viscoelastic forces as where k b and c b are, respectively the stiffness and damping associated to the bouncing ball.It is worth noting that k b = kb f nl may eventually be, in turn, nonlinear function of s(t) by introducing a penalty function f nl that avoids the ball to go out by the limits of the tank.The penalty function coefficient α is one of the design parameters of the bouncing ball, which defines the intensity of this nonlinear function.The overall sloshing force exchanged between the tank and fluid mass inside consists of two terms yielded by bouncing ball force and frozen mass inertia:

Simulation environment and bouncing ball tuning
Numerical modeling of fluid-structure interaction problems is performed in Sim-ulink , casting the bouncing ball model as representative of vertical sloshing dynamics (as it was done in Ref. [12]). Figure 6 displays the flowchart of the complete system comprising the structure and sloshing.In this figure, the sloshing forces denote the input of the structural subsystem, while the input of the sloshing subsystem includes the vertical rigid displacement and the tank velocity.This numerical model ideally reproduces the dynamic behaviour of the SDOF experimental (or numerical) system presented in Sect. 2. Note that, the bouncing ball based sloshing block can also be simulated in isolation from the structure, with an appropriate seismic excitation assigned as input.
The identification of the bouncing ball model is performed following the same procedure as in Ref.  [12], in which the five design parameters (r 0 , kb , c b , α, β) of the model are estimated through an optimisation process.Specifically, this procedure uses the numerical data of vertical tank acceleration calculated in the SPH simulations of the FSI problem.From the envelope of the acceleration signal, the instantaneous damping ratio is derived by means of the logarithmic decay.By reporting this quantity as a function of the amplitude of acceleration, the influence of viscous and frictional (Coulomb) structural damping is also highlighted.The instantaneous damping ratio defines the objective function of the optimisation procedure.In other words, the bouncing ball model is identified by finding a combination of parameters that when simulated in the FSI environment (see Fig. 6), returns the closest instantaneous damping ratio to the one retrieved by the decaying SPH acceleration data.The objective function that is minimized in the optimisation process , where ζ sph i and ζ bb i represent respectively the instantaneous damping obtained from the SPH acceleration data and that obtained from the simulations with the ball.Note that the summation in the objective function is performed on the damping values as the acceleration varies (N values corresponding to the acceleration values considered).

Hyper-surface ROM of the bouncing ball parameters
The bouncing ball presented in the previous subsections is used to create a dataset aimed at identifying a reduced-order model capable of estimating a new bouncing ball model, for a generic liquid type and without performing the optimisation process.This ROM consists of an interpolating hypersurface that maps the physical parameters of the liquid (together with the fill The logic of what is proposed in this work is schematized in Fig. 7. For the development of this map, many SPH simulations were carried out by varying the physical parameters in order to have a sufficient of points in the three-dimensional space of physical quantities -which then form the design of experiment (DoE).Specifically, 54 randomly selected free response SPH simulations were performed at the same operating natural frequency, initial amplitude, tank dimensions and structural parameters.The physical parameters that were varied are the following ones: -Filling level fill: [10-90] % -Reynolds number Re: [62-1257480] -Weber number We: [2640-26400] (corresponding to γ : [10-90]) All these simulations started at the same initial conditions, equal to z T (0) = −0.0584m and żT (0) = 0.Among the 54 simulation cases we define the reference case for which we have fill = 50%, Re = 140000 and We = 3668.The Re and We values given for the reference case correspond to a test with water.As each SPH simulation is performed at a constant frequency value, the identification process must be repeated if the system is operated at other frequencies.In fact, the bouncing ball is unable to adequately extrapolate the dissipative behaviour of sloshing at frequencies other than the identification frequency.The identification procedure presented in Sect.3.2 is used to derive the optimal bouncing ball design parameters for each of the FSI fluid-dynamic simulations.Each identified bouncing ball will correspond to a set of physical parameters, associated with the SPH simulation used as the target of the optimisation.By employing this method, a data set of parameters can be generated to enable the development of the interpolating response model.The resulting model provide estimation of ball parameters for a different set of physical liquid parameters, distinct from those utilized in SPH simulations.The notable advantage of this model is its ability to remove the need for the optimisation process to identify a new ball, resulting in a significant reduction in computational cost.

Results
This section presents a brief overview of the obtained results, illustrating the comparisons of achievements in terms of acceleration, damping, and sloshing force.First, the design of the experiment on which the reduced-order model is built is presented, followed by the simulation results for both the reference and test cases.A further analysis of the performance of the identified model with respect to important aspects related to the physics of vertical sloshing is also conducted.

Reference case
Following the identification of the bouncing ball model for each of the 54 triplets of physical parameters, a set of ball design parameters is obtained.
Figure 8 shows how the 54 points are distributed in the three-dimensional space defined by the parameters Re, We and fill.Each of them is associated with a colour indicative of the value of the objective function J (defined in Sect.3.2), obtained after the identification process of the bouncing balls.With the design parameter sets associated with each of the bouncing ball models identified at the 54 points, an interpolation process is implemented in MATLAB , aimed at obtaining a response surface model (hypersurface) capable of providing the prediction of the five ball parameters based on a set of physical parameters different from those selected for the starting simulations.For the sake of brevity, of all the SPH simulations considered, only the results obtained for the reference case mentioned in Sect.3.3 are shown in Fig. 9.The rest of the simulation cases and bouncing ball responses can be found and downloaded at http://canal.etsin.upm.es/ftp/SLOWD_DATABASE/ROM_data/SPH_BB/. The optimisation process presented in Sect.3.2, results in a bouncing ball having the constitutive parameters shown in Table 1.
They lead the definition of a bouncing ball model capable of exhibiting dissipative behavior comparable to that observed in the reference SPH free response simulation.In fact, as can be seen from Fig. 9a and b, the identified bouncing ball model is able to return an acceleration and an instantaneous damping ratio that present a good match with the numerical ones.The force predicted by the model can also be compared to that obtained with the SPH test.As can be seen from Fig. 9c, the equivalent mechanical model is able to capture the nonlinear trend of the load exerted by the liquid.However, the force predicted by the ball turns out to be more impulsive, and this is due to the way it was designed, namely, as a viscoelastic impact force.

Blind test case
Similar results were also obtained for the remaining simulations.It should also be noted that, although not included in the objective function to be optimised, the sloshing force that the bouncing ball is able to provide is, for each of the tests, very similar to the SPH force (again leaving aside the presence of somewhat more impulsive peaks present in the first cycles of the response, where the fluid impacts are more violent).
To validate the process of identifying the interpolating surface, an additional test with physical parameters different from those used for the other simulations is considered.However, these always fall within the ranges defined in Sect.3.3, and correspond specifically to fill = 90%, Re = 251496 and We = 4400.
Figure 10 shows the results obtained, again comparing acceleration, damping and sloshing force.In particular, the curves indicated in the legend with opt, are those obtained by performing the direct optimisation aimed at identifying the ball parameters (as presented in Sect.3.2).While, the curves indicated with int, are those obtained by performing free response simulations with a ball having the parameters predicted by the identified interpolating surface (in correspondence with the physical quantities of the additional test).The above parameters are compared with those obtained from direct optimisation in Table 2.
The parameters kb and α shown in Table 2 have a considerable relative error of 13% and 25% respectively, while c b and β are more close (both have an error of 7%).The radius of the ball r 0 , on the other hand, presents a negligible relative error when compared to the identified ball model.However, these parameters make it possible to obtain a bouncing ball model that, once simulated, is able to return accurate responses.This is also seen in the comparisons in Fig. 10, where for the acceleration, instantaneous damping ratio and sloshing force curves there is almost an overlap for the two bouncing ball models (direct optimisation and interpolating response surface).In particular, the distance between the damping curves obtained by direct optimisation and by interpolation (see Fig. 10b) is again quantified by the J function, which in this case is equal to 0.022.This value is comparable to the lower values obtained for the simulations shown in Fig. 8.
Taking advantage of the test case presented in this section, it is also possible to perform a sensitivity analysis as the number of physical parameter triplets that are considered to generate the interpolating surface changes.Specifically, our objective is to evaluate whether an increase in the number of considered simulations leads to improved accuracy in predicting In terms of computational costs, the direct optimisation process aimed at identifying the parameters of the bouncing ball takes 12.3 minutes, while the construction of the interpolating surface with the full number of simulations available takes 2.5 minutes.Once available, the interpolating surface makes it possible to obtain the parameters of the bouncing ball in considerably less time than is necessary for the direct optimisation.It is important to note that free response simulation in Simulink is of the order of a few seconds so it was not considered when evaluating the computational cost.In comparison, an SPH simulation run with AQUAgpusph took approximately 6 h to complete using stateof-the-art GPU technology which translates in 95.86 % reduction of computing time.

Comparison of sloshing magnitudes
This section explores other aspects related to the physics of vertical sloshing, and how the bouncing ball model is able to describe them.Fig. 12 shows three snapshots where the ball is impacting the tank ceiling (t = 0.4 s), floating (t = 0.44 s) or impacting the tank floor (t = 0.47 s).In the δ-SPH conunterparts we see that the bulk liquid is located at the ceiling or floor of the tank when the ball is impacting these walls and we find that the fluid is scattered across the tank space when the ball is floating.Similar observations were made in Ref. [12] when the bouncing ball model was compared with experimental snapshots.Building upon this rationale, a comparative analysis is proposed between the position of the liquid's center of gravity (COG) (from the SPH simulation) and the vertical displacement of the bouncing ball mass.The COG of the liquid in the SPH simulation z COG , is calculated in the same way as a discrete system of particles: where z i and m i are the vertical coordinate and mass of particle i, N is the total number of particles and m l is the total liquid mass.Fig. 13 shows the comparison between the velocity of liquid COG and the velocity of the bouncing ball, for the reference case introduced in Sect.3.3.Overall, there is a good match between the bouncing ball model and SPH simulation except for the first 4 cycles when the tank amplitudes and sloshing forces are the highest.For these first cycles, the flow is highly fragmented and the bouncing ball model tends to underestimate the COG amplitude.The good match between ball velocity and liquid COG velocity is confirmed in the analytical study conducted in Ref. [15], which was further corroborated through numerical validation in Ref. [33].These references emphasise the direct correlation between the dissipative component of the vertical sloshing force, which includes liquidinduced dissipation, and the relative motion of the liquid with respect to the tank.Therefore, if the bouncing ball model matches the dissipation of the problem by optimising the objective function J , it should match the COG of the liquid.
Another way of checking that the bouncing ball model is properly reproducing the sloshing physics in the vertical sloshing problem is by performing a sensitivity analysis with respect to filling level variation.The influence of the filling level in the liquid induced power dissipation and the phase-shift between the tank motion and the vertical sloshing force has been experimentally confirmed in several studies (see Refs. [3][4][5]).In these experiments, the maximum dissipation, computed as the global damping ratio of the experiment ξ , is found close to half of the filling level fill ∈ [40,60].On the other hand, the maximum sloshing force F increases with the filling level and the phase-shift ϕ decreases as the fill level increases.To see how these magnitudes are computed the reader is encouraged to see Ref. [3].These trends have been reproduced by the SPH and bouncing ball simulations and shown in Fig. 14, where they have been normalized with respect to the reference case (corresponding to magnitudes with subscript 0).While the bouncing ball model may not precisely reproduce the exact values of the magnitudes, it successfully captures the overall trends.Hence, despite the simplifications inherent in the model and considering the predictive capability offered by the response surface, it can be concluded that the bouncing ball model is a valuable tool for investigating vertical and violent sloshing flows within the defined parameter space.

Conclusions
The approach presented in this paper makes it possible to identify a reduced-order model for vertical sloshing on the basis of numerical data from SPH fluid dynamic simulations.The novelty lies in the possibility of modelling the sloshing phenomenon for a wide range of possible fluids with different intrinsic properties, i.e. different viscosity and surface tension.But also for different tank levels.Therefore, by assigning the desired physical parameters, it is possible to obtain an equivalent mechanical model, represented by a bouncing ball, capable of emulating the dynamic behaviour of a liquid contained in a tank subjected to vertical excitation.In particular, this model will reproduce the same dissipative behaviour that the fluid induces in the system by impact.
By building up a database of numerical simulations in which the fluid properties are varied, hypersurfaces can be constructed that map the physical properties to the bouncing ball parameters.This mapping procedure gives the ROM a predictive capability which has been tested against a blind test.The predicted response of the ball, when validated with the blind test, agrees well with the model identified with the optimisation procedure and the numerical SPH data.In addition, a sensitivity analysis is performed to see how the objective function is affected by the number of simulations used to build the interpolation hypersurfaces.It is seen that the accuracy of the procedure increases as the number of randomly selected cases increases.
Having validated the predictive capability of the ball model, the ROM is used to investigate other aspects of the vertical sloshing problem.It is seen that the relative motion of the ball is consistent with the liquid COG motion of the SPH simulation.In addition, a level variation study is carried out and it is observed that the bouncing ball model captures the maximum sloshing force and phase-shift trends of the study.This observation indicates that the bouncing ball model could be used as a tool to predict some key characteristics of the vertical and violent sloshing problem while consuming very low computational resources compared to other CFD tools.
Despite the strong performance observed with the bouncing ball, it is important to note that its estimation accuracy decreases when the frequency deviates from the specific value for which it was identified.Reducedorder models based on neural networks are able to bypass this problem, allowing the accurate estimation of forces in a larger operational domain.However, the model identified in this work may be preferable, as the use of neural networks carries the risk of obtaining an unstable response outside the identification input space.In such cases, a physics-based model like the bouncing ball can be advantageous as it provides a response that, while potentially less accurate, remains stable and does not lead to divergence.

Fig. 6
Fig. 6 Diagram of the fluid-structure interaction problem implemented in Simulink

Fig. 7
Fig. 7 Outline of the mapping procedure between fluid parameters and bouncing ball parameters

Fig. 9
Fig. 9 Comparisons between SPH data and bouncing ball predictions in terms of tank acceleration, instantaneous damping ratio and vertical sloshing force for the reference case (Re = 140000, We = 3668 and fill = 0.5)

Fig. 10 Table 2
Fig. 10 Comparisons between SPH data and bouncing ball predictions in terms of tank acceleration, instantaneous damping ratio and vertical sloshing force for the blind test case (Re = 2514960, We = 4400 and fill = 0.9)

Fig. 11
Fig. 11 Trend in objective function J as the number of simulations considered in identifying the interpolating surface varies

Fig. 12 Fig. 13
Fig. 12 of snapshots coming from a δ-SPH simulation (top) and a BB model simulation (bottom) for instants t = 0.4 (floating), 0.44 (ceiling impact) and 0.47 (floor impact) seconds.The sloshing forces predicted by SPH F sph slosh and the bouncing ball model F bb slosh are included for comparison

Fig. 14
Fig. 14 Global damping ratio, sloshing force maximum and phase-shift as a function of the filling level.Comparison between the SPH simulation and the BB model for water tests at Re = 140000 and We = 3668.All magnitudes have been normalized with respect to the fill = 0.5 case