Thermally induced secondary atomization (TISA) is a complex phenomenon that accelerates phase change in the combustion chamber. It occurs if multi-component fuels, having a wide boiling range, are exposed to high temperatures. Several airlines are recently experimenting with bio- and fossil fuels blends. However, the characteristics of droplet TISA are primarily unknown because of the challenges associated with experimental activities like suspended or falling droplets. In this scenario, numerical models become essential to study TISA. That is why a new multi-component, multi-phase volume of fluid computational fluid dynamics solver was developed to simulate droplets TISA. The solver takes advantage of the OpenFOAM framework and uses the *isoAdvector* methodology. The bio- and fossil fuels were represented by n-heptane and n-hexadecane, respectively, to simplify the problem. Evaporation was implemented by assuming that the mixture could only boil at that temperature. Surface tension and other relevant mixture properties were considered as a function of species concentration and temperature to replicate all phenomena comprehensively. An analysis of bubble expansion based on the Rayleigh–Plesset equation preceded the breakup tests. The test cases consisted of a droplet suspended in microgravity having a bubble initialized at the interface. The bubble eventually expanded, and the bubble cap collapsed, leading to the micro-explosion. A parametric study of breakup cases under different pressures and at a fixed temperature of 1200 K was performed. The atomization mechanism was tested at 1, 3, 10, and 20 bar and compared. It was observed that while high pressure slows down the process, it finally leads to a higher surface area. This behavior was confirmed by testing two different bubble sizes. Together with the atomization intensity, also the morphology of the particles changed. At atmospheric pressure, the maximum surface area was reached when the droplet had a disk-like shape, while at higher pressures, it evolved in an elongated shape.

## I. INTRODUCTION

The continuous reliance on fossil fuel reservoirs and rising carbon emissions have motivated the adoption of renewable fuels. This trend is reflected in the aviation sector, where companies, such as Air France/KLM, Lufthansa, Japan Airlines, and Air New Zealand, are already testing commercial passenger flights using blends of conventional fuel and bio-derived sustainable aviation fuels (SAF).^{1} Nevertheless, exploiting renewable fuels in the aviation sector is challenging since the new and current powers must be fully compatible with maintaining existing engines and minimizing changes to logistics and infrastructures. Together with compatibility issues, the new fuel must satisfy safety and quality standards and be characterized by low carbon emissions throughout its lifecycle.^{2} The International Air Transport Association (IATA) has approved three classes of fuels: Fischer–Tropsch Hydroprocessed Synthesized Paraffinic Kerosene (FT-SPK), produced from any carbonaceous raw material, the Synthesized Paraffinic Kerosene from Hydroprocessed Esters and Fatty Acids (HEFA-SPK), created from animal or vegetable oil, and the Synthesized Iso-Paraffins from Hydroprocessed Fermented Sugars (SIP-SPK), from sugar through fermentation. The HEFA can be blended up to 50% volume with conventional fuel oil.^{3}

A crucial step in the combustion of aviation fuel is the atomization process. HEFAs, for instance, are more viscous than conventional jet fuel; consequently, they produce larger droplets when atomized, lowering the surface area and decreasing combustion efficiency. On the other hand, studies from Pacheco *et al.*^{4} showed that droplets of 25% hydroprocessed vegetable oil (HVO), with jet A-1 fuel, presented micro-explosions, which significantly increased the surface area of the droplet, decreasing the time required for its consumption.

Micro-explosions occur when a multicomponent droplet is exposed to high temperatures; the more volatile component vaporizes, forming a bubble that later causes breakup phenomena.^{5} The thermally induced secondary atomization (TISA) process depends mainly on the Lewis number and the boiling range. The Lewis number, in particular (ratio between thermal conductivity and mass diffusivity), results in larger than one when heat penetrates the droplet at a faster rate than mass diffusion to the external environment.

Dryer^{6} and others extensively studied micro-explosions for water-in-fuel emulsions. This study adopted the suspended droplet experiment; however, since quartz suspension fiber may lower the limit of superheat in the mixture,^{7} the experiment was replicated in the form of falling droplets by Lasheras *et al.*, and micro-explosive behavior was observed again. In later experiments, Lasheras *et al.* studied the combustion of paraffinic mixtures. They observed that the combinations yielded secondary atomization when the components had different normal boiling points.^{8}

^{9}The effect of pressure on this phenomenon was studied by Wang

*et al.*, who showed the occurrence of micro-explosions for both emulsions of water and n-alkanes and solutions of n-alkanes and alcohols, concluding that an increase in pressure increased the probability of micro-explosions.

^{10}Guida

*et al.*

^{11}demonstrated how substantial the deviation from the “

*d*

^{2}” law could be when micro-explosions and puffing occur. Many studies have conducted secondary atomization induced by air shear and fluid-dynamics instabilities of the liquid jet.

^{12,13}The above-mentioned phenomena were governed mainly by the Weber number. At the same time, a thermally induced breakup is a function of the thermodynamics of the system and the Ohnesorge number (

*Oh*), which relates viscous forces and surface tension,

*μ*is the dynamic viscosity of the fluid,

*ρ*is the density,

*σ*is the surface tension, and

*R*is the radius of the bubble. Exploring a broader pressure range is germane because of the operative conditions of gas turbines. For this reason, the thermally induced secondary atomization of fuel droplets at different pressures was studied here with particular focus on the evolution of the surface area of the liquid phase. The latter is a relevant parameter because it dictates the area exposed to the gas phase involved in the vaporization process and combustion. While most studies focused on the description of the conditions necessary for the atomization phenomenon to take place, the approach detailed in this work focuses on the breakup itself in order to attempt the estimation of the additional surface area obtained after the collapse of the parent droplet. The numerical study proposed in this work was performed using an in-house computational fluid dynamics (CFD) geometric volume-of-fluid-based (VoF-based) solver. The solver was implemented in OpenFOAM and validated with analytical and experimental data in the previous work.

^{11,14}Interface reconstruction was performed geometrically using the

*isoAdvector*concept.

^{15,16}Energy and species equations were solved separately for the two phases to estimate the phase change terms better.

^{17}The coupling of the energy jump condition and thermodynamic pseudo-equilibrium at the liquid–vapor interface was not implicitly resolved for simplicity. The model is briefly explained, the case setup is introduced, and assumptions and simplifications are justified in Secs. II–V. Based on the Rayleigh–Plesset equation, a simple calculation demonstrates how the expansion is significantly faster than the breakup. The results are discussed, focusing on the peculiar morphology of atomizing droplets and how pressure affects their evolution. Conclusions and some guidelines for future development are presented at the end of this work.

## II. ALGORITHM

^{18}The dynamics of two-phase compressible fluids consisting of liquid and gas were described by a unified velocity $ u ( x , t )$, pressure $ p ( x , t )$, temperature $ T ( x , t )$, and the mass fraction of the

*j*-th species $ Y j ( x , t ) , j = 1 , \u2026 , N S$, where

*NS*is the total number of species in the system. Letters in bold represent vectorial quantities, and non-bold letters represent scalar quantities. In the following chapter, the computational domain is defined as $ \Omega \u2208 \mathbb{R} 3$, while the two subsets $ \Omega 1$ and $ \Omega 2$ represent the liquid and gas phase, respectively. The intersection of the two subsets is the interface $ \Gamma \u2208 \mathbb{R} 2$. The interface Γ is the boundary between two adjacent phases, $ \Omega 1$ and $ \Omega 2$. The transport equation for the volume fraction of liquid

*α*of a compressible fluid experiencing phase change is expressed as

**g**represents the gravity force,

*μ*is the dynamic viscosity, and $ f s$ is the surface tension force. Surface tension force, $ f s$, is important and generally modeled using the CSF method introduced by Ref. 19,

*σ*is the surface tension coefficient,

**n**is the normal to the interface, and

*q*is the curvature expressed as the divergence of the normalized gradient of the volume fraction,

*j*represents the

*j*-th species, and

*D*is the mass diffusivity. Finally, the energy equation was also solved to reconstruct the temperature field across the computational domain,

_{j}*k*is the thermal conductivity,

*c*is the heat capacity, and $ \Delta H v , j$ is the heat of vaporization of the

_{p}*j*-th species. The parameter NS is the total number of species. Scalar variables (temperature and species mass fraction) are transported using a two-fluid approach. The liquid and gas phases are solved independently and coupled through a moving boundary condition at the interface. Discretized transport equations were solved using the PIMPLE algorithm, a hybrid methodology combining the pressure of the iterative

*pressure-implicit with the splitting of operators*(PISO) and the

*semi-implicit method for pressure linked equations*(SIMPLE). Momentum and pressure equations were solved sequentially until convergence was achieved. The pressure equation was solved using the

*geometric agglomerated algebraic multigrid*(GAMG) method. Time derivatives were discretized with an implicit Euler scheme; convective terms in the volume fraction equation and the momentum equation were discretized with a van Leer scheme, while a Gauss linear scheme was applied to the remainder. Simulations were performed on the supercomputer Shaheen II (Cray XC40), operating one node (32 cores) for each simulation.

## III. BUBBLE NUCLEATION AND EXPANSION

*Le*), defined as

*h*is the thermal diffusivity, and

*D*is the mass diffusivity. The fuel blend considered had a Lewis number notably higher than one. The mass diffusivity of heptane in hexadecane was, in fact, proportional to $ 1 \xd7 10 \u2212 9 \u2009 m 2 s \u2212 1$, while the thermal conductivity was proportional to $ 1 \xd7 10 \u2212 7 \u2009 m 2 s \u2212 1$. The boiling point of the two pure components was also different, because n-hexadecane boils at 560 K and n-heptane at 371 K. This meant that nucleation would likely occur when the mixture was exposed to high temperatures. However, nucleation is a non-deterministic process, and its description is somewhat complex. Works by other authors dealt with the complexity of the nucleation process.

^{22–24}The formation of a critical nuclei proceeds the enlargement of a bubble. The expansion of a spherical bubble is described by the Rayleigh–Plesset equation,

*p*and $ p \u221e$ are the pressure inside the bubble and in the liquid, respectively. The physical properties refer to the liquid phase. The equation can be solved to obtain the evolution of the bubble radius in time,

_{B}*R*(

*t*). According to Mikic

*et al.,*

^{25}two limiting cases can be identified: an inertially controlled growth and a diffusion-controlled growth. In the first scenario, the solution is $ R ( t ) \u223c A t$, while in the second one, $ R ( t ) \u223c B t 1 / 2$.

*A*and

*B*are two proportionality parameters. The authors

^{25}identified the two as

*l*refers to the liquid phase,

*v*to the gas phase, and

*sat*to the saturation properties. The two combined describe the bubble expansion in both regimes,

Figure 1 shows the evolution of the bubble radius in time [obtained using Eq. (18)] for different values of the liquid pressure. The solution showed that the timescale for bubble growth was shorter than the one observed for the micro-explosion, which was of the order of a few milliseconds.^{8,26} Consequently, it was possible to consider the bubble growth as instantaneous, thus focusing only on the subsequent breakup.

## IV. CASE STUDY

*μ*m from the interface. Because of the generally small influence of the gravity on droplet dynamics, the assumption of micro-gravity conditions was chosen. The influence of gravity is related to the Eötvös number (

*Eo*), defined as

*μ*m were also simulated to provide a comparison and depict a scenario closer to actual applications. External temperature was set to 1200 K to simplify the problem although larger gradients are required to better model real-world applications. The pressure changed depending on the case. Ambient pressures of 1, 3, 10, and 20 bar were analyzed, and two bubble sizes of 0.6 and 0.8 mm were also considered for the case with the droplet size of 1 mm, while the same bubble to droplet ratio was also set for the 200

*μ*m test. A reasonable surrogate for the mixture of this study combines n-heptane and n-hexadecane, the first simulating the fossil fuel and the second a heavier, bio-derived component. The combination was chosen since the dominant phenomena involved in thermally induced secondary atomization are controlled by the thermodynamic properties of the two components, which reasonably well represent the mixtures of interest. The composition of the particle was 50:50 wt. % of n-heptane and n-hexadecane. Those two species were sufficiently different in boiling temperature to experience thermally induced atomization, particularly under pressure. The internal bubble was assumed to be created solely by the lighter component, while the liquid was a mixture of n-heptane and n-hexadecane. The external environment was pure air, although no reactivity was considered. The initial droplet temperature was 400 K, so heat was derived from the external to the internal environment. Bubble pressure was calculated from the Young–Laplace equilibrium equation at the interface,

*i*-th component. The mixing rule for the surface tension was

*x*and

_{i}*x*are the molar fraction of the

_{j}*i*-th and

*j*-th, respectively, and

*ρ*is the average mixture density. Thermal conductivity of the liquid mixture was calculated with the correlation proposed by Ref. 27:

*V*is the molar volume of the

_{i}*i*-th component. The parameter $ k i , j$ is calculated as

*μ*is the viscosity of the

_{i}*i*-th component, and

*x*is its molar fraction. If not specified, a simple molar or mass average was used.

_{i}The equation of state used in the work was the Soave Redlich Kwong.^{29} The computational domain was a cubic box of a size 3R×3R×3R, where R is the radius of the droplet. The mesh was uniform and orthogonal and composed of a total of $ 1 \xd7 10 6$ cells. The boundaries are open; therefore, a Neumann boundary condition is imposed by assuming the derivative of the variables with respect to the space equal to zero across the interface. Mesh convergence was demonstrated in a previous work from the group.^{14} The time step was set at 1e-6 s.

## V. RESULTS

The atomization process progressed in three steps: pinching, crown formation, and ejection. The first step occurred when the liquid meniscus above the droplet became thinner, rupturing into a hole. The two leading causes of pinching were expansion of the bubbles and boiling at the surface caused by the hot external environment. Another reason for the bubble motion was the Marangoni effect, originated by surface tension gradients in the liquid. A recent work demonstrated that in most cases, boiling is not a major cause of pinching although it contributes to its occurrence.^{14} When pinching occurs, instabilities propagate along the crater crest; those instabilities form crown-like structures and eventually evolve into liquid filaments. The liquid filaments anticipate the ejection of droplets that occur through amplification of Plateau–Rayleigh instabilities. Micro- and nano-droplets formation occurs because of Rayleigh–Taylor instability in the liquid filaments. Figure 2 shows the filament formation in an azimuthal view of droplets at atmospheric pressure. The crown formation generates the droplets; however, the formation of those secondary droplets was not the main cause of surface area increase; instead, the surface area increased mainly because of the shape assumed by the parent droplet during the process. Secondary droplets spread at velocities in the order of $ 0.01 \u2009 ms \u2212 1$. The secondary particles rapidly changed phase as their temperature rose while projected to the hot external atmosphere. The temperature of secondary droplets quickly reached the boiling point of the mixture. Figure 3 shows the liquid temperature of the droplet, highlighting the much higher temperature reached by the small structures, which basically boiled immediately as they left the parent droplet. Conversely, the temperature of the parent droplet increased at a lower rate. The crown formed from instabilities propagating at its crest and enlarged following the ejection of the higher pressure internal vapor until it folded on itself. In the case in Fig. 4, the particle became almost a disk, while ejection took place. Figure 4 demonstrates the crown folding after 0.002 s, when the crater reached its maximum diameter and the ejection was initiated. Because of computational power requirements, the ejection step was not completely simulated.

### A. Pressure effect

Pressure was observed to have a major effect on the evolution of the droplet after ejection. Pinching occurred later in time, as also noted by Law *et al.*^{5} and reported in their work. Slower pinching could be attributed to a lower relative gradient between the internal vapor and the external environment, which eventually reduced the stress on the bubble cap. Another important effect of high pressure was the limited expansion of the droplet, which maintained a more rounded but elongated shape throughout the breakup.

Figures 5 and 6 report the evolution in time of an iso-contour of the volume fraction of liquid for two different experiments at 10 and 20 bar.

_{max}on pressure. The maximum normalized droplet size was defined as

_{0}is the initial droplet size. The numerical experiments were performed for two bubbles to droplet volume ratio,

Breakup time was defined as the time between pinching and the occurrence of the maximum surface area; in both cases, this time gradually increased with pressure, as expected. In fact, pressure acted as a binding force, first, compressing the liquid and later slowing the ejection of gas from the internal bubble to the external environment. Larger bubbles also slowed the atomization process at low pressures, while they performed very much the same when the pressure was higher than 10 bar. The atomization process of a large bubble was intrinsically slower as multiple pinching zones and ejection occurred simultaneously.

Finally, a last numerical experiment was performed to verify the influence of droplet size. A case, where the droplet had a diameter of 200 *μ*m and the bubble a diameter of 120 *μ*m, was tested at 10 bar and at the usual temperature of 1200 K (Fig. 9). The outcome was a faster atomization compared to the larger particle with the same surface area ratio between bubble and droplet. This result was expected as the pressure inside a smaller bubble is higher compared to other cases. Also, the temperature of the liquid increases faster given the shorter characteristic length heat has to travel.

## VI. CONCLUSIONS

Numerical simulations were used to demonstrate and quantify the effect of pressure on thermally induced secondary atomization. A binary mixture of n-heptane and n-hexadecane was considered as a surrogate. It was assumed that the significant difference in the boiling temperature would trigger bubble formation. The evolution of the critical nuclei was then studied by solving the Rayleigh–Plesset equation, but at all pressures considered, the characteristic time of bubble expansion was substantially shorter than the time required for breakup.

As pressure increased, the onset of breakup occurred later, and the entire process was substantially slower. The onset of breakup, in particular, required a time higher than one order of magnitude in the cases of breakup at 20 bar, compared to the simulation at 1 bar. Compared to an almost instantaneous breakup (order 0.01 s) in the case at atmospheric pressure, in cases at 10 and 20 bar, the droplets reached their maximum surface area after almost 0.7 s. It was also observed that the maximum surface area reached under pressure was substantially higher in the case with a smaller (0.6 mm) bubble initialized, reaching almost two times the initial surface. The effect of pressure on the extent of the atomization process was greatly reduced when the bubble was larger. Breakup time with the larger bubble was slower at a low pressure and slightly faster at higher pressure. The latter behavior could have been caused by the lower pressure excess of the inner bubble, given the larger size. To conclude, it was demonstrated that pressure has a potentially beneficial effect in the atomization of fuel blends. However, further analysis is required to generalize the observations made in this work, as, for instance, the size of the analyzed droplets was larger than those actually generated in combustion devices such as aircraft turbines.

## ACKNOWLEDGMENTS

This work was sponsored by the Clean Combustion Research Center at King Abdullah University of Science and Technology (KAUST). Computational resources were provided by the KAUST Supercomputing Laboratory (KSL).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Paolo Guida:** Conceptualization (lead); Software (lead); Validation (lead); Writing – original draft (equal). **Alberto Ceschin:** Conceptualization (supporting); Methodology (equal); Software (supporting); Validation (equal); Visualization (lead). **Chiara Canciani:** Conceptualization (equal); Data curation (supporting); Investigation (equal); Visualization (supporting); Writing – original draft (equal). **Hong G. Im:** Conceptualization (supporting); Supervision (supporting); Writing – review & editing (supporting). **William L. Roberts:** Conceptualization (supporting); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

_{2}emissions evolution towards 2050

^{®}isoAdvector solvers using single bubble benchmarks

*Classical Nucleation Theory in Multicomponent Systems*

*n*-paraffin solutions and emulsions