Application of Constrained Optimization Techniques in Optimal Shape Design of a Freezer to Dosing Line Splitter for Ice Cream Production

Design of multiple branches splitting of equal mass flow rate in complex rheological flows like ice cream near melting point temperature can be a challenging task. Pulsations in flow rate due to air pumping process and small fluctuations in temperature affecting flow rheology can determine a consistent difference in internal pipe velocity distribution, resulting in a significant difference in the distribution of ice cream dosage. Computational sciences and engineering techniques have allowed a major change in the way products and equipment can be engineered, as a computational model simulating physical processes can be more easily obtained, rather than making prototypes and performing multiple experiments. Among such techniques, optimal shape design (OSD) represents an interesting approach. In OSD, the essential element with respect to classical numerical simulations in fixed geometrical configurations relays on the introduction a certain amount of geometrical degrees of freedom as a part of the unknowns. This implies that the geometry is not completely defined, but part of it is allowed to move dynamically in order to minimize or maximize an objective function. From a mathematical point of view, OSD is a branch of differentiable optimization and more precisely the application of optimal control for distributed systems. OSD is still today numerically difficult to implement, because it relies on a computer intensive activity and moreover because the concept of “optimal” is a compromise between shapes that are good with respect to several criteria. In this work, the applications of a multivariate constrained optimization algorithm is proposed in the case of a mechanical ice cream 1 to 5 splitting system, required to distribute in an evenly way from one freezer into five dosing valves. Results allowed to design a retro-fitting system on an existing production plant reducing the dosing error down to 3% on the average.


Introduction
Mass flow rate splitter design in complex rheological flows like ice cream near melting point temperature can be a challenging task. Pulsations in mass flow rate, due to air/ frozen mixture pumping process and small fluctuations in temperature, strongly affect the flow rheology, and they can determine a consistent nonuniformity in final dosing valves, resulting in a significant difference in the distribution of the dosage. Moreover, shortcuts in piping design like the adoption of fully symmetric 3D configurations cannot be easily arranged in existing pilot plants, for example, because an odd number of dosing valves or a direction change in tubing configuration due to a freezer relative position does not allow such configuration. Application of the heuristic approach in splitter design in strongly nonlinear problems results almost impossible to adopt, especially when small dosing errors (i.e., less than 3%) are required.
On the other hand, computational science and engineering techniques [1] have allowed a major change in the way that products and equipment can be engineered, as a computational model that simulates the physical processes can be built rather than building real-world prototypes and performing experiments. Among such techniques, optimal shape design (OSD) [2] represents an interesting approach. In OSD, the essential element with respect to classical numerical simulations in fixed geometrical configurations is to introduce a certain amount of geometrical degrees of freedom as a part of the unknowns, which means that the geometry is not completely defined, but part of it is allowed to move dynamically in order to minimize or maximize the objective function. The applications of optimal shape design (OSD) are numerous. For systems governed by partial differential equations, they range from structural mechanics [3], heat transfer [4], acoustic [5], extrusion [6], electromagnetism [7], and fluid mechanics [8]. Almost none has been done in food engineering, apart from an application in aseptic processing [9] and another in hopper discharge design [10].
In food industry, optimum design is not a once and for all solution tool because engineering design is made of compromises owing to the multidisciplinary aspects of the problems and the necessity of doing multipoint constrained design [1]. From a mathematical point of view, OSD is a branch of differentiable optimization and more precisely the application of optimal control for distributed systems [11], where gradient and Newton-based methods are the natural numerical tools. The problem is that OSD is still numerically challenging, because it is computer intensive and moreover because the "optimal" is a compromise between shapes that are good with respect to several criteria. In this work, the applications of a multivariate constrained optimization algorithm are proposed in the case of a mechanical ice cream splitter, required to distribute in an evenly way a volumetric flow rate from one freezer into five dosing valves. In this case, the optimal configuration is defined considering the possibility to perform a final fine tuning in a real implementation using a regulation valve, taking into account the possibility to handle possible oscillations of process parameters leading to an uneven distribution of the ice cream.

Ice Cream Rheology
Ice cream can be considered a foam composed by a matrix (a solution of sugars, milk proteins, stabilizers, salts) with ice crystals, fat, and air inclusions. The rheological behavior depends on several different parameter like composition, percentage of air inclusions (overrun), and temperature, and they can significantly change over each stage of the continuous production process [12].
In ice cream production process, the preliminary step requires ingredients to be blended, pasteurized, homogenized, cooled down to 4 °C, and then aged to form a premix, which is later introduced into a scraped surface heat exchanger (the freezer) where it is frozen using a process temperature ranging from − 4 to − 8 °C, and small bubbles of air are entrained into the mixture (D air < 30 µm).
After freezing, approximately 65% of the ice cream is formed by a highly viscous solution of polysaccharides, milk proteins, sugar, lipids, and a small amount of emulsifiers, while the remaining 35% is water present as ice crystals (D ice < 50 µm). The air bubbles raise up to 50% of the volume at atmospheric pressure, and the ice crystals around 10%. A subsequent high-shear processing can induce D air to break down to less than 15 µm in diameter, resulting in a more acceptable final product [13].
After the product is inserted into appropriate packaging, it is cooled down to below − 18 °C, hardened, and an additional 40% of the water is incorporated into the existing texture, resulting into an ice cream with the familiar texture and rheology [13,14]. Rheological studies have demonstrated a viscoelastic behavior strongly related to microscale characteristics [15][16][17][18], including an investigation on the effects of microstructure on an ice cream's oscillatory rheology [19].
In this work, a base vanilla flavored ice cream mixing was used in the experimental part with the following ingredients (Table 1): During the process for final experimental test, the overrun OR was set to 100%, defined as where W l_mix is the weight of unit volume of liquid mix before freezing and W ice cream is the weight of unit volume of ice cream.
Rheological properties used in numerical simulations for OSD design were obtained from the work of Elhweg et al. [20] for industrial ice cream.
Following their work, a power-law model for the ice cream viscosity expressed in (Pa s) was adopted: A flow index n equal to 0.5 was assumed to apply across all temperatures, and the apparent viscosity data yielded an approximate consistency dependency given by where a = 2 (Pa s n ), b = − 0.7 (1/°C) and T is the temperature in degree Celsius. The thermal conductivity k, expressed in W/(m °C), below the freezing point, was estimated using the data published by Cogne et al. [21], resulting in where T is the temperature in degree Celsius. The enthalpy expressed in Joules per kilogram was obtained by fitting data reported by [13] obtaining the following: And the specific heat capacity Cp, which is temperature dependent and expressed in J/(kg °C), was obtained from differentiation of Eq. (5) (Cp = dH/dT), while density was set to 540 kg/m 3 . As it will be detailed in the following sections, rheological and physical characterizations of ice cream are important as in this work; the complete set of momentum and energy conservation equations are solved in the optimization problem. The Geometrical Constrains of the Problem Target of this work was to develop a dosing system with a theoretical maximum error among 3 to 5% of the final weight of an ice cream (cornetto ice cream typology) in a linear dosing system with 5 cones per line (one to five linear splitting system) to replace an existing production system.
The original setup in the industrial production line is shown in Fig. 1, together with the experimental average dosing error for production rate equal to 600 kg/h: Although the original setup is symmetric, the experimental data are not, due to small irregularities in internal pipe flow, small temperature changes due to viscous dissipation effects, and external insulation (made by the ice accumulating on the external piping surfaces after a transitory time) characteristics. Moreover, viscous dissipation heating plays an important role, inducing an internal asymmetric flow field into a symmetrical piping configuration, a sort of "memory" in terms of internal velocity distribution retained by the fluid at each bending.
In liquid food production, dosage control is a critical issue: in order to fulfill legal requirements, in this case, a Dosing error for the original production splitting system maximum difference between real and declared weights equal to − 5%, an error up to − 20% in the line 3 valve ( Fig. 1) implies that the average dosing weight must be increased up to + 15%, with a significant economic loss with respect to the declared ice cream weight.
In our case, the optimization of production lines introduces further constrains, due to the topology of existing production plant: in this case, a second constrain was the request to limit dimensions of the new system in a box with a maximum depth of 40 cm and 3 m height.
Several approaches could be adopted to fulfill all requirements, for example a fully symmetric geometrical splitting as in [22], but while this approach is feasible for a one-toeven number of branches, it is not for one-to-odd branches configuration. For this reason, a different approach was adopted (see Fig. 2): the main mass flow rate (Q 1 ) is split into 3 secondary branches (splitter S 1 ), while two secondary splitters S 2 provide the final splitting to the five dosing valves V 1 -V 5 .
The target relations among all fluxes are the following: The rationale for this partially symmetric choice is the following: a "perfect" split from one to five streams has been considered not achievable in practical conditions due to irregularities like pressure oscillations during ice cream production, so that a final control system for tuning was introduced in line 3. In order to minimize the number of possible control systems and avoid nonlinear instantaneous feedback effects of the other dosing lines when one of them was perturbed (for example by using a control valve for each line), a minimum ("optimal") number of final control systems equal to one was chosen. For this reason, the mass flow rate of the streams (Q 3 ) was designed to exceed the other four of a minimum amount (5%) and a regulation valve V R1 was introduced on this line for fine tuning. The other four lines were designed to self-adjust mass flow rate automatically once the regulating valve is active.
The achievement of the final configuration preserving such constrains was obtained by using the OSD approach.

The OSD Approach
Traditionally, optimal shape design has been considered as a branch of the calculus of variations and more specifically of optimal control. This subject interfaces with no less than four fields: optimization, optimal control, partial differential equations (PDEs), and their numerical solutions-this is the most difficult aspect of the subject.
Quoting Pironneau [23], we can say that "it represents a special blend of mathematics and engineering which some engineers may find too theoretical and applied mathematicians too computational." From a technical point of view, if a computer program is used to obtain the numerical solution of a PDE describing an optimal shape design problem, an optimization algorithm is required and it will have to be written, while optimal control theory provides the basic techniques for computing the derivatives of the target functions with respect to the boundaries. So, in essence, optimal control and optimization may be applied when the "control" becomes associated with the shape of the domain. However, problems are encountered with numerical discretization as they need to be integrated to a finite element or control volume methods.
As previously mentioned, OSD is a branch of differentiable optimization for distributed systems [11], where often gradient and Newton-based methods are natural numerical tools. Despite the development of the available

Fig. 2
One to five mass flow rate splitting system topology computational power, OSD is still today a technique of not easy application, first of all because the concept of optimal shape varies in relation to the objective function adopted and is often a compromise between various optimal solutions. One possible approach is the use of Pareto optimality; adopting a mathematical theorem stating that says in simple situations Pareto optimal points are minimizers of some convex combination of all the criteria. The trouble is that such linear combinations can lead to stiff problems with many suboptimal nodes, requiring a global optimization algorithm such as genetic ones. Another problem is linked to the existence of solutions and differentiability of the criteria with respect to shape deformation, and it was the focus of most of the 1980s studies [24][25][26], after it became clear [27] that oscillations of shapes could lead to nonphysical solutions of the optimization problem in the limit, a phenomenon known as homogenization, leading to a new class of problems called topological optimization. Among different approaches, descent algorithms are generally used to find local or global minimum, based on the information given by the variation of the gradients numerically obtained. The sensitivity analysis allows the calculation of the objective function gradients, but it is a complex job, fraught of possible errors. Sometimes, it is not even possible, since the use of industrial solver with closed source codes does not allow the possibility to access the required extensive code rewriting. In some cases, it is preferable to the use of algorithms that do not use the computation of the gradients, such as genetic algorithms (GA).
They offer a certain robustness addition, since it does not require the objective function conditions of regularity. They also allow the development of innovative design compared with conventional methods because they are not based on local minimum criteria.
The possibility to obtain a solution for the constrained multivariate minimization/maximization problem using gradients methods depends on the topology of the solution, on the starting point and on the stability of the solution itself. To better explain this concept, in Fig. 3, a hypothetical configuration of the solution of the target function F(x) to be minimized (in a univariate problem for simplicity) is shown. Using a steepest descend method, if the starting point is A, following the gradient, no solution can be found for M 1 or M 2 ; if the starting point is B, we will land in M 1 , while in the case of C, we will obtain the minimum M 2 . Note that in some cases, multiple minimum can exist, so that the solution minimizing F(x) can be a local minimum (M 1 ) and not the global minimum (M 3 ). For complex set of PDE governing the problem, an a priori knowledge of the topology of the solutions is impossible to obtain, and the situation is even more complicate in the case of time-varying solution for unsteady problems.
In some cases, genetic algorithms are adopted, but gradient-based algorithms are quite efficient when we start from a reasonable configuration; in this work, the constrained optimization algorithm adopted is the Broyden-Fletcher-Goldfarb-Shanno (BFGS; [28]) algorithm, an iterative approach for solving unconstrained nonlinear optimization problems.
The BFGS method belongs to quasi-Newton methods, a class of hill-climbing optimization techniques that seek a stationary point of a (preferably twice continuously differentiable) function. For such problems, a necessary condition for optimality is that the gradient be zero. Newton's method and the BFGS methods are not guaranteed to converge unless the function has a quadratic Taylor expansion near an optimum. However, BFGS has proven to have acceptable performance even for nonsmooth optimization instances [2]. In quasi-Newton methods, the Hessian matrix of second derivatives is not computed. Instead, the Hessian matrix is approximated using updates specified by gradient evaluations (or approximate gradient evaluations). Quasi-Newton methods are generalizations of the secant method to find the root of the first derivative for multidimensional problems. In multidimensional problems, the secant equation does not specify a unique solution, and quasi-Newton methods differ in how they constrain the solution. The BFGS method is one of the most popular members of this class for solving constrained nonlinear optimization problems of the Quasi-Newton class.
As any of Newton-like methods, BFGS uses quadratic Taylor approximation of the objective function in a d-vicinity of x: where g(x) is the gradient vector and H(x) is the Hessian matrix. The necessary condition for a local minimum of q(d) with respect to d results in the linear system: which, in turn, gives the Newton direction d for line search: The exact Newton direction (which is subject to define in Newton-type methods) is reliable when • The Hessian matrix exists, and it is positive definite. • The difference between the true objective function and its quadratic approximation is not large.
In quasi-Newton methods, the idea is to use matrices which approximate the Hessian matrix and/or its inverse, instead of exact computing of the Hessian matrix (as in Newton-type methods). If we define matrices, B ≈ H and D ≈ H − 1.
The matrices are adjusted on each iteration and can be produced in many different ways ranging from very simple techniques to highly advanced schemes.
The BFGS method uses the BFGS updating formula which converges to H (x*): Where As a starting point, B 0 can be set to any symmetric positive definite matrix, for example and very often, the identity matrix.
The BFGS method exposes superlinear convergence; resource intensivity is estimated as O(n 2 ) per iteration for n-component argument vector.
The application of such technique to optimal shape problems in this work has been implemented in the following way. The methodological approach requires the definition of a set of physical variables x p (temperature, velocity, etc.) and a set of geometrical variables (points on the mesh free to move) x k n , where max(n) = N is the number of geometrical degrees of freedom (DOF) and k is an iterative index. To reduce the number of controlling parameters in the discretized mesh, not all the points on the moving boundaries are considered as new degrees of freedom. This is possible by using a dynamic CAD representation of the boundary (B-spline or Bezier patchwork) setting the smoothness of the geometrical representation in a mathematical sense. Figure 4 provides a graphical representation of the concept, where points P x are let free to move into a constrained angle, and all others computational nodes are moving inside a CAD surface representation, to preserve the smoothness of the solution.
Only the control nodes on the surface are considered as new DOF, while all the others will move according to the displacement of the previous ones.
The set of geometrical DOF is constrained between the envelope defined by two functions s 1 (x k n ) and s 2 (x k n ) which delimit the maximum and the minimum geometrical variation allowed and consequently the envelope in space that the new configuration can occupy: This constrains can be assigned depending on the problem, for example, linked to the mechanical position of the part or process considerations.
For a given set of partial differential equations (PDEs), where I(x p , x k n ) represents a required property of the system on the domain. Depending on the choice of the weighting function (x) , we can have an integral or a pointwise condition. The target function represents the optimal performance required to the system, defined by solving the associated minimum problem: As already mentioned, the existence of a local or a global minimum in a mathematical sense is strongly linked to the chosen target function Θ , to the search envelope defined by the geometrical constrains (Eq. (15) where is a defined error, then x k+1 n = x optimal else go to b.
Numerically, the minimization code was integrated inside a moving mesh approach for problems with a stationary solution. The moving mesh approach allows to use the previous solution as starting point to compute the new one based on the deformed mesh.
Numerical simulations were performed using Open-FOAM 3.0, using the moving mesh approach, with a number of computational cells ranging from 500 k up to 1.2 M control volumes. Grid independence was checked on a partial configuration (single 3-way splitter), and the same computational grid density was retained for the complete setup. The optimization code was written in C++ and linked to the CFD solver by using a recursive call procedure from one code to the other.
Conservation equations of mass, momentum, and energy solved numerically in 3D Cartesian coordinates are described in the following lines: where e = h − p/ρ is the internal energy, h is the sensible enthalpy, U is the velocity vector, and τ is the stress tensor related to the strain rate by the following: neglecting compressibility effects, and τ:∇U is the viscous dissipation term.

Results
Numerical results are always required to be verified in experimental test, and for this purpose, the proposed framework was designed to be tested experimentally on a pilot plant with reduced capability with respect to the full scale one, but with the same inline geometrical dosing configuration and spacing, with one freezer and five dosing valves. The mass flow rate was set to 400 kg/h, with a target weight per ice cream equal to 33 g and approximately 40 strokes per minute.
As a matter of fact, the easiest geometrical co-planar and symmetrical configuration for splitter S 2 and S 3 (Fig. 5a) does not preserve the symmetry for the internal velocity flow field, as radial symmetry of the incoming flow is broken into 3 different tubes.
A simulation of the velocity field for splitter S 1 (Fig. 6) shows that the internal flow field symmetry planes for splitter S 2 and S 3 are defined by an angle of ± 30° with respect to the linear dosing valve system (Fig. 5b). As a consequence, a rotation of the splitter S 2 and S 3 is required, with different pipe lengths to connect the final aligned dosing valves.
Two optimization problems needed to be solved: in the first, the protrusion P 1 in the picture order to increase the pressure drop into Q 3 line, to the purpose of setting a theoretical 5% higher mass flow rate in valve V 3 .
The second optimization problem solved by using an OSD approach is related to the correct angle of splitter 2 and 3 to recover correct equal mass flow rate among lines 1 and 2 and symmetrically on 4 and 5.
In both problems, optimal solutions are obtained for the minimum value of the target functions.

OSD Problem 1
Relative orientation (30° with respect to the 5 dosing valve lines) of secondary splitters S 2 and S 3 introduces an asymmetric final connection between the dosing valve and the vertical pipes at the exit of the 2-way splitter (see Fig. 7), with an unbalance of mass flow rate due to different pressure drops connected to wall shear.
Equalization of mass flow rate in valve V 1 and V 2 was obtained by changing the relative angle of the 2-way splitter with respect to the theoretical 30°.
The splitter S 2 (mechanical part in red) was allowed to rotate changing the angle (Fig. 7a), while connection from the final bend to the valves (green parts) was allowed to slide in order to accommodate the variable inclination angles.
The control variable is the rotating angle of the splitter with respect to the symmetrical configuration at 60°. The geometry will follow the rotation preserving the position of the dosing valves (the system was going to be implemented in an existing production plant), and the tube circular diameter was set constant.
The quadratic objective target function was set as Outflow velocity field from the 3-way splitter simulation (Fig. 6), with a partially closed valve on central branch (almost 60%) to reduce Q 3 mass flow rate down to Q tot /5 × 1.05%, was used as inflow boundary conditions, and results after 360 iterations with a maximum increment of 0.005° per iteration show an error < 1% for = 56.247 •

OSD Problem 2
The second optimization problem was required to release the necessity to close the regulating valve VR up to 60% to setup the initial ratio between mass flow rate Q 3 with respect to Q 2 and Q 4 , so that the valve could be used for the final fine tuning to compensate pressure oscillation and small mechanical irregularities always present when a mechanical part is produced.
In Fig. 8, the initial geometry (a) and the final one (b) are presented. Required pressure drop was obtained by increasing the curvature of the bend (red part) rather than closing the regulation valve, by allowing 10 moving nodes ( X n k , k = 1 to 10) free to move in y direction. The moving nodes were geometrically linked by using a pretensioned B-spline, while the cross section was always perpendicular to this curve to preserve the tube diameter (35 mm).
The blue part (bend) geometry was fixed while the green one was allowed to slide to compensate the P 1 elongation.
The quadratic objective target function was set in the following way: Fig. 7 Optimal configuration (a) and convergence graph (b) for OSD problem 1 Fig. 8 Optimal shaper design problem 2: a) Initial geometry; b) Optimal elongation of "duck neck" protrusion P 1 to control mass flow rate Q 3 using wall shear-related pressure drop meaning that the central dosing valve is designed to provide 5% more ice cream when the regulating valve is open. The control variable in this case was the sliding distance L (Fig. 8).
Following the theoretically developed design, experimental validation was performed at the former Unilever R&D Center in Caivano, Italy, now no longer operational (see Fig. 9). An experimental pilot plant was used with a process mass flow rate equal to 400 kg/h, temperature T = − 5.7 °C, and an operative pressure equal to 4 bar. All pipes were insulated using an external foam, and tests were performed after 30 min of preliminary transient adjustment where the pilot plant was running up to reach a steady-state processing configuration, evaluated by weighting ice creams. During the steady-state ice cream production, target maximum dosing error was set to ± 5%. Fig. 9 Final optimal configuration of one-to-five branches splitting system In Table 2, results obtained in steady-state configuration in two different days are showed for a significative sequence of strokes. The average dose per stroke is slightly different due to pressure oscillation present in the freezing system, recorded by inserting a pressure sensor collecting data in the final pipelines before the dosing valves. The average % error, referred to the average instantaneous dose, is in mostly cases less than 3%, showing the capability of such system to self-adjust when small oscillations in pressure/mass flow rate are present.

Conclusions
Optimal shape design represents a powerful tool to design in an automatic way processing system in food engineering where the shape is strongly connected to the function. In this case, peculiar rheological characteristics of ice cream strongly affect flow behavior inside tubing systems, making very difficult to provide an optimal design using a heuristic approach. Geometrical constrains, often linked to industrial plant topologies, introduce further difficulties in the design process. In this work, the application of OSD allowed to obtain a splitting system configuration ready to be implemented into existing industrial lines, with a maximum dosing error of less than 5%.
Acknowledgments Former Unilever Research Center in Caivano, Italy is gratefully acknowledged for providing financial and experimental support.
Funding Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in Table 2 Dosing errors in ice cream production for two different test sequences. A, mass flow rate at V 1 ; B, mass flow rate at V 2 ; C, mass flow rate at V 3 ; D, mass flow rate at V 4 ; E, mass flow rate at V 5 . Italicized entries are dosing error < 3%; boldface entries are error > 3% and < 5% the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.