Numerical Investigation of Spray Impingement Heat Transfer in the Film Boiling Regime

1. Introduction

Spray cooling of hot surfaces plays a significant role in numerous industrial processes. Notable applications include the quenching of metal components [1], the use of cryogenic cutting fluids in machining operations [2], cryogenic spray systems, and the thermal management of high-power electronic devices. In all these contexts, the rate of heat transfer is of paramount importance: the resulting microstructure of a material is governed by how rapidly it cools, tool longevity is closely tied to cutting temperatures, and the effectiveness of cryogenic machining hinges largely on how efficiently heat is dissipated from the cutting zone ([3,4]). The subject is of considerable interest also in the context of internal combustion engines, encompassing fuel injection [5] as well as SCR after-treatment systems, in which the interaction between urea injection and the heated surfaces of mixers [6,7] or wall pipes [8] plays a critical role.

The heat exchanged by the spray depends on the surface temperature and droplet momentum [9]. Below the boiling point, a fraction of the droplet mass typically deposits on the wall, forming a liquid film. Conversely, when the surface temperature significantly exceeds the boiling point, reaching the Leidenfrost temperature, a vapour layer forms between the droplet and the surface [10], preventing wall contact and liquid film formation. This regime, known as film boiling, exhibits a lower heat transfer coefficient compared to conditions where a liquid film is present, such as conduction and nucleate boiling ([11], p. 472). Additionally, droplet momentum governs the post-impingement behaviour, with high droplet velocities promoting droplet break-up.

A substantial body of literature exists on the experimental characterisation of individual droplet kinematics [12] and heat exchange [13] in the film boiling regime, and a range of heat transfer correlations have been proposed for both single droplets and sprays impinging on flat surfaces [14]. Nevertheless, acquiring detailed experimental data for complex, real-world systems remains a considerable challenge. Numerical modelling offers a valuable alternative in this regard, enhancing the understanding of such processes through spray simulation and its capacity to reproduce impinging spray behaviour [15]. To date, only a limited number of studies have addressed the thermal aspects of spray impingement through simulation. A CFD model capable of simultaneously capturing spray dynamics, heat transfer at the solid surface, and heat propagation within the solid would enable the comprehensive simulation of an entire spray cooling system. Such a framework enables the systematic investigation of spray characteristics, geometry, and operating conditions, along with their mutual interactions, thereby supporting the optimisation of the technology under consideration.

For a complete numerical description of sprays in cooling applications, accurately capturing the spray evolution must be followed by an appropriate modelling of its interaction with the wall surface. Accordingly, this work aims to investigate this aspect.

In this work, a computational conjugate heat transfer (CHT) CFD model for spray–wall interaction above the Leidenfrost point is proposed. The spray is represented using an Eulerian–Lagrangian approach. The CHT framework enables real-time coupling between the fluid and solid domains by solving the heat transfer problem across the interface surfaces.

Lagrangian models to simulate the thermodynamic interaction between the droplets and the wall have been implemented. In particular, two different models for heat exchange in the film boiling regime have been adopted: those proposed by Deb [14] and Breitenbach [16]. Whereas the dynamic interaction between droplets and the wall was modelled following the approach proposed by Kuhnke [17]. A new coupling condition based on the Dirichlet-Neumann coupling for thermally coupling the solid, the spray and the air regions is also proposed.

The consistency of the developed model with respect to mass and energy conservation is first discussed, followed by an experimental validation of the heat transfer coefficient against the experimental data reported by Wendelstorf et al. [18]. This experimental database is based on a series of tests involving orthogonal water spray impingement. The choice of this work for model validation stems from the large number of tests performed, which enables an extensive assessment of the model under varying spray operating conditions. The numerical model was found to accurately estimate the total heat transfer. The variation of the heat transfer coefficient (HTC) with mass flux is well captured, with only a minor deviation observed for very dense sprays. Furthermore, both the Deb and Breitenbach models exhibit sensitivity to wall surface temperature and successfully predict the increase in heat transfer within the transition region.

2. Materials and Methods

The spray impingement model combines an Eulerian–Lagrangian approach for simulating the spray evolution in air with a Conjugate Heat Transfer (CHT) methodology to capture the heat exchange between the solid and fluid regions [19], along with an impingement model to assess the thermodynamic phenomena occurring during droplet impact on the wall surface. The model was implemented within the OpenFOAM^® framework

The adoption of a Lagrangian approach to modelling sprays eliminates the need to fully resolve the nozzle region within the computational mesh. In this framework, the spray is represented by discrete points, commonly referred to as parcels, which are representative of a group of droplets with the same thermodynamic properties. Each point can be assigned a variety of properties, including location, velocity, diameter, mass, temperature, and fuel composition. These points purely act as markers; hence, they are zero-dimensional, and they occupy no physical space in the computational domain. The Lagrangian parcels are tracked as they move through the domain, transitioning from cell to cell according to the governing physical laws. The coupling with the continuous phase is achieved by introducing source terms in the Eulerian transport equation to conserve the global system energy. The flow diagram shown in Figure 1 illustrates the steps of this model. This section describes the different aspects of this modelling strategy.

2.1. Droplet Thermo-Dynamical Impingement

An impingement model was developed for the droplets that impact the solid wall. If a droplet is close to the wall, this model is applied to determine the post-interaction status (dynamic model) and the heat transfer (heat transfer model). The implemented methodology for the dynamic wall interaction is based on the Kuhnke model [17,20]. Two alternative formulations for the heat exchanged between the droplet and the wall were implemented. A surface field stores the heat exchanged by all the impacting droplets in the considered time step, which will be used by the coupling boundary condition to actually perform the heat transfer between the regions.

2.1.1. Droplet-Wall Dynamic Impingement

The Kuhnke model is a semi-empirical framework for droplet–wall interactions, which categorises the possible outcomes into four main regimes: deposition, splash, rebound, and thermal breakup.

• Deposition occurs at low wall temperatures and low impact velocities and corresponds to the condition in which the droplet adheres to the surface, contributing to the formation of a liquid film;
• Splash occurs at higher impact velocities but relatively low wall temperatures. In this regime, part of the droplet mass is deposited onto the liquid film, while another portion, due to the high kinetic energy of impact, escapes from the film in the form of secondary small droplets;
• Rebound happens at low velocities and high temperatures and is basically a quasi-elastic interaction between the droplet and the wall;
• Thermal break-up happens at high velocities and temperature and represents a disintegration of the impingement droplet into secondary small droplets without film formation.

The classification of the impingement regimes is therefore based on two non-dimensional parameters: one representing the thermal characteristics ($T^{*}$) and the other the dynamic characteristics (K) of the specific impingement phenomenon, as defined in Equation (1):

$$\begin{array}{cc}\hfill T^{*}& =T_w{T}_{sat}^{-1},\hfill \\ \hfill K& =W{e}_n^{1/8}L{a}^{5/8}.\hfill \end{array}$$

(1)

Where $T_w$ is the surface temperature, whereas ($W{e}_n$), is the normal Weber number at the impingement state, hence defined with respect the droplet velocity component normal to the surface and ($La$) the Laplace number, representing the ratio of surface tension to the momentum-transport inside the droplet ($La=\sigma \rho D/\mu$). Kuhnke defines a deposition limit ($T_{crit}^{*}$) beyond which no liquid film formation occurs, and consequently, neither deposition nor splashing can take place. Following the assumption proposed by Birkhold et al., this limit corresponds to the Leidenfrost temperature [17]. Its value typically ranges from 1.1 [21] to 1.43–1.47 [17] times the liquid saturation temperature. Since the regime of interest in this work is above the Leidenfrost temperature, the developed model does not include a liquid film sub-model, since it would never be formed, and its inclusion would complicate the simulation with no benefit. Its validity is therefore limited to $T^{*}>T_{crit}^{*}$ and only rebound and thermal break-up interactions are considered.

The two impingement regimes of interest are illustrated in Figure 2.

Since both regimes occur above the Leidenfrost temperature, the determining property for their classification is the droplet kinetic energy. A critical value, $K_{crit}$, distinguishes between these two regimes. This parameter is assumed to be independent of surface roughness and is defined as a uniform probability function in the range of 20 to 40, to account for the inherent complexity and variability of the phenomenon. In order to fully describe the post-impingement state of the droplets, it is necessary to define both the velocity vector and the droplet diameter. Since no deposition occurs, the droplet mass is conserved during the interaction (apart from the effect of heat transfer with the wall). The following impingement properties must therefore be specified to define the dynamic state of the post-impingement droplet ($drople{t}_1$) relative to the incident droplet ($drople{t}_0$):

• Size rate:

  $$\gamma ={\displaystyle \frac{d_1}{d_0}}$$
• Velocity vector, hence magnitude ($|{\mathbf{U}}_1|$), ejection angle ($\beta$) and deviation angle ($\phi$).

2.1.1.1. Rebound

For the rebound interaction, the post-impingement state is straightforward to define, as no secondary droplets are generated; consequently, $\gamma =1$. A factor $U_f<1$ is introduced to account for the non-perfectly elastic nature of the impact and the resulting momentum loss during impingement. This factor depends on the impinging Weber number and is defined as follows:

$$U_f=0.678exp(-4.415·{10}^{-2}W{e}_0).$$

(2)

The exit velocity is hence:

$${\mathbf{U}}_1=(1-2·U_f)U_{n,0}{\mathbf{n}}_{\mathbf{w}}+U_{t,0}{\mathbf{t}}_{\mathbf{w}},$$

(3)

where ${\mathbf{n}}_{\mathbf{w}}$ and ${\mathbf{t}}_{\mathbf{w}}$ are the wall normal and tangential unit vectors and $U_{n,0}$ and $U_{t,0}$ the corresponding normal and tangential components of the impinging droplet.

2.1.1.2. Thermal break-up

In the thermal break-up interaction, the size ratio is correlated to the Weber number and the impingement angle ($\alpha =atan\left({\displaystyle \frac{U_{n,0}}{U_{t,0}}}\right)$) as follows:

$$\overline{\gamma }=3.3e^{3.6{\left({\displaystyle \frac{\alpha }{\pi }}\right)}^2}W{e}^{-0.65}$$

(4)

Three parcels are initialised for each impingement event. To assign the size ratio to each of them, a Weibull distribution is employed, centred around the average value $\overline{\gamma }$. The two parameters defining the distribution, the scale and shape parameters, are set to $\overline{\gamma }$ and 2, respectively. The Weber number of the secondary droplets ($W{e}_1$) depends on both the impingement Weber number and the impingement angle. For its determination, the following correlation is employed:

$$W{e}_1=\gamma \left(W{e}_0(1-0.85{sin}^2\alpha )+12\right)-6$$

(5)

The magnitude of the ejected velocity is therefore derived from the computed Weber number:

$$U_1=\sqrt{{\displaystyle \frac{\sigma W{e}_1}{\rho d_{i,1}}}}$$

(6)

The average ejection angle ($\overline{\beta }$) is correlated as follows to the injection angle and the Weber number:

$$\overline{\beta }=0.96\alpha exp(-0.0045W{e}_0)$$

(7)

A logistic distribution is then used to define the angle for each secondary droplet. The deviation angle ($\phi$) is evaluated by the following correlation:

$$\psi =-{\displaystyle \frac{\pi }{\omega }}ln(1-p{e}^{-\omega })\mathrm{with}\omega =\left\{\begin{array}{cc}\sqrt{{\displaystyle \frac{1+8.872cos1.152\alpha }{1-cos\alpha }}},\hfill & \alpha \le 80\mathrm{deg}\hfill \\ {\displaystyle \frac{{\pi }^2}{2}}cos\alpha ,\hfill & \alpha >80\mathrm{deg}.\hfill \end{array}\right..$$

(8)

The mass conservation is achieved by relating the number of particles of the secondary droplet to the cube of the size ratio

$$n_1={\displaystyle \frac{n_0}{{\gamma }^3}}$$

(9)

2.1.2. Droplet-Wall Heat Transfer

Once the dynamic behaviour of the impinging droplets has been defined, it is necessary to quantify the heat exchange occurring during impingement. Two alternative models were implemented to evaluate the heat transfer between the wall and the droplets. The Breitenbach model is based on a theoretical description of the orthogonal impingement of a single droplet, whereas the Deb model is an empirically derived correlation. The Breitenbach model is discussed in greater detail below to provide a clearer understanding of the underlying physical mechanisms governing this phenomenon.

2.1.2.1. Breitenbach Model

The model proposed by Jan Breitenbach et al. [16] is based on a theoretical analysis of the emergence and expansion of the thermal boundary layers in an orthogonal impingement of a single droplet. The heat transfer during the spreading of the drop is analysed considering the thermal interaction with the solid substrate and the creation of the vapour layer.

The model is developed based on the hypothesis of a droplet impacting a flat solid surface in the fully developed film boiling regime (Leidenfrost interaction). Under this assumption, a vapour film forms between the liquid and the solid surface due to rapid evaporation. A one-dimensional approximation in the direction orthogonal to the surface is employed, since the film thickness is much smaller than the droplet diameter.

The energy balance of the liquid-vapour layer-solid system, represented in Figure 3, is written as:

$${\dot{q}}_1={\dot{q}}_2+{\rho }_l{h}_{lv}{\displaystyle \frac{dh}{dt}},$$

(10)

where ${\dot{q}}_1$ is the heat exchanged from the solid to the vapour, which equals the one exchanged from the vapour to the liquid ${\dot{q}}_2$ plus the phase change contribution, defined as the product of the time derivative of the vapour layer height (${\displaystyle \frac{dh}{dt}}$) and the latent heat of evaporation ${\rho }_l{h}_{lv}{\displaystyle \frac{dh}{dt}}$.

The solution of the energy balance equation in the spreading droplet is written as:

$$T_l(\theta ,t)=T_{\mathrm{sat}}+(T_{D_0}-T_{\mathrm{sat}})\mathrm{erf}\left({\displaystyle \frac{\sqrt{5}\theta }{2\sqrt{\alpha t}}}\right),$$

(11)

where z is the orthogonal coordinate to the surface and $\alpha$ is the liquid thermal diffusivity. Note that at $z=0$ the liquid temperature is the saturation temperature since it is at the vapour interface, and hence it is where the phase change occurs. It is now possible to define the heat flux at the liquid/vapour interface as:

$${\dot{q}}_2\left(t\right)\equiv {\lambda }_f{\left.{\displaystyle \frac{\partial T_f}{\partial \theta }}\right|}_{\theta =h}={\displaystyle \frac{\sqrt{5}{e}_f(T_{\mathrm{sat}}-T_{D_0})}{\sqrt{\pi t}}},$$

(12)

where $e_f$ is the liquid thermal effusivity. Similarly, the heat flux at the solid-vapour interface is:

$${\dot{q}}_l\left(t\right)\equiv {\lambda }_w{\left.{\displaystyle \frac{\partial T_w}{\partial \theta }}\right|}_{\theta =0}={\displaystyle \frac{e_w(T_0-T_c)}{\sqrt{\pi t}}},$$

(13)

During the initial spreading of the droplet, the pressure within the liquid is dominated by inertial forces. In this stage, the formation of the vapour film is primarily governed by heat transfer from the liquid and solid substrates.

Since heat transfer between the wall and the droplet is inversely proportional to the vapour layer thickness, most heat is transferred during the early stage of impact, when the vapour thickness is minimum. The heat flux density of an impacting single droplet in the film boiling regime can be calculated as:

$${\dot{q}}_1\left(t\right)={\displaystyle \frac{2G{e}_w(T_0-T_{\mathrm{sat}})}{\sqrt{\pi }(K+2G)\sqrt{t}}},$$

(14)

where $e_w$ is the solid effusivity and K and G are two parameters function of the system’s physical properties. The total heat removed by a single impacting drop during the initial stage $t_i$ is hence the time integral of the product between mono-dimensional heat flux and the impingement area, and, after various mathematical steps, takes the following form Equation (15):

$$\begin{array}{c}\hfill Q_d={\displaystyle \frac{4.63G{D}_d^{5/2}{e}_s(T_w-T_{sat})}{\sqrt{U_{d,n}}(K+2G)}},\end{array}$$

(15)

The hypothesis of an orthogonal impingement of a single droplet is hardly ever the case for spray droplets. To generalise the model to a random impingement, only the magnitude of the velocity component normal to the surface has been considered in the formulation. This is a common choice in literature ([12,13,22,23]), and it is made since the heat transfer is determined by the spreading of the droplet on the surface and by the contact time. Both these aspects are driven by the normal momentum.

2.1.2.2. Deb model

The semi-empirical correlation proposed by Deb [14] is based on experimental data available in the literature for the impingement of a stream of droplets. This model does not directly define the heat exchanged, it is instead based on the definition of the droplet heat transfer effectiveness ($ϵ$). This parameter represents the ratio between the actual heat transferred to the wall and the total heat stored in the droplet (Equation (16)):

$$\begin{array}{c}\hfill ϵ={\displaystyle \frac{Q_d}{m_d[\Delta h+C_{p,l}(T_{sat}-T_l)]}},\end{array}$$

(16)

where $m_d$ is the droplet mass, $\Delta h$ is the latent heat of vaporization and $Cp$ the heat capacity. To quantify $ϵ$, dimensional analysis was applied, and the following non-dimensional parameters are defined:

$$\begin{array}{cc}\hfill W{e}_n& =\rho d{U}_{d,n}/\sigma B=C_{p,v}(T_w-T_{sat})/\Delta h\hfill \\ \hfill K_d& ={\displaystyle \frac{k_v}{{\mu }_{v}C{p}_v}}SF={\displaystyle \frac{e_s}{e_{st}}}-1.\hfill \end{array}$$

(17)

A statistical analysis was carried out, and the correlation presented in Equation (18) was obtained:

$$\begin{array}{c}\hfill ϵ=0.027exp\left[{\displaystyle \frac{0.08\sqrt{ln(W{e}_n/35+1)}}{{(B+SF/60.5)}^{1.5}}}\right]+0.21K_{d}Bexp\left[{\displaystyle \frac{-90}{W{e}_n+1}}\right]\end{array}.$$

(18)

Where $W{e}_n$ is the normal Weber number, B is the wall super-heat parameter, $K_d$ is the vapour parameter,$SF$ is the surface material parameter, which is referenced to the steel thermal effusivity $e_{st}=\sqrt{k_{st}C{p}_{st}{\rho }_{st}}$.

The final model for the droplet heat transfer in the spray condition is shown in Equation (19):

$$\begin{array}{c}\hfill Q_{d,spray}=X·Q_d.\end{array}$$

(19)

Here, X is a fitting parameter introduced to align the model predictions with experimental data. Breitenbach [16] introduced this coefficient because the model incorporates several parameters derived from empirical correlations and does not explicitly account for heat exchange during the final stage of droplet spreading. A coefficient value close to unity is considered acceptable by the authors. Conversely, Deb [14] did not include such a coefficient; however, due to the inherent complexity of the impingement phenomenon, significant scatter is observed among different experimental studies. Consequently, a correction coefficient is also applied to the Deb model to fit the predicted heat transfer to experimental data. The magnitude of this coefficient is then used as a reference parameter to assess the validity of both models for the considered cases.

2.2. Multiphase-Solid CHT Coupling

The coupling between the solid and fluid regions is performed within the Conjugate Heat Transfer (CHT) framework. CHT is a simulation methodology that enables the simultaneous analysis of regions governed by different physical phenomena, in this case, fluid and solid domains. The coupling between these regions is achieved through a boundary condition that resolves the heat transfer at their interface. In the case of spray impingement, heat exchange occurs across both the droplet–solid and air–solid interfaces; therefore, a dedicated boundary condition was implemented to accurately account for these combined effects.

At the interface, in the case without spray cooling, two conditions must be satisfied, namely the continuity of temperature and the energy conservation (Fourier Law for heat conduction). The following system of equations can be written:

$$\left\{\begin{array}{cc}{T}_{a,w}\hfill & =T_{s,w}\hfill \\ -k_a{\displaystyle \frac{\partial T_a}{\partial n}}{|}_w\hfill & =-k_s{\displaystyle \frac{\partial T_s}{\partial n}}{|}_w\hfill \end{array}\right..$$

(20)

In which $T_{a,w}$ is the wall temperature on the fluid side and $T_{s,w}$ the wall temperature on the solid side. This system can be solved through two approaches. The first consists of applying the Fourier law on both sides of the interface. The second relies on a Dirichlet-Neumann coupling, where the temperature is imposed on one side and the heat flux on the other. The former represents the standard approach in OpenFOAM, as it requires no user specification regarding the boundary condition assignment on each side. During spray impingement, the heat flux on the surface of each solid boundary face ($q_s^{\prime \prime }$) consists of two contributes: the convective heat transfer with the gaseous region ($q_a^{\prime \prime }$) and the heat exchange between the solid and the impacting droplets ($q_d^{\prime \prime }$) (Figure 4).

Under the hypothesis of a dimensionless dispersed phase, consistent with the Lagrangian framework, the system of equations is modified. While temperature continuity between the Eulerian gaseous phase and the solid region is still enforced, an additional heat transfer contribution on the solid side arising from the impinging droplets must be accounted for:

$$\left\{\begin{array}{cc}{T}_{a,w}\hfill & =T_{s,w}\hfill \\ -k_a{\displaystyle \frac{\partial T_a}{\partial n}}{|}_w+q_d^{\prime \prime }\hfill & =-k_s{\displaystyle \frac{\partial T_s}{\partial n}}{|}_w\hfill \end{array}\right..$$

(21)

This condition renders the standard OpenFOAM procedure of solving the Fourier equation on both sides of the interface unsuitable, since the heat transfer is no longer balanced between the Eulerian gaseous phase and the solid region. As a result, the boundary condition introduced in this work is based on a Dirichlet-Neumann coupling approach, in which the heat flux is enforced on the solid side of the interface while the temperature is imposed on the gaseous side.

2.3. Validation Case

The validation of the impingement model was carried out using a water orthogonal spray impingement case, based on the experimental study conducted by Wendelstorf et al. [18]. This analysis was performed to evaluate the performance of the droplet–wall heat transfer models under varying conditions of mass flux and wall temperature.

The main validation parameter is the surface heat transfer coefficient, defined in Equation (22):

$$\begin{array}{c}\hfill HTC={\displaystyle \frac{q_w^{\prime \prime }}{T_w-T_{spray}}},\end{array}$$

(22)

where $q_w^{\prime \prime }$ and $T_w$ are the surface heat flux and temperature, while $T_{spray}$ is the droplet temperature.

2.3.1. Experimental Set-Up

According to the authors [18], the experimental testes where conduced employing three full-cone nozzles (Spraying Systems VKE6/60, VKE6/90, and VKE8/60) at nozzle-to-sample distances ranging from 62 to 105 mm. The test specimens consisted of 1.0 mm thick discs made from commercially pure nickel with a diameter of 70 mm. To measure the Heat transfer coefficient (HTC), five thermocouple pairs were welded to the underside of each disc, equally spaced starting from the centre and moving radially every 10 mm. The droplet velocity and size distributions were estimated from empirical correlations [24] and, for the investigated spray conditions, the average droplet velocity was found between 13 and 15 m/s, with mean droplet diameters in the range of 300–400 µm. 28 experiments were performed for water mass flux densities ($V_s$) varying from 3 to 30 $\mathrm{kg}/\left({\mathrm{m}}^2\mathrm{s}\right)$ and temperature differences between the wall and the droplet ranging from 150 to 1100 K. The mass-flux density is proportional to the ratio between the mass flow rate and the impinging area ($V_s={\displaystyle \frac{\dot{m_s}}{A_{imp}}}$).

The following analytic correlation for heat transfer dependency from the wall temperature and the mass flux was hence proposed (Equation (23)):

$$\begin{array}{c}\hfill HTC=190+tanh\left({\displaystyle \frac{V_s}{8}}\right)\left(140V_s\left[1-{\displaystyle \frac{V_s\Delta T}{72000}}\right]+3.26\Delta T^2\left[1-tanh\left({\displaystyle \frac{\Delta T}{128}}\right)\right]\right),\end{array}$$

(23)

where $\Delta T=T_w-T_l$ is the difference between the wall surface temperature and spray temperature.

2.3.2. Simulation Set-Up

Nine simulations were performed for validation purposes, varying the wall temperature and mass flow rate. The corresponding operating conditions are summarised in table:validation1.

The spray has been simulated, resorting to a cone injection model. The parcels are introduced in the domain at the nozzle position, with a fixed velocity (14 m/s) and a random direction limited by the characteristic cone angle of the spray (60 deg). Droplet diameters are initialised based on the Rosin-Rammler mass distribution [25] with a minimum of 300 $\mu \mathrm{m}$, a maximum of 400 $\mu \mathrm{m}$ and an average of 350 $\mu \mathrm{m}$.

The domain analysed is a box of 150mmX150mm and 100 mm in the dimension along the spray axis, while the solid domain is 150mmX150mm and 1 mm, as the real disk thickness. The computational mesh is fully hexahedral. The cell size was chosen as a compromise between computational effort and accuracy of the results. Since the validation with experimental value is performed in steady-state conditions, the temperature profile along the normal direction is linear, and therefore substantially insensitive to mesh resolution of the solid region. The mesh configuration consists of 225000 cells for the solid region and 112500 for the fluid region, the mesh resolution is 0.1 mm for the solid region and 2 mm for the fluid one. The dimension of the wall layer in the fluid region is less critical since the majority of the energy transferred is due to the droplet impingement, which is modelled as a marker and hence their properties do not depend directly on the dimension of the mesh.

As concerns the simulation set-up, PIMPLE resolution algorithm has been used, with two internal pressure correctors and two outer correctors. The Euler implicit time scheme was used for time discretisation, and secondary order limited schemes are used for divergence and Laplacian discretisation. Turbulence is treated in RANS fashion, and the well-known k-$\epsilon$ model was employed [15].

Adiabatic boundary conditions were applied to the lateral patches of the solid disk, while the bottom face was maintained at a constant temperature equal to the initial condition. This configuration was adopted to enable direct comparison with experimental data during post-processing. Fixing the bottom temperature drives the system towards a steady-state solution, at which point the HTC is computed. Figure 5 reports the temporal evolution of the average solid surface properties in the impingement region for the case ${\mathrm{V}}_{\mathrm{s}}=15\mathrm{Kg}/{\mathrm{m}}^2\mathrm{s}$ and $\Delta \mathrm{T}=700\mathrm{K}$. The surface temperature reaches steady-state conditions after 0.15 s. Accordingly, the HTC is computed as a time-average from this instant to the end of the simulation, which was fixed at 0.3 s for all investigated cases. The observed fluctuations in the HTC are ascribed to the stochastic variation in the number of droplets impinging on the surface at each time step.

3. Results

The results section is structured into two distinct parts. The first focuses on a detailed characterisation of the thermodynamic state of the system under a reference operating condition, defined as ${\mathrm{V}}_{\mathrm{s}}=15\mathrm{Kg}/{\mathrm{m}}^2\mathrm{s}$ and $\Delta \mathrm{T}=700\mathrm{K}$. The second addresses the validation of the proposed methodology over the entire range of operating conditions considered.

3.1. Thermo-Dynamical Analysis of the Impingement Process

Figure 6 shows the evolution of the dispersed phase during the impingement process.

As is visible, the droplets undergo thermal break-up, as the diameters of the post-impingement droplets are significantly smaller than those of the incoming droplets. This occurs because, for this specific spray configuration, the droplets are characterised by a kinetic non-dimensional number (K) greater than the critical value ($K_{crit}$). Also, the velocity of the post-impingement droplets drastically decreases, both as a direct effect of the energy dissipation at the impingement and due to the lower momentum of the smaller droplets, which leads to a stronger deceleration due to the drag with the air.

Focusing on the heat transfer at the solid-fluid interface, Figure 7 illustrates the evolution of the heat flux within the impingement area under steady-state conditions. With reference to Equation (20), three heat transfer contributions can be distinguished:

• On the solid side, heat is transferred by conduction: $WH{F}_{Solid}=-k_s{\displaystyle \frac{\partial T_s}{\partial n}}{|}_w$;
• on the gaseous side, within the Eulerian domain, heat exchange with the solid takes place through convection: $WH{F}_{Gas}=-k_a{\displaystyle \frac{\partial T_a}{\partial n}}{|}_w$;
• on the Lagrangian side, the heat transfer corresponds to the cumulative contribution of the droplets impinging on the surface at each time step: $WH{F}_{Liquid}=\sum Q_{d,i}/\Delta T$.

As anticipated, direct droplet impingement constitutes the dominant heat transfer mechanism, corroborating the point-particle assumption. Moreover, the sum of $WH{F}_{Gas}$ and $WH{F}_{Liquid}$ is found to equal $WH{F}_{Solid}$, thereby ensuring energy conservation at the solid-fluid interface.

Concerning the phase change of the dispersed phase, evaporation occurs predominantly at the impingement zone, since the spray propagates through a nearly isothermal environment before reaching the wall. Moreover, as a consequence of the Leidenfrost phenomenon, the heat exchanged between the impinging droplets and the wall directly induces phase change, driven by the formation of a vapour cushion at the liquid-solid interface.

Figure 8 confirms this behaviour, showing that the region of highest water mass fraction — corresponding to the evaporated dispersed phase — is concentrated in the cells closest to the wall, while the water mass fraction is nearly negligible further away from the impact zone.

Figure 9 provides a quantitative description of the phase change process, evaluated both in the Eulerian domain — by integrating the product of density and water mass fraction — and in the Lagrangian domain — by computing the difference between the injected and current Lagrangian mass in the domain. The mass transfer rate increases sharply when the first droplets contact the wall, approximately 0.008 s after the start of injection, and continues to rise thereafter. This is partly attributed to post-impingement droplets, which exhibit a smaller diameter and thus a larger surface-to-volume ratio, resulting in an enhanced phase change rate. Additionally, the agreement between the mass transfer rates computed in both domains confirms the mass conservation of the proposed approach.

3.2. Validation

This section presents the results of the simulation campaign conducted across the operating conditions defined in Table 1. For each case, the heat transfer coefficient is computed and its trend is discussed.

The influence of the mass flux on the heat transfer is shown in Figure 10. The simulations demonstrate a satisfactory agreement with the experimental data across most of the investigated conditions. The heat transfer increases sub-linearly with ${\mathrm{V}}_{\mathrm{s}}$, due to the non-lenear dependency of the thermal impingement models (Equations (15) and (18)) on the impact velocity. The largest discrepancies are observed for the densest spray conditions, where coalescence effects become significant, thereby reducing the accuracy of the Lagrangian approximation in the solid interface region.

The effect of wall surface temperature on heat transfer is shown in Figure 11. Both the Deb and Britenbach models are sensitive to wall temperature and show a significant increase in HTC in the transient regime (smaller $\Delta \mathrm{T}$). The Breitenbach model fits the experimental data better in the fully developed film boiling regime, the Deb model in the transient one.

The fitting parameter (X in Equation (19)) was set to 3.2 for the Breitenbach model and 0.3 for the Deb model across all simulated cases. The value adopted for the Breitenbach model is very close to that proposed by the original authors under similar conditions, indicating that the present simulations are consistent with their findings. Conversely, the coefficient required for the Deb model deviates significantly from unity. Nevertheless, this fitting parameter was kept constant across the different operating conditions to highlight the model’s flexibility and robustness under varying system configurations.

Figure 12 compares the distributions of the heat transfer coefficient on the surface for $V_s=15$ kg/(${\mathrm{m}}^2$s) and $\Delta T=700$ K, as obtained from the two models. The resulting profiles are very similar, indicating that both models capture the same overall behaviour. The small local discrepancies can be attributed to the inherent stochasticity in the droplet spatial distribution within the spray.

To assess the validity of the CHT implementation, Figure 13 presents the cooling rate measured 1 mm below the surface along the nozzle axis for the case with a mass flux of 12 kg/(${\mathrm{m}}^2$s), using the Breitenbach model. The agreement between the simulated and experimental curves is good, demonstrating that the CHT impingement model can accurately predict the transient cooling behaviour of a workpiece.

4. Discussion

Spray cooling is used in different industrial applications, and the possibility to have numerical methods to simulate real industrial applications is valuable for their optimization. This work hence focuses on the development of a simulation methodology for estimating the heat transfer of sprays above the Leidenfrost temperature. In particular, a CHT model for the heat transfer of an impinging spray above the Leidenfrost temperature is presented and validated under various operating conditions. For this purpose, a dynamic interaction model [17], along with two models for quantifying the heat exchanged by the droplets (Breitenbach [16] and Deb [14]), are implemented. A coupled boundary condition for the solid–multiphase interface was developed to account for heat exchange between the solid and both the droplets and the surrounding air. For validation, simulations with impinging water were conducted to evaluate model flexibility across varying spray mass fluxes and surface temperatures, and to compare the performance of the implemented droplet heat exchange models.

The experimental data used to validate the model were obtained from the work of Wendelstorf et al. [18]. The thermodynamic analysis confirms that droplet impingement constitutes the dominant heat transfer mechanism, with evaporation occurring predominantly in the near-wall region. Energy conservation at the solid-fluid interface was verified through the consistency between the Eulerian and Lagrangian mass transfer rates. The validation against the experimental correlation of Wendelstorf et al. demonstrates satisfactory agreement across most operating conditions. The heat transfer coefficient increases sub-linearly with mass flux, and both the Deb and Breitenbach models correctly capture the sensitivity to wall temperature, including the enhancement observed in the transition regime. The Breitenbach model exhibits superior performance in the fully developed film boiling regime, while the Deb model provides better predictions in the transitional region. The transient cooling behaviour of the solid is also accurately reproduced, confirming the reliability of the CHT implementation.