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. Reference Dauxois, Peacock, Bauer, Caulfield, Cenedese, Gorlé, Haller, Ivey, Linden and Meiburg2021). 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 Reference Birouk and Gökalp2006). 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 (Reference Mittal, Ni and Seo2020), Balachandar et al. (Reference Balachandar, Zaleski, Soldati, Ahmadi and Bourouiba2020), Bourouiba (Reference Bourouiba2021). Therefore, several studies have been conducted to understand how relative humidity affects evaporation and condensation (Rosti et al. Reference Rosti, Olivieri, Cavaiola, Seminara and Mazzino2020; Chong et al. Reference Chong, Ng, Hori, Yang, Verzicco and Lohse2021; Ng et al. Reference Ng, Chong, Yang, Li, Verzicco and Lohse2021) as well as on the interaction between droplets and turbulent flows (Bourouiba Reference Bourouiba2020; Rosti et al. Reference Rosti, Cavaiola, Olivieri, Seminara and Mazzino2021) 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 Reference Wong and Lin1992; Mandal & Bakshi Reference Mandal and Bakshi2012), on the interaction between droplet dispersion and vapour clouds structure (De Rivas & Villermaux Reference De Rivas and Villermaux2016; Villermaux et al. Reference Villermaux, Moutte, Amielh and Meunier2017; Sahu, Hardalupas & Taylor Reference Sahu, Hardalupas and Taylor2018) and on evaporation enhancement of super-Kolmogorov (i.e. finite-size) droplets (Marti et al. Reference Marti, Martinez, Mazo, Garman and Dunn-Rankin2017; Verwey & Birouk Reference Verwey and Birouk2018, Reference Verwey and Birouk2020). On the computational side, the most common approach has been the point particle method (Kuerten Reference Kuerten2016; Maxey Reference Maxey2017) for the studies of sub-Kolmogorov droplets in different flow configurations, i.e. forced homogeneous isotropic turbulence (Mashayek Reference Mashayek1998a; Weiss, Meyer & Jenny Reference Weiss, Meyer and Jenny2018), homogeneous shear turbulence (Mashayek Reference Mashayek1998b; Weiss et al. Reference Weiss, Giddey, Meyer and Jenny2020), turbulent channel flow (Bukhvostova et al. Reference Bukhvostova, Russo, Kuerten and Geurts2014; Russo et al. Reference Russo, Kuerten, Van Der Geld and Geurts2014) and turbulent jets (Reveillon & Demoulin Reference Reveillon and Demoulin2007; Wang, Dalla Barba & Picano Reference Wang, Dalla Barba and Picano2021). 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 Reference Elghobashi2019) and hinders the possibility to investigate directly aspects such as the reciprocal influence of nearby droplets (Lupo et al. Reference Lupo, Gruber, Brandt and Duwig2020; Chong et al. Reference Chong, Ng, Hori, Yang, Verzicco and Lohse2021) or the evaporation enhancement due to turbulence (which typically occurs for droplets of the order of the Taylor length, Birouk & Gökalp Reference Birouk and Gökalp2006). 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. (Reference Albernaz, Do-Quang, Hermanson and Amberg2017) 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. (Reference Dodd, Mohaddes, Ferrante and Ihme2021) 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. (Reference Kasbaoui, Patel, Koch and Desjardins2017), 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 Reference Pumir1996; Sekimoto, Dong & Jiménez Reference Sekimoto, Dong and Jiménez2016). The SS-HST has been recently considered in multiphase flows, e.g. in Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b) for emulsions and in Tanaka (Reference Tanaka2017), Yousefi, Ardekani & Brandt (Reference Yousefi, Ardekani and Brandt2020) 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 Reference Verwey and Birouk2020) 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. Reference Lupo, Ardekani, Brandt and Duwig2019, Reference Lupo, Gruber, Brandt and Duwig2020). 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.
(a) What is the level of approximation of the estimated evaporation rates when the gas thermophysical properties are assumed constant?
(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?(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 (Reference Scapin, Costa and Brandt2020). 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. (Reference Kasbaoui, Patel, Koch and Desjardins2017), 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 Reference Maxey1982) 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. Reference Wang, Chen, Wang and Yang2019; Ni et al. Reference Ni, Hespel, Han and Foucher2021), 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,
\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 Reference Majda and Sethian1985). 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
$$\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}$$
$$\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}$$
\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}
\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 Reference Reid, Prausnitz and Poling1987). 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.
$$\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}$$
$$\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$),
$$\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}$$
$$\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}$$
$$\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
\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}'$,
\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}'$,
\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
\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}
\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}
\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 Reference Leonard1979), 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. (Reference Scapin, Costa and Brandt2020), while the details of the weakly compressible two-phase code are reported in Dalla Barba et al. (Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021). Both implementations are based on the solver CaNS (Costa Reference Costa2018), 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. (Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012), Rosti, De Vita & Brandt (Reference Rosti, De Vita and Brandt2019a). We start from the topological equation for the phase indicator function (2.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. (Reference Scapin, Costa and Brandt2020) 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
\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. (Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012) for more details. Next, the convection by mean shear is included,
\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.
\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. Reference Scapin, Costa and Brandt2020). 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}$,
\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. (Reference Scapin, Costa and Brandt2020), 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,
\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. (Reference Scapin, Costa and Brandt2020) for details). Next, the coefficients 
$(\rho _gD_{lg})_{i\pm 1/2}$ are obtained as the arithmetic mean,
\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
\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}$,
\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 (Reference Gerz, Schumann and Elghobashi1989), Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017),
\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 Reference Tanaka2017), 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.
\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. (Reference Scapin, Costa and Brandt2020), with additional care paid to the time discretization. The adopted approach follows that proposed in Gerz et al. (Reference Gerz, Schumann and Elghobashi1989) 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,
\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,
\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,
\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. (Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021) 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 Reference Demou, Frantzis and Grigoriadis2019; Dalla Barba et al. Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021),
\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,
\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),
\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,
$$\begin{gather} \boldsymbol{RU}^{n-1} = \boldsymbol{\widetilde{RU}}^{n-1}(\boldsymbol{x} -\Delta t^{n}\mathcal{S}z\boldsymbol{e}_x), \end{gather}$$
$$\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}$$
$$\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}$$
$$\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}$$
$$\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. Reference Rosti, Ge, Jain, Dodd and Brandt2019b).
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 Reference Dodd and Ferrante2014). 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 (Reference Motheau and Abraham2016), Dalla Barba et al. (Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021), 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 (Reference Tanaka2017), Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b).
3.2. Time-step restriction
The time step 
$\Delta t^{n+1}$ is estimated from the stability constraints of the overall system,
\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. (Reference Scapin, Costa and Brandt2020), Dalla Barba et al. (Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021),
\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 Reference Hubbard, Denny and Mills1975) (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. (Reference Sekimoto, Dong and Jiménez2016), 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. Reference Scapin, Costa and Brandt2020), 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 Reference Abramzon and Sirignano1989; Lupo et al. Reference Lupo, Ardekani, Brandt and Duwig2019, Reference Lupo, Gruber, Brandt and Duwig2020) 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. (Reference Hubbard, Denny and Mills1975) 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 Reference Langmuir1918). 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. Reference Piedra, Lu, Ramos and Tryggvason2015). 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 (Reference Birouk and Gökalp2006). First, we obtain the dimensionless evaporation rate in laminar conditions 
$K_{L}$ (Abramzon & Sirignano Reference Abramzon and Sirignano1989) as
\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 Reference Ranz and Marshall1952; Birouk & Gökalp Reference Birouk and Gökalp2006),
\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 Reference Ranz and Marshall1952) and more recent works (Ng et al. Reference Ng, Chong, Yang, Li, Verzicco and Lohse2021; Wang et al. Reference Wang, Dalla Barba and Picano2021), 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. Reference Ng, Chong, Yang, Li, Verzicco and Lohse2021). 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. Reference Gökalp, Chauveau, Simon and Chesneau1992),
\begin{equation} \frac{K}{K_L} = 0.74Da_v^{{-}0.115}.\end{equation}In the above,
\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 Reference Abramzon and Sirignano1989; Birouk & Gökalp Reference Birouk and Gökalp2006). The characteristic time scale of evaporation 
$\tau _{v,L}$, is computed as
\begin{equation} \tau_{v,L} = \frac{\delta_M}{V_r},\end{equation}
where 
$\delta _M=d_0/(Sh_{RM}-2)$ (Abramzon & Sirignano Reference Abramzon and Sirignano1989) 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 Reference Lupo and Duwig2020).
Next, we examine the temporal evolution of the Sherwood number, 
$Sh$, for the three different models, see figure 6, using directly the definition
\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. (Reference Méès, Grosjean, Marié and Fournier2020) 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. Reference Lupo, Gruber, Brandt and Duwig2020).
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 Reference Verwey and Birouk2018, Reference Verwey and Birouk2020) and in the presence of a mean flow (Marti et al. Reference Marti, Martinez, Mazo, Garman and Dunn-Rankin2017). 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 (Reference Birouk and Gökalp2006) 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}$,
\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. (Reference Ni, Hespel, Han and Foucher2021), 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 Reference Wong and Lin1992; Pinheiro et al. Reference Pinheiro, Vedovoto, da Silveira Neto and van Wachem2019), 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}$,
\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}$,
\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 Reference Tonini and Cossali2013), and experimentally for sessile droplets (Sáenz et al. Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017). 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.
\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. (Reference Kasbaoui, Patel, Koch and Desjardins2017). 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 Reference Ranz and Marshall1952), 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. Reference Méès, Grosjean, Marié and Fournier2020).
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 Reference Tonini and Cossali2013) and the experimental observations for sessile droplets (Sáenz et al. Reference Sáenz, Wray, Che, Matar, Valluri, Kim and Sefiane2017).
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.
Acknowledgements
F. Picano and M.C. Esposito are thanked for the useful discussions and, the latter, also for providing the droplet tagging algorithm. K. Khamis is acknowledged for pointing out some typos in the dimensionless formulation of the mathematical model.
Funding
N.S. and L.B. acknowledge the support from the Swedish Research Council via the multidisciplinary research environment INTERFACE, Hybrid multiscale modelling of transport phenomena for energy efficient processes, grant no. 2016-06119. M.E.R. was supported by the JSPS KAKENHI, grant no. JP20K22402. The computer time was provided by SNIC (Swedish National Infrastructure for Computing) and by the National Infrastructure for High Performance Computing and Data Storage in Norway (project no. NN9561K).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Equations of state
In this appendix we report the equations of state employed for the gas density, the gas thermal diffusivity, the gas viscosity, thermal conductivity, heat capacity at constant pressure and the Span-Wagner equation of state.
A.1. Gas density
The gas density variations as a function of 
$p_{th}, T_g$ and 
$Y_{l}^v$ follow the ideal gas model
\begin{equation} \rho_g = \varPi_{p,2}\frac{p_{th}\bar{M}_{m,av}}{T_g},\end{equation}
where 
$\bar {M}_{m,av}$ is the mixture molar mass computed using a harmonic average between the gas and liquid molar masses and 
$\varPi _{p,2}=(p_{th,r}M_g)/(R_uT_{g,r})$, where 
$p_{th,r}$ and 
$T_{g,r}$ are the reference thermodynamic pressure and gas temperature, taken equal to the corresponding initial value and 
$R_u$ is the universal gas constant.
A.2. Liquid–gas diffusion coefficient
The liquid gas diffusion coefficient is computed using the relation (Reid et al. Reference Reid, Prausnitz and Poling1987)
\begin{equation} D_{lg} = \frac{T_g^{3/2}}{p_{th}}, \end{equation}
where 
$T_g$ is the gas temperature.
A.3. Gas viscosity
The gas viscosity varies with the temperature according to the simplified Sutherland's law (Reid et al. Reference Reid, Prausnitz and Poling1987),
\begin{equation} \mu_g = T_{g}^{2/3}. \end{equation}A.4. Specific heat capacity at constant pressure
The specific heat capacity at constant pressure is evaluated as a function of the temperature using the virial polynomials (Reid et al. Reference Reid, Prausnitz and Poling1987)
\begin{equation} c_{p} = A_1+A_2T_g+A_3T_g^2+A_4T_g^3+A_5T_g^4, \end{equation}
where the coefficients 
$A_1=1.012, A_2=0.0553, A_3=0.006, A_4=2\boldsymbol {\cdot } 10^{-3}$ and 
$A_5=5\boldsymbol {\cdot } 10^{-4}$ are adequate to estimate the heat capacity at constant pressure for dry air over a wide range of temperatures.
A.5. Gas thermal conductivity
As for the gas viscosity, the gas thermal conductivity is a function of the temperature only and computed as (Reid et al. Reference Reid, Prausnitz and Poling1987)
\begin{equation} k_g = \frac{\mu_gc_{pg}}{Pr}, \end{equation}
where 
$Pr$ is the Prandtl number.
A.6. Span–Wagner equation of state
In order to compute the value of 
$Y_{l}^{v}$ at the interface, see § 3.1.2, we assume that the gas mixture is ideal and composed of ideal components. Hence, 
$Y_{l,\varGamma }^{v}$ can be computed using Rault's law,
\begin{equation} Y_{l,\varGamma}^{v} = \frac{\dfrac{p_{s,\varGamma}}{p_{th}}}{\dfrac{p_{s,\varGamma}}{p_{th}}+ \left(1-\dfrac{p_{s,\varGamma}}{p_{th}}\right)\dfrac{1}{\lambda_M}},\end{equation}
where 
$\lambda _M=M_l/M_g$ is the molar mass ratio, while 
$p_{th}$ and 
$p_{s,\varGamma }$ are the thermodynamic pressure and the partial pressure of the vaporized species at the interface. Since the interface is assumed at saturation, 
$p_{s,\varGamma }$ is computed as a function of the liquid temperature at the interface. For most of the substances, the Span–Wagner equation of state represents a valid approximation over a wide range of thermodynamic pressures and temperatures (Span & Wagner Reference Span and Wagner1996),
\begin{equation} \frac{p_{s,\varGamma}}{p_{th}} = \frac{\varPi_{p,3}}{p_{th}}\exp\left[\frac{B_1\eta_{sw}+B_2\eta_{sw}^{1.5} +B_3\eta_{sw}^{2.5}+B_4\eta_{sw}^5}{1-\eta_{sw}}\right],\end{equation}
where 
$\eta _{sw}=1-T_{\varGamma }(T_{g,0}/T_c)$ and 
$\varPi _{p,3}=p_c/p_{th,r}$, with 
$T_{g,0}, T_c$ and 
$p_c$ the initial gas temperature, the critical temperature and critical pressure, respectively. The coefficients 
$B_{i=1,4}=\{-7.32714,1.82365,-2.272744,-2.711929\}$ are experimentally determined and correspond in the current study to those of pentane; see Span & Wagner (Reference Span and Wagner1996). Note that the use of the Span–Wagner model is also convenient since it does not involve the saturation temperature (like the Clausius–Clapeyron relation), which in a weakly compressible phase-changing system is not a constant anymore, but should be computed as a function of the time-varying thermodynamic pressure.
As a final remark, to justify the assumption of incompressible liquid with constant properties, we evaluate a posteriori the variations of 
$\rho _l, c_{pl}, k_l$ and 
$\mu _l$ given the maximum variation of liquid temperature discussed in § 4.2. Using the data reported in Reid et al. (Reference Reid, Prausnitz and Poling1987) and taking as reference the cases at 
$T_{g,0}/T_c=1.5$ (when the liquid temperature variation is more pronounced), these amount to 
$\Delta \rho _l/\rho _{l,0}=6\,\%, \Delta \mu _l/\mu _{l,0}=8\,\%, \Delta c_{pl}/c_{pl,0}=4\,\%, \Delta k_l /k_{l,0}=5\,\%$.
Appendix B. Derivation of the velocity divergence
The system composed of (2.2), (2.3), (2.5) and (2.4) reported in § 2 is not closed and an additional equation should be included. To derive it, we adapt the approach proposed in Majda & Sethian (Reference Majda and Sethian1985) in the context of reacting flows and recently employed for phase change in Dodd et al. (Reference Dodd, Mohaddes, Ferrante and Ihme2021). Nevertheless, since a detailed derivation is missing in literature, it is provided here for completeness. The idea is to compute the missing relation from the divergence of the velocity field, 
$\boldsymbol {u}$, which can be computed as 
$H\boldsymbol {u}_l+(1-H)\boldsymbol {u}_g$. Taking the divergence yields
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}=(\boldsymbol{u}_g-\boldsymbol{u}_l)\boldsymbol{\cdot}\boldsymbol{\nabla} H+(1-H)\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_g+H\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_l.\end{equation}
The first term on the right-hand side represents the volume variation at the interface due to phase change, while the second and the third come from the density variations in the gas and liquid phases. We start with the phase-change contribution and consider a reference frame moving with the interface. Accordingly, we decompose the vector 
$\boldsymbol {u}_g-\boldsymbol {u}_l$ along the normal 
$\boldsymbol {n}_{\varGamma }$, tangential 
$\boldsymbol {t}_{\varGamma }$ and bi-normal directions 
$\boldsymbol {n}_{\varGamma }$,
\begin{equation} \boldsymbol{u}_g-\boldsymbol{u}_l=(u_{g,n}-u_{l,n}) \boldsymbol{n}_{\varGamma}+(u_{g,t}-u_{l,t})\boldsymbol{t}_{\varGamma} +(u_{g,b}-u_{l,b})\boldsymbol{b}_{\varGamma}.\end{equation}
If the interface has zero thickness, the mass balance imposes that the velocity is continuous along the tangential and bi-normal direction (i.e. 
$u_{g,t}=u_{l,t}$ and 
$u_{g,b}=u_{l,b}$). Along the normal direction, conversely, the velocity has a discontinuity proportional to the mass flux 
$\dot {m}_{\varGamma }$ between the two phases, i.e.
\begin{gather} \left.\begin{gathered} \dot{m}_{\varGamma} = \rho_{g,\varGamma}(u_{g,n}-u_{\varGamma,n}), \\ \dot{m}_{\varGamma} = \rho_{l,\varGamma}(u_{l,n}-u_{\varGamma,n}), \end{gathered}\right\} \end{gather}
where 
$u_{\varGamma,n}$ represents the normal component of the interface velocity and 
$\rho _{i=l,g,\varGamma }$ the phase density at the interface location. Combining the previous two equations yields
\begin{equation} u_{g,n}-u_{l,n}=\dot{m}_{\varGamma}\left(\frac{1}{\rho_{g,\varGamma}} -\frac{1}{\rho_{l,\varGamma}}\right).\end{equation}Inserting (B2) and (B4) into (B1), we obtain
\begin{equation} (\boldsymbol{u}_g-\boldsymbol{u}_l)\boldsymbol{\cdot}\boldsymbol{\nabla} H=\dot{m}_{\varGamma} \left(\frac{1}{\rho_{g,\varGamma}}- \frac{1}{\rho_{l,\varGamma}}\right)\delta_{\varGamma},\end{equation}
where 
$\delta _{\varGamma }=\boldsymbol {\nabla } H\boldsymbol {\cdot }\boldsymbol {n}_{\varGamma }$.
The second term of (B1) is computed from the continuity in the gas region. If chemical reactions are absent, this reads as
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_g ={-}\frac{1}{\rho_g}\frac{{\rm D}\rho_g}{{\rm D}t}. \end{equation}
Using the equation of state for the gas density (A1), the total derivative of 
$\rho _g$ can be expanded as
\begin{equation} \frac{1}{\rho_g}\frac{{\rm D}\rho_g}{{\rm D}t} = \frac{1}{p_{th}}\frac{{\rm D}p_{th}}{{\rm D}t}-\frac{1}{T}\frac{{\rm D}T}{{\rm D}t}-\bar{M}_{m,av}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\bar{M}_{m,av}}\right). \end{equation}Combining equations (B5) and (B7) finally provides the velocity divergence in the gas region (Dodd et al. Reference Dodd, Mohaddes, Ferrante and Ihme2021),
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_g={-}\frac{1}{p_{th}}\frac{{\rm d}p_{th}}{{\rm d}t}+\frac{1}{T}\frac{{\rm D}T}{{\rm D}t}+\bar{M}_{m,av}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\bar{M}_{m,av}}\right). \end{equation}
Note that, given the low-Mach assumption, the total derivative of 
$p_{th}$ in (B8) has been replaced by the time derivative of 
$p_{th}$. A similar strategy can be employed to compute 
$\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {u}_l$; however, we consider here constant and uniform liquid density and, therefore, 
$\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {u}_l=0$ in (B1) and 
$\rho _{l,\varGamma }=\rho _l$.
Differently from what it is commonly done in literature (Motheau & Abraham Reference Motheau and Abraham2016), we found preferable to replace the total derivative of the molar mass and of the temperature with the corresponding spatial derivatives. By doing so, this solution completely removes the time discretization errors when computing the terms in (B8). For the contribution due to composition, we employ (2.3) and we treat the mixture of inert gas and vapour as ideal, so that the mixture molar mass is given by the harmonic average of the molar mass of each component (Reid et al. Reference Reid, Prausnitz and Poling1987). Accordingly, using the relation 
$Y_g=Y_N=1-\sum _{j=1}^{N-1}Y_{l,j}^v$ the third term on the right-hand side of (B8) becomes
\begin{equation} \bar{M}_{m,av}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\bar{M}_{m,av}}\right) = \bar{M}_{m,av}\sum_{j=1}^{N-1}\left(\frac{1}{\lambda_{M,j}}-1\right)\frac{{\rm D}Y_{l,j}^v}{{\rm D}t}, \end{equation}
where 
$\lambda _{M,j}=M_{l,j}/M_g$. For the temperature 
$T$, we start from the generic enthalpy equation for the gas phase,
\begin{equation} \rho\frac{{\rm D}h}{{\rm D}t} = \frac{1}{RePr}\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla} T) + \varPi_{p,1}\frac{{\rm d}p_{th}}{{\rm d}t}+\frac{1}{ReSc}\sum_{j=1}^N h_{j}^v\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg}\boldsymbol{\nabla} Y_{j}).\end{equation}In a weakly compressible system, being the enthalpy of the gas phase function of temperature, thermodynamic pressure and vapour composition, the left-hand side of (B10) can be expanded as
\begin{align} \rho\frac{{\rm D}h}{{\rm D}t} &=\rho\left(\frac{\partial h}{\partial T}\right)_{Y_{l,j}^v,p_{th}}\frac{{\rm D}T}{{\rm D}t} + \rho\left(\frac{\partial h}{\partial p_{th}}\right)_{Y_{l,j}^v,T}\frac{{\rm d}p_{th}}{{\rm d}t}+\rho\sum_{j=1}^{N-1}\left(\frac{\partial h}{\partial Y_{l,j}^v}\right)_{T,Y_{i\neq j}}\frac{{\rm D}Y_{l,j}^v}{{\rm D}t} \nonumber\\ &=\rho c_p\frac{{\rm D}T}{{\rm D}t} + \rho(1-\beta T)\frac{{\rm d}p_{th}}{{\rm d}t} + \rho\sum_{j=1}^{N-1}\left(\frac{\partial h}{\partial Y_{l,j}^v}\right)_{T,p_{th},Y_{l,i\neq j}^v}\frac{{\rm D}Y_{l,j}^v}{{\rm D}t}. \end{align}
Note that for an ideal gas, the isothermal compressibility coefficient 
$\beta =1/T$ and, thus, 
$(1-\beta T)=0$, which can be omitted. Considering the term varying with the composition 
$Y_{l,j}^v$, this can be recast as (see Lupo et al. (Reference Lupo, Ardekani, Brandt and Duwig2019) for details)
\begin{equation} \sum_{j=1}^{N-1}\left(\frac{\partial h}{\partial Y_{l,j}^v}\right)_{T,p_{th},Y_{l,i\neq j}^v}\frac{{\rm D}Y_{l,j}^v}{{\rm D}t} = \sum_{j=1}^{N-1}(h_{v,j}-h_g)\frac{{\rm D}Y_{l,j}^v}{{\rm D}t}. \end{equation}Inserting (B12) and (B11) into (B10) and using (2.3), we get
\begin{equation} \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_g}{ReSc}\sum_{j=1}^{N}D_{lg,j}\boldsymbol{\nabla} h_{j}\boldsymbol{\cdot}\boldsymbol{\nabla} Y_{j}\right)(1-H). \end{equation}
Note that (B13) is the same as (2.4) except for the source term due to phase change which acts at the interface location. Finally, by combining (B13), (B9) (B8) in expression (B1), the velocity divergence 
$\boldsymbol {u}$ can be expressed as
\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} &= \underbrace{\dot{m}_{\varGamma}\left(\frac{1}{\rho_{g,\varGamma}}- \frac{1}{\lambda_{\rho}}\right)\delta_{\varGamma}}_{f_{\varGamma} (\boldsymbol{x}_{\varGamma},t)} \nonumber\\ &\quad +\frac{1}{p_{th}}\underbrace{\frac{1}{ReSc} \frac{\bar{M}_{m,av}}{\rho_g}\sum_{j=1}^{N-1} \left(\frac{1}{\lambda_{M,j}}-1\right)\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_gD_{lg,j}\boldsymbol{\nabla} Y_{l,j}^v)}_{f_Y(\boldsymbol{x},t)}(1-H) \nonumber\\ &\quad +\frac{1}{p_{th}}\underbrace{\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_g}{Sc}\sum_{j=1}^ND_{lg,j}\boldsymbol{\nabla} h_j\boldsymbol{\cdot}\boldsymbol{\nabla} Y_j\right]}_{f_T(\boldsymbol{x},t)}(1-H) \nonumber\\ &\quad -\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}(1-H), \end{align}
where 
$f_{\varGamma }, f_Y$ and 
$f_T$ are functions representing the different contributions to the total velocity divergence, i.e. phase change, (
$f_{\varGamma }$), change of the gas density due to composition (
$f_Y$) and due to temperature (
$f_T$).
The last step is to derive an expression for the thermodynamic pressure, 
$p_{th}$. For an open domain, 
$p_{th}$ is constant, whereas for a closed or triperiodic domain, it can be obtained by imposing the volume conservation on the global domain 
$V$. This can formally be expressed as a constrain on the velocity divergence
\begin{equation} \int_V\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u} \,{\rm d} V = 0.\end{equation}
By employing (B15), we can easily rearrange the rate of change of 
$p_{th}$ as
\begin{equation} \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{equation}Appendix C. Verification/validation of the low-Mach solver with phase change
In this appendix we provide a verification and two validation cases of the weakly compressible code employed to perform the numerical simulations in the present work.
C.1. Droplet evaporation due to a prescribed, constant mass flux
The set-up of this verification case consists of a two-dimensional circular droplet with initial diameter 
$d_0$, which evaporates due to a prescribed, constant mass flux 
$\dot {m}_{\varGamma,A}$. In this case, the momentum equation is decoupled from the transport of energy and vapour mass fraction equations and it is straightforward to show that the droplet diameter evolves as
\begin{equation} \frac{{\rm d}(t)}{d_0} = 1-\left(\frac{2\dot{m}_{\varGamma,A}}{d_0\rho_l}\right)t. \end{equation}
In our previous work (Scapin et al. Reference Scapin, Costa and Brandt2020), we reproduce this test case using a zero-pressure outflow boundary condition, whereas here we prescribe periodic boundary conditions. Accordingly, as evaporation starts, the thermodynamic pressure builds up and eventually reaches a stationary value when the droplet is completely evaporated. By setting 
$f_Y$ and 
$f_T$ to zero in (B16) and after some manipulations a nonlinear ordinary differential equation can be derived for 
$p_{th}$,
\begin{equation} \frac{{\rm d}p_{th}(t)}{{\rm d}t} - \frac{42\dot{m}_{\varGamma,A}{\rm \pi} d^2(t)}{30G_{tot}-5\rho_l{\rm \pi} d^3(t)}\left(1-\frac{p_{th}(t)}{\rho_l(R_u/M_g)T_0}\right)p_{th}(t)=0,\end{equation}
where 
$d(t)$ is computed using (C1), 
$G_{tot}$ is the mass inside the system, taken equal to its initial value (i.e. 
$G_{tot}=\rho _{g,0}V_{g,0}+\rho _lV_{l,0}$), 
$T_0$ and 
$(R_u/M_g)$ are chosen so that the group 
$p_{th,0}/(\rho _{g,0}(R_u/M_g)T_0)=1$. Equation (C2) is here solved with the fourth-order Runge–Kutta scheme. The test has been repeated for three density ratios 
$\lambda _{\rho }=10-50-100$, with the remaining dimensionless physical parameters 
$Re=25, We=0.10, \lambda _{\mu }=50$; gravity is absent. The governing equations are solved in a square domain 
$[-2d_0;2d_0]^2$, discretized with 
$128\times 128$ grid points. The droplet is held at the centre of the domain. Results are reported in figure 18: an excellent agreement between the numerical and the analytical solutions is found for all the cases.
Figure 18. (a) Temporal evolution of the instantaneous droplet diameter over the initial value, 
$d/d_0$. (b) Temporal evolution of the thermodynamic pressure over the initial value, 
$p_{th}/p_{th,0}$.
C.2. Hexadecane static droplet evaporation in a hot gas
As a validation case, we reproduce the numerical test proposed in Ni et al. (Reference Ni, Hespel, Han and Foucher2021) consisting of a static hexadecane droplet of initial diameter 
$d_0=550\ \mathrm {\mu }\textrm {m}$, which evaporates in dry air at 
$T_g=673$ K. For this test case, we set 
$Re=25, We=0.1, Pr=Sc=1, \lambda _{\rho }=770, \lambda _{\mu }=202, \lambda _k=20, \lambda _{c_p}=2.1$. Zero-pressure outflow boundary conditions are prescribed (therefore, 
$p_{th}=p_{th,0}$) and the assumption of constant liquid bulk density is relaxed (see Ni et al. (Reference Ni, Hespel, Han and Foucher2021) for more details). The whole set of governing equations are solved in a square domain 
$[-10d_0;10d_0]^2$ using 
$256\times 256$ grid points. Once evaporation starts, the liquid droplet undergoes a local initial expansion (i.e. 
$(d/d_0)^2>1$) until 
$t/d_0^2\approx 2$, after which the 
$d^2$ regime is approached (see figure 19a). Overall, good agreement between our simulations and the reference data is observed, confirming once more the validity of our numerical algorithm.
C.3. Single droplet evaporating in HIT
As a final validation case, we reproduce the experimental results in Verwey & Birouk (Reference Verwey and Birouk2018), reproduced also in Dodd et al. (Reference Dodd, Mohaddes, Ferrante and Ihme2021), consisting of a single 
$n$-heptane droplet of initial diameter 
$d_0=200\ \mathrm {\mu }\textrm {m}$, which evaporates in nitrogen at 
$T_g=348$ K and 
$p_{th,0}=10$ bar. For this test case, we consider HIT sustained by an artificial forcing (Podvigina & Pouquet Reference Podvigina and Pouquet1994) at a Reynolds number based on the Taylor's microscale 
$Re_{\lambda }=u_{rms}'\lambda /\nu _{g,r}\approx 32$. The initial droplet Weber number 
$We_{rms}=\rho _{g,r}u_{rms}'^2d_0/\sigma$ is set equal to 
$0.0044, Pr=0.715, Sc=2.81, Ste=1.05, \lambda _{\rho }=25.27, \lambda _{\mu }=11.82, \lambda _k=3.72$ and 
$\lambda _{c_p}=2.45$. The whole set of governing equations are solved in a cubic domain 
$[-16d_0;16d_0]^3$ using 
$768^3$ grid points corresponding to about 
$48$ grid points per initial diameter. Figure 19(b) reports 
$(d/d_0)^2$ as a function of 
$t/d_0^2$ for both the experimental and numerical results. A good agreement is observed with a deviation of evaporation rate 
$K_{DNS}$ of less then 
$7\,\%$ with respect to the experimental data, which is considered a satisfactory result and within the uncertainty of the measurement (Verwey & Birouk Reference Verwey and Birouk2018).
Appendix D. Improved Adams–Bashforth scheme for HST simulations
Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017) show that the classical Adams–Bashforth scheme (AB2) employed in Gerz et al. (Reference Gerz, Schumann and Elghobashi1989) is not suitable for DNS of HST. Using as a benchmark the Kelvin modes derived in the framework of the RDT (Maxey Reference Maxey1982; Isaza & Collins Reference Isaza and Collins2009), they show that the method fails to reproduce the analytical solution with an unbounded growth of the error. In this appendix we propose a modification of the classical Adams–Bashforth scheme able to reproduce the Kelvin modes and, more generally, stable and accurate simulations.
D.1. Proof of the numerical stability
We follow the approach proposed in Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017) and compute the time discretization error introduced by the classical Adams–Bashforth scheme using the simplified equation
\begin{equation} \frac{\partial u}{\partial t} + \mathcal{S}z\frac{\partial u}{\partial x} = Au, \end{equation}
where 
$A$ is a positive constant and 
$\mathcal {S}$ is the applied shear. As noted in Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017), (D1) allows us to remove the complications due to the pressure gradient while preserving some important key features of homogeneously sheared flows, i.e. the exponential growth of the turbulent kinetic energy (Maxey Reference Maxey1982). The solution of (D1), 
$u=\hat {u}(t)\exp (\textrm {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x})$, needs to satisfy two conditions for the amplitude 
$\hat {u}(t)$ and for the wave vector 
$\boldsymbol {k}$,
$$\begin{gather} \frac{{\rm d}\hat{u}}{{\rm d}t} = A\hat{u}, \end{gather}$$
$$\begin{gather}\frac{{\rm d}\boldsymbol{k}}{{\rm d}t}\boldsymbol{\cdot}\boldsymbol{x} ={-}\mathcal{S}z\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{e}_x. \end{gather}$$
Assuming that the discretization error evolves as a wave of amplitude 
$\hat {\varepsilon }$ and wave vector 
$\boldsymbol {k}$, at the 
$n$th time step
\begin{equation} \varepsilon^n=\hat{\varepsilon}^n\exp({\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}\boldsymbol{x}). \end{equation}Employing the classical Adams–Bashforth integration scheme, we obtain a first intermediate error neglecting the mean shear contribution,
\begin{equation} \bar{\varepsilon}^{n+1}=\varepsilon^n+A\Delta t^{n+1} [\gamma_1\boldsymbol{RU}(\hat{\varepsilon}^n,\boldsymbol{k}^n)- \gamma_2\boldsymbol{RU}(\hat{\varepsilon}^{n-1},\boldsymbol{k}^{n-1})],\end{equation}
where 
$\boldsymbol {RU}$ represents the discrete operator for the spatial terms in the governing equation of the error and 
$\gamma _1=(1+0.5\Delta t^{n+1}/\Delta t^n)$ and 
$\gamma _2=0.5\Delta t^{n+1}/\Delta t^n$ are the two coefficients of the Adams–Bashforth scheme for a variable time-step size. Note again that, at this stage, the mean shear contribution is not included in 
$\boldsymbol {RU}$.
Following the approach proposed by Gerz et al. (Reference Gerz, Schumann and Elghobashi1989), we set 
$\boldsymbol {RU}(\hat {\epsilon }^n,\boldsymbol {k}^n)=\hat {\varepsilon }^n \exp (\textrm {i}\boldsymbol {k}^n\boldsymbol {\cdot }\boldsymbol {x})$ and 
$\boldsymbol {RU}(\hat {\epsilon }^{n-1},\boldsymbol {k}^{n-1})= \hat {\varepsilon }^{n-1}\exp (\textrm {i}\boldsymbol {k}^{n-1}\boldsymbol {\cdot }\boldsymbol {x})$. Moreover, by using the identity 
$\exp (\textrm {i}\boldsymbol {k}^{n-1}\boldsymbol {\cdot }\boldsymbol {x})=\exp (-\textrm {i}\Delta \boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x})\exp (\textrm {i}\boldsymbol {k}^n\boldsymbol {\cdot }\boldsymbol {x})$, (D4) becomes
\begin{equation} \bar{\varepsilon}^{n+1}=\{\varepsilon^n+A\Delta t^{n+1}[\gamma_1\hat{\varepsilon}^n-\gamma_2\hat{\varepsilon}^{n-1}\exp(-{\rm i}\Delta \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x})]\}\exp({\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}\boldsymbol{x}),\end{equation}
where 
$\Delta \boldsymbol {k}=\boldsymbol {k}^n-\boldsymbol {k}^{n-1}$. Next, the mean shear contribution is included and this results in an error at time 
$n+1$,
\begin{align} \varepsilon^{n+1}&=\overbrace{\{\varepsilon^n+A\Delta t^{n+1} [\gamma_1\hat{\varepsilon}^n-\gamma_2\hat{\varepsilon}^{n-1}\exp(-{\rm i}\Delta \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{x}-\Delta t^{n+1}\mathcal{S}z \boldsymbol{e}_x))]\}}^{\hat{\varepsilon}^{n+1}}\nonumber\\ & \quad \boldsymbol{\cdot}\exp[\underbrace{{\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}(\boldsymbol{x}-\Delta t^{n+1}\mathcal{S}z\boldsymbol{e}_x)}_{{\rm i}\boldsymbol{k}^{n+1}}]. \end{align}
If we analyse (D6), we see that the new wave vector 
$\boldsymbol {k}^{n+1}$ satisfies the first condition, i.e. (D2a). In fact,
\begin{equation} \boldsymbol{k}^n\boldsymbol{\cdot}(\boldsymbol{x}-\Delta t^{n+1}\mathcal{S}z\boldsymbol{e}_x) = \boldsymbol{k}^{n+1}\boldsymbol{\cdot}\boldsymbol{x}\rightarrow \frac{\boldsymbol{k}^{n+1}-\boldsymbol{k}^n}{\Delta t^{n+1}}\boldsymbol{\cdot}\boldsymbol{x} ={-}\mathcal{S}z\boldsymbol{k}^n\boldsymbol{\cdot}\boldsymbol{e}_x. \end{equation}
Conversely, the amplitude at the new time level, 
$\hat {\varepsilon }^{n+1}$, does not satisfy the condition (D2b) since it contains a spatial-dependent term, 
$\Delta \boldsymbol {k}$, which is caused by the difference in orientation between the wave vector at the current and old time level, 
$n$ and 
$n-1$.
The approach proposed in this work improves the method by Gerz et al. (Reference Gerz, Schumann and Elghobashi1989) by removing the spatial-dependent term, 
$\Delta \boldsymbol {k}$. To this end, the third term on the right-hand side of (D4) is first shear mapped to the current time level,
\begin{equation} \hat{\varepsilon}^{n-1}\exp({\rm i}\boldsymbol{k}^{n-1}\boldsymbol{\cdot}(\boldsymbol{x}-\Delta t^{n}\mathcal{S}z\boldsymbol{e}_x))=\hat{\varepsilon}^{n-1}\exp({\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}\boldsymbol{x}). \end{equation}
Equation (D8) is key to obtain two wave errors of different amplitudes (i.e. 
$\hat {\varepsilon }^n$ and 
$\hat {\varepsilon }^{n-1}$) but equal wave vector, i.e. 
$\boldsymbol {k}^n$. Inserting (D8) in (D4) gives
\begin{equation} \bar{\varepsilon}^{n+1}=\{\varepsilon^n+A\Delta t^{n+1}[\gamma_1\hat{\varepsilon}^n-\gamma_2\hat{\varepsilon}^{n-1}]\}\exp({\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}\boldsymbol{x}). \end{equation}
Next, the shear-mapping step is applied to (D9), providing an error at 
$n+1$,
\begin{equation} \varepsilon^{n+1}=\overbrace{[\varepsilon^n+A\Delta t^{n+1} (\gamma_1\hat{\varepsilon}^n-\gamma_2 \hat{\varepsilon}^{n-1})]}^{\hat{\varepsilon}^{n+1}} \exp[\underbrace{{\rm i}\boldsymbol{k}^n\boldsymbol{\cdot}(\boldsymbol{x}-\Delta t^{n+1} \mathcal{S}z\boldsymbol{e}_x)}_{{\rm i}\boldsymbol{k}^{n+1}}].\end{equation}
As expected, and similarly to what is observed for (D5), (D10) satisfies the condition for 
$\boldsymbol {k}^{n+1}$ given by (D2b). This time, however, also the condition (D2a) is met since the term 
$\hat {\varepsilon }^{n+1}$ does not contain any spatial-dependent term, i.e. 
$\boldsymbol {k}$. The effectiveness of the proposed approach is shown in the next subsection with an analytical benchmark.
D.2. Validation
To validate the modified algorithm proposed here, we employ the two-dimensional analytical Kelvin modes derived in the framework of the RDT,
\begin{gather} \left.\begin{gathered} v^K(\boldsymbol{x},t) = v_0^K\left(\dfrac{k_0}{k}\right)^2\exp\left[-\dfrac{t\mathcal{S}}{Re_{\lambda}} \left(k_0^2-k_{x,0}k_{y,0}\mathcal{S}t+\dfrac{k_{x,0}^2 \mathcal{S}^2t^2}{3}\right)\right]\exp({\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}), \\ u^K(\boldsymbol{x},t) ={-}\dfrac{k_y(t)}{k_x}v^K(\boldsymbol{x},t), \\ k_y = k_{y,0}-k_x\mathcal{S}t, \end{gathered}\right\} \end{gather}
where 
$Re_{\lambda }$, defined as 
$\mathcal {S}\lambda _0^2/\nu$, is the Reynolds number based on the perturbation wavelength 
$\lambda _0=\sqrt {(2{\rm \pi} /k_x)^2+(2{\rm \pi} /k_{y,0})^2}$, taken as a reference length. As done by Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017), we consider a two-dimensional square box 
$l_x=l_y=2{\rm \pi}$, discretized with 
$N=128$ points along each direction. The analytical solution (D11) is employed as the initial condition, with initial vertical amplitude 
$v_0^K=\sqrt {2}/2$, initial wave vector 
$\boldsymbol {k}_0=(4,1)$ and 
$Re_{\lambda }=208$. We compare three time-integration methods: (1) the Adams–Bashforth (AB2) in the two variants (classical and modified), (2) the third-order Runge–Kutta (RK3) method proposed in Tanaka (Reference Tanaka2017), Yousefi et al. (Reference Yousefi, Ardekani and Brandt2020) and (3) the explicit Crank–Nicholson (CN2) method proposed in Kasbaoui et al. (Reference Kasbaoui, Patel, Koch and Desjardins2017). In all cases, the momentum equation is advanced with a constant time step, 
$\Delta t=0.25\,(\mathcal {S}N)^{-1}$. Figures 20–22 report the horizontal mode obtained with the three different time-integration schemes, which according to the RDT should go to zero at 
$t\mathcal {S}=1$. Using the modified AB2 scheme and CN2, we obtain almost identical results and zero horizontal mode up to machine precision at 
$t\mathcal {S}=1$. Conversely, using the original approach by Gerz (figure 21), the normal mode does not go to zero for 
$t\mathcal {S}=1$ and, for 
$t\mathcal {S}=5$, the numerical errors largely affect the solution.
Figure 20. Temporal evolution of the horizontal mode for 
$t\mathcal {S}=1, 3$ and 
$5$ using the modified Adams–Bashforth scheme.
Figure 21. Temporal evolution of the horizontal mode for 
$t\mathcal {S}=1, 3$ and 
$5$ using the standard Adams–Bashforth scheme.
Figure 22. Temporal evolution of the horizontal mode for 
$t\mathcal {S}=1, 3$ and 
$5$ using the explicit Crank–Nicolson method.
The excellent agreement between the results obtained with our approach and the analytical solution from RDT are reported in figure 23. A more detailed analysis considers the evolution in time of the relative error, reported in figure 24 for the horizontal and vertical Kelvin modes. The modified AB2, RK3 and CN2 provide similar errors, especially for 
$t\mathcal {S}>2$, while for 
$t\mathcal {S}\leq 2$, the CN2 method appears to be slightly more accurate. On the contrary, the error of the classical AB2 is well above that from the other three methods, with a local peak at 
$t\mathcal {S}=1$ in the horizontal mode 
$u^K$, as already observed in the first panel of figure 21. Furthermore, the error evolves at a much faster pace, which is the main reason for the poor performance of the classic AB2 in more demanding simulations.
Figure 23. Comparison between the numerical solution obtained with AB2 (proposed method) and the analytical results from RDT.
Figure 24. Temporal evolution of the discretization of error (in log scale) for (a) the horizontal mode (i.e. 
$u$), (b) the vertical mode (i.e. 
$v$) using the three time-integration methods, AB2 - Gerz et al. and AB2 - Modified, RK3 and CN2.
Open access
