Description of a Eulerian–Lagrangian Approach for the Modeling of Cooling Water Droplets

: The present paper describes a tool developed in-house for the modeling of free-falling water droplet cooling processes. A two-way coupling model is employed to account for the interactions between the droplets and the carrier ﬂuid, following a Eulerian–Lagrangian approach. In addition, a stochastic separated ﬂow technique is employed, involving random sampling of the ﬂuctuating ﬂuid velocity. In physical modeling, two empirical correlations are considered for determining the heat and mass transfer coefﬁcients, with the possibility of accounting for vibrations. The numerical results indicate the preponderance of the interactions between droplet and carrier ﬂuid at various humidity ratios.


Introduction
The impact of cooled water droplets and their accretion in lifting surfaces and engine intakes is of paramount importance for safe aircraft operation [1][2][3]. Other areas where cooling phenomena play a preponderant role include refrigeration in the food industry [4], the performance of wind turbines [5] and generally adverse effects related to infrastructures [6].
In order to gain insight into the physical aspects regarding cooling and freezing mechanisms, the present paper describes a Eulerian-Lagrangian approach to modeling the cooling rate of free-falling water droplets. Particles are described following a Lagrangian reference frame, while a finite volume method (FVM) is used for the carrier's fluid description, corresponding to a Eulerian frame. The interactions between carrier fluid and particles are described following a two-way coupling model, accounting for the increase in particle mass-load ratio.
Eulerian-Lagrangian formulations have seen widespread use in the literature for the modeling of turbulent particle dispersion in dilute flows [7]. There are many reviews dealing with either physical or mathematical aspects of the approach. While Gouesbet and Berlemont [8] and Patankar and Joseph [9] focused on particle-laden and particulate flows, respectively, Subramaniam [10] deals with numerical implementations of Eulerian-Lagrangian approaches for the description of multiphase flows.
In order to study the ice shape and aerodynamic performance of a NACA 0012 airfoil with a heater on the leading edge, Uranai et al. [11] considered the flow field following a Eulerian method, with the droplets' trajectories computed in a Lagrangian frame. A contribution to a better understanding of the thin ice layer formed behind the heater due to the well-known run-back ice phenomenon was achieved by proving that a modified extended Messinger model based on the work of Messinger [12] was more suitable for de-icing simulations for both rime and glaze ice conditions. Similarly, a computational method based on statistical theory considering a Euler-Lagrange framework was used by Peng et al. [13], in which the effects of turbulent dispersion on water droplet impingement in a NACA 0012 airfoil were considered. This study enhances the predominance of turbulent dispersion in particle-fluid interaction time scale, and is usually used in these problems. Nevertheless, other dispersion models in particle-laden flows exist [27].
Eulerian-Lagrangian formulations are commonly encountered in numerical studies dealing with the study of evaporating droplets related to engine fuel consumption [28][29][30]. A three-dimensional computational tool to study turbulent two-phase flows was implemented and tested by Barata [31] and Rodrigues et al. [32] to model the dispersion and evaporation of fuel droplets in two-phase turbulent jets in a cross-flow, characteristic of the conditions encountered in the injector and combustion chambers of aero-propulsive systems. Retaining the original numerical framework, the present work is developed by recasting the thermo-physical formulation from dealing with different fuels to water and the evaporation processes into the description of cooling phenomena. The implemented models, thermodynamic relations and properties of the relevant fluids are discussed throughout the manuscript. Two distinct phases are considered, formulated in a Eulerian-Lagrangian approach, with mass, momentum and energy exchanges. Evidence is found on the role of vibrations and deformations for low relative humidity ratios while at high humidity ratios the impact of these phenomena is lower, suggesting the presence of a transition criterion set at intermediate humidity levels closely related to increasing droplet diameter.
The remainder of the manuscript is organized as follows: first, the governing equations for the gas and droplet fields are described, followed by the thermodynamic relations of the cooling process. Then, the coupling between the considered phases is described, for which the corresponding numerical algorithm is presented. Next, the numerical tests and method validation are addressed, attending first to the influence on the cooling process of a variable droplet diameter considering a constant humidity ratio and then to a variable humidity ratio and a constant droplet diameter. Lastly, conclusions are drawn, taking into account the advantages and limitations of the proposed methodology.

Governing Equations
The RANS equations for mass, momentum and energy are considered for an incompressible flow. These equations are written for a stationary, viscous and Newtonian fluid under a conservative form to include source terms that consider mass, momentum and energy exchanges. In addition to these mean flow properties, the mass fraction of the flow is also considered. Reynolds stresses, U i U j , are approximated through the eddy viscosity concept and Boussinesq's hypotheses. Finally, the system of partial differential equations (PDEs) is closed with the κ-turbulence model [33] for the determination of the turbulence kinetic energy and its dissipation.
The governing equations can be reduced into a single advective-diffusive conservation equation for a general property per unit volume, φ, represented in Equation (1), where Γ φ and S φ correspond to the effective diffusion coefficient and the source term for each property, φ, respectively.
The trajectories of representative samples are used to obtain the position and velocity of the droplets. This is achieved by solving the particle momentum equation, according to Newton's second law of motion via the Eulerian fluid velocity field, which is consequently obtained by solving the RANS equations. Finally, the particle position is represented by an equation that considers both the previous and current locations, the particle velocity and the time it took to travel that distance.
A mathematical expression for the particle motion equation, considering a spherical droplet that accounts for all four types of forces, is suggested by Shirolkar et al. [7]. Nevertheless, in the model employed for cooling free-falling water droplets, only steady-state mechanisms for particle dispersion are considered, resulting in a simplified equation for a three-dimensional Cartesian coordinates system. Equations (2) and (3) respectively represent the new particle position and velocity, with subscript d denoting a droplet property and g a gas field property. The trajectory equations depend directly on three additional parameters: the particle relaxation time, d ; the integration time step, ∆t; and the gravitational acceleration, g. The superscripts NEW and OLD correspond to a certain property about the next step and one related to the current one. The parameter X i and U i are the generalized position and velocity coordinates.
The particle relaxation time is defined as the response rate of particle acceleration to the relative velocity between the particle and the carrier fluid, depending on the particle's inertia and free-fall velocity. This parameter is represented in Equation (4). The parameter d is the droplet diameter, ρ the density, µ the dynamic viscosity and Re d the diameter-based Reynolds number.
The droplet drag coefficient, C D , is represented in Equation (5), estimated according to experimental data fitting [34].
The remaining two parameters required to solve Equations (2) and (3) are the integration time step and the instantaneous fluid velocity at the particle location differing from the averaged value. For that purpose, a stochastic separated flow (SSF) technique is employed, involving random sampling of the fluctuating fluid velocity from certain known distributions. The SSF technique estimates the fluctuating component of fluid velocity, considering a concept which states that the particle interacts with a succession of eddies as it moves along the computational domain. The fluctuating fluid velocity corresponding to a particular eddy is randomly sampled from a probability density function (PDF) obtained from local turbulence properties. This fluctuating fluid velocity PDF at each particle location is assumed to be a Gaussian, with zero mean and standard deviation equal to √ 2k/3. The eddy size and lifetime scales needed at each particle location to determine the next interaction time are obtained, resorting to local turbulence properties (κ and ) along the particle trajectory Kolmogorov time and length scales for isotropic flows. However, when using isotropic SSF models, it is possible to use different expressions for these two turbulence scales [7]. Equations (6) and (7) respectively represent a general expression for the turbulence eddy lifetime time scale, τ e , and length scale, l e . The parameters a and b are two dependent constants, determined either by scaling analysis or from experimental data [7].
In addition, the crossing trajectory effect must also be considered, being translated mathematically by the eddy transit time, υ e , represented in Equation (8). This parameter is interpreted [26] as the minimum time a particle would take to cross an eddy with a characteristic dimension, l e . It becomes necessary to account for this parameter since a particle may remain trapped inside an eddy for the entire lifetime of that eddy, or may prematurely migrate from one eddy to another before the decay due to the turbulence of the original eddy. This premature migration is usually related to the significant free-fall velocity of the particle under consideration. Thus, if the minimum crossing time calculated is smaller than the eddy lifetime, the particle would jump to another eddy.
Looking at Equation (8), it becomes clear that there is no solution when the characteristic eddy size is greater than the fluid-particle relative velocity multiplied by the relaxation time. This fact is a consequence of the linearized stopping distance of the particle being smaller than the eddy size, causing the particle to be trapped by the eddy. In this case, the interaction time will be the eddy lifetime. Therefore, the particle will interact with a particular eddy for a time, which is the minimum between the eddy lifetime and the eddy transit time. The fluctuating velocity associated with a particular eddy is assumed to be constant over the interaction time. The time step selected is the eddy-particle interaction time, t int , given by Equation (9).
The fluid velocity sampling and the interaction time make it possible to obtain the particle trajectory from Equations (2) and (3). At the end of each step, a new fluctuating fluid velocity is sampled from a new PDF, generated through local turbulence properties. Consequently, the following interaction time is determined by the local properties at the new particle location.
Concerning the source term, its contribution is divided into a component accounting for the interactions related to the gas field, S φ,g , and the other for the interactions of the droplet field, S φ,d . Equation (10) clarifies the terms that constitute the general source term. Additionally, the contribution of the dispersed phase can also be divided into two contributions calculated for each Eulerian cell of the continuous phase. S φ,ipt specifies the source term due to inter-phase transport and S φ,c accounts for the transfer phenomena caused by the cooling phenomenon. Table 1 summarizes the system of PDEs, where the subscripts t, a and wv respectively indicate properties related to turbulence, dry air and water vapor. Sc d and Pr d are the Schmidt and Prandtl numbers. Υ is the mass fraction, T the temperature, p the pressure, V the volume, m the mass, N the number of particles, Q rel the relative rate of heat transfer, L ph the latent heat of phase change and c p the specific heat capacity at constant pressure. The turbulence model constants are defined as C µ = 0.09, C 1 = 1.44, C 2 = 1.92, C 3 = 1.1, σ κ = 1.0 and σ = 1.3. Table 1. General representation of the system of PDEs, adapted from Gouesbet and Berlemont [8], Sommerfeld [30].

Cooling Process
The cooling process of water droplets has been studied by different authors over the years [35][36][37][38][39][40][41]. In the present work, a basic heat balance model is used to obtain the rate at which the temperature of a droplet varies over time, assuming that the particle integrates a thermodynamic system in which heat transfer takes place from the droplet to the surrounding environment and vice versa. The heat transfer rate is given by Equation (11), assuming that the droplets' internal motion is so vigorous that complete mixing is achieved. The effects of vibrations and deformations on the particle are neglected. The surface area is represented by A and the parameters q h , q m and q r respectively denote the heat fluxes (per unit area) due to convective heat and mass transfer, and thermal radiation. h h is the convective heat transfer coefficient, h m the convective mass transfer coefficient, e the emissivity and σ the Stefan-Boltzmann constant of radiation.
In terms of the convention used, the heat fluxes due to the convective heat transfer and thermal radiation have an outward flux to the surrounding environment. However, since we are undergoing a cooling process, the heat flux due to convective mass transfer has an inwards movement. As a result, a diffusive-convective region occurs around the droplet due to significant temperature variations between the droplet and the far-field region of the gas. The temperature in this area is determined by a combination of the droplet and gas temperatures, and is referred to as a reference temperature.
The values of the convective heat and mass transfer coefficients are usually obtained by considering correlations based on empirical data. When considering the steady-state heat transfer of solid spheres, the Ranz-Marshall correlations [42,43] Additionally, we consider the Ranz-Marshall corrected formulation [44] of Equations (13), which introduces the effects of vibrations and deformations into the computations, whose range of validity is given according to Equation (14).
3 mm ≤ d ≤ 6 mm; 10 ≤ z d ≤ 600 The thermodynamic relations needed to obtain the previous relations as well as the thermophysical properties of dry air and water are presented, respectively, in Appendices A and B.

Coupling of Continuous and Dispersed Phases
To obtain a converged solution for both phases, it is necessary to couple them. For that purpose, the iterative procedure described in Algorithm 1 is followed. In order to reach convergence, a comparison between the source terms of the previous and current iteration is necessary.
The velocities are evaluated at the cell's edge in a staggered grid configuration for pressure-velocity coupling, discretized following the Quadratic Interpolation for Convective Kinematics (QUICk) scheme [45]. The Semi Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm [46] is used for the pressure-velocity coupling. Lastly, the system of equations is solved through the Tridiagonal Matrix Algorithm (TDMA).

Algorithm 1:
Iterative procedure for the two-way coupling of the dispersed and continuous phases.
Input: Initial conditions of the spray and surrounding environment as well as computational domain dimensions; (1) Initialization of the source terms; (2) Mesh generation; (3) Initialization of the continuous phase without the dispersed phase source terms (source terms of the discrete phase equal to zero); (4) Initialization of the dispersed phase. The particles are tracked through the flow field in this phase, and their source terms are calculated; (5) Recalculate the continuous phase using the new source terms from the dispersed phase; while convergence criterion not satisfied do Recalculate the dispersed phase using the new source terms from the continuous phase; Recalculate the continuous phase using the new source terms from the dispersed phase; end (6) Recalculate the dispersed phase one last time using the source terms obtained from the continuous phase in the last run; Output: Temperature variation of water droplets along with the computational domain.

Initial and Boundary Conditions
We used the experimental conditions of Yao and Schrock [44] to validate the numerical predictions, and they are summarized in Table 2.
The experimental setup consisted of a plastic column of a squared base with a 0.170 m edge, 3 m height and a droplet generator. The droplet generator was set on the center of the upper face and, an air stream was injected downwards from the upper into the lower face with a constant velocity of 0.03 m s −1 .
In this way, four distinct droplet diameters (3, 4, 5 and 6 mm) and four different humidity ratios (0.29, 0.36, 0.52 and 1.00) were considered. For each droplet diameter, three distinct humidity ratios were used to represent a low, medium and high humidity content. The computational domain is described in Figure 1. A symmetry plane was considered at x = 0, resulting in only half the domain being considered. As seen in the figure, the Cartesian reference frame has its origin on the upper face. An inlet was considered on the top face, where the variables were specified with uniform profiles according to the experimental conditions. In the lower plane, an outlet was considered where the normal gradients of the dependent variables were set to zero, and the remaining boundaries were defined as no-slip walls. Lastly, the initial integration interval was of 5 × 10 −5 s.

Numerical Tests and Validation
To evaluate the cooling rate over time, two distinct parameters were monitored: the droplet temperature, T d , and its falling height, z. These parameters are presented in a non-dimensional way to enable a fair comparison of the phenomenon for different conditions. The droplet non-dimensional falling height was evaluated in relation to its diameter, while a temperature ratio ((T d − T d,0 )/ T g,∞ − T d,0 ) was considered for the cooling process.
However, to ensure that the numerical results were not affected by the solver's numerics, a grid independence analysis was carried out using three levels of refinement in a structured orthogonal mesh of rectangular elements, with 9660 (coarse mesh), 19, 200 (fine mesh) and 38, 000 (finest mesh) points. Figure 2 shows a comparison between these three grids for a droplet diameter of 3 mm and a relative humidity ratio of 0.29. The falling height was represented as a function of the vertical velocity component. Since no significant deviation was observed, the fine mesh was used in the computations.  Figure 2. Grid independence study for three distinct grids. Figure 3 depicts the temperature variation of free-falling water droplets with diameters of 3 and 5 mm and a humidity ratio of 0.29. In the figure, the Ranz-Marshall correlations (RM) are compared with the corrected formulations (RMcf) and validated with experimental measurements [44]. Differences between both correlations are evident in the results. While the complete mixing model, corresponding to the Ranz-Marshall correlations, led to an under-prediction of the droplets' cooling rates, the introduction of vibrations and deformations (RMcf) led to a closer agreement with the experimental data up to a distance of 80 diameters. Moreover, by including the effects of vibrations and deformations, the initial curvature of the temperature profile was retrieved. Nevertheless, as the distance increased, the experimental data behavior and the numerical predictions started to diverge.
Upon increasing the humidity ratio from 0.29 to 0.36, the cooling of water droplets with diameters of 3, 4, 5 and 6 mm was observed (Figure 4). A similar pattern was discovered concerning the previous case. While a linear evolution of the cooling rate was observed, when considering the Ranz-Marshall correlation, the corrected formulation led to a close agreement with the experimental data. Given the results, it is possible to observe that the numerical predictions (RMcf) followed the experimental data for all the considered diameters up to a distance of approximately 200 diameters, after which results started to diverge. Comparing the results for a humidity ratio of 0.29 in Figure 3 with the ones of Figure 4 for a humidity ratio of 0.36, it is possible to infer the effect of vibrations and deformations for the prediction of droplet cooling rate.    Figure 5 depicts the results for a humidity ratio of 0.52 and variable droplet diameters of 4 and 6 mm. A distinct behavior was observed concerning the previous results. Overall, the corrected Ranz-Marshall formulation slightly over-predicted the droplet cooling for small distances. In the case of a 4 mm diameter, a generally accurate representation of the experimental data was observed, while for the 6 mm case an accentuated over-prediction was evident up to a falling distance of 100 diameters. Concerning the results obtained with the RM correlation, the predicted linear growth approached the experimental data at around a distance of 400 diameters in the 6 mm droplet. Lastly, considering the humidity ratio of 1.00 and droplet diameters of 3, 4, 5 and 6 mm, an almost linear growth of the experimental cooling rate for all droplet diameters considered is evident in Figure 6. Given the high humidity ratio considered and the numerical results obtained with the Ranz-Marshall corrected correlation (RMcf), it is possible to conclude that the influence of vibrations in the cooling rate under these conditions was meager. However, the linear variation retrieved from the RM correlation led to a linear growth that was sharper than the experimental observation. Nevertheless, a good agreement with the experimental data was evident up to a distance of around 200 diameters for the 3, 4 and 5 mm droplets. The cooling rate over-prediction in the case of the 6 mm droplet started at around 100 diameters.

Concluding Remarks
In this paper, our computational tool developed in-house for the prediction of the cooling rate of free-falling water droplets was successfully implemented and compared with experimental data. This tool was a RANS-based two-way coupling model using a Eulerian-Lagrangian approach to model each of the considered phases.
The two-way coupling model, which considers the interactions between droplets and the surrounding environment, showed an overall close agreement with experimental data. Furthermore, the consideration of vibrations and deformations through the introduction of an appropriate empirical correlation allowed for the determination of initial curvature of the cooling rate, while the impact of vibrations was found to be meager at high humidity ratios. This suggests the presence of a transition criterion set at intermediate humidity levels that is closely related to increasing droplet diameter.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analysis and interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Thermodynamic Relations
Equations (A1)-(A4) represent, respectively, the reference specific density, the dynamic viscosity, the specific heat at constant pressure and the conductive heat transfer coefficient [47], considering the process of droplet evaporation. The subscripts a and wv respectively correspond to dry air and water vapor. The subscript re f indicates that the parameter was calculated using a reference temperature. This temperature combines the droplet temperature and the ambient gas temperature, aiming to represent a near-droplet region in which the temperature of the droplet has a great influence on the surrounding air (diffusive-convective region). The reference temperature is given by Equation (A5).
The water vapor-gas mass diffusivity needed for the determination of the Schmidt number, considering a temperature range between 223.15 K and 293.15 K, is given in Equation (A6) [48]. The relation between the dry air mass fraction at reference conditions and the water vapor mass fraction at reference conditions is given by Equation (A7).
The mass ratio between the water vapor and air (dry air plus water vapor) is known as specific humidity, SH (Equation (A8)), considering an approximation with the humidity ratio or mixing ratio (the amount of water vapor in a certain volume of air may be defined as the ratio between the water vapor mass and the total mass). This parameter is, by definition, the water vapor mass fraction (Equation (A8)). SH = m wv m wv + m a HR 1 + HR (A8)

Appendix B. Thermophysical Properties of Dry Air and Water
The thermodynamic and transport properties of dry air at atmospheric pressure for a temperature range between 200 K and 400 K were evaluated according to Table A1. A non-linear squares algorithm combining a Gaussian method, also known as Taylor series, and the method of steepest descent, was used. It is possible to conclude that all correlations deviated by less than 0.15 % from the tabulated values [49]. For liquid water, the thermophysical properties were obtained by resorting to a fourthdegree polynomial interpolation of data available from the National Institute of Standards (NIST) database [50]. These interpolations were evaluated according to Table A2 and are valid for a temperature range between 273.16 K and 335.16 K at 1 atm.