Experimental–theoretical analysis of cooling and freezing of a droplet in contact with a cold substrate: influence of substrate wettability

TL;DR

This paper studies how droplets cool and freeze on cold surfaces, using experiments and models to show how surface properties affect the process.

Contribution

The study introduces a new mixed lumped-differential model for heat transfer during droplet freezing, validated with experimental data.

Findings

The mixed model accurately predicts droplet surface temperatures during supercooling.

Hydrophilic and hydrophobic substrates significantly influence freezing dynamics.

Classical lumped system analysis fails to capture the observed freezing behavior.

Abstract

The physics and modelling of cooling and freezing of droplets in contact with a colder substrate are of interest in various engineering applications. This work provides experimental results of this process employing infrared thermography for temperature measurements at the droplet’s surface. Also, a high-speed camera is employed to observe the recalescence period and measure the freezing front movement and the droplet shape change. Three substrates are prepared with distinct wettability ranges, i.e. one hydrophilic and two hydrophobic surfaces. From the experimental observation of a solidification front parallel to the substrate plane, a mixed lumped-differential model of the heat transfer process based on the Coupled Integral Equations Approach is proposed, reformulating the two-dimensional partial differential formulation in cylindrical coordinates into a one-dimensional transient energy equation for the droplet external surface temperature. Direct comparisons of the experimental and theoretical results for the supercooling period show excellent agreement for the droplet surface temperatures at different heights and for different values of the substrate–droplet contact angle. It is also shown that the classical partial lumped system analysis does not provide adequate predictions in the present problem. Finally, the dynamics of the recalescence and freezing stages are experimentally evaluated and physically interpreted.

This article is part of the theme issue ‘Heat and mass transfer in frost and ice’.

Figures

Sequence of side view images showing the spatial and temporal evolution of the temperature field of a 20 µl water droplet. (a) A super cooled droplet, (b) onset of recalescence (dashed circle) approximately 5 ms later, (c) completion of recalescence approximately 40 ms later, (d) approximately 7 s post-recalescence, the second stage of freezing is shown by the moving ice front, (e) end of freezing (14 s post-recalescence) of the water droplet characterized by a sharp tip at the apex.

Schematic of the experimental set-up. The thermoelectric stage 1 is CP-065, thermoelectric stage 2 is a generic Peltier module. TC 1 and TC 2 are used to monitor substrate and air temperature, respectively. The camera used is either ×901sc (FLIR) with a MIC 1× lens for thermal imaging or v411 (Phantom) with a K2 Distamax lens. The distance between the lens and droplet is denoted as x and equals 30 mm for the IR imaging and 98 mm for optical imaging.

Configuration of a droplet supported on a cold surface. Mathematical model in cylindrical coordinates–spheroidal solids.

Temperature variation at the sessile droplet surface on SHB substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). L1=2.11mm, L2=2.08mm,  Tb=−19.7, T∞=−2.8 ∘C and Tw(Z,t) at the droplet surface.

Temperature variation at the sessile droplet surface on HYB substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). . L1=2.87mm, L2=1.9mm,  Tb=−18.46∘C, T∞=−2.9 ∘C and Tw(Z,t) at the droplet surface.

Temperature variation at the sessile droplet surface on HYL substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). L1=3.63mm, L2=1.55mm,  Tb=−15.2∘C, T∞=−1.04 ∘C and Tw(Z,t) at the droplet surface.

Evolution of the droplet surface temperature, Tw, on three different substrates where SHB is superhydrophobic (plot marker circle), HYB is hydrophobic (plot marker triangle) and HYL is hydrophilic (plot marker diamond). Three locations are chosen for each of the substrates. They are (a) top plane (T1exp), (b) middle plane (T2exp) and (c) bottom plane (T3exp). The subscript ‘exp’ is for the experimental origin of the plotted data.

Surface temperature variation for the sessile droplet on SHB substrate. (a) Temperature variation at three different locations (T1, T2, T3) for the entire duration of cooling and freezing stages (recalescence is indicated in the dashed bo); (b) dashed box is magnified and replotted.

Snapshots showing solidification for substrates with different wettability: (a) HYL (1, 5, and 7 s), (b) HYB (1, 6, and 10 s), and (c) SHB (1, 7 and 15 s). The arrow indicates the position and direction of ice front development. Time stamp is inset in the top-left corner. Scale bar equals 1 mm.

Position of the ice front (S) relative to the substrate as a function of time for experiments on three different substrates.

Tables

Table 1.. Nomenclature.

B(θ)	auxiliary parameter of the nonlinear boundary condition
Bic	Biot number for contact heat transfer between droplet and substrate surface
Bi	Biot number for convective heat transfer between droplet and air
Bim	dimensionless group for mass transfer
Bir	dimensionless group for radiative heat transfer
c	specific heat (J kg−1 K−1)
D	droplet diameter (m)
Da	diffusivity (m2 s−1)
g	gravity acceleration (m s−2)
h	heat transfer coefficient (W m−2 K−1)
hm	mass transfer coefficient (m s−1)
k	thermal conductivity (W m−1 K−1)
L1	liquid droplet height (mm)
L2	liquid droplet base diameter (mm)
Le	latent heat of evaporation (J kg−1)
Nu	Nusselt number
nz	direction cosines of the unit normal vector with respect to the z coordinate
nr	direction cosines of the unit normal vector with respect to the r coordinate
Pr	Prandtl number
q	heat flux (W m−2)
Rd	initial spherical droplet radius (m)
Re	Reynolds number
RH	relative humidity
r	radial coordinate (m)
rw	droplet surface radius along axial variable (m)
R	dimensionless radial variable
S(t)	freezing front position (mm)
Sc	Schmidt number
Sh	Sherwood number
t	time variable (s)
T(r,z,t)	temperature distribution (K)
Tn	nucleation temperature (K)
T0	initial temperature (K)
Tb	substrate temperature (K)
Tw	temperature distribution on the droplet surface (K)
T∞	ambient air temperature (K)
Z	dimensionless axial variable
z	axial variable
Greek symbols
α	thermal diffusivity (m2 s−1)
β∞	volumetric expansion coefficient (1 K−1)
ε	emissivity
θ(R,Z,τ)	dimensionless temperature in cylindrical coordinates
θ¯(Z,τ)	average dimensionless temperature
θ~(Z,τ)	auxiliary dependent variable
μ∞	dynamic viscosity of air (N s m−²)
ρ	density (kg m−3)
ρv,0	water vapour density at 273 K (kg m−3)
ρv,l	water vapour density at liquid droplet surface (kg m−3)
ρv,∞	water vapour density at surrounding air temperature (kg m−3)
σ	Stefan–Boltzmann constant (W m−2 K−4)
τ	dimensionless time variable
Λ	dimensionless parameter (related to the front position)
ν	kinematic viscosity (m2 s−1)
γ	ratio between the radius of the base and the height of the droplet
Subscripts and superscripts
∞	air
av	average
b	base-substrate
c	related to droplet–substrate contact
m	related to convective mass transfer
r	related to radiative heat transfer
ice	solid phase
l	liquid phase
w	wall (substrate)
1 to 4	stage 1 (supercooling), 2 (recalescence), 3 (freezing), 4 (cooling)

Table 2.. Physical properties and input data.

property	value			source
	SHB	HYB	HYL
substrate temperature, Tc (°c)	−19.7	−18.46	−15.2	experiments
air temperature, T∞ (°c)	−2.8	−2.9	−1.04	experiments
initial droplet temperature, t0 (°c)	0 ±1	0 ±1	0 ±1	experiments
liquid droplet height, L2 (mm)	2.08	1.9	1.55	experiments
liquid droplet base diameter, L1 (mm)	2.11	2.87	3.63	experiments
contact angle (°)	130	110	82	experiments
air thermal conductivity, k∞ (W m−1 K−1)	0.0234			Hindmarsh et al. [53]
air diffusivity, Da	2.060 × 10⁻⁵			Hindmarsh et al. [53]
air dynamic viscosity, μ∞ (Ns m−2)	1.663 × 10⁻⁵			Incropera et al. [65]
air specific mass, ρ∞ (kg m−3)	1.3317			Hindmarsh et al. [53]
volumetric expansion coefficient, β∞(1 K−1)	0.034			Incropera et al. [65]
liquid specific heat, cl (kJ kg−1 K−1)	4.2			Hindmarsh et al. [53]
liquid thermal conductivity, kl (W m−1 K−1)	0.607			Incropera et al. [65]
latent heat of evaporation, Le (J kg−1)	2.502 × 10⁶			Hindmarsh et al. [53]
emissivity, ε	0.96			Hindmarsh et al. [53]
liquid specific mass, ρl(kg m−3)	1000			Hindmarsh et al. [53]
specific mass of sat. water vapour, ρv,o (kg m−3)	4.8473 × 10⁻³			Hindmarsh et al. [53]
specific mass of sat. water vapour at the surface, ρv,l,s (kg m−3)	equation (3.2b)			Hindmarsh et al. [53]
Stefan–Boltzmann constant, σ(W m−2 K−4)	5.670 × 10⁻⁸			Incropera et al. [65]

Equations

Full text

Introduction

Water exists as a liquid far below its thermodynamic freezing temperature (super cooled) but readily freezes when in contact with a sufficiently cold solid surface. This is due to the lowered barrier of heterogeneous nucleation compared with homogeneous nucleation [1]. While suppression and removal of frost in closed systems, such as in HVAC applications, has been achieved to some extent [2], atmospheric ice accretion remains a matter of grave concern for aviation [3–6], maritime transport [7–9], wind turbine [10–12], power transmission [13–15] and road transport [16,17] sectors. De-icing by heating and/or vibrating are effective [18] but energy and time consuming [19–21]. At the same time, the use of freezing point depressants are effective passive strategies albeit with an environmental penalty [22,23]. Where applicable, the development of liquid repellent and/or photo thermal [20] coatings are promising but yet to be commercialized. Ice accretes through impact of super cooled droplets or solidification of condensed vapour. A better understanding of atmospheric icing, through observing a super cooled droplet freezing process, is quite desirable towards designing effective icephobic strategies.

Solidification of a liquid sessile droplet is a complex interplay of several physical factors but may be described in four stages [24] as shown in figure 1: (a) *supercooling—*the sessile droplet is cooled below its thermodynamic freezing temperature (T < Te) as shown in figure 1, depending on the substrate’s surface roughness, wettability [25] and ambient conditions [26–28]; (b) *nucleation—*a stable ice embryo forms, either at the liquid–substrate interface or the liquid–air interface [29] by overcoming the Gibbs’ energy barrier to heterogeneous nucleation [1]. The temperature at which nucleation occurs is an important characteristic of icephobic surface; (c) *recalescence—*proliferation of ice nuclei within the droplet accompanied by a rapid rise in its bulk temperature to Te. The droplet is now a frozen mixture of ice and water where the ice concentration is estimated by balancing the release of latent heat to the sensible heat required to raise the droplet’s temperature. The entire process is assumed near adiabatic since the duration of recalescence ( $[eqn]$ is considerably shorter than the thermal diffusion timescale $[eqn]$ . Here Rdrop is the nominal size of a liquid droplet (~O(1) mm), thermal conductivity of water (kwater) is considered as 0.607 W/m.K, density as 997 kg m^−3^ and specific heat capacity (Cp) as 4.2 J/g.K; (d) *solidification—*starting at the droplet’s base and moving in the direction imposed by the temperature gradient [29–31]; (e) finally, solidification culminates (in case of water) into a conical tip at the apex [32,33]. Solidification is limited by the rate of heat transfer to the ambient [34], after which the droplet cools down to the substrate temperature.

Sequence of side view images showing the spatial and temporal evolution of the temperature field of a 20 µl water droplet. (a) A super cooled droplet, (b) onset of recalescence (dashed circle) approximately 5 ms later, (c) completion of recalescence approximately 40 ms later, (d) approximately 7 s post-recalescence, the second stage of freezing is shown by the moving ice front, (e) end of freezing (14 s post-recalescence) of the water droplet characterized by a sharp tip at the apex.

Icephobic surfaces usually focus on delaying the onset of nucleation (freezing delay). Freezing delays up to 25 h [35,36] have been reported. However, atmospheric conditions can trigger nucleation at the liquid–air interface [29] and render icephobic surfaces useless. Solidification or the freezing time should be equally important; a frozen ice droplet is harder to remove, likely to anchor other droplets impinging on the surface and accelerate frost formation [37]. The freezing time can be in principle estimated from the one-dimensional Stefan’s model of solidification [38]. Attempts of theoretically analysing the ice formation process include classical numerical techniques such as the enthalpy method [24,39], enthalpy-porosity technique [34], front tracking technique [40] and heat capacity method [31]. A fairly complete review was provided by Akhtar et al. [41], discussing the various mathematical approaches undertaken to model all stages of droplet freezing. Two less known hybrid numerical–analytical approaches in the droplet freezing field have been recently adopted to simulate the four stages of a suspended droplet freezing process [42–44]. One of them is the Generalized Integral Transform Technique or GITT [45–52], which is a hybrid numerical–analytical approach based on the classical integral transform method for linear partial differential equations, but generalized to handle nonlinear formulations, including moving boundary problems. The GITT approach is demonstrated [42] for one-dimensional transient droplet freezing, adopting a nonlinear eigenvalue problem proposal [51] with markedly improved convergence rates. Accurate benchmark results were then generated and used to co-validate experimental results for a suspended droplet from Hindmarsh et al. [53]. The same hybrid approach was employed in the analysis of a two-dimensional heat conduction formulation for a droplet in contact with a cold superhydrophobic substrate [43]. The second approach considered for simulating freezing of a suspended droplet [44] is in fact a problem reformulation strategy, known as the Coupled Integral Equations Approach (CIEA) [52,54–58], which provides improved lumped-differential formulations by appropriate averaging the full partial differential problem in one or more spatial coordinates, and using information from the boundary conditions in the Hermite formulae for the approximated averaged quantities, such as temperatures and heat fluxes. Such improved lumped formulations were then tested on the same suspended droplet freezing problem [43,44] and compared against the same set of experimental results [53]. Experimental co-validation with Stefan-type models is based on optical imaging of the moving ice front [28,31] and measuring the droplet temperature through both intrusive [24,26,59] and non-intrusive techniques [24,34,60]. Carefully calibrated infra-red measurements [34] are usually better for being non-intrusive, thus avoiding nucleation and heat conduction from within the droplet and providing better spatio-temporal resolution.

In the present work, first the cooling and freezing of a supercooled droplet in contact with a cold substrate is analysed experimentally. Non-intrusive infrared thermography is employed throughout the process to accurately measure the droplet surface temperatures, while a high-speed camera is used to determine the freezing front evolution. Different substrate surfaces have been prepared to provide a range of wettability behaviours, spanning from hydrophilic to superhydrophobic characteristics. Based on the observed almost vertical advance of the freezing front, a lumped-differential model has been proposed, using the CIEA [52,54–58]. Considering the droplet as a deformed cylinder, with a variable external radius as a function of the vertical coordinate, a lumping process is proposed in the radial coordinate only, thus transforming the axisymmetric two-dimensional partial differential formulation for the entire droplet, into a one-dimensional one for the droplet temperatures at the external surface. Then, the theoretical predictions are critically compared with the experimentally determined surface temperatures along the droplet cooling phase, at different vertical positions and for the three different substrate preparations. Finally, the experimental behaviour of the entire freezing process, including recalescence, is illustrated and compared for the three types of substrates. Please see table 1 for nomenclature.

Experimental methods

Aluminium sheets (grade 6061) 1.5 mm thick are cut into 2.5 cm × 2.5 cm coupons. They are degreased by sequentially washing in acetone, iso-propanol, ethanol and water baths. The degreased samples are kept aside to be used as hydrophilic substrates (denoted as HYL in the text). The rest is dipped in 3 M KOH solution to remove the native oxide layer before being anodized in 0.3 M oxalic acid solution. See Grizen et al. [61] for details. The anodized samples are then immersed in 1 wt% solution of octadecyltricholosilane in hexane for 1 h and subsequently heated at 120°C in a temperature-controlled oven for 2 h. The static contact angle on the HYL sample is approximately 80 ± 5°. For the silane-treated anodized sample (referred to in the text as HYB), the static contact angle is approximately 113 ± 5°. In order to decrease the area of substrate–droplet contact, some of the anodized samples are further immersed in 0.3 M phosphoric acid solution for 2 h with constant stirring. This leads to wider anodic pores and combined with the silane treatment gives a static contact angle of approximately 135 ± 5°. This substrate is referred to as SHB.

The substrate (HYL or HYB or SHB) is placed in a three-dimensional printed enclosure (40 mm × 40 mm × 40 mm), in the experimental schematic shown in figure 2. An air cooled thermo-electric cooler (CP-065, TE Technology) is used for all experiments. The topside of the enclosure is fitted with the thermoelectric module (with a heat spreader) to further cool the air inside the enclosure to −1.5 ± 1°C. Two thermocouples are used to measure the substrate temperature (positioned beside the actual substrate) and air temperature (position 10 mm from the substrate). Thermocouple readings are logged using a TC-08 (USB temperature logger, Picolog) at 1 Hz. Prior to the experiments, the thermocouples are calibrated in a re-circulating chiller bath at different temperatures in the range of interest.

Schematic of the experimental set-up. The thermoelectric stage 1 is CP-065, thermoelectric stage 2 is a generic Peltier module. TC 1 and TC 2 are used to monitor substrate and air temperature, respectively. The camera used is either ×901sc (FLIR) with a MIC 1× lens for thermal imaging or v411 (Phantom) with a K2 Distamax lens. The distance between the lens and droplet is denoted as x and equals 30 mm for the IR imaging and 98 mm for optical imaging.

Thermal infra-red (IR) images are acquired using the ×901sc camera (FLIR) fitted with a MIC 1× microscope lens (pixel resolution of 40 µm). The distance between the lens and the droplet is x = 30 mm. In order to avoid spurious background effects, the inside of the enclosure is covered with a black masking tape (Thorlabs). A circular window covered with IR transparent sheet of polypropylene is used to view the droplets. The acquisition frame rate is maintained at 1 frame per second (fps) till the substrate temperature is below −10°C and then continued at 400 fps till the end of freezing. Error due to emissivity and view factor was found to be between −1 and −2°C, which is within 10% of the temperature range measured in this current work. A constant emissivity of 1 was used for all IR image acquisitions.

To better visualize the icing front, separate freezing experiments (under the same conditions) are imaged using a high-speed camera (Phantom v411) fitted with a long-distance microscope (K2 Distamax + CF2 lens, Infinity Photo-Optical, USA) with spatial resolution of 7 µm. The distance between the lens and the droplet is x = 98 mm. Two light sources are used, one opposite to the camera to outline the droplet profile and the other kept perpendicular to illuminate the droplet interior. Images are acquired at 100 fps.

The following procedure is followed for all experiments. The droplet is dispensed using a handheld pipette at room temperature and the enclosure is secured. The thermoelectric modules are switched on and thermal images are continuously acquired as described above. For optical images, the camera is triggered 2 s before recalescence (image-based trigger available on Phantom cameras).

Problem formulation and solution

The supercooling stage (first stage) for the sessile droplet in contact with a cold substrate is discussed here. The transient two-dimensional heat conduction equation is used, with the mathematical formulation written in cylindrical coordinates (for a deformed cylinder with external radius variable with the axial coordinate; see figure 3) with proper initial and boundary conditions given by:

$$[eqn]$$

Configuration of a droplet supported on a cold surface. Mathematical model in cylindrical coordinates–spheroidal solids.

Initial condition:

$$[eqn]$$

Boundary conditions:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where L1 is the liquid droplet base diameter, L2 is the liquid droplet height, $[eqn]$ (m² s^−1^) is the thermal diffusivity of water, T(r,z,t) is the droplet temperature, t is the time, r is the radial coordinate, z is the axial coordinate, T0 is the initial droplet temperature, k (W m^−1^K^−1^) is the liquid thermal conductivity, hc (W m^−^²K^−1^) is the heat transfer coefficient at the substrate–droplet contact, Tb (K) is the substrate temperature, h (W m^−^²K^−1^) is the heat transfer coefficient at the air–droplet interface, hm (m s^−1^) is the mass transfer coefficient at the air–droplet interface, Le (J kg^−1^) is the latent heat of evaporation, $[eqn]$ (kg m^−1^) is the specific mass of saturated water vapour at the droplet surface, $[eqn]$ (kg m^−1^) is the specific mass of saturated water vapour at the surrounding air temperature, $[eqn]$ is the droplet surface emissivity, $[eqn]$ (W m^−1^K^−1^) is the Stefan–Boltzmann constant, $[eqn]$ (K) is the air temperature, nr and nz are the direction cosines of the unit normal vector with the r and z coordinate axes. Average external heat transfer and mass transfer coefficients were assumed at the air–droplet interface. The typically low values of Biot number (Bi) < 0.1 and the small sizes of the droplet favour the employment of the lumping procedure here proposed, due to the mild temperature variations expected within the droplet. For the same reason, spatial variations inherent to the local heat and mass transfer coefficients are not expected to markedly affect the model’s accuracy. Nevertheless, the correlations for the average coefficients here adopted carry information on the droplet surface temperature, which varies along the droplet height.

Murphy & Kopp [62] provide a comprehensive literature review of correlations for saturated water vapour pressure as a function of temperature. Among the reviewed correlations, the one proposed by Bohren & Albrecht [63] was selected for this study. Consequently, the equations for $[eqn]$ and $[eqn]$ , in terms of the temperature for liquid water are presented below, where RH represents the relative humidity in the air:

$$[eqn]$$ $$[eqn]$$

For the calculation of the convective heat $[eqn]$ and mass $[eqn]$ transfer coefficients, correlations for the Nusselt number were taken from Ozisik [64] and for the Sherwood number the equivalent expression was adopted using the Schmidt number instead of Prandtl, as follows.

$$[eqn]$$ $$[eqn]$$

The following dimensionless parameters are used:

$$[eqn]$$

The resulting dimensionless formulation is then written as:

$$[eqn]$$

Initial condition

$$[eqn]$$

Boundary conditions

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

For the supercooling stage, to simplify the manipulation that follows, the radially averaged dimensionless temperature, $[eqn]$ is written in terms of the auxiliary dependent variable, $[eqn]$ , as:

$$[eqn]$$

Thus, the radially averaging operator $[eqn]$ is applied over equation (3.5a):

$$[eqn]$$

Each term in equation (3.7) is evaluated as follows:

• first term

$$[eqn]$$

• second term

$$[eqn]$$

• third term

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

So,

$$[eqn]$$

Replacing equations (3.8a), (3.8b) and (3.8f), into equation (3.7), it results in:

$$[eqn]$$

which is reorganized as

$$[eqn]$$

The last term of equation (3.10) is obtained through the boundary condition, equation (3.5f), which can be rearranged as follows:

$$[eqn]$$

Applying equation (3.11) into equation (3.10):

$$[eqn]$$

and manipulating the above expression,

$$[eqn]$$

However, equation (3.12a) has two dependent variables $[eqn]$ and $[eqn]$ , besides the derivative of $[eqn]$ in the Z variable. Therefore, it is required to relate the external surface temperatures to the radially integrated temperatures, so that one of them can be eliminated from equation (3.12b). To find a relation between $[eqn]$ and $[eqn]$ , we will then use Hermite integration formulae to approximate the radially integrated temperatures and heat fluxes, thus employing the so-called CIEA [52,54–58]. Two Hermite approximations shall be here adopted, the H_0,0_∣H_0,0_ and H_1,1_∣H_0,0_ formulae, where H_0,0_ and H_1,1_ correspond, respectively, to the trapezoidal and corrected trapezoidal integration formulae [52,54–58].

For the radially integrated temperature, the Hermite integration formulae, $[eqn]$ and $[eqn]$ , yield:

$$[eqn]$$ $$[eqn]$$

And, for the radially integrated flux, only the $[eqn]$ formula will be employed:

$$[eqn]$$

The boundary conditions in the R-coordinate are represented by

$$[eqn]$$ $$[eqn]$$

Applying the boundary conditions equations (3.14a) and (3.14b) into equations (3.13a) and (3.13b), one finds:

$$[eqn]$$ $$[eqn]$$

To find the dimensionless temperature at R = 0, we use the boundary condition, equations (3.14a) and (3.14b), into equation (3.13b), to find the expression for the central temperature in terms of the surface temperature:

$$[eqn]$$

In this way, one obtains the following lumped-differential formulation using the H_0,0_∣H_0,0_ approximation:

$$[eqn]$$ $$[eqn]$$

with the following initial condition and boundary conditions:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Similarly, the following lumped-differential formulation is achieved using the H_1,1_∣H_0,0_ approximation:

$$[eqn]$$ $$[eqn]$$

with the following initial and boundary conditions:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The derivation of the reduced models and the numerical solution of the resulting PDEs as obtained by the CIEA, for both the H_1,1_∣H_0,0_ and H_0,0_∣H_0,0_ approximations, are obtained through a symbolic-numerical code built on the Wolfram Mathematica® v.13.3 platform.

Results and discussion

Comparison of experimental and theoretical results

(a)

The reduced models, using the H_1,1_∣H_0,0_ and H_0,0_∣H_0,0_ approximations, were symbolically derived and numerically solved with the NDSolve function in the Wolfram Mathematica® v.13.3 platform and compared with the experimental data collected in this study during the supercooling stage of the sessile droplets. To perform this comparison, the input data were obtained from the experiments, complemented by physical data from the literature, as presented in table 2.

The analysis was conducted for the three different substrates (HYB, SHB and HYL), detailed in §2 and at three distinct positions along the droplet height: T1 (located at 3/4 of the initial droplet height, h0), T2 (at 1/2 h0) and T3 (at 1/4 h0), as shown in figures 4–6.

Temperature variation at the sessile droplet surface on SHB substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). L1=2.11mm, L2=2.08mm,  Tb=−19.7, T∞=−2.8 ∘C and Tw(Z,t) at the droplet surface.

Temperature variation at the sessile droplet surface on HYB substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). . L1=2.87mm, L2=1.9mm,  Tb=−18.46∘C, T∞=−2.9 ∘C and Tw(Z,t) at the droplet surface.

Temperature variation at the sessile droplet surface on HYL substrate for cooling stage: comparison between theoretical (CIEA: H0,0∣H0,0 and H1,1∣H0,0) and experimental results: (a) T1 (at 3/4 h0), (b) T2 (at 1/2 h0), and (c) T3 (at 1/4 h0). L1=3.63mm, L2=1.55mm,  Tb=−15.2∘C, T∞=−1.04 ∘C and Tw(Z,t) at the droplet surface.

It is observed that the results obtained using the H_0,0_∣H_0,0_ approximation, which essentially reduces to the classical lumped system analysis, show a significant deviation from the experimental data. This occurs because this approximation assumes that the droplet temperature along the radial coordinate is practically uniform, or that the surface temperature is essentially equal to the radial average, which appears to be a poor approximation in this problem. On the other hand, the H_1,1_∣H_0,0_, approximation, which incorporates more information from the boundary condition at the droplet surface (see equation (3.15b)), thus rewriting the surface temperature as a function of the radially averaged temperature, shows excellent agreement with the experimental results. The differences observed between the simulated and measured values remain within the 1°C uncertainty of the measurements (infrared camera), validating the theoretical model in predicting the thermal behaviour of the droplet, regardless of the wettability of the substrate (hydrophilic HYL, hydrophobic HYB and superhydrophobic SHB substrates).

Additionally, figures 4–6 show that the droplets start the supercooling stage at an initial temperature close to 0°C and, over time, cool down until they reach the nucleation temperature. This point marks the onset of the recalescence stage, when latent heat is released and then freezing starts, temporarily halting the temperature drop. The nucleation times, which vary depending on the substrate, were used as a stopping criterion in the simulations: 716 s for the SHB substrate, 626 s for HYB and 422 s for HYL. The above comparisons provide a co-validation of the experimental and theoretical results and encourage pursuing further analysis towards the simulation of the freezing and final cooling stages.

Temperature variation during nucleation and freezing of the sessile droplet

(b)

Using the IR thermography videos, the surface temperature history (Tw(Z,t)) at three different planes (T1_exp_, T2_exp_ and T3_exp_) of the droplet is plotted for different substrate wettability, based on a single experiment. Repeating this experiment is highly challenging due to difficulties in controlling the initial droplet and ambient conditions figure 7, in addition to the probabilistic nature of the freezing onset itself. Therefore, the different experimental runs cannot be considered actual replicas of the experiment. Nonetheless, three experimental runs were conducted, and deviations in the temperature readings in individual runs with respect to the average value during the cooling stage were observed, varying between 0.2 and 0.8°C, depending on the droplet shape. Observed variation in Tw for the three substrates (figure 7a–c) is expected due to dependence of heat transfer on the droplet geometry and the relative distance of the temperature plane from the substrate. This variation is accurately modelled, for the supercooling stage of the droplet as shown in the preceding section (see figures 4–6). Additional factors also play a role, such as variations in radiative transfer view factor [66] due to the droplet’s curvature (see §2), being different for hydrophobic and hydrophilic droplets. Variations in the ambient temperature also contribute to the temperature-dependent convection and radiation exchange with the environment. These factors eventually affect the onset of the recalescence stage (figure 7) when the droplet’s temperature rises sharply to 0.6–0.8°C, higher than the expected equilibrium temperature of 0°C. Clearly, the onset of freezing is significantly earlier for HYL than either HYB or SHB due to the reduced energy barrier of nucleation and is consistent with the literature [36].

Evolution of the droplet surface temperature, Tw, on three different substrates where SHB is superhydrophobic (plot marker circle), HYB is hydrophobic (plot marker triangle) and HYL is hydrophilic (plot marker diamond). Three locations are chosen for each of the substrates. They are (a) top plane (T1exp), (b) middle plane (T2exp) and (c) bottom plane (T3exp). The subscript ‘exp’ is for the experimental origin of the plotted data.

Next, the dynamics of recalescence and solidification are better elucidated by measuring the temperature variation at different heights of the same droplet (SHB), as shown in figure 8a. The rise in droplet temperature is first observed in T3_exp_ (bottom most location). This is expected due to the higher probability of heterogeneous nucleation triggering at the droplet–substrate interface [29]. The temperature rise for T2_exp_ and T1_exp_ are 5 and 10 ms later (not visible on the plot), respectively. Given the known distance between T1_exp_, T2_exp_ and T3_exp_, a rough estimate of the average recalescence speed is found to be O (10^2^–10^1^) mm s^−1^, slightly higher than previously reported values [25,34]. The overestimation is possibly due to the limited temporal resolution of the IR imaging.

Surface temperature variation for the sessile droplet on SHB substrate. (a) Temperature variation at three different locations (T1, T2, T3) for the entire duration of cooling and freezing stages (recalescence is indicated in the dashed bo); (b) dashed box is magnified and replotted.

In the second stage of freezing, the temperature at any given location in the droplet remains constant till the freezing front moves past this location, as shown in figure 8b. In case of the SHB surface shown above, the onset of temperature fall at T1_exp_, T2_exp_ and T3_exp_ is nearly 3, 7 and 11 s post-recalescence, respectively. The average speed of the icing front is approximately O (10^−1^) mm s^−1^. However, the solidification rate is not constant throughout the stage as discussed in the next section.

Dynamics of ice growth

(c)

High-speed optical images are used to track the ice front (see §2). The images are spatially band-pass filtered using the FFT plugin in ImageJ (open source image processing software) to remove regions of extreme contrast and improve clarity. Figure 9a–c show the growing ice front at different time instants. The final morphology of the droplet is characterized by a pointy tip at the apex position, irrespective of the substrate wettability. This is caused by expansion of water on freezing and the curving of the ice front near the droplet’s edge, not visible here due to limitations of the imaging system. Marin et al. [32,33] used a Hele-Shaw arrangement to visually confirm the front’s curving and argued that the front itself is a part of a larger circle centred at the apex. As the front progresses, the radius of the front (ice–liquid contact) reduces and finally vanishes into a singular tip at the apex.

Snapshots showing solidification for substrates with different wettability: (a) HYL (1, 5, and 7 s), (b) HYB (1, 6, and 10 s), and (c) SHB (1, 7 and 15 s). The arrow indicates the position and direction of ice front development. Time stamp is inset in the top-left corner. Scale bar equals 1 mm.

The freezing front development for the three substrates is plotted in figure 10, for a single experiment, since it is extremely difficult to perform repetitions of this experiment due to the challenges in controlling the shape of the droplet and initial and control conditions, as previously mentioned (§4b). However, three experimental runs were conducted, and the front position deviations of the individual runs vary from 0.02 to 0.2 mm, calculated relative to the mean value of the experiments conducted. Since the temperature difference between the freezing front and solid substrate remains roughly constant as the height of solid mass increases, the temperature gradient decreases due to the increase in thermal resistance offered by the thicker solid layer, which is however counterbalanced by the decreasing solid–liquid contact area after the largest radius cross section. Since the solidification rate is governed mainly by heat dissipation to the substrate, the front speed slows down in the first stages, as observed in figure 10, followed by a more or less constant value in the intermediate stage, and then markedly accelerating towards the final stages, due to the vanishing water–ice interface [24] and thus vanishing liquid water volume. The hydrophilic droplets have a slightly higher rate of ice progression compared with the hydrophobic cases. This is consistent with the numerical model developed by Zhang et al. [31].

Position of the ice front (S) relative to the substrate as a function of time for experiments on three different substrates.

Conclusions

Supercooling and freezing of sessile droplets deposited over colder substrates have been studied both experimentally and theoretically. Non-intrusive temperature measurements at the droplet surface by IR thermography, together with high-speed camera images, have provided valuable data for the analysis of the heat transfer process and the detailed dynamics of the recalescence phenomena, of the solidification front movement, and of the droplet shape changes. The study also explored the three distinct ranges of substrate wettability, namely hydrophilic, hydrophobic and superhydrophobic, to investigate the differences in the phase change process. The onset of freezing was observed earlier on the hydrophilic substrate due to lower energy barrier to nucleation. The speed of recalescence is O (10^2^–10^1^) mm s^−1^ while the speed of the solidification is O (10^−1^) mm s^−1^. The analysis also provided a co-validation of measurements and simulations from a proposed reduced model for the supercooling stage. A problem reformulation strategy was employed known as the CIEA, which is essentially an improved lumping procedure $[eqn]$ that approximates averaged temperatures and heat fluxes using Hermite formulae for integrals, instead of just the classical lumped system assumption of uniform radial temperature distributions $[eqn]$ . The latter model leads to significant error since the assumption of a uniform radial temperature is physically inadequate. The present work shall form the basis for the proposition of more complete models that account for the recalescence and freezing stages, as well as for the droplet deformation along the process. Experiments should now proceed to different liquids, suspensions and droplet sizes and shapes, as well as to further understanding the interactions with the treated substrate.