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

(3.1a–d)\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)

(3.2)\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,

(3.3)\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.