Multibody Analysis of Sloshing Effect in a Glass Cylinder Container for Visual Inspection Activities

: This paper addresses the phenomenon of sloshing and the issues that arise during liquid handling at visual inspection stations. The pharmaceutical industry, recently put under pressure by the pandemic, has long adopted modular solutions consisting mainly of robotic islands. This work focuses on a visual inspection island for glass vials and ampules called VRU. This machine uses robotic arms to optimize the inspection process and enables automated control of a wide range of products using image recognition techniques and AI algorithms. However, the handling of containers in the presence of liquids requires special precautions to avoid the occurrence of bubbles inside the fluid that can prevent the cameras from correctly capturing any defects present. The banal solution involves a drastic reduction in the speeds and accelerations to which the liquids are subjected. However, using appropriate techniques makes it possible to achieve performance values similar to those obtainable when manipulating solid materials. The developed algorithms were tested using multibody simulations in the Mathworks Simscape environment and then validated using a six-axis Fanuc robot. In this study, however, the analysis conducted aimed to determine the correlations between trajectories, laws of motion, and sloshing in containers handled at high speed in industrial applications. In this study a multibody model was developed using a CFD analysis. The container consisted of a glass vial for pharmaceutical uses containing a liquid inside. The results obtained from the CFD analysis allowed us to calibrate the multibody model for the next phase of optimization of the laws of motion to be followed by the manipulator.


Introduction
The first studies on sloshing began in the 1960s when NASA began to use liquid fuels in their rockets.The inertial effects of a fluid and the gradual emptying of the tanks, with the consequent increase in the wave motion, created severe problems of stability and control of the trajectories during the launches of the rockets [1,2].Given the low capacity of computing systems, simplified models were developed, derived mainly from experimental activities, able to capture the phenomenon for specific fields of application [3].Today, however, most studies present in the literature consist of sloshing applications for vibration mitigation activities in structures [4][5][6].In this work, the authors defined a simplified multibody model developed through a CFD analysis to determine the optimal operating conditions for serial manipulators used in visual control stations for glass containers [7][8][9][10][11].The SimScape multibody multidomain simulation environment is increasingly used in complex systems analysis due to its ability to model the different subdomains that characterize real systems in a single environment [12][13][14].In this way, the estimation of the dynamic behavior of new systems becomes straightforward, allowing the dimensioning of mechanical, electrical, hydraulic, or other subsystems [15,16].This work, which represents the first phase of an extensive research activity, is divided as follows: Section 2 describes the mathematical model underlying the simplified system and the corresponding multibody model developed in Mathworks' SimScape multidomain simulation environment and developed using a CFD analysis; Section 3 describes the numerical activity undertaken with the estimation of equivalent parameters for the simplified system and their validation.In Section 4, a discussion of the results is presented while in Section 5, the conclusions are finally reported.

Materials and Methods
The analysis of the dynamic behavior of a fluid inside a container necessarily passes through the resolution of the Navier-Stokes equations [17].
The accurate estimation of the vibration characteristics of the sloshing motion is crucial for analyzing the behavior of tanks or containers experiencing base motion.The virtual mass force, drag and stagnation pressure force, lift force, and others should all be evaluated to analyze the behavior of a fluid in a container [18].Such modeling can be simplified when the oscillations are confined and non-linear sloshing phenomena are avoided [19].Under these conditions, it is possible to use the equivalent mechanical model shown in Figure 1a, which takes into account the motion of the fluid directly proportional to the acceleration of the rigid container, described by a mass m 0 integral to the container, and the superposition of j different modes of oscillation described using a series of simple pendulums in which the remaining mass of fluid is concentrated, linearized to the vertical position [20].In (1), from the linear wave theory obtained by the linearization of the Navier-Stokes equations, the equation of the natural frequency ω n of the j-th oscillation mode has been derived: The dependence of the damping coefficient σ j of the j-th mode of oscillation on the internal forces of the system means that σ j can be calculated by observing how the wave motion decreases according to a logarithmic decay due to the dissipation of energy.In addition, only the first asymmetric mode was considered because it is dominant over higher modes.The damping factor from Ibrahim's theory was obtained by an empirical law reported in (2) that considered the kinematic viscosity ν, the depth of the liquid h, and the radius of the container R.
From Ibrahim's theory, it was possible to derive the relationships between the different physical parameters, where it was found that the values of L j were reduced by approaching the points where they are bound to make the j-th in pendulum oscillate around the free surface.Furthermore, there is an inverse proportionality between the indices j for the lengths L j and masses m j .Such a result allowed us to simplify the problem further, as shown in Figure 1b, and model the system as a single pendulum of mass m 1 and length l 1 , bounded to oscillate around a point at the center of the free and undisturbed surface of the liquid [21,22].The plane of the free surface was always orthogonal to the pendulum and planar during the motion.Further simplifying hypotheses were as follows: • The incompressibility of the liquid; With this model, only the fundamental oscillation mode (mode 1) is adequately described, which, among the various j-th modes, is the one with the lowest frequency and characterized by the most significant amplitude (see Figure 2a).In the literature, it has been demonstrated experimentally that the influence of the other higher frequency modes is not particularly significant, and therefore, their inclusion in the model would lead to excessive complications in the control system [23,24].What has been described is valid exclusively for symmetrical containers concerning a single axis.Otherwise, for each main axis of the section of the container, there is a fundamental mode [25,26].Based on the hypotheses made, a one-dimensional model was developed consisting of a trolley subject to motion that translated along the horizontal direction, visible in Figure 2b, in which a linear damper was connected to the pendulum to add to the model the dissipative characteristics given by the viscosity and by the friction created between the fluid and the container [27].For the definition of the equations of motion that governed the chosen 2-DOF model, the Lagrange energy method, reported in (3), was used.
where q indicates the following vector of Lagrangian coordinates x θ , where x the onedimensional motion of the container is indicated, and where θ, the oscillatory motion of the equivalent pendulum which represents the motion of the fluid, is indicated.Furthermore, E c is the total kinetic energy of the system, D is the dissipative function equal to half of the power dissipated by the damper during its operation, V instead describes the potential energy of the system.Finally, Q represents the active component of forces and torques applied to the system.
The two equations of motion are presented in a non-linear form, but for small oscillations, it is possible to linearize them [28].In the following expression the linearized equations are reported in matrix form: The container's acceleration becomes the system's control variable and is consequently the variable we focused on in this study.Specifically, by rearranging the terms known to the second member, and dividing the expression by the inertia of the pendulum, the following relation was obtained: Therefore, it is clear from (7) that to calibrate the model, it is necessary to define the equivalent mass of the system m, the equivalent length of the pendulum L, and the equivalent damping σ of the liquid.
The mass of the equivalent pendulum corresponding to the sloshing mass of the fluid m n can be estimated by using Equation ( 8), while the equivalent length of the pendulum can be estimated by using relation ( 9) For both relationships, the integer n indicates the mass and equivalent length of the pendulum representative of the nth generic mode of vibration.
Experimental testing or a CFD analysis may be used for estimating equivalent damping [29].In this work, it was decided to test the numerical approach.A CFD analysis makes it possible to find a solution for discretized Navier-Stokes partial differential equations because they do not have an analytical solution.These equations describe processes of momentum, energy transfer, and mass by using a strong discretization of the geometry to be analyzed (see Figure 3a).In sloshing problems, a further complication involves tracing the free surface and defining the boundary of the different phase fluids.The Volume of Fluid (VOF) model tracks the free surface's speed, position, and shape between the different phases of the fluids [30].The finite volume method solves the Navier-Stokes equations by storing the values of all fluid properties, such as velocity, pressure, density, and volume fraction, at the center of each control volume.The VOF model then extracts the volume fraction data in each control volume to determine the shape and position of the free surface.For example, in the case of a multiphase system consisting of air and water, if the water volume fraction is 1 with respect to the control volume, it means that the control volume is completely filled with water.The solution of the volume fraction conservation equation defined by Hirt and Nichols is used to trace the free surface throughout the volume [31]: The F function represents the volume fraction corresponding to each control volume.Such function is also used to describe the default domains.As an example, near the container mouth, where only air is present, the air volume fraction is 1, while the water volume fraction is 0.More than the information given by F, there is still the need to solve the problem, because it does not define the tilting of the free surface.To do this, the VOF uses the gradient of the volume fraction at each control volume across the free surface to determine the slope of the free surface, as shown in Figure 3b.The materials' properties and the parameters used in the CFD analysis are reported in Tables 1 and 2, respectively.

Numerical Activity
The numerical activity was divided into three linked phases.A flowchart of the study has been created to help the reader follow the process.Such a scheme is reported in Figure 4, in which the inputs of the different processes are indicated in red, while in white, the CFD simulation and the simulation of the multibody model with the respective results are reported in the lower part.This activity's primary goal was to evaluate the damping ratio, one of the main characteristic factors that characterize the motion of the fluid.A CFD VOF analysis was performed for extrapolating the value of ξ through the gradient of the long volume fraction z, the liquid height variation considered the system response.
As for the mesh, the model was designed with a CAD modeler.and developed within the program Ansys Meshing-2023 R1 Free Student Software, whose reference system was positioned on the edge of the bottle at the vertical axis of symmetry in y, in the direction of motion in z.Two mesh domains were created: one corresponding to the upper part of the vial and the other corresponding to the container walls.The same boundary conditions were used by considering a harmonic law for the displacement with a frequency of 1 Hz along the z direction, as reported in Figure 5a.Starting from the variation in the height reached by the liquid as a function of time, it was possible to calculate the logarithmic decrement by implementing the Matlab 2023b function findpeaks to find the peaks of the function.This value was estimated at 0.163.
Equation (11) shows the expression used to evaluate the logarithmic decrement of the discrete signal, while Equation ( 12) describes the relation that expresses the damping ratio ξ 1 in the case of an underdamped vibrating system: The graph of the logarithmic decrement is shown in Figure 5b.By estimating the damping factor, we were able to begin the study's second phase, the optimization of the law of motion to be imposed on the slide.To summarize the first phase, all the values used for modeling the simplified system realized in Mathworks' SimScape multidomain environment have been reported in Table 3.  Figure 6 shows the rigid multibody model developed in the SimScape environment of a slide mounted on guides and actuated by a timing belt.A vial rack was placed on the guide.In the image, it is possible to see the vial in which the pendulum is located placed on the rack, whose properties have been shown previously.Given all the geometric and physical parameters of the model, we focused our attention on the laws of motion and how these interacted with sloshing.
A trapezoidal, cubic polynomial, V-degree motion profile, and finally, a polynomial that optimized the time derivative of acceleration (jerk) were generated as input data for a defined linear trajectory of 0.5 m.
The system's response was analyzed by comparing the different amplitudes of oscillation of the pendulum, determining which law of motion could minimize sloshing.The most common law of motion is the trapezoidal one, typical of mechanisms driven by electric motors, considering a maximum speed equal to 0.5 m/s, setting starting and ending points as waypoints.Third-and fifth-degree polynomial laws and optimized jerk polynomial laws were also used to test the pendulum's behavior.In Figure 7, the velocity and jerk laws for the four laws of motion tested are reported.From the literature, the influence of the discontinuities in the accelerations on the onset of sloshing in manipulating fluids is evident.From the graph reported in Figure 7b, it is evident that unlike the jerk-optimized law of motion, the other laws of motion always present a discontinuity.To verify this deduction, it was decided to test the behavior of the fluid by subjecting the cart to the optimized law and the fifth-degree polynomial law of motion, where the initial and final conditions of the acceleration were controllable.

Results
The first test carried out between the simplified model and the CFD analysis carried out in Ansys Fluent concerned the comparison between the oscillation of the pendulum θ MB (t) and the oscillating motion of the plane corresponding to the free surface of the liquid θ CFD (t) due to the motion imposed on the cart.Analyzing the free underdamped responses for both models allowed us to estimate the natural oscillation frequencies by studying the frequency response of the signals.
For evaluating the free response θ CFD (t) from the CFD analysis, the expression (13) in CFD-Post was implemented.areaInt(Water.VolumeFraction.GradientZ)@Isosur f ace (13) where Isosurface corresponds to the free surface of the liquid.Subsequently, the results were imported in CSV format into Matlab for the calculation of θ.
To calculate θCFD , instead, from CFD-Post, the values of the surface velocity as a function of time were imported, and from this, the angular velocity was compared with that of the pendulum and derived.For the multibody model, the data, on the other hand, were measured through the revolute joint that connects the pendulum rod with the rotational center appropriately calculated through the relationships described in the previous paragraph.The comparison graphs are shown in Figure 8. Table 4 reports the estimated parameters for the two case studies: the natural frequencies and the damping ratio.From the table, it is possible to notice that the natural frequencies of the system evaluated for the two case studies were encouraging.This is due to the fact that the natural frequency of the pendulum was calculated with a semi-empirical law (see Equation ( 1)).Finally, Figure 9 shows the last activity of this study, where the behavior of the simplified model, calibrated on the results of the CFD analysis, was tested with the four laws of motion considered in the previous paragraphs.As can be inferred from the study, the strong irregularities recorded for the trapezoidal jerk and the third-degree polynomial function translated into significant oscillations of the pendulum.Instead, in the case of the fifth-degree polynomial law and in particular, for the optimized jerk law of motion, the oscillations were highly reduced.

Conclusions
This work aimed to develop a simplified model with concentrated parameters to test algorithms of optimal trajectory and motion law planning.Such a model will minimize the sloshing effect of fluids within pharmaceutical glass vials.The ultimate goal is the optimal programming of anthropomorphic manipulators in visual control stations to identify defects in high-speed processes.In recent years, there has been a high demand for pharmaceutical products worldwide, especially vaccines, to increase productivity.Attention has been focused on this phenomenon because, at high speeds, detecting defects in the presence of bubbles due to sloshing by using cameras and sensors in moving fluids is problematic.A CFD analysis was used to properly calibrate the simplified model, from which it was possible to calculate the equivalent damping coefficient, a fundamental parameter for developing the multibody model.In the last part of the study, the correlation between the points of discontinuity of the acceleration, resulting in the presence of strong jerk irregularities, and the oscillation of the liquid implied an increase in the sloshing effect.The results obtained and in particular, the simplified model developed will allow testing of new three-dimensional trajectories minimizing the wave motion of the fluid, generally increasing the speed of the manipulators and consequently, the performance of the visual control station.
radius; • h s = static height of the liquid; • ϵ j = j-th root of the derivative of the Bessel function of the first kind (for the first fundamental mode, it is equal to 1.841).(a) (b) Figure 1.Free-surface motion equivalent mechanical models: (a) multi-DOF pendulum-based model; (b) 1-DOF pendulum-based model.

Figure 2 .
First three modes of the free surface of the fluid (a) and lumped-parameter model used in the numerical activity (b).

Figure 3 .
Glass vial mesh in Ansys Fluent (a) and water volume fraction during oscillations in a simulation (b).

Figure 5 .
Figure 5. Motion law imposed on the trolley (a) and vertical position of the oscillating mass with the peak values, reported with a red circle and black dot, used for estimating the damping factor (b).

Figure 6 .
Figure 6.Multibody model of the test-rig in the Matlab's SimScape environment environment.

Figure 7 .
Figure 7. Cart's linear speed for the different laws of motion considered (a) and related jerk signals (b).

Figure 8 .
Figure 8. Simplified system response for the different laws of motion considered.

Figure 9 .
Figure 9.Comparison between free response of simplified model and CFD analysis: angular position (a) and angular velocity (b).

Table 3 .
Geometric and physical parameters of the equivalent pendulum.

Table 4 .
Comparison of characteristic parameters evaluated in the CFD environment and in the multibody SimScape environment.