Sloshing dynamics estimation for liquid-filled containers performing 3-dimensional motions: modeling and experimental validation

Many industrial applications require the displacement of liquid-filled containers on planar paths (namely, paths on a horizontal plane), by means of linear transport systems or serial robots. The movement of the liquid inside the container, known as sloshing, is usually undesired, thus there is the necessity to keep under control the peaks that the liquid free-surface exhibits during motion. This paper aims at validating a model for estimating the liquid sloshing height, taking into account 2-dimensional motions of a cylindrical container occurring on a horizontal plane, with accelerations up to 9.5 m/s2. This model can be exploited for assessment or optimization purposes. Experiments performed with a robot following three paths, each one of them with different motion profiles, are described. Comparisons between experimental results and model predictions are provided and discussed. Finally, the previous formulation is extended in order to take into account the addition of a vertical acceleration, up to 5 m/s2. The resulting 3-dimensional motions are experimentally validated to prove the effectiveness of the extended technique.


Introduction
The transport of containers filled with liquids finds application in several industrial applications, e.g., in food & beverage or pharmaceutical production and packaging lines. Typically, the manipulation of such containers is assigned to linear transport systems or industrial serial robots; in many cases the required motion follows planar curves. The prediction of the liquid movement inside the container, referred to as sloshing, is important to prevent the liquid from overflowing. Additionally, a reliable sloshing prediction model can be exploited to limit the stirring of the liquid during task execution.
For this purpose, machine-learning methodologies are presented in [2] and [3], where, starting from data collection, predictive algorithms are built to inspect the behavior of discrete liquid particles inside a cylindrical container. This technique, though very powerful, requires experiments to be run in advance, together with a not negligible computational effort.
In [4] and [5], the Finite Element Method (FEM) is adopted for the analysis of sloshing in rectangular containers, requiring a preliminary generation of the mesh able to replicate the liquid behavior.
The meshless Smooth Particle Hydrodynamics (SPH) method is employed in [6] and [7] to model the sloshing by discretizing the liquid in tens of thousands particles: the simulations accurately match the experimental data, at a cost of days of computation.
In [8] the coefficients of the nonlinear sloshing dynamics model presented in [9] are provided to evaluate the sloshing height for 3-dimensional motions, leading to a complex formulation, which may be difficult to use.
A ready-to-use and fast alternative is represented by the development of equivalent discrete mechanical models. The literature considers two main discrete approaches for the modeling of sloshing dynamics inside a container subjected to 2-dimensional planar motion [10]: a spherical pendulum and a 2-DOF mass-spring-damper system. In the former case, the generalized coordinates describing the system are the angles defining the position of the pendulum mass, whereas in the latter they are the mass displacements from the reference position. Although being intuitive, the use of the angular coordinates of the pendulum mass to assess the sloshing behavior of the liquid (see [11,12]) lacks physical meaning, in particular when the knowledge of the liquid peak height is important. For this reason, in the spherical pendulum model used in [13], [14], and [15], the sloshing height is estimated by means of the tangent functions of the spherical coordinates. However, estimating the sloshing height by means of the tangent of the pendulum angles may lead to singularity conditions when the container acceleration is high, since in this case these angles can approach 90 • and the tangent tends to assume unrealistic high values.
To overcome this drawback, a novel approach, based on the mass-spring-damper model [16], is proposed in [17] for the sloshing-height estimation. This model is validated for 1-dimensional motions in [17] and it is exploited in [18] and [19] to plan antisloshing trajectories. The same technique is used in a software application presented in [20] to execute simulations of liquid sloshing in cylindrical and rectangular containers. In [17], the authors propose the possible extension to planar motions, but no experimental validation is provided. The latter is the objective of this paper, particularly referring to 2-dimensional planar motions of a cylindrical container, with accelerations up to 9.5 m/s 2 . With respect to the preliminary conference version presented in [1], this contribution also reports an extension of the formulation to 3-dimensional motions comprising a vertical acceleration up to 5 m/s 2 . Experiments are also provided for this case.
The paper is structured as follows. Section 2 presents the model parameters and the equations of motion (EOMs) in terms of the corresponding generalized coordinates. Section 3 Fig. 1 Mass-spring-damper model provides the formulation of the sloshing-height estimation. In Sect. 4, the results from an experimental campaign are described and discussed. Lastly, in Sect. 5 conclusions are drawn and suggestions for future developments are given.

Model parameters
We will consider a cylindrical container of radius R, filled with a liquid of height h and mass m F . A simplified discrete mechanical model can be used to reproduce the liquid-sloshing dynamics. In particular, the mass-spring-damper model comprises a rigid mass m 0 (whose signed vertical distance from the liquid's center of gravity G is h 0 ) that moves rigidly with the container, and a series of moving masses m n , with each one of them representing the equivalent mass of a sloshing mode (Fig. 1a). Each modal mass m n is restrained by a spring k n and a damper c n , and its signed vertical distance from G is h n .
The model parameters can be determined by imposing a number of equivalence conditions with the original system [10]: • the overall mass must be the same, (1) • the height of the center of gravity G must remain the same for small oscillations of the liquid, • the natural frequency associated with the nth mode must coincide with the one that can be obtained from the continuum model, • the sloshing force acting on the container wall must be the same as the one calculated from the continuum model, leading to the determination of the nth sloshing mass, In equations (3) and (4), ξ 1n is the root of the derivative of the Bessel function of the first kind with respect to the radial coordinate r, for the 1st circumferential mode and the nth radial mode [21], while g is the gravity acceleration. The damping ratio ζ n = c n 2 √ k n m n can be determined by using empirical formulas [10]. For a container under 2-dimensional motion on the horizontal xy plane, the excitation is provided by the container accelerations along the x and y directions, denoted withS 0 = [ẍ 0ÿ0 0] T . The motion of the nth sloshing mass is described by the generalized coordinates (x n , y n ), whose definition is illustrated in Fig. 1b. The latter are then used to compute the liquid sloshing height.

Equations of motion
In general, three dynamic regimes are possible [10]: • small oscillations in which the liquid free-surface remains planar (Fig. 2a); • relatively-large-amplitude oscillations in which the liquid free-surface is no longer planar (Fig. 2b); • strongly nonlinear motion, where the liquid free-surface exhibits instantaneous peaks characterized by swirling shapes.
While the third motion regime will not be the object of the present study, the first and second cases can be analyzed by means of a linear mass-spring-damper model (L model) and a nonlinear mass-spring-damper model (NL model), respectively. The NL model considers the sloshing mass m n as sliding on a parabolic surface, with an attached nonlinear spring of order w (Fig. 1c) [16]. The analytical expression of the parabolic surface allows writing the vertical coordinate z n as a function of x n and y n , namely where C n = ω 2 n R g . The time derivative of equation (5) iṡ z n = C n R (ẋ n x n +ẏ n y n ).
The nonlinear spring exerts the forces α n k n x 2w−1 n and α n k n y 2w−1 n , along the x and y direction, respectively. In this paper, we choose w = 2 and α n = 0.58, as suggested in [16]. If the radial generalized coordinate r n = x 2 n + y 2 n is introduced, the nonlinear-spring force in the radial direction can be written as α n k n r 2w−1 n . The EOMs, describing the time evolution of the generalized coordinates (x n , y n ), can be obtained by means of the Lagrange equations: where • the kinetic energy T of the nth sloshing mass can be computed by taking into account its velocityṡ n = [ẋ nẏnżn ] T , the container velocityṠ 0 = [ẋ 0ẏ0 0] T and by exploiting equation (6): • the potential energy V considers the contribution of gravity and nonlinear-spring forces, namely V = m n gz n + rn 0 α n k n r 2w−1 • the Rayleigh function D accounts for energy dissipation, The substitution of equations (8), (9), and (10) in the system (7) leads to the formulation of two coupled EOMs for the NL model: y n + 2ω n ζ n [ẏ n + C 2 n (y 2 nẏn + x nẋn y n )]+ +C 2 n (y nẏ 2 n + y 2 nÿn + y nẋ 2 n + y nẍn x n )+ where x n = x n /R, y n = y n /R. As far as the L model is concerned, the linearization of the EOMs in equation (11) provides two decoupled EOMs in the generalized coordinates (x n , y n ) of the nth mode: ẍ n + 2ζ n ω nẋn + ω 2 n x n +ẍ 0 = 0, y n + 2ζ n ω nẏn + ω 2 n y n +ÿ 0 = 0. (12)

Extension to 3-dimensional motion
If an excitationz 0 along the z-axis is added to that on the xy plane, the container is subject to a 3D motion, namelyS 0 = [ẍ 0ÿ0z0 ] T . As long as we assume to employ the massspring-damper model presented in Sect. 2.1 to reproduce the liquid behavior, the same model parameters can be used. This choice allows the derivation of a fast and easy model extension, without the complication inherent in the construction of a different discrete model. As a consequence, we assume that the additional motion along the z-axis only influences the kinetic energy of the nth sloshing mass: whereż n is still given by equation (6). Hence, combining equation (7) with equations (13), (9), and (10), the NL-model EOMs for the 3D motion become whereas the L-model EOMs yield ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ẍ n + 2ζ n ω nẋn + ω 2 n x n +ẍ 0 +z 0 g ω 2 n x n = 0, y n + 2ζ n ω nẏn + ω 2 n y n +ÿ 0 +z 0 g ω 2 n y n = 0.

1-Dimensional motion
If only a 1-dimensional excitation in the y direction is provided and the phenomenon of rotary sloshing is neglected [10], solely the generalized coordinate y n is different from zero.
In such a case, the conservation of the center of gravity y-coordinate, between the continuum model and the equivalent model, yields Considering a cylindrical container with cross-section S = πR 2 , filled with a liquid of height h, y G can be computed as where the function η(r, θ, η n ) describes the liquid free-surface shape, η n is the sloshing height of the nth mode, (r, θ ) are the polar coordinates, with x = r cos θ , y = r sin θ , and dS = r dθdr. As for the L model, the function η(r, θ, η n ) describes a plane (Fig. 2a), whereas, for the NL model, the nonplanar free-surface can be described by means of the first-kind Bessel function (Fig. 2b), namely Independently from the adopted function η, the expression of y G from equation (17) can be used in equation (16) to express η n as a function of the model parameters and the generalized coordinates (x n , y n ), with the latter being obtained by solving the EOMs (see Sect. 2.2). The L-model assumption of planar surface leads to Regarding the NL model, by exploiting one of the Bessel function properties, i.e., , y G can be evaluated as y n from (11) y n from (14) Table 2 MSH estimation for a 2-dimensional planar motion, without (left column) and with (right column) an excitation along the z-axis x n , y n from (12) x n , y n from (15) NL model x n , y n from (11) x n , y n from (14) Inserting the results from equations (20) and (21) into equation (16) allows the formulation of the nth sloshing height (SH) for the L and NL models, respectively, as shown in the leftmost column of Table 1.

2-Dimensional motion
When accounting for a 2-dimensional excitation, the plane , on which the maximum sloshing height (MSH) occurs, changes its orientation instantaneously, according to a rotation about the z-axis by the angle (Fig. 1b), If the liquid behavior is analyzed on the plane at every instant, equation (16) can be extended to the radial coordinate of G, remembering that r n = x 2 n + y 2 n : Equations (20) and (21) can be used to express r G in terms of η n , depending on the adopted model. The approach seen in Sect. 3.1 allows writing the formulas for the nth MSH evaluation, both for the L and NL models, considering a 2-dimensional motion of the container (the leftmost column of Table 2).

Remarks
By looking at the leftmost columns of Tables 1 and 2, one can point out that, for equal values of the generalized coordinates (x n , y n ), the ratio between η n obtained from the L model in equations (22), (28) and η n from the NL model in equations (24), (30) is always 4/ξ 2 1n . If only the 1st mode is considered, this ratio is equal to 4/ξ 2 11 ≈ 1.18 and shows that the assumption of planar free-surface always overestimates the sloshing height compared to the assumption of a nonplanar free-surface. Furthermore, while in equations (22), (24), η n has the same sign of y n , in equations (28), (30) η n is always positive. This means that equations (22), (24) express the trend of the sloshing height on only one side of the container (on the other side, the sloshing height is estimated as the opposite of η n ): in this case, we will simply talk about sloshing height (SH). Conversely, in equations (28), (30), η n indicates the maximum peak that occurs on the container wall (on a plane oriented as described in equation (26)): in this case, we will use the expression maximum sloshing height (MSH).

Extension to 3-dimensional motion
When an additional excitation along the z-axis is taken into account, the formulas for the SH and the MSH estimation are reported in the rightmost columns of Tables 1 and 2. They are seemingly identical to those employed in the case of a planar motion along the xy plane. However, the generalized coordinates x n , y n are obtained by solving the EOMs in (14) for the NL model and in (15) for the L model, instead of equations (11) and (12), respectively.
Summing up, once that the liquid properties are known (in terms of container radius R, static height h, and density ρ), the equivalent-discrete model parameters (see Sect. 2) can be computed, and the 3-dimensional formulation can be applied to study the liquid sloshing for a general motion of the container, namelyS 0 = [ẍ 0ÿ0z0 ] T , by simply employing the EOMs in (15) if the L model is adopted, or by solving the EOMs in (14) if the NL model is chosen. Once that the generalized coordinates (x n , y n ) of the nth sloshing mass are computed, if the container is not subject to an acceleration along the x-axis throughout the whole motion (ẍ 0 = 0 [m/s 2 ]), for the SH estimation, we must refer to Table 1. Conversely, if bothẍ 0 and y 0 are other than zero, the reference for the MSH estimation is Table 2. The aforementioned generalization for the sloshing-height estimation is illustrated by the flowchart in Fig. 3.
Hence, the sloshing-height estimation requires the numerical solution of the EOMs to find the generalized coordinates (x n , y n ) of the nth sloshing mass, together with the computation of the formulas in Tables 1 and 2, depending on the adopted model. Regarding our (nonoptimized) MATLAB implementation, the average computational time is 0.7 s, for either the linear and the nonlinear formulations.

Experimental setup and motion laws
The experimental setup comprises a cylindrical container with radius R = 50 mm and a liquid static height h = 70 mm. The liquid is water, which is colored by adding dark brown powder, in order to obtain a better contrast for the image processing analysis. Motions are performed by an industrial robot (Stäubli RX130L) and recorded by a GoPro Hero3 camera. The trajectories are planned so that the robot follows three geometrical paths on the xy plane (Fig. 4), each of them with different motion profiles, characterized by increasing container accelerations: • a back-and-forth linear path (indicated as l-motion); • an eight-shaped path (e-motion);    In Fig. 5, the trends of the second time derivative 1σ of the path parameter are illustrated: for every path, all three motion profiles are shown. Note that the legend refers to the maximum of the container acceleration norm S 0 max reached during the corresponding motion.
Additionally, the 2D planar motions obtained from a modified trapezoidal motion law with 6 segments ofσ are extended into 3D motions (l3-motion, e3-motion, c3-motion) through the inclusion of an excitation along the z-axis (Fig. 6a). In particular, for the e3motion and c3-motion, the accelerations along the x and y directions are kept unchanged with respect to the corresponding e-motion and c-motion of the 2D case, respectively. The same cannot be said about the l3-motion, where the acceleration along the y-axis was slightly modified to meet the robot limits. As a result of the extension, if 2D and 3D paths are observed from the top of the xy plane, they share the same shape, whereas, from a front perspective, the vertical coordinate changes along the path according to an excursion z (Fig. 6b). On each 2D path, two different trends ofz 0 are considered, thus providing two different 3D paths, characterized by two different values of z: • a profile with moderate dynamics, namely |z 0 | max ∈ [2, 3] m/s 2 ; • a profile with more prominent dynamics, namely |z 0 | max ∈ [4, 5] m/s 2 . Regarding the extraction of the experimental sloshing height from the recorded videos, once each frame is converted into black & white, the image must be processed in two different ways, depending on the excitation type: • SH detection. When the excitation is given by a 1-dimensional acceleration along the y-axis with or without the addition ofz 0 on the vertical direction, the trend of the SH can be examined on only one side of the container (Sect. 3.3); in this case, the SH is detected by identifying the black pixel with the highest z-coordinate on the rightmost side of the container (Figs. 7a, 7b); then, the SH is evaluated as the difference (converted in mm) between the vertical coordinate of the detected pixel and that of the pixel representing the liquid static height h. • MSH detection. When the excitation is a general planar or spatial motion, the MSH can occur anywhere on the container wall. The MSH is again detected as the black pixel with the highest z-coordinate, but in this case the whole surface of the container is considered. Furthermore, a distinction is made between a peak occurring on the front part of the container and a peak on the rear wall, noticing that in the former case the liquid image presents a uniform black shape (Fig. 7c), whereas in the latter case the liquid image is characterized by white regions due to the light reflection on the liquid free-surface (Fig.  7d). Through this distinction, the peak can be correctly located on the container surface and the knowledge of the image depth can be used to obtain a more realistic measure of the MSH.

Experimental results
In Fig. 8, the 2-dimensional L and NL model predictions are compared with the results from the experimental motions, only considering the 1st sloshing mode. A good adherence between the experiments and the models can be appreciated for the 1-dimensional motions (Figs. 8a, 8b, 8c), and tracking remains reliable also for 2-dimensional motions, especially when considering lower values of S 0 max (Figs. 8d, 8g). As the value of the 2-dimensional excitationS 0 is increased, the model predictions still capture the trend of the real liquid MSH, although they seem to lose accuracy in correspondence to the peaks reached by the liquid (Figs. 8e, 8f, 8h, 8i). This can be eventually attributed to two reasons: • The high dynamics given to the container causes a regime of strongly nonlinear motions, where the liquid free-surface loses the assumed shape and shows instantaneous swirly peaks, as illustrated in Fig. 9a, that cannot be tracked by the models; • The height of the frames that are employed for the image processing analysis, grants a greater field of view when the liquid peak occurs on the rear wall of the container (Fig. 9b), whereas, for a peak on the front wall (Fig. 9c), the value of the real MSH is saturated by the frame upper limit; this explains the discrepancy between the experiments and the prediction models in the red areas that are highlighted in Figs. 8e, 8f, 8h, 8i.
Regarding the former case, the L-model assumption of a planar surface during motion loses adherence with reality when the container acceleration is roughly greater than S 0 max ≈ 8 m/s 2 , with a maximum jerk of ... S0 max ≈ 60 m/s 3 . This is reflected in a less accurate correspondence between the model prediction and the experimental result, even though the model evaluation is still reliable. It is worth observing that the maximum peaks are always overestimated by the L model. As far as the NL model is concerned, the assumption of a free surface described by means of the Bessel function of the first kind finds a better correspondence with reality, compared with the L model. The model estimation begins to lose adherence with reality when the container acceleration reaches a value of roughly S 0 max ≈ 9.5 m/s 2 , with a maximum jerk of ... S0 max ≈ 74 m/s 3 . Additionally, the global maxima predicted by the NL model are always below the experimental ones. Table 3 summarizes the obtained results by reporting the accuracy index mod expressing the error between the model and the experimental maxima: where the subscripts mod and exp denote the adopted model (2-dimensional L/NL) and the experimental results, respectively. For all motions, | 2D−L | is always below 18%, and | 2D−NL | never exceeds 19%, with the NL model granting a better tracking during the whole time period. Furthermore, the positive sign of 2D−L proves that the L model always overestimates the real SH and MSH peaks, as expected (see Sect. 3.3), hence providing a more conservative estimation.  In Fig. 10, the results from the 3D motions are illustrated: for each motion, the 3dimensional L and NL model predictions are compared with the 2-dimensional L and NL ones, to show the benefit obtained by employing the extended formulation. In general, the 3-dimensional models exhibit a better correspondence with respect to the 2dimensional ones, especially when the vertical acceleration is more significant (e.g., when |z 0 | max ∈ [4, 5] m/s 2 ). The definition of the index σ mod expresses the mean absolute error between the experimental results and the model predictions, where the subscripts mod and exp denote the adopted model (2-dimensional L/NL or 3dimensional L/NL) and the experimental results, respectively, whereas t i refers to the ith time instant and N represents the number of samplings. The analysis of Table 4, in which the index σ mod is evaluated for the 3-dimensional motions, confirms the qualitative evidence, i.e.: • for equal formulations (3D or 2D), the NL model grants a better tracking of the realliquid trends with respect to the L model (for instance, in most cases, σ 3D−NL < σ 3D−L and σ 2D−NL < σ 2D−L ); • in general, the 3D formulation provides a more reliable correspondence with reality, compared to the 2D formulation; indeed, in most cases σ 3D−NL < σ 2D−NL and σ 3D−L < σ 2D−L .

Remarks
In the described experiments, the employed amount of liquid inside the container (R = 50 mm, h = 70 mm, with h/R = 7/5) has a mass of m F ≈ 0.55 kg. Considering only the first sloshing mode to model the liquid dynamics, from equation (4), we can theoretically compute the mass m 1 which is responsible of the liquid sloshing (m 1 ≈ 0.18 kg). If the case with the highest dynamics is taken into account (e.g., the c-motion with S 0 max ≈ 9.5 m/s 2 ), it can be shown, from the simulations, that the forces transmitted by the liquid on the robot end-effector are negligible, hence not influencing the robot dynamics. Conversely, as the value of m F is increased (and consequently of m 1 , if h/R is fixed), the forces acting on the robot end-effector get more important.   It is also worth mentioning the influence of the ratio h/R. In particular, for a fixed value of m F , when h/R increases, the value of the sloshing mass m 1 drops down [10], whereas m 1 represents a more important contribution as h/R decreases. As long as the sloshing mass m 1 grows, the sloshing of the liquid becomes more significant, and the free surface might not satisfy the model assumptions. Consequently, for lower values of h/R, the prediction models may lose adherence with reality. This also finds confirmation by looking at the correlation between the first natural frequency and the ratio h/R (see equation (3)): with the decreasing of h/R, ω 1 tends to lower values, meaning that the fundamentals of the liquid free-surface can be easily excited by the motion-law spectrum, in correspondence of the highest amplitudes [18]. This may result in a more nonlinear behavior and shape of the liquid free-surface that might not be detected by the models, although a higher number of sloshing modes is considered.

Conclusions and future work
A novel technique for the sloshing-height estimation of a liquid inside a cylindrical container subject to 2-dimensional planar motions was proposed, extending what was presented in [17]. The technique is based on simple discrete mechanical models, rather than machinelearning or complex fluid dynamics methodologies, thus granting a reliable and easy-tocompute estimation.
Experiments, considering three container planar paths performed by an industrial robot with different motion profiles, were presented and the relative results were discussed. An accuracy index, expressing the error between the model and the experimental sloshing-height maxima, was used to prove the effectiveness of the estimation, even for high values of the container acceleration (up to 9.5 m/s 2 ).
Additionally, the 2-dimensional formulation was extended into a 3-dimensional one, to take into account an additional excitation along the vertical direction, hence resulting in a 3D motion of the container. The planar paths of the 2D case were extended by adding a z-axis acceleration (up to 5 m/s 2 ). The mean absolute error between the experimental results and the model predictions was examined to evaluate the benefit, in terms of accuracy, obtained by employing the 3D formulation instead of the 2D one.
Future work will include a further validation study considering higher values of the vertical acceleration and the use of the proposed sloshing-height estimation for square-section containers, adapting the formulation that was presented in Sect. 3. In addition, a detailed sensitivity analysis will be addressed to assess the influence of the container dimensions on the accuracy of the model predictions (Sect. 4.3).