1. Introduction
Aerodynamic breakup of drops is a classical multiphase flow problem commonly encountered in various spray applications such as liquid fuel injection and spray cooling. When the relative velocity between the drop and the surrounding gas is sufficiently high, the drop deforms and may even break as it accelerates. In many atomization problems, bulk liquids first break into larger drops, which then undergo secondary atomization into even smaller child droplets. In the present study, the term aerobreakup refers to any scenario in which the drop undergoes a topology change, such as forming a hole. For cases near the critical breakup condition, the drop may not fully atomize into a large number of smaller droplets. For theoretical and numerical studies, drop aerobreakup is often formulated in an idealized configuration, i.e. a stationary drop suddenly exposed to a uniform gaseous stream (O'Rourke & Amsden Reference O'Rourke and Amsden1987; Jain et al. Reference Jain, Prakash, Tomar and Ravikrishna2015, Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Marcotte & Zaleski Reference Marcotte and Zaleski2019; Rimbert et al. Reference Rimbert, Castrillon Escobar, Meignen, Hadj-Achour and Gradeck2020). Experimentally, this sudden change in relative velocity is achieved using a shock wave (Hsiang & Faeth Reference Hsiang and Faeth1995; Theofanous, Li & Dinh Reference Theofanous, Li and Dinh2004) or a jet (Flock et al. Reference Flock, Guildenbecher, Chen, Sojka and Bauer2012; Opfer et al. Reference Opfer, Roisman, Venzmer, Klostermann and Tropea2014; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). The shape evolution and acceleration of the drop are governed by key dimensionless parameters, including the Weber (
${We}$), Ohnesorge (
${Oh}$) and Reynolds (
${Re}$) numbers, and the liquid-to-gas density ratio (
$\eta$). Extensive experimental, theoretical and numerical studies have been dedicated to investigating drop aerobreakup (Guildenbecher, López-Rivera & Sojka Reference Guildenbecher, López-Rivera and Sojka2009; Theofanous Reference Theofanous2011). In recent years, high-fidelity experimental diagnostics and large-scale numerical simulations have emerged (Jackiw & Ashgriz Reference Jackiw and Ashgriz2021; Chirco et al. Reference Chirco, Maarek, Popinet and Zaleski2022; Ling & Mahmood Reference Ling and Mahmood2023; Tang, Adcock & Mostert Reference Tang, Adcock and Mostert2023), providing important details to understand the underlying physics of interfacial multiphase flow. While 
${We}$ is often used to characterize different breakup modes, the effects of 
$\eta$ (Jain et al. Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Marcotte & Zaleski Reference Marcotte and Zaleski2019), 
${Re}$ (Aalburg, Van Leer & Faeth Reference Aalburg, Van Leer and Faeth2003; Strotos et al. Reference Strotos, Malgarinos, Nikolopoulos and Gavaises2016b) and heating (Strotos et al. Reference Strotos, Malgarinos, Nikolopoulos and Gavaises2016a) may also be important for low-
${We}$ scenarios.
In applications where drop aerobreakup occurs in high-temperature environments, such as injection of liquid fuel into combustion chambers, drop vaporization occurs simultaneously as the drop accelerates and deforms (Strotos et al. Reference Strotos, Malgarinos, Nikolopoulos and Gavaises2016a). Phase change is known to significantly influence canonical interfacial instabilities, such as Rayleigh–Taylor and Kelvin–Helmholtz instabilities (Hsieh Reference Hsieh1978; Pillai & Narayanan Reference Pillai and Narayanan2018). Since interfacial instability plays a significant role in drop deformation and breakup, it is expected that vaporization can also have a substantial influence on drop aerobreakup. However, a detailed characterization of the vaporization effect is lacking in the literature. The vaporization of drop liquid on the surface can be driven by heat transfer as a result of the temperature gradient or by mass transfer due to the vapour concentration gradient. This study focuses on the former scenario, where the rate of vaporization is determined by the imbalance of heat fluxes across the interface. This introduces additional parameters, such as the Prandtl (
${Pr}$) and Stefan (
${St}$) numbers. Experimentally, the vaporization of a spherical drop in high-temperature environments has been studied, resulting in commonly used empirical relations for the drop vaporization rate (Renksizbulut & Yuen Reference Renksizbulut and Yuen1983; Abramzon & Sirignano Reference Abramzon and Sirignano1989). These correlations are strictly valid in the zero-
${We}$ limit and will significantly underestimate the vaporization rate for deformable drops with finite 
${We}$ (Boyd, Becker & Ling Reference Boyd, Becker and Ling2023; Setiya & Palmore Reference Setiya and Palmore2023). For drops with finite 
${We}$ and 
${St}$, the strong coupling between drop shape deformation and vaporization complicates the problem. On the one hand, significant drop deformation will change the gas flow around the drop and the temperature distribution in the thermal boundary layer near the surface, which leads to the time-varying and non-uniform distribution of the local vaporization rate on the drop surface. However, vaporization, in particular, the asymmetric vaporization rate on the windward and leeward sides of the drop, modulates the interfacial dynamics and drop acceleration, which in turn change the drop deformation and the breakup criteria. The goal of this study is to comprehensively investigate these complex interactions between drop aerobreakup and vaporization through detailed parametric numerical simulations. Beyond the commonly considered 
${We}$, we will systematically vary two other key dimensionless parameters, i.e. 
${St}$ and 
$\eta$, to fully characterize the vaporization effect.
2. Methodology
To accurately predict the aerobreakup of vaporizing drops, two-phase interfacial flows with phase change must be resolved faithfully with the proper physical models and numerical methods. The simulation approach used for the present simulations has been presented in recent work of the authors (Boyd & Ling Reference Boyd and Ling2023; Boyd et al. Reference Boyd, Becker and Ling2023), therefore only a brief summary will be given here. The Navier–Stokes equations with surface tension are solved for two-phase flows using the one-fluid approach. The momentum and continuity equations are given as
$$\begin{gather} \rho \left( \frac{\partial \boldsymbol{u}}{\partial t}+ \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} \right) ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} (2\mu \boldsymbol{\mathsf{D}}) + \rho \boldsymbol{g}+ \sigma \kappa \delta_{\gamma} \boldsymbol{n}_{\gamma}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = s_{\gamma} \left( \frac{1}{\rho_g}-\frac{1}{\rho_l} \right), \end{gather}$$
where 
$\boldsymbol {u}$, 
$p$, 
$\mu$, 
$\rho$, 
$\sigma$ and 
$\kappa$ are the velocity, pressure, dynamic viscosity, density, surface tension coefficient and interfacial curvature. The deformation tensor is defined as 
$\boldsymbol{\mathsf{D}} = (\boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^{\rm T})/2$. The interface normal and the interface Dirac distribution function are denoted by 
$\boldsymbol {n}_{\gamma }$ and 
$\delta _{\gamma }$, respectively. The subscripts 
$l$ and 
$g$ refer to the liquid and gas properties, respectively, while the subscript 
$\gamma$ denotes properties related to the surface.
The two different phases are distinguished by the liquid volume fraction 
$f$, which follows the advection equation with the phase-change-induced source term, i.e.
\begin{equation} \frac{\partial f}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\, f \boldsymbol{u}\right ) = \frac{-s_{\gamma}}{\rho_l}. \end{equation}
In (2.2) and (2.3), 
$s_{\gamma }$ is the mass flow rate per unit volume, which in turn depends on the mass flux at the interface (
$j_{\gamma }$) and the interfacial area density (
$\phi _{\gamma }$) as
\begin{equation} s_{\gamma} = j_{\gamma} \phi_{\gamma}, \end{equation}
where 
$\phi _{\gamma }=A_{\gamma }/V_c$, with 
$A_{\gamma }$ the liquid–gas interface area in a cell with volume 
$V_c$. The source term on the right-hand side of (2.2) is responsible for the generation of the Stefan flow; while the counterpart on the right-hand side of (2.3) accounts for the additional change in the interface location due to phase change.
The phase change at the interface is driven by heat transfer, so 
$j_{\gamma }$ is determined by the gas temperature gradient at the interface,
\begin{equation} j_{\gamma} =\frac{1}{h_{l,g}}(k_l\,\boldsymbol{\nabla} T|_{l,\gamma} \boldsymbol{\cdot} \boldsymbol{n}_{\gamma} - k_g\,\boldsymbol{\nabla} T|_{g,\gamma} \boldsymbol{\cdot} \boldsymbol{n}_{\gamma}), \end{equation}
where 
$T$, 
$k$ and 
$h_{l,g}$ are the temperature, thermal conductivity and specific latent heat of vaporization, respectively. The gas and liquid temperature fields required to calculate 
$j_\gamma$ are obtained by solving the energy conservation equation for both the liquid and gas phases:
$$\begin{gather} \rho_g C_{p,g} \left( \frac{\partial{T_g}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} T_g \right) = \boldsymbol{\nabla}\boldsymbol{\cdot} (k_g\,\boldsymbol{\nabla} T_g) , \, \end{gather}$$
$$\begin{gather}\rho_l C_{p,l}\left(\frac{\partial{T_l}}{\partial{t}} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} T_l \right) = \boldsymbol{\nabla}\boldsymbol{\cdot} (k_l\,\boldsymbol{\nabla} T_l) , \end{gather}$$
where 
$C_p$ is the isobaric-specific heat. Since 
$T_g$ and 
$T_l$ are solved only in the gas and liquid regions, the gas–liquid interface is treated as an embedded boundary where the temperature remains as the saturation temperature 
$T_{sat}$.
To rigorously resolve the sharp and vaporizing interface, a novel volume-of-fluid (VOF) method was used (Boyd & Ling Reference Boyd and Ling2023; Boyd et al. Reference Boyd, Becker and Ling2023). As the projection method is used to incorporate the continuity equation, the pressure Poisson equation is solved. The volumetric source due to vaporization was added to the pressure equation to account for the non-zero divergence of the velocity near the interface and the resulting Stefan flow. To avoid contaminating the velocity at the interfacial cells, the vaporization-induced volumetric source is distributed to the pure gas and liquid cells adjacent to the interface in a conservative manner. This treatment will guarantee that the velocities at the interface cells are correctly represented and can be used directly in VOF advection. The energy equations for both phases are solved with the Dirichlet boundary condition at the interface. To avoid calculating the gradient across the interface, we estimate the temperature gradient for each phase by extrapolating the neighbouring pure cells for the corresponding phase. The contribution of vaporization to the motion of the interface towards the liquid side is accounted for by geometrically shifting the planar VOF interface.
The physical model and numerical methods were implemented in the open source Basilisk solver, and the solver has been thoroughly validated to accurately simulate two-phase interfacial flows, including the aerobreakup of drops without and with phase change (Zhang, Popinet & Ling Reference Zhang, Popinet and Ling2020; Boyd & Ling Reference Boyd and Ling2023; Boyd et al. Reference Boyd, Becker and Ling2023; Ling & Mahmood Reference Ling and Mahmood2023). A key feature of the Basilisk solver is that the quadtree/octree mesh is used to discretize the domains, providing important flexibility to dynamically refine the mesh in user-defined regions. The adaptation criterion for the present study is based on the wavelet estimate of the discretization errors of the volume fraction, temperature and velocity.
3. Results
3.1. Problem set-up
We consider a freely moving drop with diameter 
$D_0$, initially stationary and at saturation temperature in an unbounded domain, suddenly exposed to a uniform superheated stream of vapour of the drop liquid with temperature 
$T_\infty$ and velocity 
$U_\infty$; see figure 1(a). Using the free-stream velocity 
$U_\infty$, the initial drop radius 
$R_0=D_0/2$, where 
$D_0$ is the initial drop diameter, the liquid density 
$\rho _l$, and the difference between the free-stream and saturation temperatures (
$T_\infty -T_{sat}$) as the characteristic scales, the dimensionless variables can be defined as
\begin{equation} u^*=u/U_\infty, \quad x^*=x/R_0, \quad \rho^*=\rho/\rho_l, \quad T^*=(T-T_{sat})/(T_\infty-T_{sat}). \end{equation}Following previous works on aerobreakup, time is non-dimensionalized by the drop breakup time (Ranger & Nicholls Reference Ranger and Nicholls1969)
\begin{equation} t^*=t U_{\infty} /(D_0 \sqrt{\rho_l/\rho_g}). \end{equation}
The definitions and values of the key dimensionless parameters, including the Weber (
${We}$), Stefan (
${St}$), Ohnesorge (
${Oh}$) and Reynolds (
${Re}$) numbers, and the liquid-to-gas density (
$\eta$) and viscosity (
$m$) ratios, are provided in table 1. The fluid properties are similar to those of acetone: the saturation temperature is 
$T_{sat}=359$ K, the surface tension is 
$\sigma =0.0153\ {\rm N}\ {\rm m}^{-1}$, the latent heat is 
$h_{l,g}=4.88\times 10^{5}\ {\rm J}\ {\rm kg}^{-1}$, and the gas density, viscosity and thermal conductivity are 
$\rho _g=5.11\ {\rm kg}\ {\rm m}^{-3}$, 
$\mu _g=9.59\times 10^{-6}\ {\rm Pa}\ {\rm s}$ and 
$k_g=0.0166\ {\rm W}\ {\rm m}^{-1}\ {\rm K}^{-1}$, respectively. In the present study, 
${Re}=1000$, 
${Pr}=0.84$ and 
$m=19.29$ are kept constant, while 
$\eta$, 
${We}$ and 
${St}$ are varied for parametric studies. We vary 
$\eta$ by changing 
$\rho _l$, and three values, i.e. 
$\eta =139$, 453 and 766, are considered, which correspond to the density ratios from acetone (
$\eta =139$) to ammonia (
$\eta =766$). By varying the free-stream temperature 
$T_{\infty }$, 
${St}$ changes from 0 to 2. The range of 
${We}$ considered is from 1 to 40, which covers the non-breakup and bag-breakup regimes, and is sufficient to identify the critical Weber number 
${We}_{cr}$. To vary 
${We}$ but keep 
${Re}$ fixed, 
$U_\infty$ and 
$D_0$ are varied simultaneously; as a result, 
${Oh}$ will vary from 
$7\times 10^{-4}$ to 
$0.01$, which is generally small and is not expected to influence the critical Weber number (Hsiang & Faeth Reference Hsiang and Faeth1995).
Figure 1. (a) Computational domain for the 3-D simulations. Morphological evolutions for (b) non-vaporizing (
${St}=0$) and (c) vaporizing (
${St}=2$) drops. The results are for 
$\eta =139$ and 
${We}=10.5$.
Table 1. Key dimensionless parameters for the present simulations.
Under the non-zero relative velocity between the drop and the vapour free stream, the drop is accelerated along the streamwise direction, while in the meantime it deforms and vaporizes. Both two-dimensional (2-D) axisymmetric and three-dimensional (3-D) simulations were performed, in which adaptive quadtree/octree meshes were used. The resolution of the mesh is controlled by the level of refinement, 
$\mathcal {L}$, which corresponds to 
$2^\mathcal {L}$ cells in a coordinate direction. For 2-D axisymmetric simulations, we have used 
$\mathcal {L}=13$, and the length of the computational domain is 
$L_0=32D_0$, which gives a minimum cell size equivalent to 
$\varDelta =D_0/512$. For 3-D simulations, we have used a coarser mesh, i.e. 
$\varDelta =D_0/256$, due to the high computational cost. Grid refinement studies were conducted by varying 
$L=11$, 12 and 13 (
$\varDelta =D_0/128$, 
$D_0/256$ and 
$D_0/512$). The results, presented in the Appendix, confirm that the mesh resolution used is sufficient to accurately resolve the vaporization of a drop undergoing significant deformation.
3.2. Vaporization-induced breakup suppression
The key observation from the simulation results is that drop vaporization tends to stabilize aerodynamic deformation and can suppress drop breakup. The strength of vaporization is characterized by 
${St}$, and the vaporization rate increases with 
${St}$. The 3-D simulation results for non-vaporizing (
${St}=0$) and vaporizing (
${St}=2$) drops are shown in figures 1(b) and 1(c), respectively, with all other parameters, including 
$\eta =139$ and 
${We}=10.5$, kept unchanged. It can be observed that the non-vaporizing drop breaks at the centre, while the vaporizing counterpart does not. The early-time shape evolutions for both drops are quite similar. The initially spherical drop expands radially into a disk shape (
$t^*=0.8$) due to the pressure difference between the windward pole and the periphery of the drop (Rimbert et al. Reference Rimbert, Castrillon Escobar, Meignen, Hadj-Achour and Gradeck2020). The internal radial flow within the liquid causes the axial thickness of the disk to reduce in time. The difference between the two cases arises at approximately 
$t^* = 1.2$, when the periphery edge of the disk starts to retract towards the centre, driven by surface tension. The radial retraction of the vaporizing drop is faster, and due to mass conservation, the thinning at the disk's centre slows down, preventing the disk thickness from diminishing to zero. As a consequence, the two surfaces of the disk do not pinch to form a hole. In contrast, the non-vaporizing drop eventually undergoes a bag-mode breakup, characterized by the formation of a central hole, whereas the vaporizing drop maintains its initial topology without undergoing topological change.
We have performed the same cases using 2-D axisymmetric simulations, which yield similar results for the present case, mainly due to the moderate 
${Re}=1000$. The less expensive 2-D simulations will then be performed for parametric simulations varying 
${We}$, 
${St}$, and 
$\eta$.
3.3. Modulation of drop surface
To confirm that the vaporization-induced suppression of bag breakup is not a coincidence for the specific density ratio 
$\eta =139$ considered, additional cases 
$\eta =453$ and 766 are simulated. The time evolutions of the drop surfaces for different 
${St}$ and 
$\eta$ are shown in figure 2(a). When 
$\eta$ increases, the dynamics of drop deformation changes, and the critical Weber number that characterizes the onset of drop breakup, 
${We}_{cr}$, decreases. The decrease in 
${We}_{cr}$ with 
$\eta$ is significant when 
$\eta \gtrsim 100$ (Marcotte & Zaleski Reference Marcotte and Zaleski2019). For the three different 
$\eta$ considered here, we have used 
${We}=12$ for 
$\eta =453$ and 766, and 
${We}=10.5$ for 
$\eta =139$, which are slightly over the corresponding 
${We}_{cr}$. If there is no vaporization (
${St}=0$), then drops of all 
$\eta$ will break in the bag mode.
Figure 2. (a) Temporal evolution of the drop surfaces for different 
$\eta$ and 
${St}$. The results for 
$\eta =139$, 453 and 766 correspond to 
${We}=10.5$, 
$12$ and 12, respectively. (b) The flow field and vaporization mass flux 
$j_\gamma ^*=j_\gamma /(\rho _l U_\infty )$ (upper half) and temperature field 
$T^*=(T-T_{sat})/(T_\infty -T_{sat})$ (lower half) at 
$t^*=1.2$ for 
$\eta =139$, 
${We}=10.5$ and 
${St}=1$.
Figure 2(a) clearly shows that the vaporization-induced breakup suppression is present for all 
$\eta$. The difference between the interface shapes for different 
${St}$ is actually larger when 
$\eta$ increases. As a result, the detailed mechanisms behind the stabilizing effect of vaporization are better revealed in the cases with higher 
$\eta =453$ and 766. After the drop deforms from a sphere to a disk, a rim is formed at the periphery; see 
$t^*=1.08$ in figure 2(a). It is observed that the rim moves downstream faster, causing the disk to bend forwards, and the level of bending increases with 
${St}$ and 
$\eta$. This modulation of the drop shape is due to the asymmetric vaporization rates on the windward and leeward sides of the rim; see figure 2(b). Since the local vaporization mass flux 
$j_\gamma$ is non-zero only in the interfacial cells (cells with fractional values of VOF), we have graphically ‘thickened’ the non-zero 
$j_\gamma$ region in the top half of figure 2(b) for better visualization, from which different vaporization rates on the two sides of the rim can be clearly seen. The vaporization rate depends on the temperature gradient near the interface, which is in turn inversely proportional to the thickness of the thermal boundary layer near the interface, given that 
$T_{\infty }-T_{sat}$ is fixed. Due to the stagnation flow on the windward side, the boundary layer thickness is much smaller than that on the leeward side where the flow separates; see figure 2(b). As a result, the vaporization is significantly stronger on the windward side of the rim. Due to the density difference between the liquid and the vapour, Stefan flow is induced at the interface, which repels the interface in the opposite normal direction. The higher vaporization rate on the windward side of the rim leads to a stronger repelling force, which causes the disk to bend forwards.
The excessive forward bending of the disk hinders bag formation and development. As shown in previous studies on bag breakup (Ling & Mahmood Reference Ling and Mahmood2023; Tang et al. Reference Tang, Adcock and Mostert2023), the bag forms when the centre of the disk becomes sufficiently thin, as seen at 
$t^*=1.44$ for 
$\eta =766$ and 
${St=1}$ in figure 2(a). The centre of the disk moves downstream faster than the periphery, forming a bag with an upstream-facing opening, as observed at 
$t^*=1.80$ for 
$\eta =766$ and 
${St}=1$. Vaporization, however, modulates the drop shape in the opposite manner, causing the edge to move downstream faster than the centre, forming a ‘reverse bag’ with a downstream-facing opening, as seen at 
$t^*=1.08$ for 
$\eta =766$ and 
${St}=2$. This modulation hinders bag development. Near 
${We}_{cr}$, the force driving bag growth (the pressure difference across the drop) is only slightly higher than the resisting force (the surface tension on the periphery). This unfavourable modulation by vaporization is sufficient to suppress bag development and prevent hole formation in the bag.
As the Stefan flow is related to the density difference between vapour and liquid, Stefan flow and the resultant repelling force are enhanced when 
$\eta$ increases. For the highest 
$\eta =766$, vaporization results in a more significant modulation of the shape and drag of the drop. As a result, even a smaller 
${St}=1$ is sufficient to suppress drop breakup, in contrast to the cases for 
$\eta =139$ where suppression is not observed until 
${St}=2$.
3.4. Modulation on gas flow and drag
Another important effect of vaporization is the enhancement of drag and drop acceleration, which can be observed qualitatively in figure 2(a). Detailed quantitative measurements of the drag coefficient (
$C_d$) are given in figure 3(d). Here, the drag coefficient is defined based on the instantaneous drop frontal area,
\begin{equation} C_d=\frac{F_d}{\dfrac{1}{2}\rho_g (U_\infty-u_d)^2 {\rm \pi}R^2} = \frac{F_d^*}{\dfrac{1}{2\eta}(1-u_d^*)^2 {\rm \pi}(R^*)^2}, \end{equation}
where 
$F_d^*=F_d/(\rho _l U_\infty ^2 R_0^2)$. While 
$C_d$ varies over time, it is observed that 
$C_d$ generally increases with 
${St}$; however, the detailed mechanisms behind this trend remain to be investigated.
Figure 3. The gas (a) vorticity (
$\omega ^*_g=\omega _g R_0/U_\infty$) and (b) pressure (
$p^*_g=p_g /(\rho _g U_\infty ^2)$) fields at 
$t^*=1.2$ for different density ratios 
$\eta =139, 453, {766}$, and Stefan numbers 
${St=0}$ and 2. The temporal evolutions of the drop lateral radius 
$R^*$ and the drag coefficient 
$C_d$ for different 
$\eta$ and 
${St}$ are shown in (c,d), respectively. The dashed lines in (a,b) indicate the lateral radius of the wake (
$r_w^*$) estimated based on vorticity.
Since vaporization on the two sides of the rim is asymmetric, the Stefan flow and the resulting repelling force are stronger on the upstream side of the rim, leading to a net force in the streamwise direction, which contributes to an increase in drag in the rim region. However, this contribution to the overall drag is expected to be minor, as vaporization on the drop surfaces away from the edge rim is approximately symmetric; see figure 2(b). Another possible reason for the higher drag coefficient is the larger lateral drop radius when vaporization is present, as seen in the results for 
$\eta =453$ in figure 3(c). While the drag on the drop will increase with 
$R^*$ and the frontal area, the drag coefficient 
$C_d$ (see (3.3)) will not. Interestingly, for 
$\eta =139$, 
$R^*$ actually decreases slightly with 
${St}$ (see 
$t^*=1$), while 
$C_d$ remains higher for larger 
${St}$. Therefore, the increase in 
$C_d$ when 
${St}$ increases, as observed in figure 3(d), cannot be explained by the modulation of 
$R^*$ alone.
The key mechanism for drag enhancement is the modulation of the wake structure by the Stefan flow. It can be clearly observed from the vorticity fields (figure 3a) that the maximum lateral radius of the wake, 
$r_w^*$, increases when vaporization is present. Near the flow separation point, the Stefan flow pushes the vorticity layer away from the interface, resulting in a larger slope of the wake boundary and eventually in a larger 
$r_w^*$. The pressure in the wake is generally low, as shown in figure 3(b), thus the enlargement of the wake radius contributes to the increase of the low-pressure region area and the net drag acting on the drop. This mechanism is similar to the drag crisis for a solid sphere, in which the drag is reduced when the wake lateral radius decreases due to the delay of the turbulent boundary layer separation (Tiwari et al. Reference Tiwari, Pal, Bale, Minocha, Patwardhan, Nandakumar and Joshi2020). Since the strength of the Stefan flow is enhanced as 
$\eta$ increases, its modulation of the wake becomes more significant for higher 
$\eta$, as shown in figures 3(a) and 3(b). As depicted in figure 3(a), when 
${St}$ increases from 0 to 2, 
$r_w^*$ for 
$\eta =139$ increases from 2.4 to 2.5, while for 
$\eta =453$ and 766, 
$r_w^*$ increases from 2.7 to 3.0 and from 2.8 to 3.3, respectively. Correspondingly, the vaporization-induced enlargement of the low-pressure region for 
$\eta =766$ is more pronounced than for 
$\eta =139$; see figure 3(b). Consequently, the increase in 
$C_d$ with 
${St}$ is more significant for higher 
$\eta$. When 
${St}$ increases from 0 to 2, 
$C_d$ at 
$t^*=1$ for 
$\eta =766$ increases by approximately 67 %, whereas the increases for 
$\eta =453$ and 139 are approximately 40 % and 5 %, respectively; see figure 3(d).
3.5. Effects of Stefan flow and interface receding
The drop vaporization results in two key modulations in interfacial dynamics: (1) the interface receding towards the liquid, reducing the drop volume; and (2) the Stefan flow induced by the fluid volume generated near the interface that occurs when a high-density liquid turns into a low-density vapour. While the previous discussions are focused on the effect of the Stefan flow, it remains to affirm that it indeed dominates the other contribution in the breakup suppression and drag enhancement. For this purpose, additional numerical ‘experiments’ are performed to investigate the individual effects on these two features.
Four numerical tests are performed, and the results are shown in figure 4: (1) 
${St=0}$, where there is no vaporization; (2) 
${St=2}$-IR (interface receding), where the interface recedes but Stefan flow is turned off; (3) 
${St=2}$-SF (Stefan flow), where Stefan flow is produced but the interface receding is turned off, so the drop volume remains unchanged; and finally, (4) 
${St=2}$, where the full vaporization effect is incorporated. In the simulation for 
${St=2}$-IR, the source term on the right-hand side of (2.2) is set to zero, while the source term on the right-hand side of (2.3) remains unchanged. As a result, the VOF interface recedes towards the liquid side due to vaporization, but Stefan flow is not produced. For 
${St=2}$-SF, the opposite is done: the source term on the right-hand side of (2.3) is set to zero, while the counterpart for (2.2) remains unchanged. Consequently, Stefan flow is induced, but the interface does not recede due to vaporization, and the drop volume remains constant. It should be noted that these two cases are ‘unphysical’ and can be examined in a numerical study only by switching vaporization-related source terms on and off in different equations. However, these cases are very helpful for understanding the individual effects of vaporization and identifying the dominant one.
Figure 4. Temporal evolutions of (a) drop surfaces, (b) drop lateral radius 
$R^*$, and (c) drag coefficient 
$C_d$ for different numerical tests: 
${St}=0$ represents cases without vaporization; 
${St=2}$-IR represents the numerical tests where the interface recedes due to vaporization but the source term on the right-hand side of (2.2) is turned off to avoid the production of the Stefan flow; 
${St=2}$-SF represents the numerical tests where the Stefan flow is generated due to vaporization, but the source term on the right-hand side of (2.3) is set to zero to avoid interface receding towards the liquid phase due to vaporization; 
${St}=2$ represents a full simulation with both interface-receding and Stefan flow effects accounted for.
For both 
$\eta$, the drop surfaces for 
${St=2}$-IR are very similar to those for 
${St}=0$ throughout the drop evolution; see figure 4(a). Holes are formed for both cases, leading to a change in the drop topology. This indicates that vaporization-induced recession of the interface and resulting volume change have minimal impact on drop deformation. In contrast, the drop surfaces for the case 
${St=2}$-SF, with Stefan flow turned on but with the interface receding turned off, are close to the full simulation results for 
${St}=2$ for all times, and eventually the breakup is suppressed for both cases. The results for the drop lateral radius shown in figure 4(b) are consistent, and the evolutions of 
$R^*$ for 
${St}=0$ and 
${St=2}$-IR are similar, while the results for 
${St=2}$-SF agree well with those for 
${St}=2$. Therefore, we can conclude that the Stefan flow is the primary effect of vaporization that influences the drop deformation.
The drag coefficient for the numerical tests is presented in figure 4(c). The good agreement between the results for 
${St=2}$-SF and 
${St}=2$ confirms that the Stefan flow is responsible for the drag enhancement. Without incorporating the Stefan flow, as in the test 
${St=2}$-IR, one will significantly underestimate the drag and drop acceleration.
3.6. Impact on breakup criteria
The drop breakup criteria are often defined based on the critical Weber number 
${We}_{cr}$. When 
${We}$ is lower than 
${We}_{cr}$, the drop will remain unbroken throughout the acceleration. Parametric 2-D axisymmetric simulations have been performed to identify 
${We}_{cr}$ for different 
$\eta$ and 
${St}$. We have considered three values of 
$\eta$, and five values of 
${St}$, so there are 15 cases with different 
$\eta$–
${St}$ combinations. For each case, simulations of different 
${We}$ are run to find 
${We}_{cr}$ iteratively, using a bisection method with minimum and maximum limits 
${We}=1$ and 40, respectively. The tolerance for iteration convergence is 0.2, that is, the root finding error for 
${We}_{cr}$ is less than 0.2. For each case, approximately 5–7 simulations are needed to reach the desired tolerance. For all 15 cases, more than 100 runs were performed.
For some cases just above 
${We}_{cr}$, the drop lateral radius 
$R^*$ is decreasing when the drop breaks at the centre. Those drops will eventually return to the ellipsoidal shape and will not atomize into a large number of child droplets. Such a near-
${We}_{cr}$ behaviour has been identified in previous studies as well (Ling & Mahmood Reference Ling and Mahmood2023). We considered the formation of a hole in the drop bag to also be a breakup case to give a robust criterion for determining 
${We}_{cr}$.
The results of 
${We}_{cr}$ for different 
$\eta$ and 
${St}$ are summarized in figure 5. For non-vaporizing drops (
${St}=0$), 
${We}_{cr}=9.8$, 10.9 and 11.4 for liquid-to-gas density ratios 
$\eta =139$, 
$453$ and 766, respectively. Experimental measurements for 
${We}_{cr}$ available in the literature are mainly for non-vaporizing drops with large 
$\eta$, and 
${We}_{cr}\approx 11$ for water drops in air (
$\eta =831$) (Guildenbecher et al. Reference Guildenbecher, López-Rivera and Sojka2009). Therefore, our results for 
${We}_{cr}$ agree well with the experimental data in the literature. Furthermore, the present results show that 
${We}_{cr}$ decreases slightly with 
$\eta$, which is also consistent with previous studies on the density ratio effect (Jain et al. Reference Jain, Tyagi, Prakash, Ravikrishna and Tomar2019; Marcotte & Zaleski Reference Marcotte and Zaleski2019). To the best knowledge of the authors, there are no experimental measurements for 
${We}_{cr}$ for vaporizing drops. Since vaporization tends to stabilize the drop, 
${We}_{cr}$ generally increases with 
${St}$, meaning that a higher gas dynamic pressure is required to break the drop. It was found that 
${We}_{cr}$ increases by approximately 13 % for 
$\eta =766$ when 
${St}$ increases from 0 to 2.
Figure 5. The critical Weber number 
${We}_{cr}$ for the aerodynamic breakup of a vaporizing drop as a function of the Stefan number 
${St}$, for different liquid-to-vapour density ratios (
$\eta =139$, 453 and 766). The corresponding drop surfaces upon breakup are also shown.
The drop shapes upon breakup at the corresponding 
${We}_{cr}$ are also shown in figure 5, where it can be observed that vaporization modifies the drop shape at the point of breakup. As expected, the breakup outcomes, including the velocity and size distributions of the resulting child droplets, will also change. However, a detailed study of the child droplet statistics is beyond the scope of the present work.
Although the general increasing trend of 
${We}_{cr}$ with 
${St}$ holds for all density ratios, a closer examination reveals that the variations in 
${We}_{cr}$ depend on 
$\eta$. For 
$\eta =139$, 
${We}_{cr}$ increases approximately linearly as 
${St}$ increases, and the drop shapes remain similar. In contrast, the variation of 
${We}_{cr}$ for higher 
$\eta$ is more complex. Specifically, for 
$\eta =453$, 
${We}_{cr}$ shows little change within the searching tolerance 0.2 for 
${St}$ values between 0.5 and 1.5. In these cases, vaporization delays the breakup and modifies the drop shape at the point of breakup, but it is not sufficient to suppress the breakup. For 
$\eta =453$, the breakup time for 
${St}=1.0$ is 
$t^*=1.907$, which increases to 1.925 for 
${St}=1.5$, while 
${We}_{cr}$ for these two cases remains almost the same.
For 
$\eta =766$, although 
${We}_{cr}$ increases monotonically with 
${St}$, similar to the trend observed for 
$\eta =139$, the drop shape upon breakup varies more significantly with changes in 
${St}$. While the shape for the non-vaporizing case (
${St}=0$) at 
$\eta =766$ closely resembles that of its lower-
$\eta$ counterparts, the drop experiencing strong vaporization, such as at 
${St}=2$, displays a very different breakup shape, deforming into a long bag along the streamwise direction. For low 
${St}$ values, the drop at 
$\eta =766$ breaks in a typical bag mode by pinching the two interfaces and forming a hole near the axis, similar to what occurs for lower 
$\eta$ values, which is expected when 
${We} \approx {We}_{cr}$. However, as 
${St}$ increases, the pinching location shifts away from the axis, and the breakup morphology appears to shift towards a bag-stem mode, though a more detailed investigation would be required to confirm this observation.
Finally, it is noted that the variations of 
${We}_{cr}$ discussed above are strictly valid for the current range of parameters considered, namely 
$139\leqslant \eta \leqslant 766$ and 
$0\leqslant {St} \leqslant 2$. Extension of the present study to higher 
${St}$ and ambient temperature will be relegated to our future work.
4. Conclusions
Parametric 3-D and 2-D axisymmetric interface-resolved simulations are conducted in the present study to investigate the aerobreakup of vaporizing drops. The Weber and Stefan numbers, and the liquid-to-gas density ratio, are varied to systematically study the effect of vaporization on the drop deformation/breakup and drag. For a wide range of liquid-to-density ratios, vaporization contributes to hindering the drop deformation, and if the Stefan number is sufficiently high, can suppress bag breakup. The vaporization also contributes to enhancement of the drag: it is shown that the drag coefficient for a vaporizing drop can be 84 % higher than its non-vaporizing counterpart. The Stefan flow induced by the interface vaporization plays the dominant role in modulations of the drop deformation stabilization and drag enhancement. The asymmetric vaporization rates on the windward and leeward sides of the edge rim and the resulting different strength of Stefan flow cause the disk to deform in an opposite manner compared to the bag, which hinders the bag development and breakup. The Stefan flow modifies the boundary layer dynamics near the separation point, leading to a wider wake and a large drag. More than 100 simulations were run to identify the critical Weber number for different Stefan numbers and density ratios. When the Stefan number increases from 0 to 2, the critical Weber number can increase up to 13 %.
Funding
This research was supported by ACS-PRF (no. 62481-ND9) and NSF (no. 1942324). We also acknowledge the ACCESS program for providing the computational resources that have contributed to the research results reported in this paper. The code used for the present simulations is an extension of the open-source multiphase flow solver Basilisk, which is made available by S. Popinet and other collaborators.
Declaration of interests
The authors report no conflict of interest.
Appendix. Grid refinement studies
Grid refinement studies for non-vaporizing drops can be found in our previous paper (Ling & Mahmood Reference Ling and Mahmood2023). The results of grid refinement studies for a vaporizing drop (
$\eta =139$, 
${We}=10.5$ and 
${St}=2$) are shown in figure 6. This case represents a drop with strong vaporization and high 
${St}$ in the present study. Three different refinement levels, 
$L=11$, 12 and 13, corresponding to 128, 256 and 512 cells per initial drop diameter, are used. The time evolution of the drop surface for different meshes is shown in figure 6(a), where it can be observed that the results for different meshes are almost indistinguishable until 
$t^*=1.08$. At later times, the results for 
$L=12$ and 13 remain very close. The distribution of vaporization mass flux 
$j_\gamma ^*$ on the drop surface at 
$t^*=0.72$ and 1.08 is shown in figures 6(b) and 6(c). Here, 
$j_\gamma ^*$ is plotted as a function of 
$\cos \theta$, where 
$\theta$ is the colatitude on the drop surface with respect to the origin, which is in the middle between the windward and leeward poles (see 
$t^*=0.72$ in figure 6a). Again, the results for 
$L=12$ and 13 agree very well, confirming that the mesh resolution (
$L=13$ or 
$D_0/\varDelta =512$) used in the present study is sufficient to resolve the interfacial dynamics when vaporization is present.
Figure 6. (a) Temporal evolution of the drop surfaces for different mesh refinement levels 
$L=11$, 12 and 13, for the case 
$\eta =139$, 
${We}=10.5$ and 
${St}=2$. (b,c) Distribution of the vaporization mass flux 
$j_\gamma$ on the drop surface for 
$t^*=0.72$ and 1.08, respectively.
