Simulation Tests of a Cow Milking Machine—Analysis of Design Parameters

: The objective of this paper was to create a mathematical model of vacuum drops in a form that enables the testing of the impact of design parameters of a milking cluster on the values of vacuum drops in the claw. Simulation tests of the milking cluster were conducted, with the use of a simpliﬁed model of vacuum drops in the form of a fourth-degree polynomial. Sensitivity analysis and a simulation of a model with a simpliﬁed structure of vacuum drops in the claw were carried out. As a result, the impact of the milking machine’s design parameters on the milking process could be analysed. The results showed that a change in the local loss and linear drag coefﬁcient in the long milk duct will have a lower impact on vacuum drops if a smaller ﬂux of inlet air, a higher head of the air/liquid mix, and a higher diameter of the long milk tube are used.


Introduction
A milking machine is a device designed for the efficient yielding of milk from a cow and for transporting milk to the milk tube. The milking process should be effective; therefore, milking machine settings are aimed at achieving high milking efficiency, that is, a short machine operation time combined with the maximum quantity of milk [1]. The function of the milking parlour is one of the factors which affect the efficiency of milk production on the farm. The increased dairy herds also require the modernization of milking equipment. [2,3]. It is important to find the appropriate criteria and technical parameters for milking parlours that would allow for optimal type of milking parlour to be chosen [4,5]. In large dairy farms robotic milking, also known as Automatic Milking Systems (AMS), are used [6]. AMS are internationally accepted as a valid alternative to conventional milking parlours, and as an advanced means of dairy farm management [7,8]. The conditions for the mechanical milking of cows are analyzed, taking various design parameters of the milking machine into account [9][10][11]. The goal of these studies is to make the machine-based milking of cows more efficient [12,13] and improve the efficiency of the milking process [14]. Appropriate milking parameters have a significant influence on cow health [15][16][17][18], and a key role in yielding milk of an appropriate level of value [7,19,20].
Among these studies, a special position is occupied by tests of vacuum changes in the milking equipment [21], which is the main technical parameter of a mechanical milking machine. Vacuum drops and oscillations in the milking cluster reduce milking efficiency and constitute the factors that are mostly responsible for udders diseases in mechanically milked cows. These cause return flows of the air/milk mix, which can achieve significant speeds with certain design solutions for the milking machine [22][23][24][25][26].
Irregular vacuum oscillations during milking, combined with cyclical vacuum oscillations and a very high value of vacuum at the end of the teat, seem to increase the frequency of udder infections and reduce milk flow [27]. A higher vacuum at the end of the teat decreases the peak flow rate of milk but may also increase hyperaemia or oedema of the teat [28]. If the working vacuum in the milking machine is too high, milking cups climb up the teats more intensively, possibly causing teat tip injuries. A high vacuum in the milking machine may cause udder infections, fluid accumulation, and teat tip tissue clogging [1,29,30]. Too high a pressure can lead to increased teat wall thickness [31]. Several authors analysed the factors affecting the vacuum inside a milking cluster with the use of simulated milking equipment and by measuring milking claw vacuum when adjusting the flow rate [32,33]. Others tested how a high vacuum affects the efficiency of mechanical milking, the condition of teats, and the health of milking cows [34][35][36]. Ambord and Bruckmaier [37] assessed the impact of various milking machine settings, causing different values of vacuum loss, on milking properties and changes in teat issue. Wiercioch et al. [38] found that a milking cluster with ventilated teat rubbers, despite creating the lowest vacuum, reduces vacuum oscillations in a cycle to the greatest extent.
The vacuum drop is increased with high milk flow rates [37,39] and small or blocked diameters of air inlets in the claw [39,40]. The operating parameters of the milking cluster should be selected to ensure the lowest possible vacuum drops.
The parameters of the milking cluster can be tested using real equipment or mathematical models and simulations. Mathematical modelling is an essential tool for analysing the impact of selected design components of the milking cluster on the course and extent of these phenomena. A model tying all milking cluster components makes it possible to analyse, by simulation, the qualitative and quantitative variations in the described values. A properly completed model may be used to control the process or for optimisation.
Demba [41] studied the effect of different milking settings on the teat-liner interface with the help of a pressure-indicating film. The experiment was performed with a conventional milking cluster equipped with round silicone liners. Ipema and Hogewerf [42] presented the design of a milking stall with special functions for monitoring and control, and tested the effect of technical milking machine settings (pulsation or milking vacuum, pulsation rate or ratio and detachment level) on quarter milking performance. Upton et al. [43] described a novel quarter milking analysis device for use in the field of bovine milking research. Wiercioch [44] benchmarked the operating parameters of pulsators in a milking process, simulated in a laboratory during disturbances in the milking machine's vacuum system.
The milking cluster was tested with the use of various methods and models, such as linear regression [45,46], non-linear regression [22,25], or square dependencies [47]. To control the milking procedure, fuzzy logic models were used [48,49].
Tan et al. [50][51][52] presented dynamic deterministic models that define pressure levels in a milking machine. The change in air mass in milking cluster chambers was also defined by differential equations, and the existing pressure levels existing were calculated based on the ideal gas law. Kupczyk [23] presented a dynamic model of pressure variations in under-teat chambers of milking cups in the claw's milk chamber.
Skalska et al. [53] prepared models to describe changes in the selected milking parameters, i.e., the average suction vacuum, suction vacuum amplitude, and massage vacuum amplitude in measuring milking clusters. The operation of a milking cluster was also modelled with computer simulations using Matlab ® Simulink. The authors modelled vacuum in the claw's milk chamber for 0, 1, 2, and 3 mm changes in air flux diameters [54,55] and developed a simulation model for a system controlling milk flow in the collector column of an autonomous milking cluster [56]. A control system that will adjust the area of the aerating hole in such a way as to maintain a minimum vacuum drop throughout the milking procedure when the milk flow rate is changing was proposed [57].
The literature also mentions applications that support the operation of a milking cluster: one that controls an autonomous milking cluster to regulate vacuum levels in the under-teat chamber of milking cups [58], and one that supports the calculation of selected pressure parameters of milking, which is based on mathematical formulas available in the source literature [59].
The examination of physical phenomena occurring in the milking cluster of a milking machine and their analysis using mathematical models make it possible to determine the factors that are most responsible for ensuring a proper milking procedure and identifying the factors that cause injuries and diseases in cow udders. The available literature examining the milking process does not include any optimisation models. Mathematical models based on a physical description of phenomena occurring in a milking machine make it possible to simulate the operation of a milking machine and build optimisation models. By using the sensitivity analysis in research using mathematical models, it is also possible to assess the influence of the input variables on the value of the output variable [60,61]. One example is Bernoulli's equation for the head loss, which defines the vacuum drop in a long milk tube [23].
The objective of this paper was to create a mathematical model of vacuum drops in a form that enables the testing of the impact of design parameters of a milking cluster on the values of vacuum drops in the claw.

Materials and Methods
Simulation tests and a sensitivity analysis of a cow milking cluster were conducted with the use of a simplified model of vacuum drops in the form of a fourth-degree polynomial.

Simplified Structure of the Vacuum Drop Model
Vacuum drops in a milking machine are defined by Bernoulli's equation which is a physical model of the flow of a mixture of milk and air in the long milk tube. In this equation, the reduced velocity of mixture u M is the sum of the reduced velocities of milk and air. The paper includes the following assumptions [62]: • Linear velocity is the sum of the velocity of the milk/air mixture with coefficient characterising a turbulent flow (1.2 u M ) and rise velocity of a single bubble ∞ in a quiescent liquid; • The reduced air velocity requires the actual conditions to be lowered to normal conditions (a change in pressure causes a change in gas volume); • The rise velocity of a single bubble in a quiescent liquid is the function of Archimedes (Ar) and Eötvös (Eo) numbers, which, due to the air density, depend on pressure p in the long milk tube; • Volumetric air coefficient α pdpm in the long milk tube is the quotient of the reduced velocity and linear velocity (calculated relative to the tube walls); • Pressure in the milking claw is equal to operating pressure p o plus drop ∆p kol .
In addition, considering the calculation's assumption that p = (p o + ∆p kol )/2, and taking the assumed relationships between the reduced velocity of the (milk/air) mixture and volumetric air coefficient α pdpm into account, Majkowska [63] determined a simplified Bernoulli's equation in the following form: Velocity v ∞ is calculated using the following formula: where: where: The discussed relationships made it possible to assume the bubble velocity in the vacuum drop model to be constant for a given milk tube diameter. After substitutions ∆p kol = x and v ∞ = c, and taking into account that p N = 100 kPa = 10 5 Pa, expression (2) takes the following form: Transforming this, we get: Multiplying by the factor πD 2 (2p o +x) 3 and denoting v ∞ = c, polynomial is obtained: After transformations and carrying the terms to one side, Equation (8) takes the following form: Giving Bernoulli's equation as a polynomial enables the use of more effective methods of calculating x from Equation (9), simplifying the analysis of vacuum drops in the long milk tube.
The searched vacuum drops are real roots of the following function: Absolute term w w of polynomial (10) has the form: Given that all its factors are positive, it can be concluded that, irrespective of the design parameters of a milking cluster, the vacuum drop will not be equal to zero (the number zero cannot be a zero of a polynomial).

Verification of a Model with a Simplified Structure of Vacuum Drops in the Claw
The vacuum drops calculated by the model (1) and the simplified model (10) were compared. The comparison of the models concerned different milking speeds and different air streams admitted in milking machines with a high-lying milk tube and canning machines for the diameters of a long milk tube in the range from 0.01 to 0.02 m. Relative errors of the simplified model (10) and model (1), regardless of the type of milking machine, milking speed and admitted air stream, do not exceed 3%. The small discrepancy in the models allows for Bernoulli Equation (1) to be replaced with its approximation, with a fourth degree polynomial (10). Figures 1 and 2 show the relative errors of the drops, calculated on the basis of model  (10), in relation to the drops calculated using model (1). The results confirm that the models are compatible.
ent air streams admitted in milking machines with a high-lying milk tube and canning machines for the diameters of a long milk tube in the range from 0.01 to 0.02 m. Relative errors of the simplified model (10) and model (1), regardless of the type of milking machine, milking speed and admitted air stream, do not exceed 3%. The small discrepancy in the models allows for Bernoulli Equation (1) to be replaced with its approximation, with a fourth degree polynomial (10). Figures 1 and 2 show the relative errors of the drops, calculated on the basis of model (10), in relation to the drops calculated using model (1). The results confirm that the models are compatible.   The model of vacuum drops (10) was also verified with the experimental data presented in the work [23]. The comparison of the simplified model with laboratory data also shows the high compatibility of the simplified model (10). The mean absolute error of this comparison is 0.03 kPa. The obtained results authorize the adoption of a simplified form of the vacuum drops model in the long milk tube for further consideration.

Sensitivity of a Model with a Simplified Structure of Vacuum Drops in the Claw
A model sensitivity test includes the examination of the value of changes in the model status variable when model parameters are changed at preset input variables. In the vacuum drop model, the division of input values into input variables and parameters arises from the underlying assumptions. When the parameters are changed from As to As +ΔAs, the status of the model will change from x to (x +Δx + o(║As║) ). The drop model equation The model of vacuum drops (10) was also verified with the experimental data presented in the work [23]. The comparison of the simplified model with laboratory data also shows the high compatibility of the simplified model (10). The mean absolute error of this comparison is 0.03 kPa. The obtained results authorize the adoption of a simplified form of the vacuum drops model in the long milk tube for further consideration.

Sensitivity of a Model with a Simplified Structure of Vacuum Drops in the Claw
A model sensitivity test includes the examination of the value of changes in the model status variable when model parameters are changed at preset input variables. In the vacuum drop model, the division of input values into input variables and parameters arises from the underlying assumptions. When the parameters are changed from A s to A s +∆A s , the status of the model will change from x to (x + ∆x + o(||A s ||) ). The drop model equation meets the assumptions of the implicit function theorem. A set of values of polynomial W(x) was analysed for coefficients that were the functions of preset parameters A s and input variables U s present in a milking cluster (vectors A s and U s ). Elements of vector U s belong to set Ω s : For input variables U s belonging to domain Ω s (12), the polynomial includes real roots. At each point, the polynomial is a function of class C 1 . Given the simple form of function (10), the following derivatives were calculated analytically: It was numerically tested that, in domain Ω and for the drops observed in reality, that is, for x∈(0, 5000) [Pa], the derivatives of polynomial (10) after x do not equal zero. At the determined point in domain P = [x 0 , A s0 , U s0 ], the derivative of the function defining an explicit form of variable x relative to the parameters composing vector A s is expressed by the following formula: The differential ∆x of the vacuum drops is expressed as a scalar product: The vacuum drop model and the related sensitivity model are shown in the block diagram in Figure 3.
The differential x of the vacuum drops is expressed as a scalar product: The vacuum drop model and the related sensitivity model are shown in the block diagram in Figure 3.  Equation (16) makes it possible to examine the sensitivity of the drop model in a simplified structure, relative to any disturbances introduced to the parameters (and input variables). It is a convenient tool for simulating variations in drops resulting from changes in the input of the drop model. Equation (15) allows for the values of partial derivatives (11) to be traced. By analysing the sensitivity of the model to changes in the value of the vector of variables and input parameters, the strength and direction of the impact of changes in input values on changes in the status of the model can be examined. In some cases, sensitivity analysis may indicate the existence of an extremum of the analysed function.

Results and Discussion
Sensitivity analysis of the vacuum drop model was performed with the use of formula 16. The values of differential ∆x(A s0 , U s0 , Q m ) referred to the values of drops x(A s0 , U s0 , Q m ), and are expressed as a percentage. The sign of the differential indicates the direction of changes, while its absolute value shows the strength of the impact of the change in parameter A si on changes in value x at point (A s0 , U s0 , Q m ). Quotient ∆x is often used as a measure of model sensitivity. In this paper, ten percent changes in the parameters were analysed. The analysis of the results presented later in the paper will concern the expression ∆x x * 100%. The sensitivity of model (9) was tested for points from the allowable set Ω s (12).
The diagrams in Figures 4 and 5 show relative values of the differentials expressed as a percentage. They correspond to the values calculated for a series of points P 0i , given in the legend below. A point in the legend is clearly identified with a colour and an index. The value of the differential (Figure 4) is calculated at the point with index "i". The diagrams in Figures 4 and 5 show relative values of the differentials expressed as a percentage. They correspond to the values calculated for a series of points P0i, given in the legend below. A point in the legend is clearly identified with a colour and an index. The value of the differential (Figure 4) is calculated at the point with index "i".
Indices i of the points given in the legend   Milk mass flux q m is a value set for the actual milking phase. In practice, this value is not constant. Similar calculations were performed for different milk mass fluxes. Change q m does not alter the opinions on the model's sensitivity to changes in the parameters. The observations point to the following relationships between the design parameters of the milking machine and the milking process: • At lower heads of the milk/air mixture, the model is more sensitive to variations in parameters ζ and λ ( Figure 5, black colour). Local mechanical energy losses in the pipes are caused by various obstacles located in the tubes. The values of local loss coefficients ζ should be determined empirically based on measurements; these values are affected by, for example, pipe curvatures, kinks in tubes, an abrupt increase or decrease in the pipe cross-section, etc. Linear loss coefficient λ depends on two parameters: Reynolds number (Re) and relative roughness of a pipe (e), which is a non-dimensional parameter. Relative roughness is defined as:

Summary and Conclusions
The preparation of a vacuum drop model in the claw, in the form of a model with a simplified structure, by expressing Bernoulli's equation as a fourth-degree polynomial, allowed the sensitivity model to be analytically developed. Examinations of the sensitivity of a simplified-structure drop model to any disturbances introduced to the parameters resulted in the formulation of a number of recommendations for the design parameters of the milking machine. Based on these examinations, certain general relationships can be presented:

•
A change in local loss coefficient ς and linear drag coefficient λ in the long milk tube will have a lower impact on vacuum drops if a smaller flux of inlet air Q p , a higher head of the air/liquid mix H, and a higher diameter of the long milk tube D are used; • Operating pressure does not affect these changes; • The model is insensitive to changes in parameter T: ten percent changes in temperature T caused changes in the drops to be lower than one percent.
The methodology applied in this paper makes it possible to carry out a broader analysis of the impact of design parameters of milking machines on the milking procedure, leading to improvements in the design of these machines.