Numerical simulation of evaporating charged sprays in spray chilling

Abstract In this research we investigate whether charged sprays offer any improvements over conventional spray in the spray-chilling of meat. A Eulerian-Lagrangian CFD modelling approach is developed and validated with results from automotive spray painting. The model is then extended to include evaporation and the transport of evaporation products. The simulation results show doubling of the spray transfer efficiency as well as a notable increase in the rate of cooling when electric charge is applied to the spray. Furthermore, the presented modelling approach yields consistent results at low computational costs, making it suitable for use in future design optimisation studies.


Introduction
Electrostatically charged spraying systems are a proven solution for improving the transfer efficiency of liquid or solid particle sprays [1,2]. Such systems are commonly used in applications where the spray material is expensive or overspray is undesirable, such as painting and pesticide application [1,3]. In the present research a new application for charged sprays is investigated, namely chill rooms in the meat industry. Spraying carcasses with water during the chilling process has two purposes; counteracting the drying effect that cold air has on warm meat, and providing additional evaporative cooling effort [4]. The main challenges with conventional sprayers in this application are ice accretion in the heat exchangers and hygiene issues from condensation and dripping, both related to excess water [5]. Electrostatic sprayers, with their higher transfer efficiency, are therefore an interesting prospective solution for the meat industry.
In this paper we propose a Computational Fluid Dynamics model that can be used to simulate the effects of charged spray in such a chill room. With this model we simulate two model cases as a proof of concept for the use of charged sprays in this application. Specifically, the potential improvement in the spray transfer and deposition pattern is investigated, as well as the effects of charged spray on evaporative cooling.
The simulation of charged sprays is not new in and of itself, dating back to the 1980's [6]. Various applications of electrostatic sprays have been simulated, with automotive spray-painting resembling our application the closest in terms of scale and spray volume. Most works in this field opt for a Eulerian-Lagrangian simulation approach, with either an unsteady [7] or steady [2,8,9] flow solution. The electric fields are solved in similar ways, and particle charge is either considered constant from the moment of injection [7,9] or variable with particles accumulating charge from the surrounding air [10].
However, evaporation of the droplets and heat transfer, which are important factors for our case, are not modelled in the previous paintspraying simulations. These effects are considered in some smallerscale works [11,12]. Such models are a significant step up in complexity compared to spray painting, both in the number of Eulerian phases as well as in the different interactions. In addition to the physics involved in uncharged droplet evaporation, the electric charge also affects droplet breakup. Moreover, evaporated droplets leave residual ions in the gas flow, which in turn interact with the electric field and the remaining particles.
Fully modelling all these mechanisms comes at the cost of increased computational costs. We aim to strike a balance between fidelity and practicality, while building a CFD model that is suitable for simulating evaporating, charged sprays in industrial scale applications. Our previous work [13], in which a simplified model was presented and validated against measurement data from literature, serves as a starting point to the present research.

Modelling approach
The purpose of the model is to simulate the spray chilling process, which involves spraying (charged) droplets at carcasses with a height of up to 2 m. Because of the requirement to model these large length scales, a fully resolved modelling approach is not feasible. Instead, the problem is divided into separate but interacting ''phases'', each of which is solved separately, with interactions represented by source terms. The main phases in this case are the spray droplets, the airflow, and the electrostatic field. In addition, the transport of heat, water vapour, and the residual charges from evaporated droplets are modelled as part of the fluid phase.
Both the validation case and the selected model case, see sections 3.1 and 4, are steady state problems. This suggests a RANS based approach for the fluid phase. Likewise, steady-state Eulerian models are used to model the transport of heat, vapour and residual charges. The droplet motion is simulated with quasi-steady Lagrangian particle tracking.

Droplets
For the Lagrangian droplet motion model the spray is treated as a sparse collection of (spherical) point-particles. Their motion is described by Newton's second law, see equation (1). The driving forces in this case are drag from the surrounding air ( F ! drag ), the electrostatic force ( F ! es ), and gravity ( F ! g ). Mutual aerodynamic interactions, electromagnetic effects and, droplet collisions are neglected, following the approaches of Shrimpton and Laoonual [11] and Arumugham et al. [12].
The following relations are used for the forces [12,14]: Note how the electrostatic interactions between droplets are computed indirectly, via the two-way coupling with the electrostatic field. Charged droplets contribute to the electric field as per equation (7), which states that the electric field has a positive divergence in regions of positive charge density. By the principle of superposition the total electric field can be seen as the sum of the electric fields of each individual charge. Therefore, the electrostatic force on each particle, as in equation (4), effectively includes the sum of all electrostatic particleparticle interactions.
To maintain conservation of momentum the drag force is coupled with the momentum equation for the fluid phase, where the reaction forces are included as momentum sources.
For the heat and mass transfer between the droplets and the continuous phase the correlations of Ranz and Marshall [15] are used, given in equations (5) and (6). Heat transfer inside the droplets is neglected, since the Biot number for our droplets is small (Bi≪ 1). Charge is not transferred between droplets and the continuous phase, except upon complete evaporation of a droplet.
Finally, the secondary breakup of droplets is neglected, considering the small droplet sizes and low Reynolds and Weber numbers involved. Contrary to Arumugham et al. [12] we also neglect droplet breakup due to electrostatic instability. The reasoning for this is the low initial droplet charge in the model-case simulations. By the time an evaporating droplet exceeds the stability limit, its mass has been reduced to less than 0.01% of its initial mass. The child droplets created by Rayleigh breakup would be even smaller. These small droplets would enter a rapid cascade of evaporation and breakup, simulation of which would take up significant computational resources. Instead, in the current work, droplets which reach the charge limit are considered evaporated, and their charge and mass are added to the continuous phase.

Electric field
Due to the steady-state nature of the simulation the electric field is assumed to be static, and can therefore be described by the Poisson equation (7).
Here ρ q denotes the volumetric charge density and ε 0 the permittivity of the vacuum (8:85 � 10 À 12 [F/m]). Both the charges carried by the droplets, as well as the residual charges contribute to this term, coupling the electric field with the discrete particles and the residual charge phase.

Gas phase
The turbulent gas flow in a steady state case can be described by the familiar 3D RANS equations, with a k-ε model for turbulence closure, following the examples set in previous literature [9,12,16]. The energy and the advection-diffusion equations govern the temperature of the gas phase and the transport of water vapour respectively. These equations and their solution methods are widely documented [14] and not the main topic of this research, and will therefore not be discussed in detail.
The residual charge left by completely evaporated droplets is treated as part of the continuous phase. Practically these residual charges will take the form of charged solid precipitates or molecular ions. Following the example of Arumugham-Achari [12] these charges are treated as mono-mobile. These charges are transported through advection, diffusion, and electrophoretic motion, as described by equation (8).
Here K denotes the electrical mobility of the residual charges and D q the diffusivity. The source term S q represents the charges being left by droplets evaporating, providing the coupling with the discrete phase.

Numerical methods
We selected the commercial Ansys Fluent code (version 18) as our CFD and DPM solver. The submodels for the electrostatic field, transport of residual charges and the electrostatic forces were appended to the base code using so-called User Defined Functions (UDFs).

Droplet tracking
Due to the steady-state nature of the simulated cases a quasi-steady Lagrangian particle tracking approach could be used. As with standard Lagrangian tracking, droplets are grouped into parcels, whose trajectories are computed by integrating their equations of motion forward in time. However, unlike with unsteady methods, each iteration of the Lagrangian model does not correspond to a certain limited time-step. Instead, each parcel's trajectory is computed from the point of injection up until the parcel exits the domain or disappears through evaporation. Along the resulting pathlines sources of momentum, heat and water vapour are created and handed back to the continuous phase solver for the next iteration. Similarly, the droplet charges passing through each element are recorded for the purpose of computing charge density, and a source of residual charge is created on droplet evaporation.
This quasi-steady method has a number of advantages over full unsteady tracking. First, it is easy to control the total number of parcels in the domain, leading to predictable computation times. In addition, the complete droplet trajectories are contained in a single iteration, aiding post-processing. A disadvantage is that parcels cannot persist in the domain between iterations, since their trajectories are integrated to completion. This makes modelling direct interactions impossible, and makes it necessary to model electrostatic forces via a separate electric field phase. In effect this means that the parcels of the next iteration will be repelled by the trajectories of the previous iteration, and care must be taken to prevent instabilities in the solution.

Electric field
The Poisson equation (7) governing the electric field can be solved by the same method as any advection-diffusion equation. It is treated as a Eulerian phase without advection, with diffusivity ε 0 and the charge density ρ q acting as a source term. It is solved in parallel with the other continuous phase equations. Dirichlet boundary conditions are used to represent conductors with a known electric potential. Fig. 1 illustrates the process flow of the solver. The solution cycle begins by running a number of iterations of the continuous phase solver to obtain a solution for the flow and electric field. Then, particles are injected and their trajectories integrated using the newest field solutions. Based on these trajectories the volumetric source terms for the continuous phases are updated, and the solution is finally checked for convergence to complete the cycle.

Solver sequence
The number of iterations of the CFD solver within each solver cycle can be varied. A large number of iterations is used on the first cycle, in order to have a converged field solution prior to the first DPM iteration. On later cycles this number is reduced to 2, to promote concurrent convergence of the discrete and continuous models.

Validation case
The simulation is validated against the results by Domnick, Scheibe and Ye [9]. They investigate the paint deposition pattern created by a rotary bell sprayer for automotive industry. Their work was chosen because it features both numerical as well as experimental data, thereby giving a better indication of accuracy than exclusively numerical studies. The validation setup and results are covered concisely, with a more comprehensive description having been presented in our preceding work [13].

Case setup
The basic geometry of the simulation consists of a rotary bell sprayer suspended above a 1 m square target plate, enclosed in a rectangular bounding box. A hybrid grid with 10 million elements is used, see Fig. 2. The sprayer assembly consists of two significant components; the stationary and electrically insulated mount (top) and the electrified rotating bell (bottom). On the bottom of the mount, surrounding the base of the bell, the so-called ''shaping air'' is injected through a 1 mm wide annular jet. This jet has a downward velocity of 75 m/s and a tangential velocity of 75 m/s opposite the rotation of the bell.
Droplets are injected along the outside perimeter of the bottom of the bell, with a velocity and direction equal to the tangential velocity of the bell perimeter. The size of the injected droplets follows the distribution shown in Fig. 3. The (initial) charge of the injected droplets is set at 5% of the Rayleigh charge limit Q R according to equation (9), where σ denotes surface tension (see Table 1) and r p the droplet radius. Because of the low volatility of paint, evaporation is not modelled by Domnick et al. [9]. For this reason the heat transfer, species transfer and residual charge transfer models are disabled in our simulation as well.
The target plate is treated as an electrically grounded wall, and absorbs the droplets that impact upon it. The top of the bounding box is treated as an inlet, creating a downwash with a velocity of 0.3 m/s. Symmetry conditions are applied to the sides of the bounding box, and the bottom of the bounding box acts as the outlet for the system. The key process parameters are summarised in Table 1.

Results
Figs. 4 and 5 show snapshots from the simulation results. To qualitatively present the droplets' behaviour three features are displayed. Fig. 4 shows a cross-section of droplet trajectories, coloured according to droplet diameter, as well as contours of the electric potential. The direction and closeness of the equipotential lines indicate the direction and magnitude of the electric field, which scales with the gradient of the potential. Fig. 5 displays the deposition rate of paint onto the target surface.
From Fig. 4 it becomes apparent that the electric field is strongest near the sprayer bell, and oriented outward. Towards the target plate the field weakens and becomes more vertically oriented. Looking at the droplet trajectories, a number of distinct regimes of behaviour can be seen. After injection at the bell edge, the droplets first follow a nearly horizontal inertial trajectory. The larger, heavier droplets are flung further from the bell edge before drag overcomes their initial velocity, while the smallest droplets are almost immediately captured by the downward airflow. After the initial straight flight, the trajectories curve down toward the target plate under the influence of the electrostatic force and gravity. The droplets follow the electric field lines closely, since the trajectories are roughly perpendicular to the equipotential lines. Finally, near the plate, the droplet trajectories bend radially outward, under the influence of the airflow resulting from the co-flow air jet impinging on the plate. The paint deposition pattern on the target plate (Fig. 5) shows two concentric rings of high deposition rate. The inner ring is the result of the smallest droplets following the airflow vertically downward, and being deposited directly underneath the bell. The exact mechanism responsible for creating the outer ring is unclear. Domnick et al. [9] do not discuss the presence of such a feature in their results, although their simulation data (see Fig. 6b) does show a bump in the corresponding location. Fig. 6a shows the accretion rate along the X-and Z-centrelines (see Fig. 5) of the target obtained from our simulation. Fig. 6b shows the results obtained by Domnick et al. [9], for comparison. For their experimental results they measured the thickness of the deposited paint layer along the middle of the target plate. This paint layer thickness should directly scale with the average deposition rate and spraying time, and depend further on the paint density and volatile fraction. Since spraying time and volatile fraction are not given by Domnick et al. [9], we qualitatively compare our deposition rate to their final paint thickness.
In general we find the same features in their results as in our simulation; a concave peaked profile with a local minimum at the centre. Moreover, in both their and our simulations some 'bulging' of the profile can be observed around 0.3 m from the centre, which is not apparent in the experimental data. The most notable difference between the results is that Domnick et al. [9] show a deeper local minimum around the target centre and a wider peak. The likely reason for this mismatch is minor differences between our and their geometry and air-jet boundary conditions.
Due to the overall similarity between the results we believe that our implementation of the electric field and charged droplet model is proven sufficiently accurate to proceed with simulations representative for the spray-chilling case.

Model case
The model case is designed to investigate whether charged spray offers any improvement over uncharged spray in a chill room. Three metrics in particular, namely the transfer efficiency of the spray, the uniformity of the spray deposition pattern, and the rate of cooling of the target, are considered.
The exact geometry of installed spraying setups varies between warehouses. Similarly, the shape of the spray target carcasses will depend on the species of animal, and the processing that has been done prior to chilling. As the goal is to study the behaviour of the spray at a general level, and in order to not introduce spurious details, we choose to define a simplified model geometry. This geometry consists of a cylinder with a radius of 10 cm, and a length of 1 m. The cylinder is enclosed in the centre of a cubic bounding box with a side of 1 m. We use a solid-cone spray, representative for an air-assisted atomiser, originating from the centre of the 'upstream' face of the bounding box. This geometry is shown in Fig. 7.  Fig. 4. Droplet trajectories and electric potential for the validation case.

Case setup
To study the differences between charged and uncharged sprays, two cases are simulated, differing only in the droplet charge. The remaining spray parameters and operating and boundary conditions are shared between cases, and explained below.
The parameters for the conical spray are given in Table 2. The sprayer body is represented by a 5 cm diameter disk centered on the upstream face of the bounding box. In the case with electrically charged droplets it is kept at a constant electric potential of 50 kV. Droplets are injected at the surface of this disk, distributed such that the mass flow rate and droplet size distributions are uniform across the injector. All droplets are injected at a velocity of 20 m/s, and with a direction vector such that the spray forms a uniform solid cone with an apex located 2.5 cm upstream of the centre of the disk.
The boundary face surrounding the sprayer is treated as an inlet, providing a 0.3 m/s co-flow. The co-flow air comes in at standard atmospheric pressure, a temperature of 293 K and 50% relative humidity. It does not carry any initial ''dissolved'' charge. The target cylinder is treated as a no-slip wall that absorbs incoming droplets. It is kept at a temperature of 293 K, and is considered a grounded conductor with a 0 Volt electric potential and absorbs residual charges. The boundary opposite the sprayer is treated as an outlet. All other boundaries are treated as planes of symmetry. Finally, the effects of gravity are neglected.
A hybrid grid with 115 � 10 3 elements (see Fig. 7) was created and used to calculate the continuous phase solutions. For the corresponding discrete particle solution 80,000 Lagrangian parcels were injected and tracked. To ensure grid-independence the simulations were repeated on a refined mesh (685 � 10 3 cells), using 320 � 10 3 Lagrangian parcels. The results obtained on the refined grid match those on the original grid, with less than 1% difference in the key metrics, confirming gridindependence.
The simulations were run on a consumer grade personal computer (4 cores, 2.8 GHz processor). Reaching convergence took a few hours when using the basic grid, and a day for the refined grid, with the computational time per iteration scaling with the number of elements as expected. The Lagrangian particle tracking was responsible for the majority of the computational effort. . 8a and b illustrate the differences between the behaviour of charged and uncharged droplets in a qualitative way. The uncharged droplets generally follow the airflow, forming a narrow plume which impinges on the centre of the target cylinder. Much of the spray plume flows around the cylinder, and is lost to the domain outlet. Conversely, the charged droplets form a wider plume due to their mutual repulsion, before curving back toward the target. A smaller but nevertheless significant portion does not get captured by the electric field and follows the airflow out of the domain. Fig. 9 shows the electric charge density and the electric potential   around the cylinder, in the mid-plane of the spray. The interaction between the droplets and the electric field is apparent by the large region of high potential near the injection, corresponding with the high charge density of the spray. The direction of the equipotential lines confirms the tendency of the charged spray to spread wider due to the mutual repulsion of the droplets. Another point of note is the large variation in the distance between equipotentials. The electric force drawing the droplets toward the target will significantly favour certain parts of the cylinder over others. Fig. 10a and b shows the rate of deposition of water on the target's surface. This data is plotted in cylindrical coordinates, with θ ¼ 0 ∘ corresponding to the spray-facing side of the cylinder. Integrating the deposition rate over the target surface we find total deposition rates of 0:92 � 10 À 3 kg/s and 1:89 � 10 À 3 kg/s for the uncharged and charged sprays respectively. The transfer efficiency of the spray, η tf , is defined as the ratio of deposited spray mass to injected spray mass. Computing this we find η tf ¼ 36:7% for the uncharged spray and η tf ¼ 75:5% in the charged case. Thus, the transfer efficiency is doubled by using charged droplets.

Water deposition
The deposition patterns follow the initial expectation, with the uncharged spray only collecting in a fairly small area. The charged spray has a weaker but still notable hotspot, and is more spread out vertically and horizontally. The charged spray also deposits more on the reverse side of the target, with a local minimum around θ ¼ 100 ∘ , corresponding to the area of low field strength shown in Fig. 9. Surprisingly, the uncharged spray also collects in a small region on the reverse of the cylinder (θ ¼ �180 ∘ ), likely due to a turbulent recirculation zone.
The deposition pattern for the charged spray does not appear fully symmetric. The amount of asymmetry is somewhat exaggerated by the logarithmically spaced isocontours, the major asymmetric features all have deposition rates less than 5 � 10 À 3 kg=m 2 s. Initially a numerical instability was suspected as the cause of this asymmetry, but it did not disappear or change when the simulation was allowed to run for a further hundred iterations. The actual cause could not be determined.
Based on these results we can conclude that the use of charged droplets improves the uniformity of the deposition pattern significantly. However, considering the electric field and the droplet pathlines it appears unlikely that a completely uniform distribution can be achieved through droplet charging alone. The single sprayer configuration used in the present model case appears to be sub-optimal in this regard. Further research is required to design the ideal sprayer setup, for which the present CFD model can be used.

Heat transfer
Since the temperature of the inlet air and the target cylinder are equal, all observed cooling can be attributed fully to the evaporation of the spray. The total convective cooling rate for the target cylinder is 48.2 W with uncharged spray, and 56.0 W with charged spray, as computed by integrating the surface heat flux over the cylinder surface. The increase in cooling rate is smaller than the increase in the transfer efficiency of the spray, η tf , discussed in section 4.2.1. This may be explained by examining the temperature along the vertical cross-section of the domain, shown in Fig. 11a and b. In both cases the evaporating spray cools the air surrounding it, resulting in a plume of colder air impinging on the cylinder. The temperature of the plume is equal in both cases, since it is limited by the amount of evaporation (and thus cooling) that can occur before the air is saturated. The main difference between the charged and uncharged case is that the wider dispersion of the charged spray creates a wider plume of cool air. This wider plume covers more of the cylinder's surface, and more heat can be convected away from it. However, the coverage difference between the narrow and wide impinging plumes is significantly smaller than the difference in liquid

Conclusions
In this work we developed a CFD modelling approach for evaporating electrically charged sprays. The results of the validation study show that the modelling approach yields results that are in line with previous literature. The results are consistent under grid refinement, while using a reasonable number of grid cells and Lagrangian parcels. While steady-state simulations of charged particles can lead to numerical instabilities, we find that these can be suppressed by applying appropriate under-relaxation factors in the interactions coupling the Eulerian and Lagrangian phases.
In terms of performance, the developed modelling approach strikes an economical balance between accuracy and computational costs. The model case simulations took several hours to run on personal computing level hardware (4 cores, 2.8 GHz processor), making them practical to use in design and optimisation studies.
With the developed CFD model we simulated a model case, representative for the application of charged sprays in a chill room. In this case the transfer efficiency of a charged spray is twice that of an uncharged spray. Additionally, the charged spray particles are distributed over a wider portion of the target. Further improvement of the spray distribution is left for a future optimisation study, which may be carried out using the developed CFD model.
Considering heat transfer, the use of charged spray results in a moderately (16%) increased rate of cooling of the target. This can be largely attributed to the wider spray cone, resulting from the mutual repulsion between spray droplets, allowing more evaporation to take place.
The improvements in the spray transfer and the increased evaporative cooling show that electrically charged sprays are a promising technology for the improvement of spray chilling systems in meat industry.