Finite-size evaporating droplets in weakly compressible homogeneous shear turbulence

Abstract

We perform interface-resolved simulations of finite-size evaporating droplets in weakly compressible homogeneous shear turbulence. The study is conducted by varying three dimensionless physical parameters: the initial gas temperature over the critical temperature $T_{g,0}/T_c$, the initial droplet diameter over the Kolmogorov scale $d_0/\eta$ and the surface tension, i.e. the shear-based Weber number, $We_{\mathcal {S}}$. For the smallest $We_{\mathcal {S}}$, we first discuss the impact on the evaporation rate of the three thermodynamic models employed to evaluate the gas thermophysical properties: a constant property model and two variable-properties approaches where either the gas density or all the gas properties are allowed to vary. Taking this last approach as reference, the model assuming constant gas properties and evaluated with the ‘1/3’ rule is shown to predict the evaporation rate better than the model where the only variable property is the gas density. Moreover, we observe that the well-known Frössling/Ranz-Marshall correlation underpredicts the Sherwood number at low temperatures, $T_{g,0}/T_c=0.75$. Next, we show that the ratio between the actual evaporation rate in turbulence and the one computed in stagnant conditions is always much higher than one for weakly deformable droplets: it decreases with $T_{g,0}/T_c$ without approaching unity at the highest $T_{g,0}/T_c$ considered. This suggests an evaporation enhancement due to turbulence also in conditions typical of combustion applications. Finally, we examine the overall evaporation rate and the local interfacial mass flux at higher $We_{\mathcal {S}}$, showing a positive correlation between evaporation rate and interfacial curvature, especially at the lowest $T_{g,0}/T_c$.

1. Introduction

The understanding of droplets evaporation in turbulent flows is a crucial aspect in many different contexts, such as geophysics and spray combustion to name a few. In the first example, evaporation is central for the formation and evolution of clouds and, more in general, in many of the grand challenges of environmental fluid mechanics (Dauxois et al. 2021). In the second example, droplet evaporation is a precursor of combustion; ensuring that all the liquid vaporizes before chemical reactions occur is fundamental to minimize the pollutants formation and to maximize the efficiency of the entire process (Birouk & Gökalp 2006). More recently, evaporation in turbulence acquired a prominent role also in our understanding of the fluid dynamics aspects of COVID-19 spreading (e.g. droplet generation due to exhalation and airborne dispersion) as described in Mittal, Ni & Seo (2020), Balachandar et al. (2020), Bourouiba (2021). Therefore, several studies have been conducted to understand how relative humidity affects evaporation and condensation (Rosti et al. 2020; Chong et al. 2021; Ng et al. 2021) as well as on the interaction between droplets and turbulent flows (Bourouiba 2020; Rosti et al. 2021) with the ultimate goal of improving social distancing guidelines.

Historically, evaporating droplets in turbulence have been the subject of numerous and extensive experimental campaigns which have mainly focused on the flow topology inside the droplets (Wong & Lin 1992; Mandal & Bakshi 2012), on the interaction between droplet dispersion and vapour clouds structure (De Rivas & Villermaux 2016; Villermaux et al. 2017; Sahu, Hardalupas & Taylor 2018) and on evaporation enhancement of super-Kolmogorov (i.e. finite-size) droplets (Marti et al. 2017; Verwey & Birouk 2018, 2020). On the computational side, the most common approach has been the point particle method (Kuerten 2016; Maxey 2017) for the studies of sub-Kolmogorov droplets in different flow configurations, i.e. forced homogeneous isotropic turbulence (Mashayek 1998a; Weiss, Meyer & Jenny 2018), homogeneous shear turbulence (Mashayek 1998b; Weiss et al. 2020), turbulent channel flow (Bukhvostova et al. 2014; Russo et al. 2014) and turbulent jets (Reveillon & Demoulin 2007; Wang, Dalla Barba & Picano 2021). The main limitation of this approach is represented by the inherent need of empirical closure equations to model the mass, momentum and energy coupling between the disperse and the continuous phase. This feature limits their rigorous applications only to sub-Kolmogorov droplets (Elghobashi 2019) and hinders the possibility to investigate directly aspects such as the reciprocal influence of nearby droplets (Lupo et al. 2020; Chong et al. 2021) or the evaporation enhancement due to turbulence (which typically occurs for droplets of the order of the Taylor length, Birouk & Gökalp 2006). Moreover, assessing the impact of droplets deformation on the evaporation rate as well as the local interfacial mass flux distribution over the droplet surface is not possible, unless further models are introduced. On the other hand, less work has been dedicated to numerical simulations of droplets larger than the Kolmogorov scale and the only few works available consider homogeneous isotropic turbulence (HIT). Specifically, in Albernaz et al. (2017) the authors studied by means of a hybrid lattice Boltzmann method the deformation and heat transfer of a single droplet with a diameter between $25\eta$ and $40\eta$, with $\eta$ the Kolmogorov length. In this set-up, however, little can be said about the mass transfer because evaporation and condensation compensate for a statistically constant droplet volume. Recently, Dodd et al. (2021) employed a geometric volume of fluid (VoF) method to study the evaporation of droplets at different volume fractions ($0.01\,\% \leq \alpha \leq 1\,\%$) with an initial droplet diameter ranging between $4\eta$ and $17\eta$, in order to highlight the limitations of point particle closures for the calculation of the evaporation rate and the semi-empirical correlations of the Sherwood number in absence of mean flow and for non-isolated droplets.

The limited number of studies and the need of further understanding of such a complex process motivates us to further investigate finite-size evaporating droplets in homogeneous shear turbulence (HST). As noted in Kasbaoui et al. (2017), this configuration can be regarded as one of intermediate complexity between homogeneous isotropic and non-homogeneous flows and it represents a particularly convenient set-up for two reasons. First, it allows us to study shear flows without the complications induced by the presence of the walls. Next, given the intrinsic production of turbulent kinetic energy by the mean shear, turbulence is self-sustained without any external forcing and the flow achieves a statistically stationary condition (SS-HST) (Pumir 1996; Sekimoto, Dong & Jiménez 2016). The SS-HST has been recently considered in multiphase flows, e.g. in Rosti et al. (2019b) for emulsions and in Tanaka (2017), Yousefi, Ardekani & Brandt (2020) for turbulence modulation by rigid particles of different shapes. In both cases, the characteristic length of the disperse phase was chosen larger than the Kolmogorov scale.

In the current study we therefore consider finite-size evaporating droplets in SS-HST, assuming an incompressible liquid surrounded by a compressible gas phase at higher temperature. The initial size, ranging between $10.5\eta$ and $21.5\eta$, is chosen to focus on evaporation enhancement by turbulence (Verwey & Birouk 2020) and to elucidate the effects of the interface deformation. This last aspect is less discussed in literature and, so far, the spherical assumption has been invoked to describe the droplet shape also in fully resolved simulations (Lupo et al. 2019, 2020). Given the large parameter space which characterizes evaporating flows, here we focus on changing the ratio between the droplet initial diameter and the Kolmogorov scale ($d_0/\eta$), the surface tension and the initial gas temperature. Moreover, in all cases, to study the isolated droplet behaviour, we consider a small initial liquid volume fraction, $\alpha _0=0.14\,\%$, corresponding to five droplets. The resulting parametric study aims at addressing the following questions.

1. (a) What is the level of approximation of the estimated evaporation rates when the gas thermophysical properties are assumed constant?

2. (b) What are the effects on the evaporation rate and on the liquid temperature of the droplet size? Moreover, how does the ratio $K/K_0$ (actual evaporation rate in turbulence over the evaporation rate in stagnant conditions) changes when increasing the gas temperature in conditions relevant for combustion applications?

3. (c) How does the interface deformation affect the evaporation rate and does the local interfacial mass flux correlate with changes in the droplet shape?

To investigate this complex phenomenon numerically, we propose a new VoF method for evaporating flows in weakly compressible homogeneously sheared turbulence, extending the algorithm in Scapin, Costa & Brandt (2020). The tool addresses the two main issues arising when performing this kind of simulation with more realistic and challenging conditions. First, as already remarked in Kasbaoui et al. (2017), numerical simulations in HST are demanding even in single phase since the commonly employed multistep time-integration schemes (e.g. Adams–Bashforth and Runge–Kutta), if employed in their classical formulation, are weakly unstable and, therefore, not adequate for long-time simulations. Next, since the HST computational domain does not possess any outflow boundary, a rigorous description of the two-phase evaporating system requires a compressible formulation that allows the thermodynamic pressure to vary with the state variables. To address the first issue, we present a modified version of the Adams–Bashforth scheme which recovers the analytical solution of Kelvin modes in the limit of the rapid distortion theory (RDT) (Maxey 1982) and, overall, ensures a stable integration over long times. The second issue is addressed by deriving and presenting a new mathematical formulation for evaporating flows with phase change in the low-Mach limit (weakly compressible formulation) with a detailed numerical implementation. Differently from other approaches available in the literature (Wang et al. 2019; Ni et al. 2021), this formulation relaxes the assumption of constant thermodynamic pressure and allows its dynamic variation according to the global expansion and contraction of the compressible gas phase, as well as, on the amount of mass flux exchanged at the interface.

This paper is organized as follows. In § 2 we introduce the mathematical model employed to describe a two-phase evaporating system, adapted to the HST configuration. In § 3 the numerical algorithm is presented, with a note on an improved Adams–Bashforth scheme for HST simulations. The results, discussed in § 4, focus on the role played by the specific thermodynamic model used to describe the weakly compressible phase and the effects on the evaporation induced by the variation of the shear-based Reynolds number and the shear-based Weber number. Finally, the main findings and conclusions are summarized in § 5.

2. Governing equations

2.1. Mathematical model for weakly compressible evaporating flows

We consider a system of two immiscible and Newtonian fluids: a single component liquid (phase $1$) and an ideal mixture of an inert gas and vaporized liquid (phase $2$). The two phases are bounded by an infinitesimally small interface, through which mass, momentum and energy are exchanged. The evaporation is driven by the partial pressure of the inert gas in phase $2$. To represent this system, a phase indicator function $H$ is defined at position $\boldsymbol {x}$ and time $t$ to distinguish between the two phases,

(2.1) \begin{equation} H(\boldsymbol{x},t) = \begin{cases} 1 & \text{if }\boldsymbol{x} \in V_l, \\ 0 & \text{if } \boldsymbol{x} \in V_g; \end{cases} \end{equation}

where $V_l$ and $V_g$ are the domains pertaining to the liquid and gas phases, divided by a zero-thickness interface $\varGamma (t)=V_l\bigcap V_g$.

Hereinafter and unless otherwise state, we assume that the liquid is incompressible with constant properties, while the gas is compressible and its properties are allowed to vary with temperature, pressure and composition. Thus, given the possible variation of the density with the state variables, we consider compressibility, which in this work is treated within the low-Mach approximation (Majda & Sethian 1985). This allows us to filter acoustic effects, while still retaining potentially large density variations in the bulk region of the compressible phase. Under this assumption, the conservation equations for species, momentum, thermal energy and mass across the interface read as

(2.2)$$\begin{gather} \rho\frac{{\rm D}\boldsymbol{u}}{{\rm D}t} ={-}\boldsymbol{\nabla} p+\frac{1}{Re}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\tau}+\frac{\boldsymbol{f}_{\sigma}}{We}, \end{gather}$$

(2.3)$$\begin{gather}\rho_g\frac{{\rm D}Y_{l}^{v}}{{\rm D}t} = \frac{1}{ReSc}\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v}) , \end{gather}$$

(2.4)\begin{gather} \rho c_p\frac{{\rm D}T}{{\rm D}t} = \frac{1}{RePr}\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla} T) +\left(\varPi_{p,1}\frac{{\rm d}p_{th}}{{\rm d}t}+\frac{\rho_gD_{lg}}{ReSc} \sum_{j=1}^{2}\boldsymbol{\nabla} h_j\boldsymbol{\cdot}\boldsymbol{\nabla} Y_j\right)(1-H)\nonumber\\ \hspace{-15pc} -\frac{(\dot{m}_{\varGamma}\delta_{\varGamma})}{Ste}, \end{gather}

(2.5)\begin{gather} \dot{m}_{\varGamma} = \frac{1}{ReSc}\frac{\rho_{g,\varGamma}D_{lg,\varGamma}}{1-Y_{l,\varGamma}^v}\boldsymbol{\nabla}_{\varGamma} Y_{l}^{v}\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}.\end{gather}

Here, $\boldsymbol {u}$ is the fluid velocity assumed to be continuous in the two phases, $p$ is the hydrodynamic pressure, $T$ the temperature, $h$ the enthalpy (with $\boldsymbol {\nabla } h=\boldsymbol {\nabla } (c_pT)$), $Y_{l}^v$ the mass fraction of the vaporized liquid in the inert gas and $\dot {m}_{\varGamma }$ is the interfacial mass flux. In (2.2), $\tau$ is the viscous stress tensor for compressible Newtonian flows and $\boldsymbol {f}_{\sigma }=\kappa _{\varGamma }\delta _{\varGamma }$ with $\kappa _{\varGamma }$ the interfacial curvature. The generic thermophysical property $\xi$ (density $\rho$, dynamic viscosity $\mu$, thermal conductivity $k$ or specific heat capacity $c_p$) is computed with an arithmetic average, i.e. $\xi =1+(\lambda _{\xi }-1)H$, where $\lambda _{\xi }=\xi _l/\xi _{g,r}$. Since $\xi _l$ is kept constant and uniform no further modelling is needed, while the generic gas property $\xi _{g}$ is computed with appropriate equation of states. The gas density $\rho _g$ is computed with the ideal gas law and the liquid diffusion coefficient with the Wilke–Lee correlation (Reid, Prausnitz & Poling 1987). More details on how the remaining gas thermophysical properties are evaluated are given in Appendix A. Unless otherwise stated, all the property ratios $\lambda _{\xi }$ are computed with respect to the reference gas property $\xi _{g,r}$, taken as the initial value.

In (2.2), (2.3), (2.5) and (2.4) different dimensionless parameters appear. By introducing a reference velocity $u_{r}$ and reference length $l_{r}$, we define the Reynolds number $Re=\rho _{g,r}u_{r}l_{r}/\mu _{g,r}$, the Weber number $We=\rho _{g,r}u_{r}l_{r}/\sigma$, with $\sigma$ the surface tension; $Sc=\mu _{g,r}/(\rho _{g,r}D_{lg,r})$ and $Pr=\mu _{g,r}c_{pg,r}/k_{g,r}$ are the Schmidt and Prandlt numbers. Note that the temperature equation (2.4) requires the definition of the Stefan number $Ste=c_{pg,r}T_{g,0}/\Delta h_{lv}$, where $T_{g,0}$ is the initial gas temperature and $\Delta h_{lv}$ is the latent heat and of the dimensionless group $\varPi _{p,1}=R_u/(c_{pg,r}M_g)$, where $M_g$ is the molar mass of the gas phase and $R_u$ the universal gas constant. To form a close set of equations, two additional equations are needed, one for the velocity divergence and one for the thermodynamic pressure $p_{th}$, i.e.

(2.6)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} = f_{\varGamma}(\boldsymbol{x}_{\varGamma},t) + \left[\frac{f_{Y}(\boldsymbol{x},t) + f_T(\boldsymbol{x},t)}{p_{th}} - \left(1-\frac{\varPi_{p,1}}{c_p\bar{M}_{m,av}}\right)\frac{1}{p_{th}}\frac{{\rm d}p_{th}}{{\rm d}t}\right](1-H) , \end{gather}$$

(2.7)$$\begin{gather}\frac{1}{p_{th}}\frac{{\rm d}p_{th}}{{\rm d}t}\int_{V_g}\left(1-\frac{\varPi_{p,1}}{c_p\bar{M}_{m,av}}\right)\, {\rm d} V_g =\int_V\left[f_{\varGamma}(\boldsymbol{x}_{\varGamma},t)+ \frac{f_T(\boldsymbol{x},t)+f_Y(\boldsymbol{x},t)}{p_{th}}(1-H)\right] \, {\rm d} V . \end{gather}$$

In (2.6) and (2.7) the functions $f_{\varGamma }(\boldsymbol {x}_{\varGamma },t), f_{Y}(\boldsymbol {x},t)$ and $f_{T}(\boldsymbol {x},t)$ represent the different contributions to the total velocity divergence from the phase change ($f_{\varGamma }$) and the change of the gas density either due to composition ($f_Y$) or to temperature ($f_T$),

(2.8a)$$\begin{gather} f_{\varGamma}(\boldsymbol{x}_{\varGamma},t) = \dot{m}_{\varGamma} \left(\frac{1}{\rho_{g,\varGamma}}-\frac{1}{\lambda_{\rho}}\right)\delta_{\varGamma}, \end{gather}$$

(2.8b)$$\begin{gather}f_{Y}(\boldsymbol{x},t) = \frac{1}{ReSc}\frac{\bar{M}_{m,av}}{\rho_g} \left(\frac{1}{\lambda_M}-1\right)\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v}), \end{gather}$$

(2.8c)$$\begin{gather}f_{T}(\boldsymbol{x},t) =\frac{1}{Re}\frac{\varPi_{p,1}}{c_p\bar{M}_{m,av}} \left[\frac{1}{Pr}\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla} T)+\frac{\rho_gD_{lg}}{Sc}\sum_{j=1}^2\boldsymbol{\nabla} h_j\boldsymbol{\cdot}\boldsymbol{\nabla} Y_j\right], \end{gather}$$

where $\lambda _M=M_l/M_g$ is the molar mass ratio. The complete derivation of (2.6), (2.7) and of relations (2.8) is provided in Appendix B.

2.2. Governing equations for the HST configuration

Equations (2.2), (2.3), (2.5), (2.4), (2.6) and (2.7) are solved in a periodic box assuming an imposed uniform mean shear $\mathcal {S}$, as depicted in figure 1. In a shear-periodic domain the streamwise $x$ and spanwise $y$ directions are periodic, whereas the so-called shear-periodic condition applies in the $z$ direction, which reads for the generic scalar quantity $g$ as

(2.9) \begin{gather} \left.\begin{gathered} g(x+l_x,y,z) = g(x,y,z), \\ g(x,y+l_y,z) = g(x-\mathcal{S}tl_y,y,z), \\ g(x,y,z+l_z) = g(x,y,z). \end{gathered}\right\} \end{gather}

Figure 1. Sketch of the computational box: the rendering represents the volume contour of the vapour mass fraction for the case at $Re_{\mathcal {S}}=6700, We_{\mathcal {S}}=0.1$ and $T_{g,0}/T_c=1.5$.

The presence of a mean velocity, i.e. $\mathcal {S}z$, suggests to decompose the velocity field $\boldsymbol {u}$ into a mean and a fluctuating component $\boldsymbol {u}'$,

(2.10)\begin{equation} \boldsymbol{u}=\boldsymbol{u}'+\mathcal{S}z\boldsymbol{e}_x,\end{equation}

where $\boldsymbol {e}_x=(1,0,0)$. Given the decomposition (2.10), the momentum equation (2.2) is written in terms of $\boldsymbol {u}'$,

(2.11)\begin{align} \rho\left(\frac{{\rm D}\boldsymbol{u}'}{{\rm D}t} + \mathcal{S}z\frac{\partial\boldsymbol{u}'}{\partial x} + \mathcal{S}w'\boldsymbol{e}_x\right) ={-}\boldsymbol{\nabla} p + \frac{1}{Re_{\mathcal{S}}}\left[\boldsymbol{\nabla} \boldsymbol{\cdot}\tau+\mathcal{S}\left(\frac{\partial\mu}{\partial z}\boldsymbol{e}_x+\frac{\partial\mu}{\partial x}\boldsymbol{e}_z\right)\right]+ \frac{\boldsymbol{f}_{\sigma}}{We_{\mathcal{S}}},\end{align}

where $\boldsymbol {e}_{z}=(0,0,1)$. Three new terms appear: the first, $\mathcal {S}z(\partial \boldsymbol {u}'/\partial x)$, is the convection of the velocity fluctuations by the mean shear; the second, $\mathcal {S}w'\boldsymbol {e}_x$, represents the production of streamwise momentum caused by the fluid parcel transport in the normal direction owing to the mean shear; the third and last term, $\mathcal {S}(\partial \mu /\partial z\boldsymbol {e}_x+\partial \mu /\partial x\boldsymbol {e}_z)$, represents the viscous dissipation due to the mean shear in the case of a fluid with variable viscosity. Note that $Re_{\mathcal {S}}$ and $We_{\mathcal {S}}$ are the shear-based Reynolds and Weber numbers, computed taking $u_{r}=\mathcal {S}l_r$.

The same decomposition (2.10) is applied to the mass fraction and temperature equations, giving

(2.12)\begin{gather} \rho_g\left(\frac{{\rm D}Y_{l}^{v}}{{\rm D}t}+\mathcal{S}z\frac{\partial Y_{l}^{v}}{\partial x}\right) = \frac{1}{Re_{\mathcal{S}}Sc}\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v}), \end{gather}

(2.13)\begin{gather} \hspace{-7pc}\rho c_p\left(\frac{{\rm D}T}{{\rm D}t} + \mathcal{S}z\frac{\partial T}{\partial x}\right) = \frac{1}{Re_{\mathcal{S}}Pr}\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla} T) \end{gather}

(2.14)\begin{gather} \quad +\left(\varPi_{p,1}\frac{{\rm d}p_{th}}{{\rm d}t}+\frac{\rho_gD_{lg}}{Re_{\mathcal{S}}Sc}\sum_{j=1}^{2}\boldsymbol{\nabla} h_j\boldsymbol{\cdot}\boldsymbol{\nabla} Y_j\right)(1-H)-\frac{\dot{m}_{\varGamma}\delta_{\varGamma}}{Ste}. \end{gather}

In (2.14) and (2.12) the two new terms, $\mathcal {S}z(\partial T/\partial x)$ and $\mathcal {S}z(\partial Y_{l}^{v}/\partial x)$, represent the convection of the temperature/vapour mass fraction by the mean shear.

3. Methodology

3.1. Numerical method for low-Mach HST simulations with phase change

The governing equations are solved on a uniform Cartesian grid of equal spacing $\Delta x=\Delta y=\Delta z$, with a staggered arrangement for the velocity while the remaining scalar fields are defined at the cell centres. The convection terms of the governing equations are discretized with the QUICK scheme (Leonard 1979), while central schemes are employed for the diffusive terms. The numerical method for phase change in the incompressible limit is presented in Scapin et al. (2020), while the details of the weakly compressible two-phase code are reported in Dalla Barba et al. (2021). Both implementations are based on the solver CaNS (Costa 2018), extended in this work to handle shear-periodic boundary conditions. In this section we therefore only describe the main modifications needed to handle weakly compressible phase-change processes in HST. The validation of the present algorithm against three benchmarks is provided in Appendix C.

3.1.1. Dispersed phase

The first step is the interface reconstruction and subsequent advection, handled in a fully Eulerian manner using an algebraic VoF method, i.e. the MTHINC by Ii et al. (2012), Rosti, De Vita & Brandt (2019a). We start from the topological equation for the phase indicator function (2.1),

(3.1)\begin{equation} \frac{\partial H}{\partial t} + \boldsymbol{u}_{\varGamma}\boldsymbol{\cdot}\boldsymbol{\nabla} H = 0, \end{equation}

where $\boldsymbol {u}_{\varGamma }$ is the interface velocity, taken as the sum of the extended liquid velocity $\boldsymbol {u}_l^{ext}$ and the contribution due to the interfacial mass flux $\dot {m}_{\varGamma }\boldsymbol {n}_{\varGamma }/\rho _l$; see Scapin et al. (2020) for more details. Note that since $\boldsymbol {u}_{l}^{ext}$ and $\boldsymbol {u}_{\varGamma }$ are derived from the one-fluid velocity $\boldsymbol {u}$, the decomposition (2.10) applies directly to the interface velocity.

Equation (3.1) is then rewritten in terms of the volume fraction, $\varPhi$, defined as the volumetric average of $H$ over a discrete computation cell of volume $V_c=\Delta x\Delta y\Delta z$. Employing the decomposition (2.10) in the colour function advection equation (3.1) yields

(3.2)\begin{equation} \frac{\partial\varPhi}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot}(H^{ht}\boldsymbol{u}'_{\varGamma}) + \mathcal{S}z\frac{\partial H^{ht}}{\partial x} = \varPhi\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}'_{\varGamma},\end{equation}

where $H^{ht}$ represents the hyperbolic tangent function approximating the indicator function $H$. From (3.2), we see that an additional term is present, $\mathcal {S}z(\partial H^{ht}/\partial x)$, which represents the convection of volume fraction by $\mathcal {S}$.

Equation (3.2) is solved in three sub-steps. First, it is advanced by $\Delta t^{n+1}$ to the new time step omitting the convection by the mean shear and the right-hand side term. By employing the classical directional splitting, we obtain a provisional $\tilde {\varPhi }$; see Ii et al. (2012) for more details. Next, the convection by mean shear is included,

(3.3)\begin{equation} \frac{\partial\tilde{\varPhi}}{\partial t}+\mathcal{S}z\frac{\partial \tilde{H}^{ht}}{\partial x} = 0. \end{equation}

Note that the time derivative in (3.3) is computed explicitly for $\tilde {\varPhi }$ while the convective term contains the phase indicator function $H^{ht}$. Since the latter is a function of the former and the mathematical form of this dependency varies according to the type of interface reconstruction method, it is possible to express the convective term as a function of $\varPhi$. Nevertheless, this would modify the advection velocity in (3.3) adding a spatial-dependent term along the mean shear direction (i.e. $x$), making more elaborated and complex the application of the method of characteristics. For this reason, we prefer to compute this extra term by an additional directional splitting along $x$, with $\mathcal {S}z$ as advection velocity.

Finally, the divergence is corrected. This consists in adding, after the four directional splittings, a correction term proportional to the discrete velocity divergence, i.e.

(3.4)\begin{equation} \varPhi_{i,j,k}^{n+1} = \tilde{\varPhi}_{i,j,k} - \Delta t^{n+1}F_{i,j,k}^n+\Delta t^{n+1}\varPhi_{i,j,k}^{n+1}(\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_{\varGamma})_{i,j,k}^n, \end{equation}

where $\varPhi _{i,j,k}^{n+1}$ is the volume fraction resulting from the directional splitting procedure, $F_{i,j,k}^n$ represents the correction used in the previous directional splitting steps, and the last term is the volume correction which ensures that the interface velocity divergence is employed to update $\varPhi _{i,j,k}^{n+1}$ (Scapin et al. 2020). The thermodynamic properties ($\rho, \mu, k$ and $c_p$) are then updated using the new value $\varPhi ^{n+1}$.

3.1.2. Vapour mass fraction equation

The conservation of the vapour mass fraction, see (2.12), is solved only in the gas domain (i.e. $V_g$), assuming saturation conditions at the interface $\varGamma$. In other words, the value $Y_{l}^{v}=Y_{l,\varGamma }^{v}$, which is a function of the thermodynamic pressure and temperature, is imposed at the interface $\varGamma$ as a Dirichlet boundary condition. To compute $Y_{l,\varGamma }^{v}$, we employ the Span–Wagner equation of state (see (A6) and (A7) in Appendix A). Equation (2.12) is advanced with the first-Euler method neglecting the mean shear contribution. This yields a provisional vapour mass fraction field $\tilde {Y}_{v,l}^{n+1}$,

(3.5)\begin{equation} \rho_g^{n+1}\frac{\tilde{Y}_{l}^{v,n+1}-Y_{l}^{v,n}}{\Delta t^{n+1}} ={-}\rho_g^{n+1}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} Y_{l}^v)^n + \frac{1}{Re_{\mathcal{S}}Sc}\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_g^{n+1}D_{lg}^{n+1}\boldsymbol{\nabla} Y_{l}^{v,n}). \end{equation}

The calculation of the gradient in the convective part of (3.5) is performed as detailed in Scapin et al. (2020), while some modifications are required for the linear term as it contains the diffusion coefficient $D_{lg}$ and the gas density $\rho _g$, which, in general, vary with the thermodynamic pressure, temperature and composition. Since the procedure can be performed dimension by dimension, we present the discretization along $x$ as example, the same approach being repeated for the other two directions,

(3.6)\begin{equation} \frac{\partial}{\partial x}\left(\rho_gD_{lg}\frac{\partial Y_{l}^{v}}{\partial x}\right)_{i} = \frac{1}{\Delta x}\left((\rho_gD_{lg})_{i+1/2}\left.\frac{\partial Y_{l}^{v}}{\partial x}\right|_{i+1/2}-\left.(\rho_gD_{lg})_{i-1/2}\frac{\partial Y_{l}^{v}}{\partial x}\right|_{i-1/2}\right)_i. \end{equation}

The evaluation of the gradients $\partial Y_{l}^{v}/\partial x_{i\pm 1/2}$ is performed on an irregular grid, employing one-sided finite difference equations for the cell cut from the interface or central scheme for cells away from the interface (see Scapin et al. (2020) for details). Next, the coefficients $(\rho _gD_{lg})_{i\pm 1/2}$ are obtained as the arithmetic mean,

(3.7)\begin{equation} (\rho_gD_{lg})_{i\pm 1/2} = 0.5[(\rho_gD_{lg})_i+(\rho_gD_{lg})^G]. \end{equation}

If the cell $i$ and its neighbours ($i\pm 1$) are not cut by the interface, $\rho _gD_{lg}^G$ is set equal to $(\rho _gD_{lg})_{i\pm 1}$. Otherwise, $(\rho _gD_{lg})^G$ is evaluated as

(3.8)\begin{equation} (\rho_gD_{lg})^G = \frac{(\rho_gD_{lg})_{\varGamma}+(\theta-1)(\rho_gD_{lg})_{i\pm 1}}{\theta}, \end{equation}

where $\theta$ represents the sub-cell resolution computed from the level-set function, reconstructed from the VoF field. Since (3.8) poorly behaves for small values of $\theta$, we set $(\rho _gD_{lg})^G=(\rho _gD_{lg})_i$ when $\theta \leq 1/4$. Note that (3.8) requires the value of the coefficient $(\rho _gD_{lg})$ at the interface location, i.e. $(\rho _gD_{lg})_{\varGamma }$. These are evaluated with the corresponding equations of state using the temperature and the vapour composition at the interface location. It is worth mentioning that we cannot access directly the liquid and gas temperatures, separately, since the temperature equation is solved over the whole domain, irrespective of the interface location. Therefore, to avoid problems of artificial heating (especially when the difference between the gas and the liquid density is high), we locally reconstruct the gas temperature in $V_l$ and the liquid temperature in $V_g$ relying on a simple constant extrapolation of $T$. The resulting two fields, $T_g$ and $T_l$, defined in a few grid cells around the interface, are then used to update the thermophysical properties; see Appendix A for more details.

Finally, the mean shear contribution is included using the method of characteristic over a time $\Delta t^{n+1}$,

(3.9)\begin{equation} \frac{\partial Y_{l}^{v}}{\partial t}+\mathcal{S}z\frac{\partial Y_{l}^{v}}{\partial x} = 0. \end{equation}

Equation (3.9) can be conveniently rewritten in more compact form as in Gerz, Schumann & Elghobashi (1989), Kasbaoui et al. (2017),

(3.10)\begin{equation} Y_{l}^{v,n+1} = \tilde{Y}_{l}^{v,n+1}(\boldsymbol{x}-\Delta t^{n+1}\mathcal{S}z\boldsymbol{e}_x). \end{equation}

Equation (3.10) is solved using a Fourier interpolation (Tanaka 2017), which ensures a spectral accuracy provided that $\Delta t^{n+1}$ is chosen lower or equal than $(\mathcal {S}N_z)^{-1}$, where $N_z$ is the number of grid points along the $z$ direction (see § 3.2 for more details).

3.1.3. Calculation of the interfacial mass flux

The interfacial mass flux (2.5) is computed only in the gas region by projecting the interfacial gradient along the normal direction and adopting a dimension by dimension approach. The interfacial vapour mass fraction $Y_{l}^{v}$, gas density $\rho _{g}$ and diffusion coefficient $D_{lg}$ should be estimated at the interface location. Similarly to the case of the vapour mass fraction described in § 3.1.2, all these quantities depend on the local gas and liquid temperatures and on the thermodynamic pressure. Therefore, we first estimate the interfacial liquid and gas temperatures in those grid cells cut by the interface. Depending on the interface position, this is done in each direction independently, providing different estimates for the gas and liquid temperatures for the three coordinates, $T_{p,\varGamma }^x, T_{p,\varGamma }^y$ and $T_{p,\varGamma }^z$, where the subscript $p$ stands for the gas and liquid phases. The resulting values are averaged using the local normal, i.e.

(3.11)\begin{equation} T_{p,\varGamma} = T_{p,\varGamma}^xn_x^2 + T_{p,\varGamma}^yn_y^2 + T_{p,\varGamma}^zn_z^2. \end{equation}

Once $T_{p,\varGamma }$ is known, the values of $\rho _{g,\varGamma }, D_{lg,\varGamma }$ and $Y_{l,\varGamma }^v$ are computed using the corresponding equations of state; see Appendix A. Note that by employing the procedure here explained, the mass flux $\dot {m}_{\varGamma }$ is available only on the grid nodes pertaining the gas region. Nevertheless, as (2.6) suggests, the values of $\dot {m}_{\varGamma }$ are needed also in some grid points inside the liquid region. Accordingly, $\dot {m}_{\varGamma }$ is extrapolated over a narrow band at the interface to populate all cells where $|\boldsymbol {\nabla } \varPhi |_{i,j,k}\neq 0$.

3.1.4. Temperature equation

The temperature equation (2.4) is solved using a whole domain approach as in Scapin et al. (2020), with additional care paid to the time discretization. The adopted approach follows that proposed in Gerz et al. (1989) and has been improved here to enhance the numerical stability and to include the additional source terms due to the gas compressibility, phase change and enthalpy diffusion. First, a prediction temperature field $\tilde {T}$ is computed using the Adams–Bashforth method,

(3.12)\begin{equation} \frac{\tilde{T}^{n+1}-T^n}{\Delta t^{n+1}} = \left(1+\frac{1}{2}\frac{\Delta t^{n+1}}{\Delta t^n}\right)RT^n-\left(\frac{1}{2}\frac{\Delta t^{n+1}}{\Delta t^n}\right)RT^{n-1}, \end{equation}

where the terms $RT^n$ and $RT^{n-1}$ include all the convective and diffusive terms at the current, $n$, and old time level, $n-1$. The term $RT^{n-1}$ is first shear mapped to the new time level $n+1$ and this step has proved to be crucial for the numerical stability of the Adams–Bashforth scheme (see Appendix D).

Next, the temperature field is shear mapped to the new time level $n+1$ by employing the same spectral interpolation described for the vapour mass fraction equation,

(3.13)\begin{equation} T^{n+1,*} = \tilde{T}(\boldsymbol{x}-\mathcal{S}z\Delta t^{n+1}\boldsymbol{e}_x). \end{equation}

Finally, the source terms in (2.4) are included using the first-order Euler scheme,

(3.14)\begin{align} \rho^{n+1}c_p^{n+1}\frac{T^{n+1}-T^{n+1,*}}{\Delta t^{n+1}} &= \varPi_{p,1}\left(\frac{{\rm d}p_{th}}{{\rm d}t}\right)^{ext}(1-\varPhi^{n+1}) \nonumber\\ &\quad +\frac{1}{Re_{\mathcal{S}}Sc}\left(\rho_gD_{lg}\sum_{j=1}^2\boldsymbol{\nabla} Y_j\boldsymbol{\cdot}\boldsymbol{\nabla} h_j^{*}\right)^{n+1}(1-\varPhi^{n+1}) \nonumber\\ &\quad - \frac{(\dot{m}_{\varGamma}|\boldsymbol{\nabla} \varPhi|)^{n+1}}{Ste}, \end{align}

where $(\textrm {d}p_{th}/\textrm {d}t)^{ext}$ represents the linear extrapolation in time of the derivative of the thermodynamic pressure. It is important to note that the second and third terms in (3.14), which include variables shear mapped to the new time level (i.e. $\varPhi ^{n+1}$ and $Y_{l,j}^{v,n+1}$), should be computed after (3.13).

3.1.5. Thermal divergence and thermodynamic pressure

The calculation of the velocity divergence is performed simply by discretizating the right-hand side of (2.6) node by node as done in Dalla Barba et al. (2021) for the case of weakly compressible two-phase solvers without phase change. The calculation of the thermodynamic pressure requires, however, more care. In the low-Mach framework, the role of $p_{th}$ is to ensure mass conservation of the compressible phase at the discrete level, since it enters the calculation of the gas density. This cannot be fulfilled simply by the advection of the colour function, which is designed to satisfy the volume conservation of the incompressible liquid, or by the pressure-correction step through the imposition of the divergence constrain (2.6) on $\boldsymbol {u}$. In fact, these ensure only the overall volume conservation of the closed and isochoric system under consideration. Therefore, we compute $p_{th}$ by integrating the equation for the gas density ((A1) in Appendix A) over the computation domain occupied by the gas phase $V_g$ (Demou, Frantzis & Grigoriadis 2019; Dalla Barba et al. 2021),

(3.15)\begin{equation} p_{th}^{n+1} = \frac{G_{g}^{n+1}\varPi_{p,2}}{\displaystyle{\int_{V} \frac{\bar{M}_{m,av}^{n+1}}{T^{n+1}}(1-\varPhi^{n+1})\,{{\rm d}\kern0.7pt x}\,{{\rm d} y}\,{\rm d} z}}. \end{equation}

The total mass of the gas $G_g$, used above, varies in time according to the following relation, derived from the material balance across the droplet surface,

(3.16)\begin{equation} \frac{{\rm d}G_g}{{\rm d}t} = \dot{m}_{t,\varGamma} = \int_{V}\dot{m}_{\varGamma}^{n+1}|\boldsymbol{\nabla} \varPhi^{n+1}|\,{{\rm d}\kern0.7pt x}\,{{\rm d} y}\,{\rm d} z. \end{equation}

Once the volume integral in (3.16) is computed, the gas mass at the new time level is computed from the time integration of (3.16),

(3.17)\begin{equation} G_{g}^{n+1} = G_{g}^{n} + \int_{t^n}^{t^{n+1}} \dot{m}_{t,\varGamma}\,{\rm d} t.\end{equation}

To evaluate numerically the integral in (3.17), we employ a trapezoidal quadrature. With this approach, we impose correctly the conservation of the compressible phase at the discrete level.

3.1.6. Momentum equation and pressure correction method

Once the thermodynamic divergence is computed, the momentum equation (2.2) is solved with a standard pressure correction method, reported below in semi-discrete form,

(3.18)$$\begin{gather} \boldsymbol{RU}^{n-1} = \boldsymbol{\widetilde{RU}}^{n-1}(\boldsymbol{x} -\Delta t^{n}\mathcal{S}z\boldsymbol{e}_x), \end{gather}$$

(3.19)$$\begin{gather}\rho^{n+1}\left(\frac{\boldsymbol{u}^{{\star}{\star}}- \boldsymbol{u}^{n}}{\Delta t^{n+1}}\right) = \left(1+\frac{1}{2}\frac{\Delta t^{n+1}}{\Delta t^n}\right) \boldsymbol{RU}^{n}-\left(\frac{1}{2}\frac{\Delta t^{n+1}}{\Delta t^n}\right) \boldsymbol{RU}^{n-1}, \end{gather}$$

(3.20)$$\begin{gather}\boldsymbol{u}^{{\star}} = \boldsymbol{u}^{{\star}{\star}}(\boldsymbol{x}-\Delta t^{n+1}\mathcal{S}z\boldsymbol{e}_x) + \frac{\Delta t^{n+1}}{\rho^{n+1}} \frac{\kappa_{\varGamma}^{n+1}\boldsymbol{\nabla} \varPhi^{n+1}}{We_{\mathcal{S}}}, \end{gather}$$

(3.21)$$\begin{gather}\nabla^2p^{n+1} = \boldsymbol{\nabla} \boldsymbol{\cdot}\left[\left(1-\frac{\rho_0^{n+1}}{\rho^{n+1}}\right) \boldsymbol{\nabla} \hat{p}\right]+\frac{\rho_0^{n+1}}{\Delta t^{n+1}} \left(\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}^\star{-}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}^{n+1}\right), \end{gather}$$

(3.22)$$\begin{gather}\boldsymbol{u}^{n+1} = \boldsymbol{u}^\star{-}\Delta t^{n+1} \left[\frac{1}{\rho_0^{n+1}}\boldsymbol{\nabla} p^{n+1}+\left(\frac{1}{\rho^{n+1}}- \frac{1}{\rho_0^{n+1}}\right)\boldsymbol{\nabla} \hat{p}\right], \end{gather}$$

where $\boldsymbol {RU}^n$ and $\boldsymbol {RU}^{n-1}$ in (3.19) include the convective and diffusive terms computed at the current and previous time level, neglecting the surface tension force and the mean shear contribution. As done for the temperature equation and explained in detail in the Appendix D, it is very important for the stability and the accuracy of the time-integration scheme that the term $\boldsymbol {RU}^{n-1}$ is shear mapped to the current time level $n$. The intermediate velocity $\boldsymbol {u}^{\star \star }$ is shear mapped to the new time level $n+1$, similarly to what is done for the temperature and vapour mass fraction, and updated with the contribution from the surface tension (Rosti et al. 2019b).

Finally, the pressure equation (3.21) is solved with a time-splitting technique which allows us to transform the variable Poisson equation into a constant coefficient one (Dodd & Ferrante 2014). Note that $\rho _0^{n+1}$ and $\hat {p}$ in (3.21) and (3.22) are the minimum density over the entire computational domain and the extrapolated hydrodynamic pressure $\hat {p}=(1+(\Delta t^{n+1}/\Delta t^{n}))p^n-(\Delta t^{n+1}/\Delta t^{n})p^{n-1}$. As already discussed in Motheau & Abraham (2016), Dalla Barba et al. (2021), the use of the pressure splitting to solve the variable density Poisson equation in a low-Mach framework is suitable when the temperature ratio between the two phases is below 2–3, which is the case of the current study. Finally, it is worth mentioning that the presence of the shear, if left untreated, would make the use of the eigenexpansion method to solve (3.21) not possible, since one direction is not periodic. For these reasons, in order to still benefit from the FFT-based solvers, (3.21) is solved in a coordinate system moving with the mean shear for which, triperiodic boundary conditions can be applied. The solution is then transformed back to the shear-periodic domain, as detailed in Tanaka (2017), Rosti et al. (2019b).

3.2. Time-step restriction

The time step $\Delta t^{n+1}$ is estimated from the stability constraints of the overall system,

(3.23)\begin{equation} \Delta t^{n+1}=\min\left[C_{\Delta t_c}\Delta t_c^{n+1},C_{\Delta t_d}\min(\Delta t_{\mu},\Delta t_{\sigma},\Delta t_m,\Delta t_k)^{n+1},C_{\Delta t_\mathcal{S}}\frac{\Delta z}{\mathcal{S}l_z}\right], \end{equation}

where $\Delta t_c, \Delta t_{\sigma }, \Delta t_{\mu }, \Delta t_m$ and $\Delta t_k$ are the maximum allowable time steps due to convection, surface tension, momentum, thermal energy and mass diffusion. These can be determined as suggested in Scapin et al. (2020), Dalla Barba et al. (2021),

(3.24) \begin{gather} \left.\begin{gathered} \Delta t_c = \dfrac{\Delta x}{|u_{max}|+|v_{max}|+|w_{max}|}, \quad \Delta t_{\mu}=\min\left(\min_{V_g}\left\{\dfrac{\mu_g}{\rho_g}\right\}, \dfrac{\lambda_{\mu}}{\lambda_{\rho}}\right)\dfrac{\Delta x^2Re_{\mathcal{S}}}{6}, \\ \Delta t_{\sigma}=\sqrt{\dfrac{We_{\mathcal{S}}(\min_{V_g}\{\rho_g\} +\lambda_{\rho})\Delta x^3}{4{\rm \pi}}}, \quad \Delta t_{m}=\min_{V_g} \{D_{lg}\}\dfrac{(\theta_m\Delta x)^2Sc}{6}, \\ \Delta t_{k}=\min\left(\min_{V_g}\left\{\dfrac{k_g}{\rho_gc_{pg}}\right\}, \dfrac{\lambda_{k}}{\lambda_{\rho}\lambda_{c_p}}\right)\dfrac{\Delta x^2Pr}{6}, \end{gathered}\right\} \end{gather}

where $|u_{i,\max }|$ is an estimate of the maximum value of the $i$th component of the flow velocity, $\theta _m=0.25$ and $\max _{V_g}\{\xi _g\}$ and $\min _{V_g}\{\xi _g\}$ denote the maximum and minimum over the computational domain, $V_g$, of the generic thermophysical property of the gas phase. For the cases presented here, the convective constrain represents the main limitation; setting $C_{\Delta t_c}=0.15, C_{\Delta t_d}=0.5$ and $C_{\Delta t_{\mathcal {S}}}=1$ was found sufficient for a stable and accurate time integration.

3.3. Computational set-up

Given the large number of dimensionless parameters that characterizes flows involving evaporation, we focus our attention on the role of the ratio between droplet initial diameter and the Kolmogorov dissipative flow scale (tuned by varying the shear-based Reynolds number, $Re_{\mathcal {S}}$), the role of surface tension (varying the shear-based Weber number $We_{\mathcal {S}}$), the ratio between the initial gas temperature and the critical temperature, $T_{g,0}/T_c$, and the type of model employed to evaluate the thermophysical properties of the gas phase during the simulations. Concerning this last aspect, we first consider all the gas thermophysical properties constant and computed with a proper averaging between the liquid and the gas temperature (i.e. with the $1/3$ rule by Hubbard, Denny & Mills 1975) (case denoted as CP); secondly, we allow only the gas density $\rho _g$ to vary (case denoted as VP$_{\rho }$) and, finally, we allow all the thermophysical properties to vary with the local thermodynamic variables and vapour composition (case VP$_{a}$). A summary of the numerical campaign is reported in table 1, together with the remaining dimensionless physical parameters, which are all kept constant. Note that the physical parameters reported in table 1 are representative of pentane evaporating droplets in dry air at high pressure ($\sim 43\ \mathrm {bar}$). In particular, we will consider three values of the ratio $T_{g,0}/T_c=0.75, 1.00, 1.50$, where $T_c$ is the critical temperature (469.69 K for pentane), which gives $T_{g,0}=354, 470$ and $705$ K. The initial liquid temperature $T_{l,0}$ is the wet bulb temperature at $T_{g,0}$ and corresponds to $T_{l,0}=334, 388$ and $432$ K.

Table 1. Left: Dimensionless parameters defining the investigated cases, the initial gas temperature over the critical temperature $T_{g,0}/T_c$, the thermodynamic model employed to evaluate the gas thermophysical property, the shear-based Reynolds number $Re_{\mathcal {S}}=\rho _{g,r}\mathcal {S}l_y^2/\mu _{g,r}$ and the shear-based Weber number $We_{\mathcal {S}}=\rho _{g,r}\mathcal {S}^2d_0^3/\sigma$ with $d_0$ the initial droplet diameter. The vaporization Damköhler number, the ratio between the turbulence time scale and the evaporation time scale in stagnant conditions, $Da_v=\tau _t/\tau _{v,L}$, is also reported (see (4.4)). Note that in the current study $d_0/l_y=0.10$. Right: Dimensionless parameters kept constant in the current study ($N_{dp,0}$ is the initial number of droplets, $\alpha _0$ is the initial liquid volume and $RH_0$ is the initial relative humidity).

The focus of the study is the behaviour of evaporating isolated droplets in HST and, therefore, we consider five isolated droplets (i.e. initial volume fraction $\alpha _0\approx 0.14\,\%$ with $N_{dp,0}=5$), injected in a single-phase statistically steady-state HST field at the desired shear-based $Re_{\mathcal {S}}$ ($2800$ or $6700$). For all the cases, the computational domain has the streamwise aspect ratio, $\mathcal {AR}_{xy}=l_x/l_y\approx 2.10$, and the cross-stream ratio, $\mathcal {AR}_{zy}=l_z/l_y\approx 1.05$. As discussed in Sekimoto et al. (2016), employing such values, the effects on the flow induced by a finite-size computational box are reduced. The domain is discretized with $1280\times 608\times 640$ grid points, thus ensuring not only an adequate resolution of the flow field (i.e. $\Delta x/\eta \approx 0.33$ for $Re_{\mathcal {S}}=2800$ and $\Delta x/\eta \approx 0.40$ for $Re_{\mathcal {S}}=6700$, where $\eta$ is the Kolmogorov length scale), but also of the droplets, whose initial resolution is $64$ points per diameter. This value is consistent with our previous study (Scapin et al. 2020), where we have assessed that a minimum of $50$ grid points per diameter is needed to fully resolve mass, momentum and energy transfer across the droplet interface. Note that the resolution is also sufficient to resolve the smallest scales of the two active scalar fields, $Y_{l}^{v}$ and $T$, since the associated Batchelor scales $\eta _{B,Y}=\eta /\sqrt {Sc}$ and $\eta _{B,T}=\eta /\sqrt {Pr}$, although smaller than the Kolmogorov scale, are always larger than $\Delta x$.

4. Results

4.1. Comparison of the three thermodynamic models

We start our analysis by assessing the influence on the evaporation dynamics of the type of thermodynamic model employed to evaluate the gas thermophysical properties. A common approach (Abramzon & Sirignano 1989; Lupo et al. 2019, 2020) is to consider the thermophysical properties of the gas phase uniform, constant and evaluate them at an intermediate temperature, $T_m=T_{g,0}+m(T_{g,0}-T_{l,0})$. Hubbard et al. (1975) show that the value $m=1/3$ guarantees a good agreement between experimental results and the theoretical predictions based on the $d^2$–law (Langmuir 1918). In the mathematical framework introduced above, this amounts to limiting the expansion and contraction at the interface, i.e. the terms $f_Y$ and $f_T$ in (2.6), (2.7) become zero. Here, we will denote the results obtained with these assumptions as CP. In many fully resolved simulations, the only thermophysical property allowed to vary is the density and this is typically done within the Oberbeck–Boussinesq approximation (Piedra et al. 2015). This model, denoted VP$_{\rho }$, is considered here to asses whether it can provide accurate results. In typical conditions of evaporating droplets, however, other thermophysical properties may play a role. In particular, $\rho _g$ and $D_{lg}$ scale differently with the gas temperature, i.e. $\rho _g\sim p_{th}/T_g$ and $D_{lg}\sim T_g^{1.5}/p_{th}$ and both appear in the expression of the interfacial mass flux. These variations may significantly affect the evaporation dynamics, especially when the difference between the liquid and gas temperatures is large. We will therefore also consider variations of all the thermophysical properties, case VP$_{a}$. To compare the different models, we consider two temperature ratios, $T_{g,0}/T_c=0.75$ and $T_{g,0}/T_c=1.50$ at $Re_{\mathcal {S}}=6700$ and the lowest Weber number under consideration in this study, $We_{\mathcal {S}}=0.02$, to limit droplet breakup and reduce the droplet deformation. Unless otherwise stated, the results refer to the averaged values over the five droplets initially in the computational domain and error bars are included to represent the droplet with the largest positive and negative deviation from the mean value. Firstly, we show the square of the normalized droplet diameter; see figure 2. From the cases at $T_{g,0}/T_c=1.50$, we observe that the complete model and the constant property model provide similar evaporation rates; when the gas density is the only varying thermophysical property, the evaporation rate is the highest. This behaviour can be attributed to the presence of colder gas around the droplets, leading to larger local gas densities (up to three times the initial gas density, see figure 3) and, thus, to increased evaporation rates, as shown in figure 4(a) where we report the time history of the surface-averaged gas temperature. Relaxing the assumption of constant liquid diffusion coefficient reduces $D_{lg} \propto T_g^{1.5}$. This counteracts the effect of the higher gas density, decreasing the overall evaporation rate, which approaches the values obtained assuming constant property values. The results at $T_{g,0}/T_c=0.75$ show a limited impact of compressibility on the evaporation dynamics, with an almost identical evaporation from the CP and VP$_{a}$ models. Once more, allowing only the gas density to vary leads to an overestimation of the evaporation rate, which is explained by a lower gas temperature and higher density at the interface; see figure 4(b). Note that in all the cases and regardless of the model, the mean evaporation rate represents a good estimation of the evaporation rate of the single droplet since the magnitude of the largest positive and the negative deviations (represented by the error bars in figure 2) is within $5\,\%$.

Figure 2. Normalized square droplet diameter, $(d/d_0)^2$, as a function of the dimensionless time (based on the diffusion time scale) for temperature ratio (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$ for $Re_{\mathcal {S}}=6700$ and $We_{\mathcal {S}}=0.02$. The dotted curves indicate the slope of the curves after the initial transient (extracted from direct numerical simulations), whereas the dashed curves indicate the estimation from the theoretical Frössling/Ranz-Marshall correlation. The length of the error bars indicates the droplet with the fastest/slowest evaporation rate among the five droplets in the simulations.

Figure 3. Volume rendering of the normalized gas density $\rho _{g}/\rho _{g,0}$ around a droplet for the case $Re_{\mathcal {S}}=6700, We_{\mathcal {S}}=0.02$ and $T_{g,0}/T_c=1.5$.

Figure 4. Normalized gas temperature at the interface as a function of the dimensionless time (based on the diffusion time scale) for temperature ratio (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$ for $Re_{\mathcal {S}}=6700$ and $We_{\mathcal {S}}=0.02$. Here $T_{g,0}$ and $T_{l,0}$ are the initial gas and liquid temperatures.

Figure 5 displays the time evolution of the ratio between the instantaneous area $A=|\boldsymbol {\nabla } \varPhi |\Delta x^3$ and the area of a spherical droplet with the same volume with error bars indicating the maximum and minimum values at each instance at $T_{g,0}/T_c=1.5$ and $T_{g,0}/T_c=0.75$. Note that the mean deviation from the spherical shape is always below $10\,\%$ (over the simulated physical time) which makes the present data suitable for comparisons with existing scaling laws, which are commonly derived under the assumption of rigid and spherical droplets. In all cases, the evaporation rate follows the linear trend predicted by the $d^2$–law since the droplets are isolated. We therefore test the closure relations available in the literature against the data from the interface-resolved simulations. To this end, we estimate the evaporation rate as proposed in Birouk & Gökalp (2006). First, we obtain the dimensionless evaporation rate in laminar conditions $K_{L}$ (Abramzon & Sirignano 1989) as

(4.1)\begin{equation} K_{L} = \frac{4 Sh_{RM}}{Re_{\mathcal{S}}Sc}\frac{\log(1+B_M)}{\lambda_{\rho}},\end{equation}

where the Spalding number $B_M=(\langle Y_{l,\varGamma }^v\rangle _{\varGamma }-Y_{l,\infty }^v)/(1-\langle Y_{l,\varGamma }^v\rangle _{\varGamma })$, with $\langle \rangle _{\varGamma }$ the average over the surface $\varGamma$. The Sherwood number in (4.1) is estimated with the Frössling/Ranz-Marshall correlation (Ranz & Marshall 1952; Birouk & Gökalp 2006),

(4.2)\begin{equation} Sh_{RM} = 2 + \frac{B_M}{(1+B_M)^{0.7}\log(1+B_M)}Sh_{0,RM},\end{equation}

with $Sh_{0,RM}=2+0.552Re_d^{0.50}Sc^{0.33}$. In this last expression, $Sc$ is the Schmidt number and $Re_d$ is the droplet Reynolds number, defined as $Re_d= \rho _{g,r}|\boldsymbol {u}_{l,d}'-\boldsymbol {u}_{g,d}'|d/\mu _{g,r}$. Following the original reference (Ranz & Marshall 1952) and more recent works (Ng et al. 2021; Wang et al. 2021), we compute $Re_d$ using the instantaneous droplet diameter, the instantaneous droplet velocity $\boldsymbol {u}_{l,d}'$ and the surrounding gas velocity $\boldsymbol {u}_{g,d}'$ at the droplet location (Ng et al. 2021). Note that this correlation is based on experiments with droplets larger than the Kolmogorov scale and in presence of a mean velocity field and is therefore suitable for the current study. Finally, we correct (4.1) for the presence of turbulence via a Damköhler number (Gökalp et al. 1992),

(4.3)\begin{equation} \frac{K}{K_L} = 0.74Da_v^{{-}0.115}.\end{equation}

In the above,

(4.4)\begin{equation} Da_v = \frac{\tau_t}{\tau_{v,L}},\end{equation}

where $\tau _t=d_0^{2/3}/\langle \varepsilon _g\rangle _{M_g}^{1/3}$ is the characteristic time scale of the turbulent eddies, $d_0$ is the initial droplet diameter and $\langle \varepsilon _g\rangle _{M_g}$ is the mass-averaged turbulent dissipation (see Abramzon & Sirignano 1989; Birouk & Gökalp 2006). The characteristic time scale of evaporation $\tau _{v,L}$, is computed as

(4.5)\begin{equation} \tau_{v,L} = \frac{\delta_M}{V_r},\end{equation}

where $\delta _M=d_0/(Sh_{RM}-2)$ (Abramzon & Sirignano 1989) is the vapour film boundary layer thickness and $V_r$ is the Stefan flow velocity computed as $V_r=8{\rm \pi} D_{lg,r}\log (1+B_M)/d_0$. Note that while (4.1), (4.2) and the expressions for $\tau _t$ and $\delta _M$ are evaluated using also the instantaneous data from the interface-resolved simulations, the Stefan velocity $V_r$ is evaluated a priori as it depends only on the initial thermodynamic conditions.

Figure 5. Ratio between the instantaneous interfacial area $A$ and the area of a spherical droplet with the same volume $V$, i.e. $A_{eq}={\rm \pi} ^{1/3}(6V)^{2/3}$, for $We_{\mathcal {S}}=0.02$ and (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$. The upper lengths of the error bar refer to the droplet with the largest $A/A_{eq}$, while the lower lengths displacement to the droplet with lowest $A/A_{eq}$. The data refers to the simulations conducted with the VP$_a$ thermodynamic models (negligible difference has been observed between the CP and VP$_a$ models).

The estimated values of $K$, dashed lines in figure 2, display a good agreement with the interface-resolved simulations at $T_{g,0}/T_c=1.50$, especially for the CP and VP$_{a}$ models, while a significant deviation is observed at $T_{g,0}/T_c=0.75$, regardless of the thermodynamic model employed. Finally, it is worth noting that the evaporation rate estimated with (4.1) depends linearly on $Sh$, whose calculation is generally affected by a higher degree of uncertainty (Lupo & Duwig 2020).

Next, we examine the temporal evolution of the Sherwood number, $Sh$, for the three different models, see figure 6, using directly the definition

(4.6)\begin{equation} Sh = \frac{\dot{m}_{t,\varGamma}d_{eq}}{\rho_{g,r}D_{lg,r}A(\langle Y_{l,\varGamma}^v\rangle_{\varGamma}-Y_{l,\infty}^v)}. \end{equation}

Here $\dot {m}_{t,\varGamma }$ is the total interfacial mass flux, $A$ the interfacial area, $d_{eq}$ the instantaneous equivalent diameter (i.e. $d_{eq}=(6V/{\rm \pi} )^3$), $\rho _{g,r}$ and $D_{lg,r}$ the reference gas density and gas–vapour diffusion coefficients, $\langle Y_{l,\varGamma }^v\rangle _{\varGamma }$ the surface-averaged vapour mass fraction at the interface and $Y_{l,\infty }^v$ the reference vapour mass fraction, i.e. its initial value. The Sherwood number defines the ratio between the mass transfer in the actual flow and the mass transfer by diffusion, i.e. evaporation in stagnant conditions, so that it can be used to quantify the effect of the flow on the evaporation. Note that since we perform a comparison among different thermodynamic models, the gas thermophysical properties in (4.6) are taken equal to the reference ones (i.e. those obtained with the ‘1/3’ rule for the CP model and those at the initial condition for the VP$_{a}$ and VP$_{\rho }$ models). Both at high and low temperatures, $Sh$ approaches an asymptotic value after the initial transient. Note that for the cases at $T_{g,0}/T_c=1.5$, the transient is faster ($\Delta tD_{lg,r}/d_0^2\approx 0.02$) than for the cases at $T_{g,0}/T_c=0.75 (\Delta tD_{lg,r}/d_0^2\approx 0.1)$ since large temperature differences increase the evaporation rates. For the high-temperature cases, the $Sh$ number differs significantly with the model employed, with the model VP$_{\rho }$ giving the largest value of $Sh$ for the reasons explained above. Conversely, at lower temperature, the three models provide a similar $Sh$, confirming the limited impact on the evaporation rate of the choice of the thermodynamic models in this regime. All other parameters being fixed, the Sherwood number is higher at the lower temperature ratio, $T_{g,0}/T_c=0.75$ than for $T_{g,0}/T_c=1.50$, regardless of the thermodynamic model used. Recalling that the Sherwood number is the ratio between the actual mass transfer and the mass transfer by diffusion, the increase of evaporation rate due to the background turbulence is thus larger at low temperatures. For the highest temperature ratio, the enhancement induced by the turbulence is lower since evaporation is mainly driven by the large temperature difference and by the larger values of $(\langle Y_{l,\varGamma }^v\rangle _{\varGamma }-Y_{l,\infty }^v)$.

Figure 6. Temporal evolution of the Sherwood number as a function of the dimensionless time (based on the diffusion time scale) for temperature ratio (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$ and $Re_{\mathcal {S}}=6700$ and $We_{\mathcal {S}}=0.02$. The dashed curves are the theoretical prediction provided from the a priori analysis based on the Frössling/Ranz-Marshall correlation. For the cases at $T_{g,0}/T_c=1.50$, only one value of $Sh$ is reported as the three models provide almost identical predictions.

As a final remark, note that the Sherwood number from direct numerical simulation (DNS) agrees well with the estimation from the Frössling/Ranz-Marshall correlation at high temperature (cf. models CP and VP$_{a}$ against the dashed line in figure 6). Conversely, at lower temperature, the estimated value is much lower than the one extracted from DNS, indicating that the Frössling/Ranz-Marshall correlation underpredicts convective effects due to turbulence at low evaporation rates, i.e. when the evaporation time becomes longer than the turbulence time scale ($Da_v \lesssim 0.2$ instead of $Da_v \approx 1$ as shown in table 1). This finding is consistent with the recent investigations in HIT performed by Méès et al. (2020) and suggests the need of including the turbulent effects (via the vaporization Damköhler number) in the available correlations for the Sherwood number and more in general in the evaporation models (Lupo et al. 2020).

4.2. Effects of the shear-based Reynolds number

We now consider the effect of the shear-based Reynolds number on the evaporation rate. Reducing from $Re_{\mathcal {S}}=6700$ to $Re_{\mathcal {S}}=2800$ results in a reduction of the ratio between the initial droplet diameter and the Kolmogorov length scale from $d_0/\eta =21.5$ to $d_0/\eta =10.5$. For this analysis, we assume that $We_{\mathcal {S}}=0.02$, three temperature ratios, $T_{g,0}/T_c=0.75-1.00-1.50$ and employ the complete model (i.e. VP$_a$). We first examine the reduction of the droplet diameter, $(d/d_0)^2$, as a function of the diffusion time scale $tD_{lg,r}/d_0^2$; see figure 7(a). In all cases, droplets with larger $d_0/\eta$ (i.e. higher $Re_{\mathcal {S}}$) evaporate faster, in agreement with previous experimental studies in HIT (Verwey & Birouk 2018, 2020) and in the presence of a mean flow (Marti et al. 2017). This is explained by the enhanced surface vapour gradient and faster dispersion of the vapour around the droplet for larger $d_0/\eta$, related to the fact that larger flow structures are more energetic. Note that a relative increase of the evaporation rate, quantified as the ratio between the evaporation rate $K$ extracted from DNS and the one computed in stagnant conditions, $K_0$, occurs at all temperatures as reported in figure 7(b) and it scales with the dimensionless gas temperature with similar exponents for the two different $d_0/\eta$ investigated, i.e. $(T_{g,0}/T_c)^{-0.76}$ for $d_0/\eta =21.5$ and $(T_{g,0}/T_c)^{-0.72}$ for $d_0/\eta =10.5$. The experimental correlations in Birouk & Gökalp (2006) suggests that, at high temperature and pressure i.e. when the evaporation is much faster than the turbulent time scale, the evaporation is only weakly dependent on the droplet size $d_0/\eta$. This is consistent with our results, where the difference in $K$ for the two droplet sizes examined decreases when increasing the temperature; however, we do not observe the rate of evaporation to become independent of the size in the parameter range considered.

Figure 7. (a) Normalized square droplet diameter $(d/d_0)^2$ as a function of the diffusion time scale $tD_{lg}/d_0^2$. (b) Ratio between the turbulent evaporation rate $K$ and the one at stagnant conditions $K_0=(8/\lambda _{\rho })\log (1+B_M)$, with $B_M=(Y_{l,\varGamma }^v-Y_{l,\infty }^v)/(1-Y_{l,\varGamma }^v)$. The data refer to the cases at $Re_{\mathcal {S}}=2800-6700$ and $T_{g,0}/T_c=1.5, 1.00$ and $0.75$. The length of the error bars indicates the droplet with the fastest/slowest evaporation rate among the five droplets in the simulations.

The variation of the evaporation rate with $d_0/\eta$ significantly affects the liquid temperature at the interface $T_{\varGamma,l}$ (computed as the surface average over the gas–liquid interface, $T_{\varGamma,l}=(1/\varGamma )\int _{\varGamma } Td\varGamma$) and the mean liquid temperature $T_{V,l}$ (computed as the volume average over the liquid, $T_{V,l}=(1/V_l)\int _{V_l} TdV_l$). This is shown in figure 8 for the three temperatures under investigation. At the highest temperature, figure 8(a), both $T_{\varGamma,l}$ and $T_{V,l}$ decrease with respect to the initial values due to the strong cooling caused by the evaporation. However, due to the heat transfer enhancement at higher $Re_{\mathcal {S}}, T_{\varGamma,l}$ remains greater than $T_{V,l}$ for droplets at larger $d_0/\eta$, whereas $T_{\varGamma,l}$ is lower than $T_{V,l}$ for the smaller $d_0/\eta$ considered. Clearly, the two temperatures eventually converge to a similar value once the thermal gradients inside the liquid reduce. For the case with intermediate temperature, see figure 8(b), the cooling due to evaporation is not sufficiently high to counteract the heat released from the gas and the two temperatures reach a regime value larger than the initial one. The final values of $T_{\varGamma,l}$ and $T_{V,l}$ are higher for the droplets with larger $d_0/\eta$ (larger Reynolds number). We also note that the interface region is warmer than the droplet bulk (i.e. $T_{\varGamma,l}>T_{V,l}$) for both droplet sizes. At the lowest temperature, figure 8(c), both $T_{\varGamma,l}$ and $T_{V,l}$ increase with time, when the evaporation is too slow to cool the droplets. Unlike the previous two cases, the values of $T_{\varGamma,l}$ and $T_{V,l}$ at $d_0/\eta =21.5$ are always lower than the corresponding ones at $d_0/\eta =10.5$, due to the higher evaporation rate at the highest $Re_{\mathcal {S}}$. For this case, the faster evaporation also changes the transient of $T_{\varGamma,l}$, which at the beginning is higher than $T_{V,l}$ while for $tD_{lg,r}/d_0>0.2$ becomes lower.

Figure 8. Dimensionless liquid temperature at the interface $(T_{l,\varGamma }-T_{l,0})/(T_{g,0}-T_{l,0})$ (dotted line) and average liquid temperature $(T_{l,V}-T_{l,0})/(T_{g,0}-T_{l,0})$ (continuous line) at (a) $T_{g,0}/T_c=1.5$ , (b) $T_{g,0}/T_c=1.0$ and (c) $T_{g,0}/T_c=0.75$.

We conclude this section by analysing the relative importance of the convective and conductive heat and mass vapour fluxes, $\mathcal {F}_{T,i=c,d}$ and $\mathcal {F}_{Y,i=c,d}$, in the gas region around the droplets and inside the liquid region. These are computed as an integral over the control surfaces $S_T$ and $S_{Y}$,

(4.7) \begin{gather} \left.\begin{gathered} \displaystyle\mathcal{F}_{T,d} ={-}\dfrac{1}{S_T}\int_{S_T} k\boldsymbol{\nabla} T\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}\, {\rm d} S_T, \\ \displaystyle\mathcal{F}_{T,c} = \dfrac{1}{S_T}\int_{S_T} \rho c_p(\boldsymbol{u}' T+\mathcal{S}z\boldsymbol{e}_xT)\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}\,{\rm d} S_T, \\ \displaystyle\mathcal{F}_{Y,d} ={-}\dfrac{1}{S_Y}\int_{S_Y} \rho_gD_{lg}\boldsymbol{\nabla} Y_{l}^{v}\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}\, {\rm d} S_Y, \\ \displaystyle\mathcal{F}_{Y,c} = \dfrac{1}{S_Y}\int_{S_Y} \rho_g(\boldsymbol{u}' Y_{l}^{v}+\mathcal{S}z\boldsymbol{e}_xY_l^v)\boldsymbol{\cdot}\boldsymbol{n}_{\varGamma}\, {\rm d} S_Y. \end{gathered}\right\} \end{gather}

In the gas phase, $S_T$ and $S_Y$ are surfaces conforming to the droplet shape (quasi-spherical at this low value of $We_{\mathcal {S}}$) at a distance approximately equal to the thermal and vapour mass fraction boundary layer thicknesses $\delta _T$ and $\delta _Y$. These are determined as the distance from the interface where $(T(\delta _T)-T_{\varGamma })/(T_{g,0}-T_{\varGamma })=0.99$ and $(Y_l^v(\delta _Y)-Y_{l,\varGamma }^v)/(Y_{l,0}^v-Y_{l,\varGamma }^v)=0.99$ as in Ni et al. (2021), moving normal to the interface thanks to a level-set function, reconstructed from the VoF field. We evaluate $\mathcal {F}$ as in (4.7) after the evaporation reaches a steady condition (i.e. $d^2$ regime) and average in time over a time interval $\Delta t=[0.20,0.50,0.60](tD_{lg,r})/d_0^2$ for the high, intermediate and low temperature ratios. To capture the dominant transport mechanism in the liquid phase, $S_T$ is defined at a distance from the interface equal to half the instantaneous droplet radius. The relative importance of the conductive and convective heat fluxes in the gas region are shown in figure 9: increasing the ratio $T_{g,0}/T_c$, the main transport mechanism changes from convection to conduction, as already anticipated in § 4.1 when discussing the Sherwood number $Sh$. The same applies to lower $Re_{\mathcal {S}}$, where we observe a slightly higher conductive contribution for $T_{g,0}/T_c=1.5$. In the liquid region, figure 10, the heat transport is mainly driven by convection in all cases. In agreement with previous experimental results (Wong & Lin 1992; Pinheiro et al. 2019), the dominant transport mechanism is weakly affected by a change of the droplet Reynolds number, here varied by changing $Re_{\mathcal {S}}$, since the motion of the liquid inside the droplet is not significantly affected by the intensity of the gas turbulence. Finally, the vapour fluxes reported in figure 11 show that $Y_l^v$ is predominantly transported by convection for the two Reynolds numbers investigated, mass diffusion being relevant only at the highest temperatures.

Figure 9. Conductive and convective heat fluxes $\mathcal {F}_T$ (over the sum of the two, $\varSigma _{\mathcal {F}}$), $\mathcal {F}_{T,i=c,d}$, in the gas region for $T_{g,0}/T_{c}=0.75-1.00-1.50$ and (a) $Re_{\mathcal {S}}=6700$, (b) $Re_{\mathcal {S}}=2800$. At $Re_{\mathcal {S}}=6700, \varSigma _{\mathcal {F}, 0.75}=0.08 \varSigma _{\mathcal {F},1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.14\varSigma _{\mathcal {F},1.50}$. At $Re_{\mathcal {S}}=2800, \varSigma _{\mathcal {F},0.75}=0.04 \varSigma _{\mathcal {F},1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.13\varSigma _{\mathcal {F},1.50}$ (in $\varSigma _{\mathcal {F},x}, x$ is the temperature ratio).

Figure 10. Conductive and convective heat fluxes $\mathcal {F}_T$ (over the sum of the two, $\varSigma _{\mathcal {F}}$), $\mathcal {F}_{T,i=c,d}$, in the liquid region for $T_{g,0}/T_{c}=0.75-1.00-1.50$ and (a) $Re_{\mathcal {S}}=6700$, (b) $Re_{\mathcal {S}}=2800$. At $Re_{\mathcal {S}}=6700, \varSigma _{\mathcal {F},0.75}=0.12\varSigma _{\mathcal {F},1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.24\varSigma _{\mathcal {F},1.50}$. At $Re_{\mathcal {S}}=2800, \varSigma _{\mathcal {F},0.75}=0.11\varSigma _{\mathcal {F},1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.48\varSigma _{\mathcal {F},1.50}$ (in $\varSigma _{\mathcal {F},x}, x$ is the temperature ratio).

Figure 11. Conductive and convective vapour mass fluxes $\mathcal {F}_Y$ (over the sum of the two, $\varSigma _{\mathcal {F}}$), $\mathcal {F}_{Y,i=c,d}$, in the gas region for $T_{g,0}/T_{c}=0.75-1.00-1.50$ and (a) $Re_{\mathcal {S}}=6700$, (b) $Re_{\mathcal {S}}=2800$. At $Re_{\mathcal {S}}=6700, \varSigma _{\mathcal {F},0.75}=0.58\varSigma _{\mathcal {F},1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.89\varSigma _{\mathcal {F},T_{g,0}/T_c=1.50}$. At $Re_{\mathcal {S}}=2800, \varSigma _{\mathcal {F},0.75}=0.56\varSigma _{\mathcal {F},T_{g,0}/T_c=1.50}$ and $\varSigma _{\mathcal {F},1.00}=0.87\varSigma _{\mathcal {F},1.50}$ (in $\varSigma _{\mathcal {F},x}, x$ is the temperature ratio).

4.3. Effects of the shear-based Weber number

At last, we consider the effect of the shear-based Weber number on the evaporation rate. For this analysis, three Weber numbers $We_{\mathcal {S}}=0.02-0.06-0.10$ and two temperature ratios $T_{g,0}/T_c=0.75 - 1.5$ are considered and, as for the previous section, the complete thermodynamic model (i.e. VP$_a$) is employed to evaluate the gas thermophysical properties. Figure 12(a,b) reports the time evolution of $d^2$ for the three values of $We_{\mathcal {S}}$ under investigation at high and low temperatures, respectively. Increasing $We_{\mathcal {S}}$ leads to a higher evaporation rate, due to the increase of the surface area available for mass transfer as a consequence of deformation. Note, also, that the increase in the evaporation rate with $We_{\mathcal {S}}$ is more pronounced at lower temperature (see figure 12c), which can be explained by the fact that the evaporation is faster at higher temperatures, the droplets are therefore smaller, which leads to a limited deformation with respect to the cases at the same $We_{\mathcal {S}}$ and $T_{g,0}/T_c=0.75$. In other words, the effective Weber number based on the droplet diameter becomes quickly smaller at higher temperatures, so differences in nominal surface tension are compensated by the reduced droplet size. To quantify deformation while accounting for the decrease of the liquid volume, we measure the deviation from the initial spherical shape as the ratio between the instantaneous interfacial area $A$ (computed numerically as $A=|\boldsymbol {\nabla } \varPhi |\Delta x^3$) and the area of a spherical droplet with the same volume $V$, i.e. $A_{eq}={\rm \pi} ^{1/3}(6V)^{2/3}$. This calculation is performed for each droplet $q$ separately and each contribution $A_q/A_{eq,q}$ averaged according to the instantaneous number of droplets, $N_{dp}$,

(4.8)\begin{equation} \frac{A}{A_{eq}}=\frac{1}{N_{dp}}\displaystyle{\sum_{q=1}^{N_{dp}}\frac{A_q}{A_{eq,q}}}. \end{equation}

The time evolution of $A/A_{eq}$ is reported in figure 13: as expected, the ratio $A/A_{eq}$ increases with $We_{\mathcal {S}}$. For $T_{g,0}/T_c=1.5$, the data in the figure also display two peaks: at $tD_{lg,r}/d_0^2\approx 0.075-0.10$ for $We_{\mathcal {S}}=0.10$ and $tD_{lg,r}/d_0^2\approx 0.085-0.10$ for $We_{\mathcal {S}}=0.06$. These are associated to the droplet breakup as we have $N_{dp}=7$ droplets for $We_{\mathcal {S}}=0.06$ and $N_{dp}=9$ for $We_{\mathcal {S}}=0.10$ at the end of the simulations. For $We_{\mathcal {S}}=0.02$, no breakup events have been observed within the simulation time. For $T_{g,0}/T_c=0.75$, given the lower evaporation rate and consequently the higher effective Weber number, breakup events have been observed for all the cases in the interval $0.05< tD_{lg,r}/d_0^2<0.125$, finally yielding $N_{dp}=7-9-12$ for $We_{\mathcal {S}}=0.02-0.06-0.10$.

Figure 12. Time history of $(d/d_0)^2$ when varying the Weber number, $We_{\mathcal {S}}=0.02-0.06-0.10$ and for temperature ratios (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$. The black dashed curves represent the theoretical prediction of the actual evaporation rate obtained with the procedure explained in § 4.1. The length of the error bars (included for $tD_{lg,r}/d_0^2>0.025$ and $tD_{lg,r}/d_0^2>0.04$, respectively) indicates the droplet with the fastest/slowest evaporation rate among the five droplets in the simulations. (c) Evaporation enhancement (i.e. $K/K_0$) as a function of $We_{\mathcal {S}}$ for $Re_{\mathcal {S}}=6700$ and $T_{g,0}/T_c=0.75-1.50$.

Figure 13. Ratio between the instantaneous interfacial area $A$ and the area of a spherical droplet with the same volume $V$, i.e. $A_{eq}={\rm \pi} ^{1/3}(6V)^{2/3}$, for $We_{\mathcal {S}}=0.02-0.06-0.10$ and (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$.

Interestingly, for the case $T_{g,0}/T_c=0.75, A/A_{eq}$ approaches a statistically stationary value for $tD_{lg}/d_0^2>0.125$, whereas a regime value is not reached in the investigated time window for $T_{g,0}/T_c=1.50$. This behaviour is interpreted by comparing the time scale of deformation as dictated by surface tension and applied shear, $\tau _{\sigma }$, with the time scale of evaporation extracted from the simulations, $\tau _{v}$,

(4.9)\begin{equation} \varPi_{\sigma v}=\frac{\tau_{\sigma}}{\tau_v}, \end{equation}

where $\tau _{\sigma }=\rho _{g,r}(\mathcal {S}d_0)d_0^2/\sigma$ and $\tau _v=d_0^2/K$, with $K$ the evaporation at statically steady state; see figure 14. At the high temperature ratio, evaporation is faster than deformation (i.e. $\tau _{\sigma }>\tau _v$) and the surface tension does not have time to adjust the droplet shape after the mass losses due to the differential evaporation across the interface. An increasing deviation from the spherical shape is therefore observed in time. Conversely, at low temperature, the deformation time is comparable or lower than the evaporation time (i.e. $\tau _{\sigma }<\tau _v$) and the droplet shape can compensate for the local deformations induced by the evaporation mass flux (note that the droplet tends to a constant deformation, increasing with $We_{\mathcal {S}}$ and not to the spherical shape, i.e. $A/A_{eq}>1$, as we have an imposed shear). From the data in figure 14, we also notice that for a fixed $T_{g,0}/T_c$, increasing $We_{\mathcal {S}}$ accelerates the evaporation time scale (i.e. higher $\varPi _{\sigma v}$), consistently with what is observed in figure 12.

Figure 14. Ratio between the deformation time scale and the evaporation time scale (in turbulent condition) as a function of $We_{\mathcal {S}}$ for $T_{g,0}/T_c=0.75$ and $1.50$.

Given the variation of the evaporation rate with the $We_{\mathcal {S}}$, it is worth investigating whether the evaporation flux at the droplet surface is correlated to the local curvature. This has been already assessed in laminar flows analytically for evaporating droplets of spheroidal shapes (Tonini & Cossali 2013), and experimentally for sessile droplets (Sáenz et al. 2017). In turbulent flows this analysis is missing and is performed here by computing the joint probability density function (p.d.f.) of the normalized interfacial mean curvature $\kappa _{\varGamma }/\kappa _{\varGamma,eq}$ and the dimensionless interfacial mass flux $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$, where $\kappa _{\varGamma,eq}$ is the curvature of a spherical droplet with the same volume $V$ (i.e. $\kappa _{\varGamma,eq}=4/d_{eq}$, with $d_{eq}=(6V/{\rm \pi} )^{1/3}$) and $\dot {m}_{\varGamma,0}$ is the interfacial mass flux for a purely diffusion-dominated evaporation process, i.e.

(4.10)\begin{equation} \dot{m}_{\varGamma,0}=8{\rm \pi}\rho_{g,r}D_{lg,r}\log(1+B_M)/d_0. \end{equation}

Note that $\dot {m}_{\varGamma }$ is evaluated using the expression (2.5) in all the nodal points cut by the interface, whereas $\kappa _{\varGamma }$ is computed directly from its definition, i.e. $\kappa _{\varGamma }=\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {n}_{\varGamma }$. The reconstructed level-set function is employed for a more accurate estimation of $\kappa _{\varGamma }$.

From the results at high temperature (i.e. $T_{g,0}/T_c=1.5$) reported in figure 15, we do not note a clear relation between $\dot {m}_{\varGamma }$ and $\kappa _{\varGamma }$, yet there is a relatively broad distribution of mass fluxes over the droplet surface, which can be attributed to the local variations of vapour concentrations in the gas phase induced by the turbulence and by the flow anisotropy. Indeed, the regions whose local outward normal vector is parallel to the mean flow direction (i.e. droplet front) experience higher surface vapour gradients, whereas larger concentrations are found in those where the normal vector is opposite to the flow (i.e.  droplet rear). Note also that the increase of $We_{\mathcal {S}}$ corresponds to a shift of the most probable values of $\dot {m}_{\varGamma }$ and $\kappa _{\varGamma }$ towards higher values. When the gas temperature is reduced (i.e. $T_{g,0}/T_c=0.75$), the correlation between $\dot {m}_{\varGamma }$ and $\kappa _{\varGamma }$ is still weak for $We_{\mathcal {S}}=0.02$, while for the case at $We_{\mathcal {S}}=0.10$, a higher interfacial curvature is associated with a higher interfacial mass flux. Also at lower temperatures, the most probable values are located at higher $\kappa _{\varGamma }$ when increasing the Weber number. Comparing with the results at higher temperatures, figure 15, we note that $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$ is larger: turbulence enhances evaporation more efficiently at low temperature, when the evaporative time scale is longer and the ratio with the turbulent time scales decreases, i.e. $Da_v$ decreases.

Figure 15. Joint p.d.f. of the dimensionless mass flux $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$ and the normalized interfacial curvature $\kappa _{\varGamma }/\kappa _{\varGamma,eq}$ for $Re_{\mathcal {S}}=6700, T_{g,0}/T_c=1.5$ and (a) $We_{\mathcal {S}}=0.02$, (b) $We_{\mathcal {S}}=0.10$.

Finally, we display in figure 17 the joint p.d.f. from the simulations at lower Reynolds number, $Re_{\mathcal {S}}=2800$, for $T_{g,0}/T_c=0.75 - 1.5$ and $We_{\mathcal {S}}=0.02$. In this case, the reduction of $T_{g,0}/T_c$ leads to more pronounced deformation, i.e. larger values of $\kappa _{\varGamma }/\kappa _{\varGamma,eq}$ (due to the slower evaporation and higher effective $We_{\mathcal {S}}$) and to higher values of $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$ (more pronounced effects of turbulence at lower absolute evaporation rates) and a broader distribution, again attributed to fluctuations induced by the turbulence. Comparing with the results at the same $We_{\mathcal {S}}$ and $T_{g,0}/T_c$ but $Re_{\mathcal {S}}=6700$ in panels (a) of figures 15 and 16, we see that, for both temperature ratios, the reduction in $Re_{\mathcal {S}}$ restricts the range of attained values for the curvature and the ratio $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$, confirming that the increase of deformation and evaporation rate is more pronounced for bigger droplets.

Figure 16. Joint p.d.f. of the dimensionless mass flux $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$ and the normalized interfacial curvature $\kappa _{\varGamma }/\kappa _{\varGamma,eq}$ for $Re_{\mathcal {S}}=6700, T_{g,0}/T_c=0.75$ and (a) $We_{\mathcal {S}}=0.02$, (b) $We_{\mathcal {S}}=0.10$.

Figure 17. Joint p.d.f. of the dimensionless mass flux $\dot {m}_{\varGamma }/\dot {m}_{\varGamma,0}$ and the normalized interfacial curvature $\kappa _{\varGamma }/\kappa _{\varGamma,eq}$ for $Re_{\mathcal {S}}=2800$ for $We_{\mathcal {S}}=0.02$ and (a) $T_{g,0}/T_c=1.50$ and (b) $T_{g,0}/T_c=0.75$.

5. Conclusions

Fully resolved simulations of finite-size evaporating droplets are performed in HST using a weakly compressible solver for evaporating flows. The new methodology, here described, combines two features. First, an improved version of the Adams–Bashforth method in order to address the limitations of the classical formulation, already highlighted in Kasbaoui et al. (2017). The improvement allows us to match the single-phase analytical solution in the rapid distortion limit and ensures a stable integration over long times. Second, a mathematical model with the details of the numerical implementation for a two-phase evaporating system composed of an incompressible liquid and a compressible gas phase. In order to remove the acoustic time-step restriction, compressible effects are here handled in the low-Mach number limit. This new methodology is applied to investigate the behaviour of finite-size evaporating droplets in HST when changing the gas temperature over the critical temperature $T_{g,0}/T_c=0.75, 1.00$ and $1.50$, the initial droplet diameter in terms of Kolmogorov scale, $d_0/\eta =10.5$ and $21.5$ and the surface tension, quantified by the shear-based Weber number $We_{\mathcal {S}}=0.02, 0.06$ and $0.10$.

First, using the data at $d_0/\eta = 21.5$ and $We_{\mathcal {S}}=0.02$, we study the differences when employing different thermodynamic models for the gas thermophysical properties. Three approaches are investigated: a constant property model where the gas properties are kept constant and initialized with the ‘1/3 rule’ (CP) and two variable-property approaches where either the gas density, VP$_{\rho }$, or all the gas properties are allowed to vary, VP$_a$. We find that the predictions by the CP and VP$_a$ models agree well, whereas the VP$_{\rho }$ model overpredicts the evaporation rate, especially at high temperature. This overestimation occurs since the local increase of the gas density (due to evaporative cooling) is captured by the VP$_{\rho }$ model, but the decrease of the diffusion coefficient with temperature, which slows down the evaporation, is not accounted for. Moreover, by extracting the Sherwood number for the three models at $T_{g,0}/T_c = 0.75, 1.50$ and comparing it with the Frössling/Ranz-Marshall correlation (Ranz & Marshall 1952), we show that the correlation provides an excellent estimation of $Sh$ at high temperature (conduction-dominated regime), whereas it substantially underestimates $Sh$ at lower temperature (convective-dominated regime), in agreement with recent experimental observations (Méès et al. 2020).

Next, reducing the ratio $d_0/\eta$ from 21.5 to 10.5 (obtained by reducing $Re_{\mathcal {S}}$ from 6700 to 2800), we show that the ratio between the actual evaporation rate, and that computed in stagnant conditions, is always much higher than 1, while decreasing with $T_{g,0}/T_c$. Interestingly, this ratio does not approach unity at the highest temperature level, suggesting an evaporation enhancement due to turbulence also in these conditions. The variation of the droplet size, $d_0/\eta$, also affects the liquid temperature at the interface and the mean liquid temperature. For $T_{g,0}/T_c=1.0, 1.5$, the regime values of both quantities are larger for droplets of $d_0/\eta =21.5$, whereas the opposite is true at $T_{g,0}/T_c=0.75$ when the regime values of the liquid temperature are lower for the largest droplets. This is explained as the result of two competing effects: cooling by evaporation and heating from the hot gas.

Finally, by varying $We_{\mathcal {S}}$ in the range $0.02$ and $0.10$, we observe an increase in the evaporation rate for higher $We_{\mathcal {S}}$ given the larger surface area available for mass transfer. At fixed $We_{\mathcal {S}}$, this increase is more pronounced at $T_{g,0}/T_c=0.75$ due to a slower evaporation rate and higher deformation (larger effective Weber number based on the instantaneous droplet diameter). By computing the joint p.d.f. of the interfacial mass flux and curvature, a weak correlation between the two is observed at high temperature regardless of the Weber number, whereas a positive correlation is recovered at $T_{g,0}/T_c=0.75$ and $We_{\mathcal {S}}=0.10$, consistently with the theoretical prediction in laminar environments (Tonini & Cossali 2013) and the experimental observations for sessile droplets (Sáenz et al. 2017).

Note that in all the cases, the initial liquid volume fraction has been kept small (i.e. $\alpha _0\approx 0.14\,\%$), the droplets do not interact with each other and, as expected, no significant deviations from the $d^2$–law have been observed. To investigate such deviations and more in general the collective dynamics of a dense suspension of evaporating droplets we are currently exploring the same configurations presented in this work at higher volume fractions.