An improved lumped model for freezing of a freely suspended supercooled water droplet in air stream

This work deals with the mathematical modeling of the transient freezing process of a supercooled water droplet in a cold air stream. The aim is to develop a simple yet accurate lumped-differential model for the energy balance for a freely suspended water droplet undergoing solidification, that allows for cost effective computations of the temperatures and freezing front evolution along the whole process. The complete freezing process was described by four distinct stages, namely, supercooling, recalescence, solidification, and cooling. At each stage, the Coupled Integral Equations Approach (CIEA) is employed, which reduces the partial differential equation for the temperature distribution within the spherical droplet into coupled ordinary differential equations for dimensionless boundary temperatures and the moving interface position. The resulting lumped-differential model is expected to offer improved accuracy with respect to the classical lumped system analysis, since boundary conditions are accounted for in the averaging process through Hermite approximations for integrals. The results of the CIEA were verified using a recently advanced accurate hybrid numerical-analytical solution through the Generalized Integral Transform Technique (GITT), for the full partial differential formulation, and comparisons with numerical and experimental results from the literature. After verification and validation of the proposed model, a parametric analysis is implemented, for different conditions of airflow velocity and droplet radius, which lead to variations in the Biot numbers that allow to inspect for their influence on the accuracy of the improved lumped-differential formulation.


Introduction
Problems involving droplets solidification find application in diverse fields such as aerospace and aeronautics, electric power transmission, meteorology, refrigeration, and cryopreservation. Besides, the freezing time of supercooled water droplets is an especially important parameter in the study of surface coatings to prevent ice and frost formation on cold solid surfaces. The aeronautical industry is particularly affected due to the extreme environmental conditions in which aircraft operate nowadays. In this sense, the mathematical modeling of the freezing mechanism is particularly useful for proper design and satisfactory performance of "icephobic" surface modifications for aeronautical sensors and components [1][2][3]. Freezing of supercooled droplets, either suspended (or flowing) in a gaseous environment or brought in contact with a cold substrate, has been studied before, and some of these previous works are here reviewed [4][5][6][7][8][9][10][11][12][13]. The freezing process of a supercooled droplet can be described in four distinct stages. Hindmarsh et al. [4] defined these stages as: (1) A supercooling stage, during which the liquid droplet is cooled from an initial temperature to below the equilibrium freezing temperature until ice nucleation occurs; (2) A recalescence stage, during which supercooling drives rapid kinetic crystal growth from the crystal nuclei. This stage results in some heat release due to partial freezing of the droplet and its return to the equilibrium freezing temperature; (3) The solidification stage, when freezing growth is governed by the heat transfer rate from the droplet to the point where the droplet liquid is completely frozen; (4) cooling or tempering stage, when the solid droplet temperature is reduced to the ambient air temperature.
A number of recent research efforts dealt with experimental investigations for supercooled droplet freezing, and concentrated on developing theoretical models for this phenomenon and applying semi-analytical or numerical techniques. For instance, Ruberto et al. [5] experimentally investigated freely suspended supercooled water droplets by using a levitation technique, when a single droplet was trapped in a test chamber. The authors investigated droplet sizes around 50 μm so as to approach the magnitude of droplet diameters appearing in clouds. The influence of the relative humidity on the evaporation of supercooled water droplets was systematically investigated. The authors observed in their experiment a linear relation between the evaporation rate and the relative humidity. Ruberto et al. [6] extended their previous work and carried out a similar experimental study, but now comparing it to the numerical solution of a proposed theoretical model. The numerical computations were performed with the in-house Free Surface Code (FS3D), which is a DNS simulator based on the volume of fluid (VOF) method to solve the incompressible Navier-Stokes equations. The simulation results for the influence of the relative humidity on the evaporation of supercooled water droplets were in good agreement with the experiments for all three temperatures considered. In Hindmarsh et al. [4], the solidification step was solved by considering it as Stefan's two-phase problem. Two different models were considered in relation to the internal temperature distribution within the droplets, one considering the transient one-dimensional heat conduction, with moving boundary in the solidification stage, and the other considering a lumped formulation in which the temperature of the droplets was assumed to be spatially uniform, the authors also carried out an experimental study, considering a thermocouple inside the droplet, analyzing this temperature variation and comparing the results with the numerical model. Among purely numerical works, Feuillebois et al. [7] analyzed the freezing of liquid droplets exposed to cold air and subjected to supercooling, based on Stefan's one-phase problem formulation. The perturbation method was used to obtain the evolution of the freezing front as a function of time, keeping the droplet surface temperature constant and equal to the external environment temperature. The solution obtained through the perturbation method was compared to the solution through a numerical method, showing good agreement between them except in the region close to the center of the droplet, where the perturbation method deviated from the numerical approach. Tabakova et al. [8] extended the work in [7] by considering heat convection at the droplet surface during the solidification stage of supercooled water droplets also as a one-dimensional one-phase Stefan problem and employing a perturbation method. The authors then suggested explicit correlations to estimate the freezing time based on their numerical results. Zhang et al. [9] experimentally and numerically analyzed the freezing process of supercooled water droplets on cold plates. Through image recognition techniques, the authors were able to inspect the behavior of the solidification front and estimate the freezing time of the droplet on both hydrophilic and hydrophobic surfaces. In their simulation, a VOF multiphase model coupled with the solidification/melting model in Fluent 14.0 was used to investigate the heat and mass transfer process during the droplet freezing. The evolution of the freezing front calculated by the proposed model agreed reasonably well with the experimental observation. In addition, through the average values of the freezing times obtained by the simulation, the authors developed a correlation to predict the freezing time of supercooled droplets, which agreed with 90% of the simulation data and all the experiments within a deviation margin of 25%. Chaudhary and Li [10] numerically and experimentally analyzed the four stages present in the freezing process of a water droplet (supercooling, recalescence, freezing, cooling) on surfaces with different wettability. The temperature evolution of the droplets was recorded using both intrusive and non-intrusive methods to identify the processes involved in the cooling and phase change within the droplets. The proposed model was written in terms of the enthalpy formulation. The numerical results of the freezing droplets temperature evolution are compared to the experimental data, showing close agreement with the experimental freezing times. Sultana et al. [11] numerically examined phase change of free-falling droplets in a sub-zero environment for droplets of fresh and salt water. The model was based on the solution of the Navier-Stokes equations coupled with the VOF methodology for tracking the droplet-air interface. They also analyzed the nucleation temperature for droplets of different sizes, concluding that large-sized droplets had higher nucleation temperature than the smaller one and the temperatures for the salty water droplet were always lower than the fresh water ones. In [12], the integral balance method (IBM) was employed in combination with a VOF multiphase model for tracking the interface position, and employed in the simulation of both freezing liquid films and droplets.
The employment of purely numerical approaches has allowed for the more computationally involved simulation of such moving boundary problems governed by an increased number of parameters, but has also confirmed the higher computational costs required for an error-controlled solution. In this context, a robust and cost-effective hybrid numerical-analytical approach known as the Generalized Integral Transform Technique -GITT [14][15][16] has been advanced in the analysis of supercooled water droplets freezing [13], based on previous hybrid implementations for such class of moving boundary heat transfer problems [17][18][19]. Specifically, Carvalho et al. [13] employed the GITT to accurately solve for the full freezing process, in all stages above described, following a more recent solution alternative through the adoption of a nonlinear eigenvalue problem for the expansion base [20], due to the associated nonlinear boundary conditions, leading to faster and more uniform convergence of the temperature distributions in the freezing droplet. This solution provides a set of reference results for this nonlinear partial differential system, with moderate computational costs.
Nevertheless, simplified reduced models are expected to be particularly useful in reducing computational costs and analytical involvement, especially in connection with very computationally intensive tasks, such as in optimization and inverse problem analysis, when the direct problem must be solved many times, or when populational dynamics analysis and stochastic simulations are undertaken, requiring a large number of samples to be simulated. Therefore, an improved lumping procedure, based on the so-called Coupled Integral Equations Approach (CIEA) [16,[21][22][23][24][25] is here employed as a formulation simplification technique for the present heat conduction problem with moving boundaries. The resulting improved lumped-differential formulation offers substantial enhancement over classical lumping schemes in terms of accuracy, without introducing additional complexity in the corresponding final simplified differential equations to be handled. The CIEA formalism approximates integrals of the temperature and heat flux profiles by a linear combination of the integrand and its derivatives at the integration limits, an idea originally developed by Hermite [26] and employed by Cotta et al. [27] in approximately solving moving boundary problems. This problem reformulation strategy has been applied to various thermal sciences and engineering problems such as in fin analysis, conjugated problems, drying, channel flow, aerospace thermal protection system, membrane metals extraction, nuclear fuel rods, heat exchangers, micro-reactors for biodiesel synthesis, nanocomposites, among others, as recently reviewed in [25]. The GITT is a hybrid numerical-analytical methodology for solving the full partial differential system for a distributed parameters formulation. As mentioned before, it has been advanced for the present problem in [13] and is particularly suitable for benchmarking purposes, due to the automatic global accuracy control that is inherent to its hybrid nature. On the other hand, the CIEA is not a solution methodology, but rather a problem reformulation tool, which allows for partial lumping in one or more spatial coordinates, thus providing model reduction alternatives to the classical lumped system analysis. It is a straightforward methodology for model reduction and these two approaches can even be employed in combination, as illustrated in [24].
The present work aims at advancing a lumped-differential reduced model for freezing of a supercooled droplet suspended in a cold air stream and subject to the three main transport phenomena at the interface between the droplet and the surroundings: convective heat transfer, convective mass transfer, and thermal radiation. The CIEA is employed to transform the partial differential equations (PDE) formulation of the energy balance, into an ordinary differential equations (ODEs) system for the boundary temperatures at each stage of the freezing process, and for finding the position of the moving freezing interface during the solidification stage. The improved reduced model is expected to extend the range of applicability of the lumping approach, in comparison to the classical lumped system analysis, in terms of the main governing parameters. The resulting nonlinear ODEs model is solved using the Mathematica ® platform, Wolfram [28]. The transient boundary temperature distributions for each stage of the process are obtained and analyzed in different scenarios. In addition, a critical analysis of the accuracy of the results achieved via CIEA is undertaken through a comparison with the reference results obtained via GITT [13] and against numerical/experimental results from the literature.

Problem formulation
In formulating the energy balance for the freezing of a suspended water droplet, the following assumptions are made: (i) The droplet is suspended in air, subject to forced convection; (ii) The droplet keeps the same volume and spherical shape all along the process; (iii) Heat transfer is assumed to be one-dimensional in the radial direction; (iv) Ice and water are isotropic and homogeneous, with constant properties; (v) Density changes at the liquid/ice interface are disregarded; (vi) In the recalescence stage, the droplet temperature will be considered uniform and equal to the equilibrium temperature for freezing (T f ); (vii) In the solidification stage, the liquid phase temperature is considered constant and equal to T f , thus leading to a one-phase Stefan problem. Then, the mathematical models for each of the four stages are written as summarized below and described in further details in the supplementary information of [13].

Supercooling (1st) and cooling (4th) stages
The supercooling stage model involves the transient one-dimensional heat conduction equation in spherical coordinates for a fixed domain, 0 < r < R, with nonlinear boundary condition that accounts for convective heat transfer, convective mass transfer, and radiative heat transfer, as in [13]. This dimensional formulation is rewritten in dimensionless form and transformed to Cartesian coordinates through a variable transformation [29], employing the following dimensionless parameters: where, in Eqs, (1) where where x is the dimensionless spacae variable, T is the dimensionless time variable, and θ * (x, t) is the dimensionless temperature in cartesian coordinates. Eq. (2), τ 1 is the dimensionless time when the droplet reaches the nucleation temperature, and in Eqs. (6)- (7), Bi c,1 , Bi m,1 , Bi r,1 , are, respectively, the characteristic Biot numbers for convective heat transfer, for mass transfer and for radiative heat transfer. Bi c,1 in particular represents a measure of the ratio of the convective and conductive heat fluxes at the droplet surface and it is an important governing parameter for this application; furthermore, its analysis is important for the application of CIEA, and this aspect is discussed in Sect. 4. The other Biot numbers similarly compare other transport mechanisms with respect to conductive (diffusive) transport. The cooling stage model is similar to the model presented above, but in this case the thermophysical properties of the liquid must be replaced by the solid phase properties. Furthermore, the initial temperature needs to be changed to the spatially varying temperature obtained at the end of the third stage and the latent heat of evaporation should be replaced by the latent heat of sublimation. The correlations for water vapor density for liquid and solid droplets, used in modelling supercooling and cooling stages, respectively, are presented in Sect. 2.4.

Recalescence (2nd) stage
For the recalescence 2nd stage, differential equations are not required in the present model. Once nucleation occurs, it is necessary to locate the ice crystals initially formed. Two hypotheses were formulated by Hindmarsh et al. [4], the first is that the nucleation initially occurs at the outer surface of the droplet, which is normally colder than the inside of the droplet and thus first reaches the nucleation temperature. This leads to the formation of a spherical shell of ice which propagates inward over time. The second hypothesis, considers that nucleation occurs homogeneously, with crystals uniformly dispersed throughout the droplet, forming a liquid-solid mixture with an opaque appearance. The recalescence model is based on the premise that the heat required to raise the droplet temperature from T n (nucleation temperature) to T f , must be equal to the latent heat released to form the volume of ice formed. This can be expressed as [4] where L[J/kg] is the latent heat of solidification, T n [K] is the nucleation temperature, V ice [m 3 ] is the volume of ice, and V dp [m 3 ] is the volume of the droplet. For the hypothesis of a spherical shell at the surface, the final position of the interface (R ini ), which is the initial position of the moving boundary in the next stage (solidification), is obtained by considering the shell volume, to yield For the second hypothesis, a water-ice mixture can be considered as a uniform phase, homogeneously distributed throughout the droplet; however, the latent heat should be substituted by a new value corresponding to the water-ice mixture:

Solidification (3rd) stage
The solidification stage model involves the transient one-dimensional heat conduction equation in spherical coordinates for a time-varying domain, is the freezing front position, again with nonlinear boundary condition at the fixed boundary that accounts for convective heat transfer, convective mass transfer, and radiative heat transfer, and prescribed temperature, T f , at the moving boundary position, as in [13]. Once more, the equations are rewritten in dimensionless form and transformed to Cartesian coordinates through a variable transformation [29]. For this 3rd stage, the following dimensionless parameters are used: where In Eq. (12), τ 1 ∼ = τ 2 , and τ 2 is the dimensionless time value after completion of the recalescence phase, and τ 3 is the dimensionless time value when the solidification stage ends. The dimensionless moving boundary differential equations, for each hypothesis of ice formation during recalescence, are given by and where ρ v,ice [kg/m 3 ] is the water vapour density at the solid droplet surface.

Correlations
Murphy and Kopp [30] present a literature review of correlations for saturated water vapor pressure as a function of temperature. Among these, the correlation of Bohren and Albrechet [31] was chosen in the present work. Thus, the equations for ρ v,l , ρ v,ice , ρ v,∞ in terms of the dimensionless temperature for water in liquid and solid states are shown below, where R H is the relative humidity in the air: For the calculation of the convective heat (h) and mass (h m ) transfer coefficients, correlations for the Nusselt and Sherwood numbers were taken from Beard [32]: where Nu is the Nusselt number, Re is the Reynolds number, Pr is the Prandtl number, Sh is the Sherwood number, and Sc is the Schmidt number, defined as is the dynamic viscosity of air, and α ∞ [m 2 s −1 ] is the thermal diffusivity of air.

Model reduction: coupled integral equations approach
The CIEA reformulation methodology [20][21][22][23][24][25] is now applied to the above dimensionless equations. Different levels of approximation in such mixed lumped-differential formulations can be used, starting from the plain and classical lumped system analysis to improved formulations, which are obtained through Hermite-type approximations for integrals [25][26][27] that are based on the values of the integrand and its derivatives at the integration limits. In the present work, we consider just the two more usual approximations, H 0,0 and H 1,1 , which correspond to the classical trapezoidal and corrected trapezoidal rules, given by These two approximations for integrals can be employed in either the average dimensionless temperature or average heat flux (temperature derivative) definitions, and lead to different final expressions with improved accuracy compared to the classical lumped system approach. The final expressions for the reduced model obtained from the CIEA implementation are consolidated in the next sections for each stage of the freezing process and for each proposed approximation, H 0,0 /H 0,0 or H 1,1 /H 0,0 , where the first symbol represents the formula employed to approximate the auxiliary average temperature, while the second one corresponds to the formula that approximates the average temperature derivative. The aim here is to reach a reduced model of similar simplicity as the classical lumped system analysis, but with an improved accuracy, thus extending the range of parameters for which the simplified analysis is still applicable. The classical lumped system would essentially be reproduced by adopting a rectangle integration rule, since the average potential would just be made equal to the boundary value in this case. Through CIEA, using either the H 0,0 /H 0,0 or the H 1,1 /H 0,0 approximations, a reduced model as simple as that for the classical lumping is achieved, but already with considerable accuracy improvement, in particular when the corrected trapezoidal rule is employed since further information on the temperature derivatives at the boundary are incorporated into the reduced model. In previous papers [20][21][22][23][24][25]27], it has been clarified that the Biot number values are essential in defining the precision of the approximations, within specified ranges of these parameters, with the H 1,1 /H 0,0 approximation offering improved results compared to the lower order H 0,0 /H 0,0 approximation. The Hermite integration formulae [26] are much more general, and provide integration rules of higher orders, which can offer further improvement in terms of accuracy to the reduced models. However, a price needs to be paid in terms of adding more dependent variables to the lumped formulation, such as additional boundary temperatures and derivatives.
3.1 1st stage: CIEA H 0,0 /H 0,0 Firstly, the spatially averaged dimensionless temperature is recalled as It should be noted that the above definition corresponds to an averaging of the temperature transformed to the Cartesian coordinates system, and thus does not correspond to the actual physical average temperature within the droplet. The same averaging process is applied to Eq. (2) through the operator 1 0 dx, and after recalling the above definition of the average dimensionless temperature, we obtain ∂θ * Equation (27) is now employed to approximate both the average temperature and the average temperature derivative based on the Hermite-type approximation for integrals [20,25], here named the H 0,0 /H 0,0 approximation, in the form The boundary conditions, Eqs. (4)-(5), are substituted into Eqs. (31), (32), to yield Substituting Eqs. (33), (34) into Eq. (30) leads to the improved lumped-differential formulation for the average dimensionless temperature within the droplet along the 1st stage: In order to avoid working with the averaged dimensionless temperature differential formulation above, which is essentially an auxiliary dependent variable, one may alternatively obtain the differential equation for the dimensionless temperature at the boundary surface θ * 1 (1, τ ), followed by the expression for the dimensionless temperature at the droplet center θ * 1 (0, τ ). Thus, substituting Eq. (33) into Eq. (35), and invoking Eq. (1a), we obtain where while the dimensionless temperature at the droplet center is obtained from with θ l,1 (1, τ ) = θ * 1 (1, τ ). For the solidification stage, the spatially averaged dimensionless temperature is defined as To simplify the manipulations that follow, the average dimensionless temperature is written in terms of the auxiliary dependent variable,θ av,3 (τ ): Following the same averaging procedure, the operator η(τ ) 0 dx, is applied over Eq. (12), and Leibniz rule for differentiation of integrals is recalled in the form which can be simplified since θ ice (η (τ ) , τ ) ≡ 0, to yield The H 0,0 /H 0,0 approximation for the solidification stage is proposed as Applying Eq. (14) into Eq. (41), the following relation between θ * 3 (0, τ ) andθ av,3 (τ ) is obtained: Substituting the boundary conditions, Eqs. (14), (15) into Eq. (42), one finds ∂θ * So, combining the equations above, the following equation is obtained: Again, to avoid working with the auxiliary variable, the same procedure as in Sect. 3.1 is applied, substituting Eq. (43) into Eq. (45), and thus obtain the differential formulation for the dimensionless temperature, θ * 3 (0, τ ), in the form: where B 3 θ * 3 (0, τ ) = Bi c,3 + 4Bi r, 3 is given in terms of the auxiliary potential by Eq. (43), while θ ice,3 (0, τ ) is the actual dimensionless temperature, as given by Eq. (11a).
The improved lumped-differential formulation is then completed with the dimensionless moving boundary position differential equations, for each hypothesis of ice formation during recalescence, Eqs. (16a, 16b) or (17a, 17b), after substitution of Eq. (44) for the temperature derivative at the moving interface, ∂θ ice In the fourth stage, the procedure is similar to that presented in the first stage, taking into account the changes in the mathematical model described at the end of Sect. 2.1.

1st stage: CIEA H 1,1 /H 0,0
Again for the first stage, but now a higher order approximation is proposed, thus the H 1,1 /H 0,0 approximation of the average dimensionless temperature and temperature derivative is written as where the corrected trapezoidal rule is adopted in approximating the average dimensionless temperature, while the plain trapezoidal rule is again used for the average dimensionless temperature derivative. The procedure is essentially the same as for the previous approximation, except for the relation between the boundary and average temperatures, θ l (1, τ ) and θ av,1 (τ ), now obtained from Eqs. (49), (50) and the boundary conditions, Eqs. (4)-(5), as Thus, the improved formulation obtained through the H 1,1 /H 0,0 approximation is given by As before, avoiding the auxiliary averaged temperature and rewritting the formulation for the boundary temperatures at the surface and at the center of the droplet, we obtain ∂θ * 1 (1, τ ) ∂τ where while the dimensionless temperature at the droplet center θ l,1 (0, τ ) is obtained similarly as presented in Eq. (36d), but now with the corrected expression for θ * 1 (1, τ ) provided by Eq. (53a), as The trapezoidal and corrected trapezoidal integration rules for the H 1,1 /H 0,0 formulation in the 3rd stage are applied to the integrated dimensionless temperature and its derivative, as before, in the form The same basic steps are followed as for the model reduction through the previously adopted H 0,0 /H 0,0 approximation. However, the new relation between θ * 3 (0, τ ) andθ av,3 (τ ) is obtained here employing the boundary conditions, Eqs. (14), (15), substituted into Eqs. (54), (55), to yield Then, the following equation is obtained: and rewritting the reduced model in terms of the droplet boundary temperature, as before, the model is redefined as where The improved lumped-differential formulation is then completed with the differential equations for the dimensionless moving boundary position, for each hypothesis of ice formation during recalescence, Eqs. (16a, 16b) or (17a, 17b), after substitution of the temperature derivative at the moving interface, ∂θ ice ∂ x x=η (τ ) . These equations are similar to Eqs. (47) and (48), where the main difference is that the information about θ * 3 (0, τ ) is now provided by Eq. (58a), in the form dη (τ ) 3.6 4th stage: CIEA H 1,1 /H 0,0 As discussed in Sect. 3.3, in the fourth stage, the procedure is similar to that presented in the first stage.

Results and discussion
The derivation of the reduced models and the numerical solution of the resulting ODEs as obtained by the CIEA, for both the H 0,0 /H 0,0 and the H 1,1 /H 0,0 approximations, are solved through a symbolic-numerical code built on the Wolfram Mathematica® platform [28]. Before presenting results of a parametric analysis for the freezing process, it is essential to verify and validate the present model reduction approach.  [13]. More information on the GITT formalism can be found in [14,15]. Second, the reduced model results were critically compared against experimental results for freezing water droplets. For suspended particles, the experiments of Hindmarsh et al. [4] provide the most relevant set of results, which has been the preferred one in validations of distributed or lumped models. As mentioned before, the experiments of Hindmarsh et al. [4] will be shown to lead to Biot number values around 0.1 or less, and the improved lumped formulations are expected to have good accuracy at such values. Nevertheless, the reduced model was also challenged to reproduce the GITT benchmark results for much higher values of Biot number, even up to 10 (two orders of magnitude larger than the experimental value). Such higher values of Biot would require larger droplet diameters and higher velocities (thus higher heat transfer coefficients) to occur, more typical of in-flight icing applications. To our knowledge, experimental values on a single droplet freezing behavior are not readily available in the open literature for such cases. The present reduced model is then subjected to a parametric analysis, for a typical droplet freezing process, by varying the radius of the droplet and the airflow velocity, which essentially affect the Biot number values.

Biot number variation
The Biot number for heat transfer represents a measure of the ratio of the convective and conductive heat fluxes at the droplet surface and it is an important governing parameter for this application. As the Biot number is defined in terms of the ratio "h R/k" the size of the droplet and the stage/phase may cause variations in its value, markedly influencing the freezing time of the droplet. A few values of Bi c,1 were chosen, covering a fairly wide range for the application, to explore the limits of applicability and to demonstrate the accuracy of the CIEA. To achieve this objective, results obtained by CIEA are compared against those for the full partial differential model, Eqs. (2)-(5), obtained by the GITT hybrid approach in [13]. Figure 1a-d presents as comparison of the dimensionless boundary temperature evolution for the first stage, θ l,1 (1, τ), for increasing values of the Biot number at the surface of the water droplet, namely Bi c,1 = 0.1, 1.0, 5.0, 10.0, as computed from the improved lumped-differential formulations (H 0,0 /H 0,0 and H 1,1 /H 0,0 ) and from the full model by GITT (M = 40), with τ ranging from 0 to τ 1 for the supercooling (1 st ) stage. It should be noted that the end of the 1 st stage (i.e. when the droplet reaches the nucleation temperature) will be different for each value of Bi c,1 and is also different for each approximation. As will be seen in Sect. 4.3, the Biot numbers for the experiment of Hindmarsh et al. [4] are around 0.1. On the other hand, from previous works on the CIEA, it is known that the classical lumped system analysis loses accuracy considerably for Biot numbers greater than about 0.1, while the improved lumping here proposed can still offer good accuracy for higher values of this parameter. Therefore, the above Biot number values of 0.1, 1, 5 and 10, up to two orders of magnitude larger than the ones typical of the considered experiment, were chosen to analyze the results and challenge the model reduction methodology. These values were defined arbitrarily, knowing that lower Biot number values lead to more uniform temperature profiles, thus favoring the lumping approach accuracy. For the lower Biot number, represented by Fig. 1a, Bi c,1 = 0.1, the two lumped formulations do not show a marked difference between them, and both are fairly accurate approximations to the partial differential formulation results (GITT). For low Biot numbers (e.g., Bi c,1 = 0.1) the water droplet has approximately uniform temperature fields along its radius, favoring the application of such lumping schemes. The classical lumped system analysis essentially equates the boundary and average temperature values and would still provide reasonable results for such lower values of Bi c,1 . On the other hand, the CIEA seeks to obtain an improved relation between the boundary and average temperatures, through the application of Eqs. (27)-(28) into the lumped form of Eqs.
(2)- (5), and the greater the order of the formulation (H 1,1 /H 0,0 is of higher order than H 0,0 /H 0,0 ), the more accurate the results are expected to be, as can be confirmed from Fig. 1. As the Biot number increases, the deviations between the H 1,1 /H 0,0 and H 0,0 /H 0,0 formulations become more evident, and their respective deviations to the GITT, become more noticeable. For the case with Bi c,1 = 1.0 (Fig. 1b), the H 0,0 /H 0,0 is still reasonable, while better accuracy of the H 1,1 /H 0,0 formulation begins to be apparent, when compared to the benchmark GITT results. In the last cases of higher Biot number values, when Bi c,1 = 5.0 and 10.0, Fig. 1c, d, the fact that the H 0,0 /H 0,0 approximation is indeed less accurate is clearly noticeable, including the marked deviations in the prediction of the recalescence stage onset. On the other hand, the H 1,1 /H 0,0 model remains reasonably accurate in relation to the GITT solution and can be used as a reliable reduced model for this problem, even within this fairly wide range of Biot number. It is clear from Eqs. (36a) and (53a) that the formulation H 1,1 /H 0,0 carries more information about the problem than the H 0,0 /H 0,0 approximation, and the correction of the average temperature approximation with the aid of the temperature derivatives at the boundaries is essential to this improved behavior. As will be seen in Sect. 4.3, the higher values of Biot adopted in these comparisons are not actually typical of this application, but allowed us to establish the applicability range in terms of this important parameter, for the enhanced lumping procedure here proposed.

Biot number and Stefan number variations
The Stefan number (St) is the characteristic dimensionless parameter encountered in phase change problems. It is defined as the ratio of the magnitudes of sensible and latent heats exchanged by the system.  Fig. 2 increases markedly as the end of the solidification is approached. On the other hand, according to Eq. (17a, 17b), this velocity is directly proportional to the spatial derivative of the temperature at the interface, which is in fact approximated by the lumped-differential formulations. Therefore, even a relatively small error in this quantity may induce a significant variation on the prediction of the final solidification time, as can be observed in Fig. 2c, e, though less noticeably in the dimensionless temperature predictions. The droplet central temperatures from the improved lumped formulation remains fairly accurate, but reach their final values at the slightly different values of final solidification time, as discussed above. As expected, the freezing process duration is noticeably linked to the Biot number and Stefan number values and it is much more rapid for both larger Bi c, 3 and St. The present methodology is now employed in a typical situation of supercooled droplet freezing. For comparison purposes, the same parameters as reported in the experimental-theoretical work in [4] are here adopted. Table 1 summarizes the input data. Figure 3 shows a comparison of the experimental results from Hindmarsh et al. [4] with the present study for the evolution of the dimensionless droplet center temperature in the super-cooling (1st) (Fig. 3a) and cooling (4th) (Fig. 3b) stages, as measured from a thermocouple that also holds the droplet. The results via both the CIEA H 1,1 /H 0,0 formulation and the GITT benchmark agree quite well with those experimentally obtained in Hindmarsh et al. [4], to within 1% relative deviation along the transient, offering an important validation of the present model.

Recalescence stage
As explained above, it was assumed that the recalescence stage occurs instantly, and therefore the application of either CIEA or GITT is not necessary. However, R ini and φ are parameters that need to be calculated to provide   Table 1 and other input data from [4], it is then obtained φ = 0.7385, R ini = 0.758 mm, for hypothesis 1, and V ice = 0.16 mm 3 , L x = 2.6 × 10 6 J kg −1 , for hypothesis 2.

Solidification stage and parametric analysis
The dimensionless position of the freezing front was computed considering the two hypotheses in Sect. 2.3, through both the full and reduced models (PDE/GITT and ODE/CIEA H 1,1 /H 0,0 ). For the two hypotheses, the dimensionless positions of the freezing front (CIEA solution) are in good agreement with the GITT solution for the solidification process (3rd stage), as can be seen in Fig. 4. Besides, the results obtained for both recalescence models were fairly close, with the freezing time for hypothesis 1 being shorter than the freezing time for hypothesis 2.
Variations on the airflow velocity ("v") and on the water droplet radius ("R"), for the same case represented by the input data in Table 1, are now analyzed. These variations directly affect important dimensionless numbers, being accounted for within the correlations for Nusselt and Sherwood numbers, and thus directly influencing the Biot numbers for heat and mass transfer. The chosen cases to describe the variations in the airflow velocity and in the water droplet radius are shown in Table 2, together with the corresponding values of the dimensionless numbers. Figure 5a-d shows the results of the proposed parametric study for the solidification (3rd) stage, through the time behaviors of the dimensionless central temperature and the dimensionless position of the freezing front, v(τ ). As Fig. 4 Comparison of GITT solution for full partial differential model and CIEA lumped-differential formulation (H 1,1 /H 0,0 ) for the dimensionless position of the freezing front along solidification stage, via the two hypotheses for the recalescence period, St = 0.11 can be seen from Table 2, an increase on the airflow velocity from 0.42 to 0.97 m/s corresponds to a significant variation in the main dimensionless parameters for the freezing model, since such an increase on the Sherwood and Nusselt numbers represents a marked effect on mass and heat transfer by convection in the system, respectively. Moreover, it can be observed that the Biot numbers remain with low values, clearly within the range previously analyzed that warrant a fairly accurate CIEA reduced model. The results show, as expected, that droplets under higher airflow velocities and of smaller radius, freeze more quickly, while the present reduced model methodology provides a fairly accurate and quite unexpensive way of estimating the total freezing time and the evolution of the moving boundary. In the typical range of parameters for this application, the CIEA results for the solidification stage also provide excellent predictions of the dimensionless central temperatures when compared to the GITT reference results.

Conclusion
A theoretical analysis was performed on the freezing of a supercooled water droplet, including all the stages of the process (supercooling, recalescence, solidification, and cooling). The energy balance for each stage was reformulated by the Coupled Integral Equations Approach (CIEA), which is a tool for generating improved lumpeddifferential formulations in diffusion and convection-diffusion problems. Thus, the original partial differential model for heat transfer with phase change was reduced to ordinary differential models at each stage of the process, in terms of the dimensionless droplet temperatures at the boundaries and, when applicable, also the dimensionless position of the freezing front. The numerical solution of the two proposed CIEA formulations (H 0,0 /H 0,0 and H 1,1 /H 0,0 approximations) of improved accuracy orders with respect to the classical lumped system analysis, were critically compared with a precision-controlled hybrid numerical-analytical solution of the full partial differential model, based on the Generalized Integral Transform Technique (GITT). Results for the dimensionless boundary temperatures and freezing front position along the transient freezing process were then analyzed in terms of the Biot numbers, exploring the limiting range for the enhanced lumping procedure here proposed. It can be concluded that at low values of the Biot number for heat transfer (Bi c < 1) both CIEA approximations, H 0,0 /H 0,0 and H 1,1 /H 0,0 , work well for the present application, but when the Biot number value increases much further, especially the lower order approximation shows noticeable deviations. In parallel, such approximations were validated against experimental results in the literature for an actual application of water droplet freezing, with excellent agreement. A comparison was provided involving two different hypotheses for modelling the recalescence stage, which define the initial conditions for the solidification stage. Finally, a parametric analysis was performed for the freezing stage, varying the values of airflow velocity and droplet radius, within typical ranges of the related application, showing that for larger droplets, the freezing time is longer and, that the greater the air flow velocity, faster will be the droplet freezing process. This analysis provides confidence in employing the advanced Coupled Integral Equations Approach as an effective model reduction tool to simulate the transient behavior along the entire freezing process of a supercooled spherical droplet. The CIEA is a fairly general problem reformulation approach and, in principle, more involved mathematical models could still be considered incorporating further physical effects. Nevertheless, the comparisons with experimental results have demonstrated that the present model is sufficiently complete to recover the actual physical behavior for the present situation and in the range of parameters considered. It is also expected to provide affordable and reasonably accurate simulations for more complex situations, such as in the analysis of droplet sizes distributions or droplets in contact with a substrate.