Hybrid integral transform analysis of supercooled droplets solidification

The freezing phenomena in supercooled liquid droplets are important for many engineering applications. For instance, a theoretical model of this phenomenon can offer insights for tailoring surface coatings and for achieving icephobicity to reduce ice adhesion and accretion. In this work, a mathematical model and hybrid numerical–analytical solutions are developed for the freezing of a supercooled droplet immersed in a cold air stream, subjected to the three main transport phenomena at the interface between the droplet and the surroundings: convective heat transfer, convective mass transfer and thermal radiation. Error-controlled hybrid solutions are obtained through the extension of the generalized integral transform technique to the transient partial differential formulation of this moving boundary heat transfer problem. The nonlinear boundary condition for the interface temperature is directly accounted for by the choice of a nonlinear eigenfunction expansion base. Also, the nonlinear equation of motion for the freezing front is solved together with the ordinary differential system for the integral transformed temperatures. After comparisons of the solution with previously reported numerical and experimental results, the influence of the related physical parameters on the droplet temperatures and freezing time is critically analysed.

for the case when the supercooled water droplets are subject to thermal convection on their surface. Again, a Stefan's one-phase problem was solved by the perturbation method, considering low Stefan numbers. The influence of convection on the freezing process was evaluated by varying the Biot number. As in Feuillebois et al. [15], the obtained perturbation solution deviated from a numerical solution based on the finite difference method under the enthalpy formulation when the droplet centre is approached. Outside this region, a good agreement with the numerical method is achieved. As expected, the results demonstrated a significant influence of the Biot number on the heat transfer within the droplet. The perturbation method was also applied to Stefan's two-phase problem for the solidification of spheres. McCue et al. [17] studied Stefan's two-phase problem for low Stefan numbers. The freezing time and temperature distribution in both phases were obtained considering the one-dimensional problem and a constant temperature condition on the outer surface of the sphere. An important conclusion of this study was that, under the conditions investigated, the influence of the liquid phase on the solidification process could be neglected. More recently, semi-analytical solutions have been used to solve Stefan's problem related to the ice crystal nucleation in the recalescence stage. Dragomirescu et al. [18] used a perturbation method as a starting solution, followed by a numerical method, to solve the one-dimensional Stefan's problem and, by combining these two methods, the authors were able to overcome the limitations of the perturbation method and to reduce the computational cost associated with using the numerical method alone. Schulte & Weigand [3] also employed a semianalytical solution to study the same class of Stefan's problem, using a similarity solution based on a sub-scale model in order to deal with the different scales involved in the initial ice growth process in the recalescence stage.
Numerical methods were also extensively used to solve Stefan's problem related to droplet solidification [19][20][21]. Hindmarsh et al. [1] carried out a numerical-experimental study for the case of solidification of water droplets immersed in a cold air stream and suspended by a thermocouple. The freezing time was obtained for various values of ambient temperature, droplet size and airflow velocity. The abrupt change in temperature measured by the thermocouple and the change in droplet opacity, allowed the authors to experimentally determine the droplet nucleation temperature for each studied scenario. The values measured for the nucleation temperatures were then employed in the computer simulation. For the numerical study, it was assumed that the droplet keeps its spherical shape and therefore the equations can be analysed in a spherically symmetric formulation (i.e. with a single spatial dimension). Each of the four steps of the process mentioned above was analysed separately, and the solidification step was solved considering Stefan's two-phase problem. Two different models were considered in relation to the internal temperature distribution within the droplets, one considering the transient onedimensional heat conduction, with moving boundary in the solidification stage, and the other considering the classical lumped system analysis (named uniform temperature method in [1]), in which the temperature of the droplets was assumed to be spatially uniform for the supercooling and cooling stages. For the solidification stage, the proposed approach involves a heat balance model, where the whole droplet freezes at a constant temperature (T f ). The model considering heat conduction in the droplet was solved using the finite difference method, while the classical lumped model was directly solved by balancing the heat loss to the external environment with the latent and sensible heat released by the droplet. The solutions obtained for both models were compared with the experimental results, showing good agreement.
With the increasing interest in the application of hydrophobic coatings to mitigate the formation of ice on different structures, several studies have also been carried out on the freezing of water droplets deposited on surfaces [22][23][24][25][26]. Song et al. [22] made a comprehensive review of the experimental data associated with the characteristics of water droplets solidification on cold surfaces, providing an overview of important concepts involving the freezing process. The discussion of significant topics, such as freezing front shape, air bubbles formation and time scales of the process, is key for future research on this subject. Regarding numerical studies, many of them had the freezing of suspended water droplets as a starting point, since several conclusions obtained in those works extend to the case of freezing of water droplets on solid surfaces. Chaudhary & Li [23], for example, numerically studied the freezing of static droplets on surfaces subjected to rapid cooling and with different contact angles. As in the work of Hindmarsh et al. [1], the solidification process was divided into four distinct stages. Two-dimensional models were solved using the finite-element method and experiments were conducted for comparison with the numerical analysis, with the satisfactory agreement.
In the context of methodologies for solving moving boundary problems, hybrid numericalanalytical techniques have been gaining special attention in the literature. The generalized integral transform technique (GITT) is one of such hybrid techniques, and has been used successfully in solving various classes of diffusion and convection-diffusion problems [27][28][29][30][31], including the important class of problems with moving boundaries [32][33][34]. The GITT as applied to transient convection-diffusion involves the proposition of an eigenfunction expansion for the original potential spatial behaviour. Then, an analytical integral transformation is performed and an ordinary differential system in the time variable results for the unknown transformed potentials. In the most general situation, the nonlinear coupled structure of the ODE system requires a numerical treatment considering a truncated version of this infinite system, which characterizes the hybrid numerical-analytical nature of the methodology. The algorithm offers automatic error control together with a fairly modest overall computational effort, since the behaviour of the solution in all space variables is reconstructed analytically, while the main computational task is concentrated on the numerical solution of an ODE's initial value problem. The relative merits of both accuracy and computational cost were more closely discussed in the realm of classical test cases, such as in [35][36][37], but it should be noted that the advantages of the hybrid approach become even more evident as the number of spatial dimensions is increased, which affects mostly the analytical involvement, and when intensive computational tasks are considered, such as in optimization, inverse problem analysis and simulation under uncertainty, when a large number of evaluations of the direct problem is required. Some recent advances in the GITT have allowed for the use of techniques to further accelerate convergence in complex problems [38], and among these advances, we can mention the methodology variant introduced by Cotta et al. [39], where nonlinearities related to the boundary conditions of convection-diffusion problems are incorporated into the eigenvalue problem adopted in the integral transformation expansion base. This approach showed a significant convergence gain when compared with the more traditional GITT approach with a linear eigenvalue problem, in addition to presenting a uniform convergence behaviour over all the solution domain, even in positions close to the boundary conditions with nonlinear formulation [39].
The goal of this work is the extension and application of the GITT in simulating the solidification of supercooled water droplets immersed in a cold air stream. The droplets are considered to be freely suspended in air and maintaining their volume and spherical shape throughout the process, since the difference between ice and liquid droplet radii is less than about 3% when considering the typical values for liquid water and ice densities. The equations were written in spherical coordinates for one spatial dimension, by exploiting the spherical-symmetric configuration of the problem allowed for this simplification. The four stages of the process were sequentially analysed. The problem formulation takes into account the nonlinear effects of convection, radiation and convective mass transfer present on the surface of the droplets. In the solidification stage, Stefan's one-phase problem was solved, considering heat conduction in the solid phase only. For each of the stages, except the recalescence stage, GITT was applied to obtain the temperature variation with position and time, and to obtain the duration of each stage. The application of GITT in this class of problem provides more flexibility, in terms of solving complex nonlinear partial differential equations, compared with the approximate analytical methods previously used (e.g. perturbation methods), pulling out restrictions such as the range of the Stefan number that allows the achievement of a reliable solution, and at the same time offering a good computational precision-cost ratio, typically encountered in purely analytical methods. The nonlinear eigenvalue problem methodology introduced by Cotta et al. [39] was extended and combined with the use of implicit filters in order to enhance the solution's convergence rates, reducing furthermore the computational costs. This work brings some advances related to the       GITT method itself, since the nonlinear eigenvalue methodology mentioned above was applied for the first time to moving boundary problems, together with the use of implicit filters, and taking into account highly nonlinear effects in the boundary conditions such as convective mass transfer and radiative heat transfer. The results obtained through the application of GITT were compared with numerical and experimental results previously reported in the literature. An analysis of the influence of the problem physical parameters on droplet temperatures and freezing times was also performed. The list and description of variables used in this work can be found in table 1.

Problem formulation
Due to space limitation, the full description of the dimensional problem formulation in all four stages is provided in the electronic supplementary material that accompanies the article, while here only the dimensionless forms are presented, without loss of content.
(a) Supercooling (stage 1) and cooling stages (stage 4) The supercooling and cooling stages are modelled using the one-dimensional heat conduction equation in spherical coordinates, considering constant thermophysical properties. The following dimensionless groups were chosen: 2) where T l is the liquid (water) temperature, k l is the liquid thermal conductivity, c l is liquid specific heat, ρ l is the liquid specific mass, r is the radial distance from the droplet centre, t is time, R is the droplet external radius, T ∞ is the air temperature and T o is the initial temperature. In addition, prior to the integral transform solution, the formulation can be rewritten similarly to the Cartesian coordinates system, by using the following classical variable transformation: Then, the supercooling stage is written in dimensionless form as and T ∞ (θ 1 (1, τ 1 )) 2 e (19.83−5417/T ∞ θ 1 (1,τ 1 )) (2.14) and where Bi c,1 represents the Biot number for convective heat transfer at stage 1, Bi m,1 is the Biot number for mass transfer, Bi r,1 is the dimensionless group related to the radiative heat transfer, h is the convection heat transfer coefficient, h m is the mass transfer coefficient, L e is the latent heat of evaporation, ε is the emissivity of the droplet surface, σ is the Stefan-Boltzmann constant, and ρ v ,o is the specific mass of the saturated water vapour at 273 K. For stage 1, when the droplet is in liquid state, ρ v ,∞ and ρ v ,l are given by the relations above, equations (2.13-2.14) [40], where R H is the relative humidity of the air. The cooling stage (stage 4) is modelled in a similar way as the supercooling stage above (stage 1). However, the thermophysical properties of the liquid phase need to be changed by the properties of the solid phase. Also, the initial temperature distribution, which for the supercooling stage (stage 1) is the uniform temperature T o , for the cooling stage (stage 4) it becomes the spatially varying temperature distribution in the droplet at the end of the solidification stage. Since the droplet is solid in this stage, the convective mass transfer occurs not through evaporation, but through sublimation, and therefore the latent heat of evaporation (L e ) also needs to be substituted by the latent heat of sublimation (L sb ). For stage 4 when the droplet is in solid state, ρ v ,s and ρ v ,∞ are given by the following relations [40]: and The recalescence stage is modelled as an adiabatic process, since the duration of this stage is considered to be very short, being much faster than the thermal diffusion through the droplet [22,23]. A global heat balance is adopted, following the relation proposed by Hindmarsh et al. [1], which is based on the assumption that the required heat to raise the droplet temperature from the nucleation temperature (T n ) to the freezing temperature (T f ), must be equal to the latent heat released to form the ice volume produced by nucleation. The relation proposed by Hindmarsh et al. [1] is given as where V s is the solid (ice) volume formed, V dp is the droplet volume, c l is the specific heat of the water in the liquid state, T f is the freezing temperature, T n is the nucleation temperature, L is the latent heat of solidification, ρ l is the liquid specific mass and ρ s is the solid specific mass.
In the present work, the ice volume formed after the recalescence stage (stage 2) will be modelled through two different paths. In the first case, the ice is formed as a spherical cask at the droplet surface and the initial position of the freezing front, R ini , is The second case considers a homogeneous distribution of the formed ice at the recalescence stage throughout the droplet. In this case, the water-ice mixture is considered as a uniform phase, and the initial interface position is at the droplet external surface, but the latent heat of solidification L is substituted by a new value for the water-ice mixture [16]. This latent heat of solidification and the liquid fraction in the mixture, ∅, are The solidification stage is a heat conduction problem with a moving boundary, also known as Stefan's problem. The equations for the temperature distribution within the droplet and for the position of the freezing front must be solved simultaneously. The equation for the position of the freezing front, also known as Stefan's condition, describes the relationship between the speed of the interface and the removal of the latent heat released at the same interface [17] and will depend on the hypothesis adopted for the position of the ice formed in the recalescence stage, as will be shown in what follows. As mentioned earlier, this step will be solved as a one-phase Stefan problem and, therefore, the temperature of the liquid phase will be considered constant and equal to the freezing temperature T f and the conduction of the heat released in the freezing front will only occur in the solid phase. Similarly to the stage 1 formulation, this stage is modelled using the one-dimensional heat conduction equation in spherical coordinates, considering constant thermophysical properties. Once more, the original partial differential equations are recast in dimensionless form and the Cartesian coordinate system transformation is implemented. The following dimensionless groups were adopted: and where T s is the solid (ice) temperature, k s is the solid thermal conductivity, c s is the solid specific heat, ρ s is the solid specific mass, r is the radial distance from the droplet centre, t is time, R is the droplet external radius, T ∞ is the air temperature, s(τ 3 ) is the position of the freezing front, L is the latent heat of solidification and St is the Stefan number. The transformation to the Cartesian coordinates system is given by and it is also convenient to introduce an additional coordinate transformation in the dimensionless system as which leads to the following dimensionless problem formulation:

Solution methodology
The above dimensionless problem formulations are now solved by the GITT for each stage governed by partial differential equations (except the recalescence stage). For convergence enhancement, the nonlinearities in boundary conditions shall be directly incorporated into the chosen eigenvalue problem, as originally proposed by Cotta et al. [39]. Implicit filters are also adopted in order to homogenize the nonlinear boundary conditions at the external surface of the droplet.
(a) Supercooling stage (stage 1) The solution of the supercooling stage is started with filtering the dimensionless form of the problem formulation. With the objective of homogenizing the boundary condition of equation (2.10), the following implicit filtering solution is proposed: where F 1 is the proposed filter and θ * 1 is the dimensionless filtered potential. A plain linear function in x 1 has been adopted as filter, given as following filtered problem is achieved: ρ v,l (a 1 (τ 1 ) + θ * 1 (1, τ 1 )) =

(b) Recalescence stage (stage 2)
As previously discussed, the recalescence stage is assumed to occur instantaneously, and the only computation required is the total ice volume formed by nucleation, once this quantity affects the values of φ and R ini . As can be seen from equation (2.20), this quantity is algebraically determined, and does not require any specific mathematical method.

(c) Freezing stage (stage 3)
Again, an implicit filter is proposed to homogenize the nonlinear boundary condition, equation (2.36), as and where F 3 is the proposed filter with a linear functional behaviour in x 3 and θ * 3 is the filtered potential, a 3 (τ 3  boundary condition source term. The filter is then applied to equations ((2.33)-(2.36)), yielding: The filtered problem is again solved using GITT, once more adopting the nonlinear eigenfunction expansion base. Thus, the integral transform pair and the eigenvalue problem, now with a moving boundary, were chosen as The eigenvalue problem is analytically solved to yield: The integral operator η(τ 3 ) 0 ψ 3i ____dx 3 is applied on both sides of equation (3.27), and after some manipulations, the transformed ODE system is obtained as The coefficients a 3 (τ 3 ) and b 3 (τ 3 ) shall be solved simultaneously with the transformed potentials.
The moving boundary equation is now considered, by substituting the filter expression, equations ((3.25)-(3.26)), and the inverse formula, equation (3.36), into equations (2.37) and (2.38), yielding: and As for the supercooling stage, an ODE system is assembled to solve for the eigenvalues, μ 3i (τ 3 ), by differentiating the corresponding transcendental equation, equation (3.42), with respect to τ 3 . The corresponding ODEs and initial conditions at τ 3 = 0, are given by Similarly, the ODEs for obtaining the coefficients a 3 (τ 3 ) and b 3 (τ 3 ) are obtained with the respective initial conditions: . (3.52) The coupled system of ODEs in equations ((3.43)-(3.44)), ((3.47)-(3.48)) and ((3.49)-(3.52)) is then numerically solved for a truncated version of the system with a finite order M. Depending on the adopted hypothesis for ice formation structure in the recalescence stage, either equation (3.45) or (3.46) should be added to the system. Then, the dimensionless temperature distribution within the droplet, θ * 3 (x 3 , τ 3 ), is obtained from the inverse formula, equation (3.36), after the transformed potentials, θ * 3 i (τ 3 ) have been numerically computed. The solidification interface position, η(τ 3 ), is directly obtained from the numerical solution of the ODE system, and the dimensional boundary position is finally determined from: (3.53)

(d) Cooling stage (stage 4)
The solution methodology for this stage is the same as for the supercooling stage (stage 1), with the necessary changes in the initial condition and other parameters, as discussed in the problem formulation section, and skipped here for brevity.

Results and discussion
This section provides validation and verification of the proposed model and integral transform solution, followed by some parametric analysis. First, a simpler model proposed by Tabakova et al. [16], which considers only convection at the external surface of the droplet and has been solved by both approximate analytical and numerical methods, will be considered and compared with results obtained from the present analysis. Only the solidification stage is then verified. Second, the present complete model shall be compared with the experimental and numerical results of Hindmarsh et al. [1], considering all the relevant effects at the droplet surface, namely, convection, radiation and convective mass transfer. Finally, a physical analysis of the most relevant parameters on the overall process is undertaken. The routine NDSolve from the Mathematica v. 11.3 system has been employed in the numerical solution of the transformed ODEs system with automatic accuracy and precision control.   [16] for St = 0.1 and Bi c,3 = 1 with the GITT solution for a fixed truncation order of M = 20. Figure 1 shows results for the dimensionless freezing front position, v(τ 3 ), with excellent agreement between the GITT and the numerical method along the entire time domain. The perturbation method solution, however, presents some deviation by the end of the freezing process. Figure 2 compares the dimensionless temperature profiles in the solid phase for different values of τ 3 . Again it can be observed that the GITT solution agrees notably well with the numerical solution based on the finite difference method with enthalpy formulation. Once more the perturbation method presents marked discrepancies as the end of the solidification stage is approached. In fact, according to Tabakova et al. [16], for a large Stefan number and close to the freezing time, the perturbation solution becomes less consistent, and the values τ 3 = 0.5 and τ 3 = 0.52 are already fairly close to the freezing time in this case.

(b) Complete model comparison with numerical method and experiments
The numerical-experimental work of Hindmarsh et al. [1] was selected for comparison with the complete model, more specifically for the following conditions: a droplet with radius R = 0.78 mm, ambient air temperature T ∞ = -19°C, nucleation temperature T n = −18.4°C and air flow velocity v = 0.42 m s −1 . Hindmarsh et al. [1] measured the freezing time and the temperatures of the water droplets suspended by a thermocouple, then compared the experimental results with an approximate numerical analysis, considering a classical lumped model (which was termed the uniform temperature method in [1]), i.e. the droplet temperature is assumed to be spatially uniform during the cooling stages and constant in the freezing stage. A finite difference solution was also implemented for the cooling stages for verification of the lumped approach. The water thermophysical properties were assumed to remain constant and evaluated at the equilibrium freezing temperature (T f ), equal to 0°C. The properties related to the air are evaluated at −19°C. Hindmarsh et al. [1] mention that the air currents had the humidity removed by silica gel, and thus we set ρ v ,∞ = 0, which corresponds to assuming that the environment is of dry air. The heat and mass transfer coefficients (h) and (h m ) were computed from the relations proposed by Beard & Pruppacher [41] for the Nusselt and Sherwood numbers. Hindmarsh et al. [1] neglected the heat conduction along the thermocouple and justified it through simulations. Table 3 Table 4. Convergence of freezing times and temperatures at the droplet surface for the solidification stage, considering two hypotheses for the spatial distribution of the formed ice after recalescence. Experiment of Hindmarsh et al.    For time values closer to the end of the stage 3, t = 23 s, convergence to six significant digits was achieved at higher orders, M = 60; however, convergence to four significant digits was achieved for truncation orders as low as M = 30. The results obtained from the two recalescence models were not markedly different, with the duration of the solidification stage for the first hypothesis being 24.11 s and 23.60 s for the other one. The duration of this stage experimentally measured by Hindmarsh et al. [1] was approximately 20 s. Figure 3 shows a comparison between the GITT solution, the classical lumped approach proposed by Hindmarsh et al. [1] and the experimental temperature evolutions at the centre of the droplet for the entire process, including all four stages. For the solidification stage, the hypothesis that considers a homogeneous distribution of the ice formed during the recalescence period has been adopted. A truncation order of M = 60 was employed in this comparison. An excellent overall agreement between theoretical predictions and experimental measurements is achieved. The deviations between the GITT solution (solid black line) and the classical lumped approach (dashed blue line) by Hindmarsh et al. [1] can be explained through the distinct modelling adopted. Hindmarsh et al. [1] model assume uniform temperature within the droplet, both for the supercooling stage and the freezing stage. In fact, Hindmarsh et al. [1] themselves also compared these two modelling approaches, distributed and lumped, pointing out deviations in the results obtained.
For additional assessment of the GITT results, a purely numerical analysis of the first stage, considering the full heat conduction equation for distributed parameters, was performed here using the NDSolve routine from the Mathematica system v. 11.3 [42]. The NDSolve algorithm implements the method of lines for the solution of partial differential equations, with variable mesh and variable difference order. This automatic solver for handling partial differential equations is the same employed in solving the nonlinear transformed ODE system, but is here also adopted in verifying the GITT solution for the first stage. The agreement between the GITT and the methods of lines numerical solution is excellent, presenting a relative deviation of around 0.07% for the duration of the supercooling stage. For improved agreement with the fully converged GITT solution, a very refined mesh is required by NDSolve through the corresponding prescribed input parameters. Excellent agreement was also observed between the GITT solution and the numerical solution via NDSolve for the cooling (fourth) stage, with a relative deviation

(c) Parametric analysis
A parametric analysis is undertaken next using the same validation case above as a reference to analyse the influence of different quantities and parameters. First, the relative importance of the three modes of energy transfer at the droplet surface is examined. The heat fluxes related to convective heat transfer (q c ), radiative heat transfer (q r ) and convective mass transfer (q m ) are evaluated from: and Figure 4 represents the evolution of all three heat fluxes during the supercooling stage, as computed from the GITT formalism. The overall contribution of each heat transfer mode to the heat exchange between the droplet and external environment at stage 1 can be obtained by integrating equations ((4.1)-(4.3)) between the limits t = 0 and t = t n , where t n is the nucleation time (i.e. the duration of the supercooling period herein). As expected, convective heat transfer accounts for most of the heat exchange with the surroundings, around 58.2%, while the convective mass transfer accounts for 39% of the overall heat exchange, and radiative heat transfer for only 2.8% of this total. From figure 4, it can be noted that, after around 19 s, the convective mass transfer becomes dominant. As the droplet surface temperature approaches the external ambient  temperature, the difference (T(R, t) − T ∞ ) approaches zero, together with the convective heat flux, while the difference (ρ v ,l (T(R, t)) − ρ v ,∞ ) present in the heat flux related to convective mass transfer approaches a fixed value, since ρ v ,∞ = 0. The importance of convective mass transfer at low-temperature differences had already been pointed out by Strub et al. [43], who observed stabilization temperatures lower than the ambient ones. Figure 5 again shows the evolution of the three modes of heat fluxes but now for the solidification stage. In this case we need to substitute the latent heat of evaporation (L e ) and ρ v ,l for the latent heat of sublimation (L sb ) and ρ v ,s in equation (4.3). The convective heat transfer remains dominant, accounting for 61% of the overall heat exchange in this period, while convective mass transfer and radiative heat transfer roughly account for 36% and 3%, respectively. Here the convective heat transfer has a higher participation in the overall heat exchange compared with the supercooling stage, because the droplet surface temperature remains practically constant during this stage, but different values of the Biot number did lead to modification of this behaviour, as will be observed in what follows.   Quantifying the pathways of heat dissipation, as in figures 4 and 5, is important in assessing the role of each transport mechanism in the freezing process and its influence on the freezing time. Castillo et al. [44] performed a similar theoretical evaluation for heat dissipation rates for the freezing of a droplet resting on a cooled substrate. They estimated the fraction of heat that is lost through the substrate and through the droplet/air interface, either by natural convection or by convective mass transfer, and their analysis has motivated the present comparison. Looking just at the droplet/air interface, the chosen parameters were the ratio between the convective heat flux and total heat flux, and the ratio between the heat flux related to convective mass transfer and total heat flux, both ratios calculated at the droplet/air interface and at the beginning of the freezing stage. These ratios were approximately 0.68 and 0.30, respectively, for the 10.1 µl droplet studied by Castillo et al. [44]. In the present work, these same ratios were nearly 0.60 and 0.36, respectively, for a 2 µl droplet. Disregarding the different boundary conditions and volumes analysed in the two studies, the ratios show plausible orders of magnitude, even more so when taking into account the fact that in the present work the air surrounding the droplet was assumed dry, increasing the influence of the convective mass transfer.
The influence of the Stefan number and the other dimensionless parameters (Bi c,3 , Bi m,3 and Bi r, 3 ) in the droplet cooling process was also analysed. The Stefan number (St) is a frequently encountered dimensionless parameter in phase change problems, being interpreted as the ratio of sensible and latent heats of the system. In the present parametric study, the external air temperature was changed in order to yield different values of the Stefan number. Table 5 summarizes the cases analysed. Case 1 corresponds to the experiments of Hindmarsh et al. [1], used in the above validation. For Case 1, the Stefan number is 0.11 and the freezing time is computed as 23.6 s.   Figure 6 shows that the droplet surface temperature, T(R, t), varies little along the solidification stage for this case. In addition to enlighten the weight of different heat transfer modes in the freezing process, as mentioned earlier, the droplet surface temperature is an important indication of the temperature gradient in the solid phase within the droplet, since the freezing front temperature remains at the equilibrium freezing temperature (T f ), i.e. 0°C. For Case 1, 22.5 s after the start of solidification, the surface temperature drops by only 2.5°C. Towards the end of the solidification period, the droplet surface temperature varies more markedly due to the reduction in the latent heat released at the freezing front.
In Case 2, the Stefan number is about 2 times higher than in the base case, after changing the external air temperature from −19°C to −40°C. In this case, the freezing time is approximately 14.7 s. From table 5 it can be observed that the values of the three parameters, Bi c,3 , Bi m,3 and Bi r,3 , also reduced with the change in the ambient temperature. The parameter Bi c,3 was the least affected one and, as a consequence, the convective heat flux contribution jumped from roughly 60% in Case 1 to around 77% in Case 2. The direct comparison of T(R, t) in Case 1 and Case 2 from figure 6 also shows that the droplet surface temperature variation is pronounced in Case 2. Around 14 s after the onset of solidification, the surface temperature has already dropped by almost 4°C. The increase in Stefan number represents a relative increase in sensible heat transfer in comparison with the latent heat, therefore, one may expect a more significant variation of temperatures during the phase change process. This observation is crucial for processes in which the temperature gradients have to be controlled, such as in the pharmaceutical, food and metallurgical industries, because the constitution of the formed solid is affected by the magnitude of this gradient. In Case 3, the values of Bi c,3 , Bi m,3 and Bi r,3 were increased with respect to Case 1, for instance through changes in the air flow velocity, but keeping the Stefan number equal to Case 1. By changing the air flow velocity from 0.42 to 2 m s −1 , Bi c, 3 and Bi m,3 values almost doubled in comparison with Case 1, although the value for Bi r,3 remained the same. The freezing time in this case was around 13.8 s, and figure 6 again shows the evolution of the droplet surface temperature along the solidification. The decrease in the droplet surface temperature in Case 3 is a little bit faster than in Case 2, i.e. in just 13 s after the start of solidification, the droplet surface temperature has already reduced by almost 4°C, illustrating that the changes in the Biot numbers have a more significant influence compared with the influence of the Stefan number on the temperature gradients along the phase change process.
With respect to the problem modelling, for sufficiently mild temperature gradients within the medium, such as in Case 1, an improved lumped-differential approach can be recalled that markedly simplifies the formulation with an acceptable loss in precision, as thoroughly discussed in [31]. Then, the original partial differential equations can be reduced to ordinary differential equations by the reformulation process. Figure 7 shows the droplet centre temperature, T(0, t), as a function of time for the base Case 1 but with three different droplet radius: (i) 0.39 mm, (ii) 0.78 mm and (iii) 1.56 mm. As expected, the larger the droplet, the greater is the freezing time, being roughly 18.7 s for the 0.39 mm radius and 183.4 s for the 1.56 mm radius, i.e. a 10 times rise for an increase of 4 times in the droplet radius. It has also been observed that with the increase in the droplet size, the contribution of the radiative heat flux became more significant in all the stages. In the supercooling stage for example, it increased from 1.7% to 2.8% when the droplet radius was doubled from 0.39 to 0.78 mm and to 4.1% when the radius was again doubled to 1.56 mm. In this case, the contribution of the convective heat flux and convective mass transfer experienced a mild decrease with the radius increase. The rise of the radiative heat flux participation is explained because of the influence of the droplet radius in the heat and mass transfer coefficients (h) and (h m ), since the value of both coefficients decreases with the increase in the droplet radius.

Conclusion
Freezing of supercooled liquid droplets suspended in air is modelled through a transient one-dimensional heat conduction formulation in spherical coordinates that accounts for the nonlinear effects of convective heat transfer, convective mass transfer and thermal radiation at the droplet surface. The freezing process is divided into four distinct stages, namely, supercooling, recalescence, freezing and cooling stages. Except for the recalescence stage, which is represented by a simplified model as an instantaneous nucleation process, the other three stages result in nonlinear partial differential equations that are solved by the GITT, including the moving boundary problem that represents the freezing stage. This hybrid numerical-analytical approach with automatic accuracy control is here employed based on a nonlinear eigenvalue problem for representing the eigenfunction expansions. The original nonlinear boundary conditions, after an implicit filtering to eliminate the associated source terms, are directly accounted for by the chosen eigenfunctions and eigenvalues. This approach results in a marked improvement on convergence rates for the desired potentials, including the temperatures at the interface and the moving interface position. The model and solution methodology are both verified with existing approximate analytical and numerical solutions for simpler formulations and validated against experimental and numerical results in the literature. Then, parametric combinations are chosen to analyse some of the relevant physical aspects in determining the freezing time, surface temperature and the evolution of the freezing front. The present approach does not introduce any approximations for the spatial dependence of the temperature field, as required in lumped system approaches, or any limitations on the magnitude of the Stefan number, which is typical in perturbation-type approaches. In addition, it can also be extended to the two-dimensional formulation in non-spherical geometries required to handle the situation of a droplet in direct contact with a solid surface, by considering the related two-dimensional eigenvalue problems.