Olivier Masbernat
email: olivier.masbernat@ensiacet.fr

Frédéric Risso
email: frisso@imft.fr

Benjamin Lalanne
email: benjamin.lalanne@ensiacet.fr

Simon Bugeat
email: simon.bugeat@sodiaal.fr

Mikaël Berton
email: m.berton@absolutcooking.com

Prediction of size distribution in dairy cream homogenization

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Homogenization of emulsions consists in reducing the average size of micrometer-size droplets down to the order of a few hundred nanometers and is commonly performed in dairy industry for more than a century (???). The knowledge of size distribution and composition of the interface of fat globules after homogenization is crucial regarding a number of issues, such as metastability (?), rheology and texture of the resulting suspension (??). High-Pressure Homogenizer (HPH) is an apparatus in which the emulsion is forced to flow at a (nearly) constant flowrate within a thin gap where it is highly accelerated, before merging in an outlet chamber or impact ring where it is highly sheared.

The high mean shear rates and induced turbulence produced in the chamber are responsible for the deformation and breakup of the emulsion droplets. Therefore, predicting the size distribution in such a flow geometry a priori requires to model the fragmentation statistics in a complex highly heterogeneous turbulent flow with very short residence times (???). In the case of concentrated emulsions with complex interfaces such as a dairy cream, this problem is even more complicated due to additional bulk and interfacial rheology issues (??). In particular, it is likely that classical Hinze-Kolmogorov theory which is often applied to model the breakup process of a drop in turbulent flows (?) is not appropriate in the case of high-viscosity milk fat globules, which can be stretched over length scales that may exceed any of the turbulent length scale present in the flow. For such systems, it seems that only macroscopic correlations relating mean or maximum stable diameter to the dissipation rate or operating pressure in the HPH can be proposed to the food processing engineers (?). Such correlations are of limited range regarding the upscaling issue of homogenization processes. In the present work, the size distributions of the fat globules resulting from the homogenization of dairy cream at various (moderate) pressures have been measured in two HPHs: a pilot one at a reduced scale, and an industrial one at practical operating conditions. The comparison of measured outlet size distributions with those of native fat globules before homogenization revealed a simple breakup mechanism that that can be modeled by considering the viscous deformation of a droplet under an unsteady shear flow, in the absence of any shape relaxation. This simple model provides an accurate prediction of the size distribution of fat globules in homogenized dairy cream.

Materials and methods

Cream properties

Raw milk is collected from farms in the area of St Etienne (France) and stored at 4°C. The cream is produced at a flowrate of 20×10 3 L/h from a centrifugal plate skimmer (Westfalia Separator Type MSA 160-01-076) operated at 4700 rpm. The fat content of the cream, measured by infrared spectroscopy (Delta Instruments, Lactoscope FTIR FTA-3.4), is close to 400 g/L. Then it is adjusted to 35% w/w by addition of skimmed milk to the cream. The cream is then stored at 4°C during 2 to 4 hours. A sample of the cream is collected for granulometric analysis. A sample of butter (82% fat) was melted at 80°C in a falcon tube during 10 minutes and the fat was collected after water demixing. The internal viscosity of the fat globules was measured on milk fat samples extracted at 80°C from melted butter produced in the site of SODIAAL in Clermont-Ferrand during the same period the cream was produced, with no significant difference in fatty acids composition between cream and butter. Rheological properties of the fat extracted from the butter and the cream are therefore similar. The viscosity of fat was determined from stress-shear rate curves measured at controlled temperatures (75, 80 and 85°C) in a cone/plate viscometer (Brookfield DV-I Prime, BE, UK). Measurements were performed with a shear rate ranging between 2 and 40 -1 . A slight shear-shinning behavior is observed in the lower range of shear rate, followed by a Newtonian plateau above 10 -1 . The fat viscosity measured at 80°C is 10.5×10 -3 .

The cream density

and viscosity at 80°C can be estimated respectively to 960 kg/m 3 and 2.9×10 -3 , from the correlations proposed by ?.

Size distribution measurement

Cream samples are prepared as follows: 1 mL of homogenized cream is diluted in 9 mL of an aqueous solution containing 1% w/w of Sodium Dodecyl Sulfate (SDS) and 35 mM of Ethylene-diamine-tetra-acetic acid (EDTA), adjusted at pH=7 at ambient temperature. SDS prevents the flocculation of fat globules and EDTA dissociates casein micelles and prevents the appearance of a parasite peak around 100 nm in the size distribution (?). Few drops of this solution are introduced in the measurement cell of a Static Light Scattering granulometer (Mastersizer 2000, Malvern

Instruments, UK). The light sources are a He/Ne laser at 633 nm and a electroluminescent diode at 466 nm. Refractive indices selected for the measurements are 1.33 for the dispersing medium (water), and 1.452 and 1.460 for the milk fat, at respectively 633 and 466 nm, as prescribed by ?. Measurements were performed with an obscuration rate of 2% after the measurement of the noise level in water. The size distributions in volume were averaged over 3 consecutive measurements.

HPH equipment

Two High Pressure Homogenizers (HPHs) of different brands were used in this study, one at pilot scale and one at industrial scale. The pilot HPH is a GEA Niro Soavi NS3015H, and the industrial HPH is a TETRAPAK Alex Homogenizer 400. The flowrate was 200 L/h in the pilot and 8800 L/h in industrial scale HPHs, which corresponds to the same residence time (ratio between the HPH volume and the flowrate) of the emulsion in both HPHs ( =0.09 s). The typical internal geometry of a HPH is schematized in Fig. 1. It's an axially symmetric valve composed of an inlet section of radius 0 (of a few millimeters) through which the dairy emulsion is injected by a volumetric pump at a flowrate = 2 0 0 and is then forced to flow through a narrow gap of thickness , which is two orders of magnitude smaller than 0 .The emulsion is highly accelerated in the gap up to a velocity forming a jet that rams the wall of the impact ring, before flowing out of the device. More on the internal geometry of the two HPHs is given in appendix B.

With this geometry, mass conservation imposes that = 0 0 ∕(2 ), so >> 0 . A simplified Bernouilli equation therefore gives an estimate of the pressure drop Δ in a HPH as Δ ∼ 2 (with the emulsion density). For a typical pressure of 100 bars in a HPH, the emulsion reaches a velocity of the order of 100 m/s. The break-up of the emulsion drops occurs in the impact ring chamber (of width ∼ (10 2 × )) where the shear rate produced is intense ( ∕ ∼ 10 5 -1 ). In this study, the operating pressure Δ used for the concentrated dairy cream homogenization ranges between 15 and 60 bars in the pilot HPH and between 8 and 23 bars in the industrial HPH. The choice of a low range of operating pressure (compared to the usual range of pressure of several hundred bars in the case of milk homogenization) is motivated by the need of producing a metastable homogenized cream (i.e. that will not recoalesce after homogenization). The reservoir of amphiphilic components (such as caseins) being smaller in concentrated cream than in milk, the increase of interfacial area in cream is limited as compared to milk, hence explaining the use of a lower range of pressure for concentrated cream homogenization. The cream is heated at 80°C prior to being processed.

At the homogenizer outlet, the cream is cooled down to 20°C in a tubular heat exchanger and a sample is collected and stored at 4°C for granulometric analysis. Size distributions of cream fat globules obtained in these conditions are presented and discussed in the next section.

Results

Pilot HPH

Size distributions of fat globules at the outlet of the pilot HPH are displayed in Fig. 2 at different pressures. While the inlet diameter ranges between 0.8 and 8 m, the outlet diameter at Δ =60 bars ranges between 0.2 and 3.6 m. Interestingly, when comparing the inlet to the outlet size distribution at a given pressure, they look self-similar, i.e. the outlet distribution can be deduced by dividing the diameters of the inlet distribution by a unique factor , keeping constant the volume fraction. This self-similarity can be checked on the CPDs: for any given diameter of the inlet CPD (blue curve in Fig. 2(b)), the corresponding diameter of the outlet distribution at same value of the CPD, is obtained by drawing a horizontal line on the graph. One then can observe that the ratio between the inlet and outlet diameters is approximately constant for all diameters. This observation is verified at all pressures investigated.

This means that the outlet CPD matches the inlet one by multiplying the abscissa by a factor . This matching is illustrated in Fig. 3(a) where all curves are collapsing with the inlet distribution curve with a pretty good level of approximation, with a factor which is an increasing function of the operating pressure. Another way to illustrate the self-similarity property is to plot outlet-to-inlet pdf ratio, noted (= ( )∕ ( )) as a function of the inlet-tooutlet diameter ratio ( ∕ ) for few characteristic diameters along the size distribution. This graph is represented in Fig. 3(b) for 20 ( ( 20 ) = 20%), (diameter class of maximum volume fraction) and 90 ( ( 90 ) = 90%), indicated in Fig. 2(b). For all pressures investigated, the evolution is very close to the first bisector ( = ), validating the self-similarity property of all outlet size distributions. As a result, the outlet pdf ( ), can be deduced from the inlet pdf ( ), applying mass conservation:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ( ) = ( ) ⇒ ( ) = ( ) = ( ) ⇒ ( ) = ( ) (1)
in which the factor is a growing function of the operating pressure Δ . Equation (1) turns to consider that in average, each diameter class of the inlet distribution will break into single size fragments of size ∕ , as schematized in Fig. 4.

Under the action of a local shear ̇ , a globule of diameter is stretched into a filament of length that eventually breaks into equal size fragments of diameter = ∕ (their number being imposed by the volume conservation).

The pdf predicted by equation ( 1) has been plotted in Fig. 5 and compared to the experimental curves for each operating pressure investigated in the pilot HPH. For each curve, has been taken equal to the ratio ∕ from the corresponding experimental size distributions (cf Fig. 2(b)). The agreement between this model and the experimental data is remarkably good. This breakup model therefore describes the size distribution of 35% cream in the pilot HPH within a significant range of operating pressure, by making use of a single parameter .

Industrial scale HPH

The pdf and CPD of cream fat globules diameters measured in the industrial HPH are reported in Fig. 6(a) and 6(b) respectively. The pdfs are bimodal in this case, with the presence of a secondary distribution of finer fat globules separated from a primary one, which tends to get closer and slightly overlap with the secondary distribution, as the pressure in the HPH is increased. As shown by Fig. 6(b), the cumulated volume fraction of this secondary distribution is less than 10%.

Due to the presence of a secondary distribution, even if of low volume fraction, the self-similarity property of the outlet distribution with the inlet one is no longer conserved. However, this property can be evaluated on the primary distributions obtained at different pressures, trying to find for each of them, a mutliplication factor of the abscissa that would make them collapse on a single curve, as for the size distributions obtained at the pilot scale HPH. The result is shown in Fig. 7(a) where the collapse of the CPDs is quite convincing. They all merge with the inlet distribution for ≥ , where is the diameter of maximum volume fraction of the inlet pdf, equal to 2.82 m (cf. Fig.

6(b))

. The corresponding multiplication factor 1 of the absicssa has been reported as a function of the inlet-to-outlet maximum volume fraction diameter, ∕ 1 for all cases in Fig. 7(b) ( 1 is the diameter at maximum volume fraction of the primary distribution of the outlet distribution, see Fig. 6(b)). Here again, the matching between these two quantities is impressive, suggesting that each diameter class of the inlet distribution breaks into two classes of fragments of contrasted sizes, the larger being composed of single size fragments, proportional to the parent globule size. The smaller size class of fragments feeds the secondary distribution displayed in Fig. 6, and the maximum volume fraction diameter of the secondary distribution, noted 2 , seems to be nearly constant in the whole pressure range (around 0.16 m). The proportionality of the smaller fragment size to that of the parent globule can also be tested, taking the ratio ∕ 2 as the proportionality factor. Doing so, each size class of the inlet distribution will feed the outlet distribution with fragments of only two different sizes:

1 = ∕ 1 and 2 = ∕ 2 , with 1 = ∕ 1 and 2 = ∕ 2
. The corresponding breakup process is described in Fig. 8. Introducing the volume fraction Φ 2 of the smaller fragments 2 , and using mass conservation, one finds:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ( 2 ) 2 = Φ 2 ( ) ( 1 ) 1 = (1 -Φ 2 ) ( ) = 1 1 = 2 2 ⟹ ⎧ ⎪ ⎨ ⎪ ⎩ ( 2 ) = 2 Φ 2 ( 2 2 ) ( 1 ) = 1 (1 -Φ 2 ) ( 1 1 ) (2)
So the outlet size pdf ( ) writes:

( ) = 2 Φ 2 ( 2 ) + 1 (1 -Φ 2 ) ( 1 ) (3)
Note that equation ( 1) is the particular case of Φ 2 =0 in equation ( 3).

Using equation ( 3), the outlet pdf ( ) can then be deduced from the inlet size distribution ( ). Corresponding results are displayed in Fig. 9. The matching between experimental pdfs and modelled pdfs deduced from equation ( 3) is quite good, with a slight underestimation of the overlapping zone of the primary and secondary distribution, which represents no more than 2-3% in volume. As the pressure is increased from 8 to 23 bars, the proportionality factor 1 corresponding to the largest fragment is increasing, in the same way as for the pilot HPH, whereas 2 remains constant ( 2 ≃ 18) and Φ 2 increases from 6.5 to 8.7%.

These results show that the homogenization of cream in HPHs at different scales can be described by a model (given by equation ( 3)) containing three independent parameters, 1 , 2 and Φ 2 . 1 describes the largest fragments resulting from the fragmentation of an initial globule, whereas 2 describes the secondary fragments of the same breakup event (see Fig. 8). When, the HPH operating pressure increases, 1 is found to increase while 2 remains almost constant. In the case of homogenization of 35% cream in the pilot HPH, the volume fraction of the smaller fragments is negligible and the model reduces to a single parameter .

Milk in pilot scale HPH

Finally, a test of milk homogenization has been performed in the pilot HPH at 80°C and 120 bars. The resulting size pdf and CPD are reported in Fig. 10(a) and (b), respectively. In this case, the pdf of the outlet size distribution is also bi-modal, but with a strong overlapping of the primary and secondary distributions between 0.5 and 1 m.

Assuming the same fragmentation mechanism as for the cream, we have used model (3) to construct the outlet size distribution, taking 1 = ∕ 1 = 2.22, 2 = ∕ 2 = 10 and Φ 2 =0.36. As can be seen in Fig. 10, the match of the calculated distribution with the experimental curve is also quite correct in this case, suggesting that this fragmentation mechanism is more related to the milk fat globule properties than to the bulk flow properties. It is also interesting to note that the identified value of 1 for milk homogenization in the pilot HPH at 120 bars is smaller than that of for the cream in the same equipment at 45 bars. This result suggests a predominance of viscous stress contribution to the deformation of fat globules in the HPH. A physical interpretation of fat globule deformation and fragmentation in the HPH is discussed in the next section.

Derivation of a predictive model

It can be concluded from the preceding section that the homogenization of cream and milk in the HPH can be described by a simple fragmentation mechanism of each fat globule, independent of its size and valid in a wide range of operating pressure. In what follows, we propose a physical interpretation of this mechanism based on a simple model of fat globule deformation in a time-dependent flow field. ? and ? developped a numerical model describing the streching, breakup and mixing of viscous filaments, in 2D-chaotic viscous flows in which the leading forcing term is due to local elongational strain rate. Here, we consider that the fat globule continuously stretches with time under the action of a local shear rate in the HPH and ends to break into a single fragment size or two fragment sizes. The largest diameter of the produced fragments ( or 1 resp. in equations ( 1) and ( 3)) is supposed to be equal to the thickness of the stretched globule at the HPH outlet. This allows us to relate the proportionality constant in equation ( 1) (or 1 in equation ( 3)) to a single hydrodynamic parameter, which is a non-dimensional shear stress. This parameter is then correlated to the operating pressure through a macroscopic scaling analysis of dissipation rate in the HPH. Therefore, the model presented below intends to provide a physical interpretation of the proportionality factor (or 1 ) in the breakup model.

Dynamic model of fat globule stretching

The model is a balance between the hydrodynamic forces exerted by the external flow on the fat globule, which tends to deform it, and the resistive forces that tend to oppose this deformation. A first important assumption is that for a fat globule of diameter at 80°C, the resistance to the deformation is supposed to be controlled by the internal viscosity of the fat globule (mixture of triglycerides), the interfacial forces due to surface tension or viscoelastic modulus (?) being supposed negligible. This is the case because the Ohnesorge number of the globule is larger than one: ℎ = ∕ √ > 1, with the fat globule density and the interfacial tension (?). Note that interfacial tension includes here the viscoelastic contributions originating from interface structure and composition (?). During deformation, the shape of the stretched globule is assumed to be a cylinder of radius ( ) and length ( ), as sketched in Fig. 4 and8.

Then it is postulated that the external force per unit surface acting on the globule in the HPH is the viscous shear stress ̇ ( ), where is the effective viscosity of the external phase (35% cream) and ̇ ( ) is the shear rate acting on the globule along its trajectory in the HPH. This viscous stress is assumed to result from the mean planar jet flow developing at the gap outlet in the ring chamber. This assumption will be discussed and validated in the next two sections. The resulting force balance along the globule trajectory through the HPH reads:

1 = ̇ ( ) (4)
Making use of mass conservation (6 2 = 3 (where is the diameter of the native fat globule), equation (4) can be integrated over the residence time of the fat globule in the HPH. The averaged radius of the filament at the HPH outlet hence reads:

≃ 2 1 + 1 ⟨ ∫ 0 ̇ ( ) d ⟩ -1∕3 = 2 1 + 1 ⟨ ̇ ⟩ -1∕3 , ( 5
)
in which the brackets denote the average over the residence time and = ∕ is the viscosity ratio. Let ′ be defined as:

′ = 1 + 1 ⟨ ̇ ⟩ 1∕3 (6)
The proposed model therefore predicts that the thickness 2 of the stretched globule at the HPH outlet is propor-tional to its diameter:

2 = ∕ ′ (7)
In equation ( 7), the parameter ′ embeds the effect of the flow hydrodynamics in the HPH on the fat globule deformation.

Breakup criterion and fragmentation model

Due to the Rayleigh Plateau instability, a fluid cylinder is unconditionally unstable (?). Thus, we assume that all fat globules passing through the HPH eventually breaks after having being stretched into a long filament. A simple fragmentation model consists in assuming that the filament will eventually break into fragments of equal size or in two fragment sizes 1 and 2 as sketched in Fig. 8, with and 1 being of the order of the filament thickness 2 :

or 1 ≃ 2 = ∕ ′ (8)
The physical reality of such a deterministic fragmentation process has been experimentally evidenced in the paper of ⟨ ̇ ⟩, the higher , the higher . This result confirms the role of external phase viscosity in the breakup mechanism of fat globules in the HPH. The assumption of a viscous stress as the major cause of globule deformation can be justified here by the fact that any inertial contribution would be globule size-dependent and would result in a non-linear relation between the fragment size and the parent globule size. Note that this observation and the deterministic character of this fragmentation process makes unlikely a turbulence-induced breakup mechanism (inertial or viscous). The next section is devoted to the establishment of the relation between ⟨ ̇ ⟩ and Δ .

?

Modeling of the relation between and Δ

As presented in section 2.3, the emulsion in the HPH experiences a radial flow in a cylindrical gap of few tens of micrometers thickness, connecting the impact ring chamber (see Fig. 1). The deformation of the globule being supposed to be due to the viscous shear stress exerted by the mean jet flow in the impact ring chamber, the averaged shear rate ⟨ ̇ ⟩ can then be estimated according to:

⟨ ̇ ⟩ ∝ ∕2 (9)
Due to mass conservation in a cylindrical channel section, is inversely proportional to . Hence equation (9) leads to:

⟨ ̇ ⟩ ∝ Δ (10)
Next step is to determine the relation between Δ and the total pressure difference in the HPH, Δ . Δ is supposed to be the sum of two main contributions, the pressure drop in the gap radial flow of thickness , noted Δ , and the pressure drop due to the discharge in the impact ring chamber, Δ :

Δ ≃ Δ + Δ (11)
The Reynolds number in the gap being independent of the gap thickness, it is straightforward to show that the pressure drop in the gap Δ is proportional to -3 , whether the flow in the gap is laminar or turbulent (see Appendix A). The pressure drop Δ is assumed to be proportional to 2 , as in a sudden expansion flow, where is the mean velocity in the gap. Due to mass conservation in a cylindrical channel section, is inversely proportional to . Hence, Δ varies as -2 and Δ is proportional to Δ 3∕2 . Equation (11) can then be written as:

Δ = Δ 3∕2 + Δ (12)
The calculation of Δ and in both laminar and turbulent regimes is developped in Appendix A. In all cases, it is possible with equation ( 12) to relate for a given fluid system and at a given flowrate, Δ to the operating pressure Δ in the HPH. For the 35% cream at 80°C, the gap flow in the pilot HPH is found to be laminar, whereas it is turbulent in the industrial HPH. In both cases, the evolution of Δ with Δ provided by equation ( 12), can be very well represented by a power law (Δ = Δ ), with the prefactor and the exponent close to 1. Hence, in the range of operating pressure Δ investigated, the discharge pressure drop Δ can be fairly identified to the pressure drop in the HPH (Δ ≃ Δ ).

From equations (10), the proportionality factor ′ given by equation ( 6) can be expressed as a function of Δ according to:

′ = 1 + ′ Δ 1∕3 ( 13
)
where ′ is a parameter which depends upon HPH geometry and flow properties. This model has been tested in Fig. 11(a) against the experimental values of and 1 obtained with the 35% cream at 80°C, respectively in the pilot (closed blue rounds) and industrial (closed red rounds) HPH. The data can be reasonably well described by the following model:

′ = 0.6 {1 + 1.2Δ } 1∕3 (14)
This result first shows that the fragmentation model proposed leads to a suitable model of the fragmentation constants and 1 as a function of Δ . The value of the prefactor is found to be smaller than 1 (0.6), probably because in the low range of values, the approximation of the stretched globule as a cylinder is not accurate (an overestimation of 22% is observed for the lowest pressure in the pilot HPH). It also validates the scaling of ⟨ ̇ ⟩ given by (10).

Prediction of ′ from numerical simulations

To complete this study, the parameter ′ given by equation ( 6) has been estimated from numerical simulations of single-phase flow in the two HPHs geometries investigated in this study. The objective was to compute the average quantities ⟨ ̇ ⟩ and (residence time of the particle in the calculation domain) seen by a fat globule along its trajectory.

In these simulations, the cream was considered as a homogeneous Newtonian fluid, and flowrates were set identical to those of the experiments ( =200 and 8800 L/h for the pilot and industrial HPH, respectively). The numerical method, the calculation domain and typical flow structure are described in appendix B. One test case was performed 14)). In Fig. 11(a), the experimental value obtained for 1 in the case of milk at 120 bars in the pilot HPH has been also reported and its value is also widely below the trend observed for cream. These results confirm that the bulk viscosity favors the breakup efficiency and increase the value of or 1 . In order to account for the effect of viscosity, a correction under the form of a power law of the viscosity ratio -has been introduced in equation ( 14). Taking =0.4, this correction seem to regroup the data on a single curve, as shown in Fig. 12, with an average discrepancy of 8.5%.

These results give support to the present model for the viscous deformation mechanism and break-up of the fat globules in the HPHs, and confirm the relevance of the hydrodynamic parameter ⟨ ̇ ⟩ ∕ to quantify homogenization efficiency of dairy cream. They also confirm the scaling proposed for its relation with the operating pressure Δ , at both pilot and industrial scales.

It must be stressed that this simple micro-scale model of the stretching of a non-relaxing viscous cylindrical fila-ment does not claim to describe exactly the flow in a HPH, but it intends to account for the main physical mechanisms in order to understand the characteristics of the shape distributions generated in HPHs. There are two underlying reasons to the success of this model. First, the fat globule resistance to deformation is controled by the internal viscosity and not by the surface stress, so its relaxation time is always much larger than the characteristic time of shearing in a HPH.

In the particular case of dairy cream, this is characterized by a value of the Ohnesorge number larger than unity. The second reason is the HPH flow geometry, which forces the each fat globule to undergo the same range of deformation within a given residence time (like in a Couette apparatus), so the description of the forcing term at the scale of an individual globule by a single hydrodynamic scale -the averaged shear rate -is relevant for the whole emulsion. We therefore believe this model builds up a valuable tool for the dairy industry engineers, regarding the understanding of dairy cream homogenization and the optimization of HPH parameters.

Secondary distribution

The fragmentation model of fat globules in HPHs considers two contrasted fragment sizes (cf. equation ( 3)), both being proportional to their parent fat globule size. The larger fragment size has been related to the operating pressure, the smaller one seems to be quasi-constant and independent of Δ , at least in the range of operating pressure in the industrial HPH (between 8 and 23 bars). This breakup mechanism hence involves two more parameters that needs to be modeled: the size reduction constant 2 and the volume fraction Φ 2 (or the number of fragments). The formation of satellite drops during the fragmentation of a viscous filament is well documented in the literature (?).

Satellite drops results from the breakage of the thin bridge developing between two separating main bulbs in the latter stage of the fragmentation process. The production of satellite drops could depend on the local hydrodynamics in the HPH (Reynolds number and turbulence level), but we have no means to propose a predictive model for these satellite drops at the moment. The prediction of their size and volume fraction requires a detailed and complex analysis involving surface stress properties, which is beyond the scope of the present paper. However, the fact that the size of the secondary fragments are found to be proportional to that of parent fat globule is in qualitative agreement with the present deformation model, since for given shear rate and residence time, a smaller drop is stretched in a thinner filament than a larger drop, leading to a thinner bridge eventually separating fragments and to the formation of smaller satellite drops. Also, the increasing volume fraction of the satellite drops as the operating pressure is increased, is consistent with the fact that the larger the fat globule extension, the longer the developing thin film, and the larger the number of satellite drops.

Conclusion

Homogenization of 35% dairy cream at pilot and industrial scales was investigated in this work. The analysis of fat globules size distributions at the pilot HPH outlet operated at different pressures shows that they are self-similar to the inlet size distribution. This means that the size distribution of the homogenized cream can be deduced by considering a single fragment size for each fat globule entering the HPH with the same volume (or mass) fraction. Size distributions measured at the outlet of the industrial HPH were found to be bimodal with two characteristic distributions centered around a large and a small fragment sizes with a small overlapping. In this case, the two distributions can be also reconstructed from the inlet distribution by considering only two fragment sizes issued from each individual inlet globule. In both cases, the size of the largest fragment produced is proportional to the inlet fat globule diameter, and the proportionality factor is independent of the inlet globule size. This simple model was also successfull at predicting the size distribution of fat globules issued from milk homogenization in the pilot HPH. The physics of this fragmentation mechanism was then interpreted by considering the dynamics of the deformation of a fat globule into a cylindrical filament, resulting from the balance between the external time dependent viscous stress and the internal viscous stress of the stretched globule. Then it was assumed that the breakup of this filament at the HPH outlet leads to the formation of a single-or two-fragment size population, with the largest fragment size scaling as the filament thickness. From that theoretical framework, the self-similarity factor deduced from the experiments could be related to a single hydrodynamic parameter, a non-dimensional stress which was modeled as a function of the operating pressure of the HPH. The evaluation of this scaling law with experimental values of the proportionality factor obtained at both HPH scales showed a good agreement, also validated by numerical simulations performed in real HPH geometries.

This elementary fragmentation mechanism of fat globules is thought to result from two specific conditions: first the relaxation time of the fat globule shape is much larger than any flow time scale in the HPH, so breakup always occurs at maximum stretching of the globule. This property is due to the high internal viscosity of the fat globule, which overcomes any viscoelastic property of the interface. The second condition is due to the flow configuration generated in a HPH, where the whole emulsion is submitted to the same level of shear rate. It results that all fat globules flowing through the thin gap of the HPH experience the same average shear rate, so the hydrodynamic stress can be reduced to a single parameter. The theoretical interpretation of this fragmentation mechanism of fat globules in HPHs, even though not an exact representation of the reality, is therefore believed to constitute a relevant guide for the optimization of process parameters and scaling up of dairy cream homogenization.

Credit author statement

The contributions of each of the authors were essential to the completion of this study, which is the result of a complex fruitful human collaboration.

A. Pressure drop contributions in the HPH

A.1. Pressure drop in the gap channel

The pressure drop in the gap channel, of length and thickness can be estimated considering that the flow is established in the gap and that the flow section is constant between the entry and the outlet of the gap. The latter assumption is justified by the condition ∕ 0 << 1 which is verified in both HPH geometries, where 0 is the radius of the cylindrical gap (cf values of 0 and in table 1). The pressure drop in the gap can be expressed through the momentum balance integrated in the gap flow section:

Δ = 2 (15)
where is the wall shear stress in the gap. The velocity in the gag, , is equal to the flowrate, , divided by the gap crosssection:

= 2 0 (16)
The Reynolds number in the gap, is equal to:

= = 2 0 ( 17
)
where and are the cream density and viscosity respectively. Therefore, the Reynolds number is independent of the gap thickness. If the flow in the gap channel is laminar, the wall shear stress can be calculated assuming a steady established parabolic velocity profile, and taking its derivative at the wall:

= 6 (18)
In laminar regime, the pressure drop in the gap reads:

Δ = 12 2 (19)
Inserting ( 16) in (19) leads to:

Δ = 12 2 0 3 (20)
For a given flowrate in a given geometry, the pressure drop in the gap varies as -3 in laminar regime. If the flow in the gap channel is turbulent and established, the wall shear stress can be expressed as a function of the friction velocity at the wall, * , which is related to the mean velocity in the gap, through the friction factor ( ):

= * 2 = 2 2 (21)
The pressure drop in the gap hence reads in turbulent regime:

Δ = 12 2 4 2 2 0 3 (22)
As is a function of gap Reynolds number which is independent of gap thickness , in turbulent regime, the pressure drop in the gap is also proportionnal to -3 . Taking =958 ∕ 3 and =2.89×10 -3 for the 35% w/w cream at 80°C, can be calculated in both HPH geometries. In the pilot HPH, at a flowrate = 200 ∕ℎ, =730, and in the industrial HPH, =3200 for a flowrate = 8800 ∕ℎ. Hence, the flow is laminar in the pilot HPH and turbulent in the larger scale HPH. With milk at 80°C and = 200 ∕ℎ in the pilot HPH, =3500, so the flow in the gap channel is turbulent.

A.2. Pressure drop in the impact ring chamber

The pressure drop in the impact ring chamber, Δ can be approximated by:

Δ ∼ 2 2 (23)
As a result, for both laminar and turbulent regimes in the gap flow, Δ scales as Δ 3∕2 :

Δ = Δ 3∕2 (24)
dissipation. Simulations of pilot and industrial HPHs have been achieved with different numerical models. At pilot scale, a transitional problem occurs in the jet while a fully turbulent flow is identified at industrial scale. It is why all the simulations have been run using LES-VMS at the pilot scale and a -model was chosen for the simulation of the industrial HPH.

In these simulations, the density of the fluid is =958 kg/m 3 and the viscosity is =1.16×10 -3 Pa.s. The operating pressure, gap (average) thickness and gap average velocity are reported in table 4 for the pilot and industrial HPH. Gap Reynolds numbers (equation ( 17)) are respectively 1800 and 8000 for the pilot and industrial HPH.

In Figs.15 and 16, typical velocity fields respectively obtained in the pilot and industrial scale HPHs are displayed.

The regions of highest shear are generated around the jet, in the impact ring wall region, and in the recirculation loops developing in the whole chamber.

Fig. 2 (

 2 Fig. 2(a) shows the volume probability density function (pdf ), noted ( ) and 2(b) the corresponding cumulative probability distribution (CPD), noted ( ). It can be observed that the pdfs are monomodal and are shifting towards smaller diameters as the operating pressure is increased, as a result of the breakup of the fat globules in the HPH.

  where images of aligned equal size fragments have been captured at the outlet of a model HPH device (flow in a duct through a planar restriction), resulting from the breakup of elongated filaments. It is important to mention that the experiments of ? have been carried out with oil-in-water emulsions of high internal viscosity, corresponding to Ohnesorge numbers larger than 1. This fragmentation model predicts that each size class of the size distribution of the fat globules at the HPH entry, of probablity density ( ), either feeds a single class of fragment or two classes of fragments 1 and 2 , the probability density of which is imposed by the mass conservation expressed in equations (1) or (2). Therefore, according to the proposed model, values of and 1 , factors of respectively the single and two classes fragment breakup models (1) and (3), are equal to ′ . Equation (6) therefore relates the breakup constants (or equivalently 1 ) to a single parameter, ⟨ ̇ ⟩ ∕ , which is a non-dimensional stress. It represents the averaged stress experienced by each fat globule along its trajectory in the HPH. The mean residence time is the HPH volume divided by the volumetric flowrate, and ⟨ ̇ ⟩ is the average shear rate exerted by the flow on each globule along its trajectory, which depends upon the HPH geometry and operating pressure Δ . The larger ⟨ ̇ ⟩, and 1∕ (i.e., the larger ), the larger , and the smaller the outlet globule diameter ( or 1 ). If one admits that at constant (or flowrate), ⟨ ̇ ⟩ is a growing function of Δ , then the present model is consistent with the observed trends of with the operating pressure. Also, the former comparison of 1 value for the milk at 120 bars with the values for the cream at much lower pressure, is consistent with the effect of the effective viscosity predicted by the model: for a given

  at a pressure of 100 bars in the pilot HPH and three test cases in the industrial HPH, at 40, 70 and 115 bars, were done. In each case, several thousands of fluid particles of the emulsion, which is treated as a homogeneous fluid, initially homogeneously distributed in the gap flow entry, were individually tracked in the impact ring chamber, and ′ was obtained by averaging the shear rate experienced by all fluid particles along their trajectory. For each particle, ⟨ ̇ ⟩ = ⟨ √ 2 ∶ ⟩ (where is the local deformation rate tensor), is directly computed and ′ is derived from equation (6). For all these simulations, the flow viscosity was arbitrarily set to 1.16×10 -3 Pa.s, which corresponds to = ∕ =9.25 (this bulk viscosity corresponds to a laminar regime of the gap flow in the pilot HPH, close to the transition, and to a turbulent gap flow in the industrial HPH). The resulting CPD ( ′ ) is displayed in Fig.11(b) for all test cases. It exhibits a multimodal shape with a principal mode corresponding to more than 50% of the distribution. The averaged value is indicated by a round symbol on each curve. These averaged values of ′ are reported in Fig.11(a) for both pilot and industrial HPHs, represented by open symbols. Their order of magnitude is close to the values of and 1 obtained for the cream, respectively in the pilot and industrial HPH, but they tend to deviate from the trend observed for the cream, with smaller values as those predicted by the model (equation (

Figure 1 :

 1 Figure 1: Schematic of the HPH. =O(10 -10 2 ) is the gap thickness, 0 is the inlet velocity of the emulsion and = O(10 2 × ) is the width of the impact ring chamber.

Figure 2 :Figure 3 :

 23 Figure 2: Experimental size distribution in the pilot scale HPH at different pressures. (a) Probability density function (pdf ) ; (b) Cumulative probability distribution (CPD). Characteristic diameters ( 20 , , 90 ) and ( 20 , , 90 ) correspond to the inlet and outlet distributions at 15 bars, respectively.

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: Scheme of fat globule breakup in the pilot HPH

Figure 8 :Figure 9 :

 89 Figure 8: Scheme of fat globule break-up in the industrial scale HPH

Figure 10 :Figure 11 :

 1011 Figure 10: Size distribution of milk fat globules in the pilot scale HPH at 120 bars and 80°C. Comparison between model (dashed line) and experiments (symbols). (a) pdf (b) CPD ( 1 =2.22, 2 =10, Φ 2 = 36%)

Figure 12 :Figure 13 :

 1213 Figure12: Fragmentation constant as a function of (1 + 2Δ ∕ 0.4 ) 1∕3 for all systems in both HPH geometries. Closed round symbols: experimental data for cream =3.69 ( =2.9×10 -3 ). Open round symbols: simulation data with =9.25 ( =1.16×10 -3 ).Triangle symbol: milk with =1.69 ( =0.63×10 -3 ). Blue color: pilot scale (PS). Red color: industrial scale (IS)

Figure 2 :

 2 Figure 2: typical RANS simulation result (velocity magnitude in m/s) obtained on a 5 degrees angle sector. Mean gap size is around 90 microns and the imposed flow rate is 8800 L/h.

Figure 14 :

 14 Figure 14: left: Domain geometry and unstructured mesh grid in the industrial HPH. right: zoom of the mesh grid in the gap outlet region.

Figure 1 :

 1 Figure 1: LES velocity magnitude (m/s) snapshots obtained at pilot scale (gap=16 microns, flow rate: 200L/h) and associated streamlines. The geometry is three-dimensional, and the computational volume is meshed with a 1 degree angular sector as the geometry is axisymmetric.

Figure 15 :

 15 Figure 15: left: Snapshot of velocity magnitude field in the pilot HPH obtained from LES simulation. right: corresponding streamlines.

Figure 2 :

 2 Figure 2: typical RANS simulation result (velocity magnitude in m/s) obtained on a 5 degrees angle sector. Mean gap size is around 90 microns and the imposed flow rate is 8800 L/h.

Table 1

 1 Inlet section radius and gap length of pilot and industrial HPHs.

	=2.9×10 -3 Pa.s	(L/h)		(10 -2 bars -1∕2 )
	Pilot	200	730	3.22	1.06 0.94
	Industrial	8800	3200	0.964	1.02 0.98

Table 2

 2 Flow rate, gap flow Reynolds number, parameter in equation (28), parameters and in relation Δ = Δ for cream in pilot and industrial HPHs.

	=0.63×10 -3 Pa.s	(L/h)		(10 -2 bars -1∕2 )
	Pilot	200	3500	1.86	1.04 0.96

Table 3

 3 Flow rate, gap flow Reynolds number, parameter in equation (28), parameters and in relation Δ = Δ for milk in pilot HPH.

		Δ (	)	( m)	(m/s)
	Pilot	100		17	129
	Industrial	40		111	88
		70		82	118
		115		67	144

Table 4

 4 Pressure, gap average thickness and velocity for the simulated cases

Present address: Absolut Cooking,

ter route de Caussidières 31560 Nailloux, France O. Masbernat et al.: Preprint submitted to Elsevier

Making use of equations ( 20) and ( 22), the factor in equation ( 24) can be determined in laminar and turbulent regime:

2. Turbulent regime

with:

The total pressure drop in HPH is the sum of Δ and Δ :

From relations (25) or (26), can be computed in both HPH geometries. At a given operating pressure Δ , Δ can then be determined using (28). In all range of values, equation (28) can be replaced by a simple power law Δ = Δ , with and close to 1, as shown in Fig. 13. As a result, in the range of pressure investigated, the discharge pressure drop is identical to the pressure drop in the HPH. The values of gap Reynolds number, parameters (in equation( 28)), and are reported in table 2 for the cream in both HPHs, and in table 3 for milk in the pilot HPH.

B. Flow simulations in the HPHs

Computations were run using ANANAS TM which solves the incompressible balance equations for mass and momentum. ANANAS TM uses tetrahedral elements and is based on a mixed finite volume/finite element method. The domain geometry and the mesh grid used for the flow computation in the industrial HPH are shown in Fig. 14. Time integration is carried out using a third order explicit scheme, while space integration is handled with a high-order scheme introduced by ?, which yields to a sixth-order accuracy for an uniform mesh spacing. Depending on the case and Reynolds number value, the turbulence can be either modelled with a -RANS model as proposed by ? or thanks to Large Eddy Simulation Variational Multi-Scale model (LES-VMS) which enables the resolution of turbulent structures scaling down to the mesh size (?). Contrarily to common LES approaches, LES-VMS reduces spurious

Boundary layer mesh at the wall