Modeling Transient Liquid Slosh Behavior at Variable Operating Speeds Induced by Intermittent Motions in Packaging Machines

Abstract

The present paper deals with the problem of modeling liquid slosh occurring in the packaging process of containers filled with liquid. Sloshing effects are induced by one-dimensional intermittent motions and are undesired due to the necessity of quality control processes, such as weighing. Therefore, motion optimizations are often applied with the intention to minimize the residual vibrations. Valid process models are required to do so. The aim of this paper is to derive models for describing the liquid slosh behavior for different motions and for common practical circumstances, e.g., different container geometries as well as machine operating speeds, and to state the model’s limits of use. Known model approaches are discussed, and their assumptions are reviewed experimentally. This leads to a set of limited ranges of operating speeds in which the applied models’ assumptions are valid. The models are derived for these sets from experimental data, and a comparison is executed that enables the determination of the models’ validity concerning their operating speed dependency. Finally, the validity of the derived models is investigated by comparing their predictive efficiency of describing the vibration for different motion profiles.

1. Introduction

The investigation of liquid slosh dates back to the 1960s and was initially motivated by aerospace research [1,2]. Later on, more practical use cases of this phenomena were investigated [3]. One practical application of the underlying physical problem occurs in packaging machines producing beverages, displayed in Figure 1. Due to high requirements in the production process of daily needs such as food and beverages, packaging machines are intended to produce at maximum productivity. Often, neither the mechanical nor the mechatronical components limit the machine’s operating speed and hence the productivity. In fact, the behavior of the materials to be processed limits the maximum achievable operating speed [4]. For instance, residual vibration induced by the applied intermittent motion profile is to be avoided, e.g., due to necessary quality control processes. Hence, the motion optimization of the conveying motion is a powerful method for minimizing the residual vibration. The optimization of the motion profile is a suitable possibility for the given problem, but it is less widespread in industrial practice [5]. Nevertheless, a valid model that predicts the sloshing behavior depending on a motion profile is needed for motion optimization.

Model approaches can principally be subdivided into two main classes, namely descriptive models and explanatory models [6]. In the case of modeling liquid slosh, explanatory approaches use the Navier–Stokes equations for describing the liquid–container interaction as a continuum problem. Therefore, material parameters, such as viscosity and density, are considered as distributed parameters. Such a general approach leads to a system of partial differential equations [7,8,9]. This system of equations is not solvable analytically, wherefore numerical solution techniques such as the finite difference method or the finite volume method have to be applied [10,11,12,13]. Models with distributed parameters can be used for motion optimization as shown in [14,15]. Such solutions are associated with enormous calculation times, although the presented solutions already represent an approximation of the partial differential equation system. Another possibility is the modeling with concentrated parameters and physical substitution models such as the pendulum model, also discussed in [7]. This approach leads to an ordinary differential equation of order two, which can be easily considered in common motion optimization algorithms in an efficient way [16,17].

A new aspect of research in the field of model-based motion optimization is the operating-speed-dependent motion profiles that are discussed in [18,19]. Applying this approach to the motion optimization for liquid slosh, the previous modeling work has to be re-investigated and discussed regarding their operating-speed-dependent usability. Hence, this paper aims to answer two main questions regarding the modeling of liquid slosh with the purpose for using it in operating-speed-dependent motion optimizations:

• Are there limitations of the substitution mechanical slosh models concerning the applicable operating speed?
• If such limitations can be stated, are the derived models applicable for the whole derived range of operating speeds and at what accuracy?

For taking into account the special case of liquid slosh occurring in packaging machines, the following limitations are assumed: only one-dimensional motions and transient excitations are considered due to considering only one machine cycle—in contrast to investigations concerning periodic behavior [20]. Thus, there will be no discussion of resonances and magnification functions, since steady-state behavior must be therefore present.

2. Materials and Methods

2.1. Test Station and Experimental Design for Investigating Liquid Slosh

To investigate the liquid slosh behavior, a test station shown in Figure 2 was build. The station consists of a direct linear drive on which different containers can be conveyed with any desired motion profile. To observe the fluid’s behavior, a high-speed camera is coupled to the mechanical setup, whereas the camera is triggered by the motion control system. Hence, the captured frames can be synchronized with the applied motion profile, which simplifies the evaluation. In addition to that, the applied frame rate (FR) is automatically determined by the observed operating speed n. Hence, all experiments supply the same amount of frames for further investigation, independently from the applied operating speed. The specifications of the used hardware are displayed in

Table 1. Specifications of used hardware elements in the test station.

Element	Specification	Remark
High speed camera	Mikrotron EoSens MC3010	Triggered by GardaSoft CC3320
Direct linear drive	Baumüller LSM 60-0916	–
LED light	Bestlight VPad-150	–
Motion control system	Beckhoff TwinCAT 3 Moog MSD	(Programmable logic controller) (Converter)

.

The captured high speed frames are evaluated automatically by detecting the region of interest, binarizing the frame, and finally detecting the fluid’s surface with the Canny algorithm [21] (see Figure 3).

To cover a wide range of practical conditions, the experiments are carried out with six container geometries at two different motion profiles and 17 different operating speeds. The geometries differ in their width w. The exact specifications are depicted in

Table 2. Specification of experimental design for creating parameter identification and validation data.

Variable Experiment Element	Specification	Remarks
Container geometries	${\mathcal{G}}_1$: $w=60$ mm	For all geometries:
	${\mathcal{G}}_2$: $w=40$ mm	Fluid: water
	${\mathcal{G}}_3$: $w=70$ mm	Filling level: 67 mm
	${\mathcal{G}}_4$: $w=47$ mm
	${\mathcal{G}}_5$: $w=79$ mm
	${\mathcal{G}}_6$: $w=54$ mm
Motion profiles	Acceleration pulse motion	Used for parameter identification ($K_1=10$)
	Quintic polynomial motion	Used for validation ($K_2=5$)
Set of operating speeds	$n_1=12$ 1/min
	$n_2=15$ 1/min
	⋮
	$n_{17}=60$ 1/min

. The stroke, velocity, and acceleration trajectories of the used motion profiles are displayed in Figure 4. Due to the high water content in everyday beverages, water is taken as the fluid to be examined.

For creating the parameter identification data, all six geometries are conveyed with the acceleration pulse motion at 17 different operating speeds, whereas every experiment is repeated 10 times, represented by the constant $K_1$. Hence, an amount of 1020 experiments is to be carried out and evaluated. The validation data is generated by applying the quintic polynomial motion profile, whereas each experiment is repeated five times ($K_2=5$). Therefore, an amount of 510 experiments is to be conducted and evaluated for creating the validation data. The idea behind this strategy is that different system inputs are used for parameter identification and validation.

2.2. Modeling Liquid Slosh

In order to be able to use the model to be created for motion optimization, as few degrees of freedom as possible should be used in order to avoid unnecessary computing time. Therefore, the already mentioned mechanical pendulum model is investigated in this paper. The real system, shown in Figure 5, can be abstracted to a model consisting of a horizontally movable mass $m_{\mathrm{M}}$ displaying the linear drive and an attached damped pendulum representing the fluid’s surface displacement with length l, mass $m_{\mathrm{P}}$, and damping constant d. For describing the system’s dynamic behavior, two generalized coordinates, namely $x_{\mathrm{M}}\left(t\right)$ as well as $\theta \left(t\right)$, are introduced, which are summarized in the vector $\mathit{q}={(x_{\mathrm{M}},\theta )}^{\top }$. The system is controlled by the time-depending force $F\left(t\right)$ affecting the mass $m_{\mathrm{M}}$. In contrast to [16], the model equations are motivated with respect to the applied drive.

According to [22], the Lagrange formalism is used for determining the equations of motion for the abstracted system. The kinetic energy T, the potential energy U of the system, and the dissipation energy D due to the damper element affecting on the pendulum are given by

$$\begin{array}{cc}\hfill T& =\frac{1}{2}(m_{\mathrm{M}}+m_{\mathrm{P}})·{\dot{x}}_{\mathrm{M}}^2+\frac{1}{2}{m}_{\mathrm{P}}{l}^2·{\dot{\theta }}^2-m_{\mathrm{P}}lcos\left(\theta \right)·\dot{\theta }·\dot{x},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill U& =-m_{\mathrm{P}}gl·cos\left(\theta \right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill D& =\frac{1}{2}d{l}^2·{cos}^2\left(\theta \right)·{\dot{\theta }}^2.\hfill \end{array}$$

(3)

Subsequently, the equations of motion can be derived with the actuating force $\mathit{Q}={[F\left(t\right),0]}^{\top }$ from evaluating

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial (T-U)}{\partial {\dot{q}}_o}\right)-\frac{\partial (T-U)}{\partial q_o}+\frac{\partial D}{\partial {\dot{q}}_o}=Q_o.$$

(4)

This leads to a system of two nonlinear differential equations:

$$\begin{array}{cc}\hfill \ddot{x}& =\frac{F+m_{\mathrm{P}}l\left(cos\left(\theta \right)·\ddot{\theta }-sin\left(\theta \right)·{\dot{\theta }}^2\right)}{m_{\mathrm{P}}+m_{\mathrm{M}}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \ddot{\theta }& =-\frac{d{cos}^2\left(\theta \right)}{m_{\mathrm{P}}}\dot{\theta }-\frac{g}{l}sin\left(\theta \right)+\frac{1}{l}cos\left(\theta \right)\ddot{x}.\hfill \end{array}$$

(6)

Considering the real system, the actuating force $F\left(t\right)$ is caused due to a direct linear drive, which is driven by a cascade control loop with a very high dynamic stiffness. Furthermore, the drive’s mass is essentially greater than the pendulum’s mass. These two statements lead to the simplification $m_{\mathrm{P}}/m_{\mathrm{M}}\approx 0$, from which

$$F\left(t\right)=m_{\mathrm{M}}·\ddot{x}\left(t\right)$$

(7)

results and, hence, the dynamical behavior of the system can be solely described with Equation (6).

Furthermore, the drive’s acceleration $\ddot{x}\left(t\right)$ is redefined as the motion’s profile acceleration a. As a last step, the equation is linearized by assuming $cos\left(\theta \right)\approx 1$ and $sin\left(\theta \right)\approx \theta$, which leads to

$$\ddot{\theta }\left(t\right)=-\frac{d}{m_{\mathrm{P}}}·\dot{\theta }\left(t\right)-\frac{g}{l}·\theta \left(t\right)+\frac{1}{l}·a\left(t\right).$$

(8)

Introducing the undamped angular frequency ${\omega }_0^2=g/l$ and the damping ratio $\zeta =d^{2}l/\left(4m_{\mathrm{P}}^{2}g\right)$, Equation (8) finally results in the model equation ${\mathcal{M}}_1$:

$${\mathcal{M}}_1(\zeta ,{\omega }_0):\ddot{\theta }\left(t\right)=-2\zeta {\omega }_0·\dot{\theta }\left(t\right)-{\omega }_0^2·\theta \left(t\right)+\frac{{\omega }_0^2}{g}·a\left(t\right).$$

(9)

This equation describes a driven harmonic oscillator and is an ordinary differential equation of order two with two free parameters, which can be used for modeling the real system’s behavior by parameter optimization. The special feature of this equation is that the factor of the external motion profile’s acceleration results from the undamped natural frequency and the acceleration due to gravity.

A generalization of this model approach is the driven harmonic oscillator of order two with three free parameters:

$${\mathcal{M}}_2(\zeta ,{\omega }_0,p):\ddot{\theta }\left(t\right)=-2\zeta {\omega }_0·\dot{\theta }\left(t\right)-{\omega }_0^2·\theta \left(t\right)+p·a\left(t\right).$$

(10)

In [23], it is asserted that the additional parameter p leads to a higher accuracy in depicting the real system’s behavior. This assumption is also claimed in [24] and will be experimentally investigated here. The model equations ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ can be solved numerically for a given set of parameters and a system input $a\left(t\right)$ with the method of Runge–Kutta [25].

2.3. Evaluating the Model’s Assumptions

The main assumption of the derivated model above is the approximation of the fluid’s surface as a straight line. In a first step, this assumption has to be investigated, because in reality the fluid surface has a characteristic waveform. For doing so, for each geometry ${\mathcal{G}}_v$, operating speed $n_u$, repetition k, and captured frame j corresponding to a discrete time value ${\widehat{t}}_j$, an approximation in the form of a line of the measured fluid’s surface has to be derived. Therefore, the measured fluid’s surface ${\eta }_{{\mathcal{G}}_v,n_u}(t,x)$ is approximated as a linear function ${\tilde{\eta }}_{{\mathcal{G}}_v,n_u}(t,x)={\theta }_{{\mathcal{G}}_v,n_u}\left(t\right)·x$ with time-dependent slope ${\theta }_{{\mathcal{G}}_v,n_u}\left(t\right)$. This slope has to be determined as part of a parameter identification. Thus, the residuals between the measured surface edge in the discrete form and the approximated line in the discrete form surface are quadratically summed:

$$\begin{array}{c}\hfill r\left({\theta }_{{\mathcal{G}}_v,n_u,k}\left({\widehat{t}}_j\right)\right):=\sum _{i=1}^I{\left({\eta }_{{\mathcal{G}}_v,n_u,k}({\widehat{t}}_j,{\widehat{x}}_i)-{\tilde{\eta }}_{{\mathcal{G}}_v,n_u,k}({\widehat{t}}_j,{\widehat{x}}_i)\right)}^2=\sum _{i=1}^I{\left({\eta }_{{\mathcal{G}}_v,n_u,k}({\widehat{t}}_j,{\widehat{x}}_i)-{\theta }_{{\mathcal{G}}_v,n_u,k}\left({\widehat{t}}_j\right)·{\widehat{x}}_i\right)}^2.\end{array}$$

(11)

The value ${\theta }_{{\mathcal{G}}_v,n_u}\left(t_j\right)$ that minimizes this sum of residuals is the value to be identified:

$$\begin{array}{c}\hfill {\theta }_{{\mathcal{G}}_v,n_u}\left(t_j\right)=argminr\left({\theta }_{{\mathcal{G}}_v,n_u}\left(t_j\right)\right).\end{array}$$

(12)

The gradient-based algorithm of Levenberg–Marquardt [26] is applied for performing the optimization, due to the convexity of the problem. Hence, short calculation times can be achieved, which are not realizable with global optimization strategies. In Figure 6, an exemplarily sketch of a surface with an approximated line is shown, whereas the difference $\Delta \eta$ is displayed over the horizontal coordinate x.

The remaining residuals are a measure for the approximation’s goodness and are therefore used to score the goodness of line approximation over all discrete time points. Due to the effect of in-phase interception discussed in [19], cancellation of residual vibration may occur at some operating speeds of the motion profile. To enable comparability, the dwell phase is not considered evaluating the goodness of line approximation, which is why only the captured frames of the motion phase $J_b$ are considered:

$${\overline{\epsilon }}_{{\mathcal{G}}_v,n_u,k}:=\sqrt{\frac{1}{J_b·I}\sum _{j=1}^{J_b}\sum _{i=1}^I{\left({\eta }_{{\mathcal{G}}_v,n_u,k}\left({\widehat{t}}_j,{\widehat{x}}_i\right)-{\tilde{\eta }}_{{\mathcal{G}}_v,n_u,k}\left({\widehat{t}}_j,{\widehat{x}}_i\right)\right)}^2}.$$

(13)

Due to the repetitions of the experiments, the mean values $\overline{\overline{\epsilon }}$ as well as the corresponding standard deviations $\Delta \overline{\epsilon }$ are examined to compare the obtained results [27]:

$$\begin{array}{cc}\hfill {\overline{\overline{\epsilon }}}_{{\mathcal{G}}_v,n_u}& :=\frac{1}{K_1}\sum _{k=1}^{K_1}{\overline{\epsilon }}_{{\mathcal{G}}_v,n_u,k}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \Delta {\overline{\epsilon }}_{{\mathcal{G}}_v,n_u}& :=\sqrt{\frac{1}{K_1-1}\sum _{k=1}^{K_1}{\left({\overline{\epsilon }}_{{\mathcal{G}}_v,n_u,k}-{\overline{\overline{\epsilon }}}_{{\mathcal{G}}_v,n_u}\right)}^2}.\hfill \end{array}$$

(15)

Setting a tolerable value for $\overline{\epsilon }$ of 0.5 mm, this evaluation allows the determination of a maximum operating speed $n_{\max }$ up to which the assumption of a fluid’s surface as a straight line is valid. This enables the answering of the first formulated question in the introduction.

2.4. Determining the Model Parameters and Evaluating the Model’s Goodness

For answering the second formulated question, it is to be examined whether the simulated values ${\theta }_{\mathrm{sim}}$ from the model approaches ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ are reproducing the approximated values of $\theta$ within the ascertained range of operating speeds in a suitable way. Therefore, two approaches are investigated.

First, for each operating speed $n_u$ and geometry ${\mathcal{G}}_v$, ${\mathcal{M}}_{1,{\mathcal{G}}_v,n_u}$ and ${\mathcal{M}}_{2,{\mathcal{G}}_v,n_u}$ models are examined by solving a least-squares problem, which considers all repetitions at once. Thereby, the least-square-optimization problem is solved with the Levenberg–Marquardt algorithm again. For calculating the objective function to be minimized, possible parameter candidates as well as the motion’s acceleration profile $a\left(t\right)$ are inserted into the model equations. The resulting differential equations are solved numerically using the Runge–Kutta method in order to finally calculate the residuals between the measured and simulated values of $\theta$.

Hence, a set of models are build to describe the measured behavior for its valid range of operating speeds. This strategy should be named the Single-Operating-Speed Model (SOSM). Second, for each geometry ${\mathcal{G}}_v$ and all valid operating speeds, the models ${\mathcal{M}}_{1,{\mathcal{G}}_v}$, and ${\mathcal{M}}_{2,{\mathcal{G}}_v}$ are examined in the same way considering all repetitions at once. Hence, for one geometry, only one model for approach ${\mathcal{M}}_1$ and one for approach ${\mathcal{M}}_2$ are derived. This strategy should be named the Multiple-Operating-Speed Model (MOSM).

Finally, the goodness of these different strategies and model approaches has to be examined. Therefore, the standard error of estimate ($SEOE$) is calculated by

$$\begin{array}{c}\hfill SEO{E}_{{\mathcal{G}}_v,n_u}:=\sqrt{\frac{1}{K·J}\sum _{k=1}^K\sum _{j=1}^J{\left(\theta \left({\widehat{t}}_j\right)-{\theta }_{\mathrm{sim}}\left({\widehat{t}}_j\right)\right)}^2}.\end{array}$$

(16)

This value is a measure of how much the simulated value of ${\theta }_{\mathrm{sim}}$ differs on average from all approximated values of $\theta$ that result from the measurements. This allows the comparison of the two model approaches ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ as well as the two strategies SOSM and MOSM in the correspondent valid range of operating speeds. Finally, it must be checked whether the identified models can make a meaningful statement regarding the process-relevant residual sloshing behavior. First of all, an appropriate quality criterion is necessary. One possibility is the effective value of the fluid surface during the dwell phase, whose definition generally has the following form:

$$\begin{array}{c}\hfill {\eta }_{\mathrm{eff},k}:=\sqrt{\frac{1}{(1-b)·T_{\mathrm{cycle}}·w}\underset{b·T_{\mathrm{cycle}}}{\overset{T_{\mathrm{cycle}}}{\int }}\underset{-w/2}{\overset{w/2}{\int }}{\eta }_k^2(t,x)\mathrm{d}x\mathrm{d}t}.\end{array}$$

(17)

Due to the discretization of all treated quantities, the integrals are calculated numerically using the trapezoidal method. The investigation of this compacted size has the advantage that it can be applied to the actually measured edge courses of the fluid as well as to the approximated straight line and the simulated courses resulting from applying the mechanical model. Thus, this parameter is universally applicable and enables an evaluation of the process success, especially with regard to the investigation of the success of optimal trajectories. Furthermore, it is possible to compare the difference between real edge progression, straight line approximation, and model prediction. Therefore, measures for these differences can be considered by the following approach:

$$\Delta {\eta }_{\mathrm{appr}}:=\sqrt{\frac{1}{(1-b)·T_{\mathrm{cycle}}·w·K}\sum _{k=1}^K\underset{b·T_{\mathrm{cycle}}}{\overset{T_{\mathrm{cycle}}}{\int }}\underset{-w/2}{\overset{w/2}{\int }}{\left({\eta }_k(t,x)-{\theta }_{\mathrm{appr},k}\left(t\right)·x\right)}^2\mathrm{d}x\mathrm{d}t}.$$

(18)

In order to describe the mean difference between the measured fluid’s edge and the approximated straight line, and the approach for describing the mean difference between the measured fluid’s edge and the simulated model response is as follows:

$$\Delta {\eta }_{\mathrm{simul}}:=\sqrt{\frac{1}{(1-b)·T_{\mathrm{cycle}}·w·K}\sum _{k=1}^K\underset{b·T_{\mathrm{cycle}}}{\overset{T_{\mathrm{cycle}}}{\int }}\underset{-w/2}{\overset{w/2}{\int }}{\left({\eta }_k(t,x)-{\theta }_{\mathrm{simul},k}\left(t\right)·x\right)}^2\mathrm{d}x\mathrm{d}t}.$$

(19)

3. Results

3.1. Determining the Validity of the Model’s Assumptions

Evaluating Equation (13) for all 1020 parameter identification experiments led to the results displayed in Figure 7. For all six investigated geometries, the mean errors of the fluid’s line approximation rise with the operating speed n. Furthermore, the corresponding standard deviations also rise with the operating speed. Observing the six different graphs, it is noticeable that every single graph can be divided into two different ranges of operating speeds. At first, an almost horizontal course can be noticed, followed by a range of a strongly rising course. The defined limit value for $\overline{\epsilon }$ of 0.5 mm forms the boundary between these two ranges and thus leads to the sought maximum operating speed $n_{\max }$ for every chosen geometry. The maximum operating speeds are stated in Figure 7 for each geometry and differ from one another. It can be noted that the container’s width and the maximum operating speed are inversely proportional to each other.

3.2. Determining the Model’s Parameters and Their Goodness

Knowing the maximum operating speeds for all investigated geometries for which assuming the fluid surface as a straight line is valid, the discussed model approaches are derived, and the related $SEOE$ values are calculated as a measure for the model’s goodness. Therefore, for each geometry and a valid operating speed, an SOSM approach is calculated. This leads to two models, each for one operating speed, with two and three parameters, respectively. Furthermore, for each geometry, an MOSM approach is calculated considering model equations ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$, resulting overall in two models with two and three parameters, respectively.

Figure 8 shows exemplarily the measured and simulated results for an operating speed of 45 1/min and the second geometry, as well as the related $SEOE$ values for the SOSM and MOSM approaches. The simulated courses appear to fit very well to the measured courses, and the differences between the different model approaches are very small.

In order to get a better impression of the $SEOE$ values determined for the different approaches and strategies for different operating speeds and geometries, Figure 9 shows the results for two exemplarily geometries. The results of the other geometries behave analogously and are therefore not shown.

It can be seen that the $SEOE$ values for model approaches 1 and 2 are essentially identical, since the data points lie approximately on the line $y=x$. Furthermore, it can be noted that the MOSM strategy provides slightly higher values than the SOSM strategy. For a comparison of all results,

Table 3. The results of model evaluations for different model approaches, strategies, and considered geometries as well as the amount of necessary parameters.

Geometry	SOSM				MOSM
	${\mathcal{M}}_{\mathit{1}}$		${\mathcal{M}}_{\mathit{2}}$		${\mathcal{M}}_{\mathit{1}}$		${\mathcal{M}}_{\mathit{2}}$
	$\mathbf{max}\left\{\mathit{SEOE}\right\}$	$\#\left\{\mathit{P}\right\}$	$\max \left\{\mathit{SEOE}\right\}$	$\#\left\{\mathit{P}\right\}$	$\max \left\{\mathit{SEOE}\right\}$	$\#\left\{\mathit{P}\right\}$	$\max \left\{\mathit{SEOE}\right\}$	$\#\left\{\mathit{P}\right\}$
1	0.785 deg	26	0.748 deg	39	0.823 deg	2	0.811 deg	3
2	1.118 deg	32	0.991 deg	48	1.265 deg	2	1.240 deg	3
3	0.570 deg	24	0.524 deg	36	0.602 deg	2	0.589 deg	3
4	0.910 deg	30	0.868 deg	45	0.965 deg	2	0.940 deg	3
5	0.444 deg	22	0.409 deg	33	0.471 deg	2	0.459 deg	3
6	0.821 deg	28	0.822 deg	42	0.852 deg	2	0.827 deg	3

summarizes the maximum $SEOE$ values for all approaches and geometries. In addition to that, the amount of necessary model parameters $\#\left\{\mathit{P}\right\}$ is stated. These quantities are identical for all geometries applying the MOSM strategy. Applying the SOSM strategy, the amount of necessary parameters depends on the maximum valid operating speed, which in turn depends on the considered geometry.

Hence, it can be summarized that, expectably, due to a higher amount of parameters, the SOSM strategy leads to better results than the MOSM strategy, and model approach ${\mathcal{M}}_2$ leads to better results than ${\mathcal{M}}_1$. Thereby, better results are understood as smaller $SEOE$ values. However, it has to be stated that the differences are nevertheless very small. Due to the very good ratio of minor model errors and a small number of model parameters, model approach ${\mathcal{M}}_1$ in combination with the MOSM strategy will be treated in the following as the preferred solution and investigated during a model validation more closely.

4. Discussion

The question posed at the beginning of the present paper concerning the maximum possible operating speed up to which the assumption of a straight line for describing the fluid surface is valid can be answered with the help of the presented evaluation method. The determined deviations increase with increasing operating speed, since the waveform of the fluid surface assumes an increasing characteristic shape and thus deviates significantly from a straight line. The fact that the standard deviations of $\overline{\epsilon }$ increase with the operating speed may indicate a higher sensitivity of the system to the initial state at higher operating speeds.

Another conclusion that can be derived from the results concerns the relationship between container width and maximum operating speed. With increasing container width, the maximum possible operating speed decreases. This is due to the fact that as the container width increases, the fluid volume to be conveyed increases as well. Due to the associated greater inertia, the characteristic surface waveform is formed at a lower motion’s profile velocity.

It has been shown that model approaches ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ both achieve good results, with the second approach having a free parameter in the form of a proportionality factor to reinforce the motion’s acceleration profile. Apparently, the factor ${\omega }_0^2/g$ given by the linearized pendulum equation is already a very good specification, which is represented almost identically by the free proportionality factor p. Therefore, no additional value is achieved in the model quality, as proposed in [23,24].

The second question formulated in the beginning concerns the applicability of the different model approaches for the whole range of valid operating speeds. Therefore, the two different strategies SOSM and MOSM were investigated. The SOSM strategy provides slightly better results, but those are hardly significant. From this it can be concluded that the model parameters can be regarded as constant within the valid operating speed range and that a model for the respective range is therefore sufficient for the considered application.

Finally, in the sense of a model validation, it should be discussed whether the models of approach ${\mathcal{M}}_1$ in combination with the MOSM strategy of the different geometries derived from the acceleration pulse motion profile also provide similar simulation results for the validation quintic polynomial motion profile. To examine this, repeat tests are performed with $K_2=5$, and the measurements are compared with the corresponding simulation results. The $SEOE$ value is used again for quantification. In Figure 10, exemplarily simulation results and a comparison to measurements are shown for geometry 2 and operating speed $n=45$ 1/min. In principle, it can be stated that the dynamic behavior of the system is also well reproduced for this modified input. Furthermore, the $SEOE$ value is similar to the $SEOE$ derived from the parameter identification measurements (see Figure 8).

For a better comparison of all considered geometries and a corresponding assessment regarding the model’s applicability,

Table 4. Results of investigations: for each geometry, the derived model parameters as well as the maximum operating speed are stated. Furthermore, the maximum $SEOE$ values are given for both the derived model applied to the identification motion profile and that applied to the validation motion profile.

Geometry	$\mathit{\zeta }$	${\mathit{\omega }}_0$	${\mathit{n}}_{\mathbf{max}}$	Identification: $\mathbf{max}\left\{\mathit{SEOE}\right\}$	Validation: $\mathbf{max}\left\{\mathit{SEOE}\right\}$
1	0.042	22.498 rad/s	48 1/min	0.823 deg	1.037 deg
2	0.050	27.347 rad/s	57 1/min	1.265 deg	1.956 deg
3	0.036	21.001 rad/s	45 1/min	0.602 deg	0.874 deg
4	0.042	25.351 rad/s	54 1/min	0.965 deg	1.718 deg
5	0.029	19.622 rad/s	42 1/min	0.471 deg	0.620 deg
6	0.039	23.729 rad/s	51 1/min	0.852 deg	1.461 deg

shows the results of the maximum $SEOE$ values recorded for the parameter identification motion profile and the validation measurements. The respective values are proportional to each other, whereby the results for the validation evaluation are higher. This is understandable for two reasons. First, these system inputs have not been used for parameter identification. Second, the quintic polynomial motion has higher acceleration values, which is why the system response also assumes higher values and therefore model inaccuracies are of greater importance.

Nevertheless, this investigation shows that the derived models are well suited to simulate realistic system responses for arbitrary motion specifications that have not been used for parameter identification. Eventually, this means that the models can also be used for automated trajectory optimization in an expedient way, as initially required.

In order to get a better impression of the derived results and model parameters, they are shown in Figure 11 as a function of the varied fluid width w. In addition, an approximation equation is given in the literature [7] for the resulting first undamped natural frequency, which is supplementary illustrated in Figure 11:

$${\omega }_0=\sqrt{\frac{\pi g}{w}tanh\left(\frac{\pi h}{w}\right)}.$$

(20)

For describing the damping ratio $\zeta$, no approximation formula is known to the author’s knowledge.

The undamped natural frequency decreases by increasing the fluid’s width. It can be seen that there is very good agreement between the derived model parameter from the experiments and the predicted parameters from Equation (20). This observation fits the vibration theory, since it is well known that, with increasing mass, which results from a larger volume, natural frequency decreases. Therefore, the results are plausible. With one exception, the damping ratio also decreases with the fluid’s width. This coherence is evident by the decreasing influence of the dissipative frictional energy between the fluid and the container walls in relation to the total kinetic and potential energy, simultaneously with the fluid width.

Whether the derived models are applicable for predicting the residual vibration had to be investigated. Therefore, the measure of ${\eta }_{\mathrm{eff}}$ introduced in Equation (17) was evaluated for the measured fluid edge, the line approximation, and the model response. An exemplarily result is displayed in Figure 12. Furthermore, the difference measures introduced in Equations (18) and (19) were calculated and are additionally shown. Overall a good accordance of the different measures can be seen, whereby expectably the difference measure $\Delta {\eta }_{\mathrm{simul}}$ is a bit greater than $\Delta {\eta }_{\mathrm{appr}}$. In

Table 5. Summary of the evaluated effective fluid edge measures for all investigated geometries and their corresponding differences from approximation and model application.

Geometry	$\mathbf{max}\left\{{\mathit{\eta }}_{\mathbf{eff},\mathbf{meas}}\right\}$	$\Delta {\mathit{\eta }}_{\mathbf{appr}}$	$\Delta {\mathit{\eta }}_{\mathbf{simul}}$
1	1.49 mm	0.06 mm	0.16 mm
2	1.36 mm	0.07 mm	0.13 mm
3	1.61 mm	0.07 mm	0.16 mm
4	1.46 mm	0.07 mm	0.16 mm
5	1.64 mm	0.06 mm	0.15 mm
6	1.56 mm	0.08 mm	0.20 mm

, these measures are shown for all investigated geometries. The evaluated differences are very small and hardly vary among the different geometries. Hence, in summary, the derived models are capable of predicting residual vibrations for motion profiles other than those used for parameter identification.

5. Conclusions

The investigation of the deviations between the measured fluid surface and the approximated line allows for the determination of the model assumption previously assumed in the literature. Thus, a method has been presented, and it can be applied generically to determine model limits in terms of the applicable operating speed range. In addition, which previous model approaches in the form of ordinary differential equations are suitable for determining sloshing behavior over a determined operating speed range was investigated. It is shown that the simplest possible approach of the linearized pendulum equation is sufficient to describe the measured data. Finally, whether the derived models are also suitable for predicting the residual vibration for other motion profiles was confirmed. The presented methodical approach thus allows for the provision of models for automated motion optimization in packaging machines. This will be treated in future works as well as increasing the trajectory’s robustness by the consideration of parameter uncertainties.