Puﬀing-enhanced fuel/air mixing of an evaporating n -decane/ethanol emulsion droplet and a droplet group under convective heating

Pufﬁng of a decane/ethanol emulsion droplet and a droplet group under convective heating and its effects on fuel/air mixing are investigated by direct numerical simulation that resolves all of the liquid/gas and liquid/liquid interfaces. With distinct differences in the boiling point between decane and ethanol, the embedded ethanol sub-droplets can be superheated and boil explosively. Pufﬁng, i.e. ejection of ethanol vapour, occurs from inside the parent decane droplet, causing secondary breakup of the droplet. The ejected ethanol vapour mixes with the outer gas mixture composed of air and vapour of the primary fuel decane, and its effects on fuel/air mixing can be characterised by the scalar dissipation rates (SDRs). For the primary-fuel SDR, the cross-scalar diffusion due to ethanol vapour pufﬁng plays a dominant role in enhancing the micromixing. When the vapour ejection direction is inclined towards the wake direction, the wake is elongated, but the shape of the stoichiometric mixture fraction isosurface is not changed much, indicating a limited effect on droplet grouping in a spray. On the other hand, when the ejection direction is inclined towards the transverse direction, the stoichiometric surface is pushed further away in the transverse direction, and its topology is changed by the pufﬁng. The trajectories of ejected ethanol vapour pockets can be predicted by the correlation obtained for a jet in cross-ﬂow, and the vapour pockets may reach a few diameters away from the droplet. Therefore, in a multiple-droplet conﬁguration, the transverse ethanol vapour ejection due to pufﬁng may transiently change the droplet grouping characteristics. In simulation cases with multiple droplets, the interaction changing the droplet grouping due to pufﬁng has been conﬁrmed, especially for droplets in the most upstream position in a spray. This implies that pufﬁng should be accurately included in the mixing and combustion modelling of such a biofuel-blended diesel spray process.

Puffing-enhanced mixing of emulsion droplet

Introduction
Liquid-fuel spray combustion is widely used in combustion engines. Reduction of carbon dioxide (CO 2 ) emissions of engines is critical to solving global environmental issues. One of the approaches for such purposes is to use fuels in a blended form, such as a blend of fossil fuels and biofuels (Agarwal 2007;Kohse-Höinghaus et al. 2010;Shahir et al. 2014), to maximise the benefits of each component fuel. In using such fuel blends, differences in physical properties may induce more complex fluid dynamic issues. Heating, evaporation, mixing and combustion will be all affected, and a thorough understanding is needed to fully utilise the benefits of such fuel blends.
In this study, addition of bioethanol to fossil oil is considered. Ethanol (C 2 H 5 OH) is oxygenated, i.e. an ethanol molecule contains an oxygen atom, which enhances the oxidation of carbon, thereby reducing particulate matter (PM). However, ethanol is not well miscible with oil (Shahir et al. 2014). With the help of a surfactant, it can be stabilised in a form of emulsion. Several experimental attempts have been conducted to demonstrate the usefulness of blended fuels. For example, for diesel/biodiesel/bioethanol emulsion mixture, the overall engine performance has been demonstrated (Satgé de Caro et al. 2001;Hansen, Zhang & Lyne 2005;Rajasekar et al. 2010;Hulwan & Joshi 2011;Pidol et al. 2012). However, the quantitative assessment sometimes tends to be combustor-specific. There have also been studies on fundamental droplet and spray processes such as evaporation, mixing and combustion of several oil/alcohol blends (Yang, Jackson & Avedisian 1990;Jackson & Avedisian 1998;Botero et al. 2012;Moon et al. 2013;Pan & Chiu 2013). Under certain conditions, microexplosion was observed. In these studies, the puffing and microexplosion dynamics was revealed and its impact on the mixing between the fuel vapour and the air was discussed. In experiments, however, detailed vapour field data are not easy to obtain. Therefore, the physical mechanisms of puffing effects on fuel/air mixing of emulsion fuel droplets remain unknown. This fact has motivated the present simulation study.
When the ethanol proportion becomes high, puffing or microexplosion will occur. Puffing is occasional ejection of boiled ethanol vapour from inside the parent diesel droplet. If the rapidly growing vapour bubble violently breaks up the parent droplet, it is usually termed as microexplosion. Puffing and microexplosion are physical processes that deserve special attention (Yang, Jackson & Avedisian 1990;Jackson & Avedisian 1998). Puffing and microexplosion are particular to blended-fuel droplets, and are caused by distinct differences in physical properties, such as boiling temperature and superheat limit (Avedisian & Glassman 1981;Avedisian & Sullivan 1984). Since puffing can break up an emulsion droplet, it enhances secondary breakup or atomisation. When it occurs, the local equivalence ratio of the reactive gas is unsteadily varied due to the ejected ethanol vapour. If puffing can be properly utilised in a blended-fuel spray, it may enhance fuel atomisation and fuel/air mixing. Therefore, it is important to better understand puffing and its effect on fuel/air mixing. In a combustor, it is generally expected that puffing (or partial microexplosion) is more likely than complete microexplosion, since the heating time for spray droplets is limited (Watanabe & Okazaki 2013;Shinjo et al. 2014). Figure 1 exemplifies the significance of puffing on fuel/air mixing and combustion. This figure shows photographs taken at different time instants in an experiment on diesel/biodiesel/bioethanol emulsion droplet combustion conducted by our research group (Avulapati et al. 2016). The blend composition was 50 % diesel, 40 % biodiesel (rapeseed oil methyl ester; RME) and 10 % bioethanol. A droplet of diameter D = 2 mm was used to enable visualisation. As shown in the images, the droplet (the transparent sphere indicated in figure 1a) was anchored on a thermocouple fibre, and ignition was achieved by a hot electric coil (the bright illuminated spot shown in figure 1a) on the right-hand side of the droplet. It was confirmed that the ratio of ethanol plays an important role in the flame characteristics. When the ratio of ethanol was low (∼5 %), smooth burning was observed. As the ratio of ethanol was raised, puffing was occasionally observed, as shown in figure 1(b). The additional bright flame region on the upper-left side of the droplet indicates that ethanol vapour was ejected and the local equivalence ratio became temporarily high. Puffing was repeated during the droplet lifetime. Similar results have been reported in Yang, Jackson & Avedisian (1990), Jackson & Avedisian (1998), Botero et al. (2012) and Pan & Chiu (2013). These results indicate that puffing can affect the gaseous fuel-vapour/air mixture around the droplet transiently and non-uniformly, thereby impacting on the droplet combustion characteristics. A further understanding of puffing and its effects on mixing is therefore needed. If multiple droplets exist in the vicinity, puffing may affect the mode of droplet group combustion, which is also important for spray combustion in the dense and transitional spray regimes. Depending on the inter-droplet distance, the combustion mode changes from the single-droplet combustion mode to the group combustion mode (Chiu, Kim & Croke 1982;Umemura 1994). A diagram for group combustion was proposed by Chiu, Kim & Croke (1982), where four modes were identified, namely single-droplet combustion, internal group combustion, external group combustion and external sheath combustion. The group combustion mode is determined by the mixing of the oxidiser with the fuel vapour, because a flame is formed around the stoichiometric gas mixture. The original concept assumed no flow convection, and the inter-droplet distance was of primary importance. Later, Umemura (1990) and Umemura & Li (1992) extended the group combustion concept to include the gas convection effect. As the droplet Reynolds number (Re) rises, the stoichiometric surface is dragged towards the wake direction. Wu & Sirignano (2011a) investigated the group combustion characteristics using direct numerical simulation (DNS), with the inter-droplet distance in the transverse direction fixed at l = 2.4D (D = 50 µm). For Re = 5.5, a large group flame was formed, while for Re > 13, separate flames were formed around each droplet and not grouped in the transverse Puffing-enhanced mixing of emulsion droplet 447 direction. This indicates that not only the geometrical distances between the droplets but also the convective flow conditions are important for group combustion (Imaoka & Sirignano 2005;Wu & Sirignano 2011a,b;Zoby et al. 2011b;Sirignano 2014). Puffing will further complicate the group combustion dynamics due to the ejection of boiled ethanol vapour (a secondary fuel) from the parent diesel droplet, which produces primary diesel fuel vapour through normal evaporation. This may temporarily change the droplet grouping characteristics. Therefore, it is interesting and important to investigate the puffing effects on mixing not only in a single droplet but also in a multiple-droplet configuration. It should also be mentioned that the reactive gas mixture under puffing would be different for the two case scenarios where a droplet in a group is already in the wake of other droplets or the droplet is located in front of other droplets (Sirignano 2014).
The main objective of the present study is to improve the understanding of puffing effects on fuel/air mixing of blended emulsion in the configurations of both a single droplet and a droplet group under realistic convective heating conditions. As described above, the puffing mode is mainly targeted to link the results to a spray in a combustor. Direct numerical simulation is used to fully resolve the evaporating, puffing and mixing dynamics of the interfacial multiphase flow. Combustion is not included in this study; it will be extended to combustion cases in the future. Decane is used as a surrogate fuel. Namely, ethanol-in-decane emulsion is considered. The boiling point temperature of decane is substantially higher than the boiling point temperature and superheat limit of ethanol, which makes the numerical set-up suitable for the study of puffing and mixing dynamics. This setting does not limit the applicability of the present results to this particular fuel combination, rather they could represent characteristics for blends of ethanol and similar hydrocarbon fuels whose boiling point temperature is higher than the boiling point temperature and superheat limit of ethanol. Results are mostly discussed from the viewpoint of the physics of fluid dynamics. It should also be pointed out that using decane, the observed phenomena were similar in our experiments (Avulapati et al. 2016).
The present study has taken advantage of the research outcome of our previous studies. In Shinjo et al. (2014), the droplet breakup dynamics due to puffing and microexplosion was directly simulated under quiescent ambient conditions, although the convective effect on heating and vapour mixing was not included in the first-stage work, where understanding the puffing dynamics was the research aim. In Shinjo et al. (2016), convective heating of an emulsion droplet prior to puffing and microexplosion was investigated. A model was proposed to approximate the temperature distribution inside an emulsion droplet under convective heating. It should be pointed out that the temperature distribution inside a parent droplet largely determines the bubble initiation location -an important initial condition for puffing simulation. In the present study, the simplifications made in the previous studies have been removed, and the model developed in Shinjo et al. (2016) is used to initialise the inner-droplet temperature distribution. Puffing-enhanced mixing will be investigated in fully three-dimensional configurations for both a single emulsion droplet and a droplet group under convective heating. By considering realistic combustor flow conditions, outcomes from the present study can be useful for practical-scale blended-fuel spray simulation, in which puffing and microexplosion effects on the spray are modelled.
The rest of this paper is organised as follows. In § 2, the formulations and numerical methods are detailed. In § 3, the case set-up is described. In § 4, the results and discussion on puffing and its effects on fuel/air mixing are presented. Finally, in § 5, the concluding remarks are given.

Formulations and numerical methods
The governing equations are the Navier-Stokes conservation equations of mass, momentum, energy and species (Shinjo & Umemura 2010Shinjo et al. 2014Shinjo et al. , 2016Shinjo, Xia & Umemura 2015). Four species, namely C 10 H 22 , C 2 H 5 OH, N 2 and O 2 , are considered, and chemical reactions are not included. The governing equations are Here, ρ is the density, u is the velocity, T is the temperature, Y i is the mass fraction of species i, p is the pressure and c v is the specific heat at constant volume. The general form for P TH is P TH = T(∂p/∂T) ρ , and for an ideal gas P TH = p. The term Q u is the acceleration due to viscous force, given as where µ is the viscosity, φ is the colour function (φ = 1 for liquid and φ = 0 for gas) and δ ij is Kronecker's delta. The subscript L denotes the liquid phase and G the gas phase. The term F k u is the acceleration due to surface tension at interface k (the interface can be between liquid and gas or liquid and liquid). It is formulated using the CSF (continuum surface force) method (Brackbill, Kothe & Zemach 1992) as where σ k is the surface tension coefficient of interface k, [ρ] k = |ρ k+ − ρ k− | is the density difference at the interface and ρ k = 1/2 · (ρ k+ + ρ k− ). The subscript k+ indicates a value on one side of the interface and k− on the other side. Here, κ k is the local curvature given by where n k is the surface-normal unit vector. Using a level-set function for interface k, F k (defined later), n k is (2.6) The term Q Y i is the mass fraction change rate due to diffusion and Q T is the temperature change rate due to heat conduction, diffusion and the work done by Puffing-enhanced mixing of emulsion droplet 449 viscous forces, given as where D dif is the diffusion coefficient, λ is the thermal conductivity, h i is the enthalpy of species i and V i is the diffusion velocity for species i modelled by Fick's law. The terms S * are the evaporation source terms described below. The pressure is obtained by solving the pressure Poisson equation derived from (2.1), where c s is the speed of sound and the superscript * denotes the updated pressure value at time t * = t + t. This method to include compressibility effects is called the CUP (combined and unified procedure) formulation (Yabe, Xiao & Utsumi 2001;Shinjo et al. 2014). For the gas phase, the perfect-gas equation of state is used. For the liquid phase, Tait's equation of state is used, where the subscript 0 denotes reference values, A w = 296.3 MPa and γ w = 7.415 (Shinjo et al. 2014). The phase change is formulated by the jump conditions at the evaporating or boiling interface between liquid and gas (Tanguy, Ménard & Berlemont 2007). The superscript k is omitted below for simplicity. At the outer surface of a parent droplet, evaporation of decane occurs since the liquid decane temperature remains below the boiling point temperature. The vapour mass fraction at the droplet surface Y V is determined by the Clapeyron-Clausius relation where p V is the vapour pressure, p ∞ is the ambient pressure, h l is the latent heat of evaporation,R is the universal gas constant, T b is the boiling temperature, W V is the vapour molecular weight and W G is the ambient gas molecular weight. The evaporation rateω is determined by the jump conditions of heat and mass transfer at the interface as (2.13) The square brackets denote the difference of a variable f between the liquid and gas phases at the interface, i.e.
[ f ] = f L − f G . In one or multiple ethanol sub-droplets inside the parent droplet, explosive boiling occurs after the ethanol is superheated. At the interface between the liquid ethanol and the gaseous ethanol vapour, the jump 450 J. Shinjo, J. Xia, L. C. Ganippa and A. Megaritis conditions are also given by the above formulation with Y V = 1. The superheat degree T = T − T b determines the mass boiling rate in (2.12). The velocity at the interface satisfies (Tanguy, Ménard & Berlemont 2007;Shinjo et al. 2014) (2.14) where the surface velocity u s is defined as the sum of the liquid velocity and the surface regression velocity, namely u s = u L + s L and s L = s L n = (ω/ρ L )n. The velocity jump is (2.15) Therefore, the source terms for phase change in (2.2b) are where δ is the delta function and is used to identify an interface, and c p is the specific heat at constant pressure (Shinjo et al. 2014).
In addition to (2.1), level-set functions are solved to capture the liquid/liquid and liquid/gas interfaces (Sussman, Smereka & Osher 1994;Sussman & Puckett 2000). The level-set function for interface k, denoted as F k , is a signed-distance function and follows (2.17) Here, F k = 0 represents the interface. For a fluid system with multiple interfaces, multiple level-set functions can be superimposed. The surface regression velocity of interface k, s k L , is zero for a liquid/liquid interface. For ethanol sub-droplets in a parent decane droplet, any neighbouring sub-droplets are captured by different level-set functions, although physically all of the sub-droplet surfaces have no difference. This is a numerical setting to prevent unphysical coalescence of ethanol sub-droplets at low temperature. This approach, called the MLE (multiple level-set functions for emulsion) method, has been validated in Shinjo et al. (2016). The volume conservation is improved by combining the level-set method with the MARS (multi-interface advection and reconstruction solver) method, a kind of VOF (volume of fluid) method (Kunugi 1997). The numerical scheme used for the advection calculation is the CIP (cubic interpolated pseudo-particle or constrained interpolation profile) method (Takewaki, Nishiguchi & Yabe 1985).
The numerical code used is MEX (MicroEXplosion) (Shinjo & Umemura 2010Shinjo et al. 2014Shinjo et al. , 2015Shinjo et al. , 2016. The code has been used for direct simulation of interfacial multiphase flows such as droplet pinch-off (Shinjo & Umemura 2010), turbulent liquid spray atomisation (Shinjo & Umemura 2010 and evaporating spray dynamics (Shinjo & Umemura 2013;Shinjo et al. 2015). In Shinjo et al. (2014), the MEX code has satisfactorily demonstrated its capability of directly simulating boiling surface dynamics (Stefan problem) and linear/nonlinear droplet oscillation, which are directly related to puffing and microexplosion. In Shinjo et al. (2016), the MEX code has demonstrated its capability of directly simulating convective heat/mass transfer of an emulsion droplet, including the internal circulation motion inside the droplet.

Case set-up
An ethanol-blended diesel spray would be the ideal configuration to study puffing effects on fuel/air mixing in a fuel spray. Due to computational cost considerations, however, convective heating of a single emulsion droplet and a droplet group has been taken as the computational configuration in the present study. It is therefore important to set up proper boundary conditions to approximate realistic flow conditions in a fuel spray process.
In a fuel spray, the relationship between the droplet scale (diameter D) and the turbulence scale (Kolmogorov scale η) is D ∼ η in the dense spray region near the nozzle, where the generation of turbulence is mainly due to the formation of boundary layers and wake flows around finite-volume droplets (Sirignano 2010;Shinjo et al. 2015). In the downstream dilute spray region where the relative velocity between the droplet and the gas is nearly zero, the scale relation is D < η. The droplet Reynolds number, Re = U r D/ν G , where U r is the relative velocity between the droplet and the gas, and ν G is the gas kinematic viscosity, is at most Re ∼ O(10) in the near-nozzle region under typical combustor conditions (Sirignano 2010). In this study, the flow conditions are set in this Reynolds number regime, with the droplet diameter D = 30 µm, which gives D ∼ η. Therefore, the free-stream gaseous flow is laminar with a fixed velocity and no fluctuations.
For a droplet group, its geometrical location in the spray also matters. Figure 2(a) shows a schematic of a near-nozzle dense spray, where some droplets are located on the periphery of the dense spray area and directly exposed to the hot air, while most other droplets are located inside a large droplet cloud. The heating conditions are different for droplet groups of these two categories. Therefore, different boundary conditions should be used to approximate the two different heating conditions for a droplet group.
Taking into account the above considerations, table 1 shows the cases simulated in this study. Cases A1-A5 are single-droplet cases illustrated in figure 2(b), where the number of nucleation is varied inside the parent decane droplet. Cases B1 and B2 are multiple-droplet cases. For cases A1-A5, the flow velocity is fixed at the upstream boundary and free outflow conditions are given at all of the other boundaries. For case B1, the incoming velocity is fixed, free outflow conditions are given at the downstream boundary and periodic conditions are imposed at the side boundaries. For case B2, periodic conditions are given at all of the boundaries. Therefore, in cases A1-A5 and B1, the upstream droplet(s) is/are always exposed to the fresh hot air flow, mimicking 452 J. Shinjo, J. Xia, L. C. Ganippa and A. Megaritis  For the upstream boundary condition, F stands for fixed-velocity convection and P for the periodic condition. It is therefore implied that the droplet group is located on the periphery of a dense spray for B1 and inside a large droplet cloud for B2.
a droplet or a droplet group on the spray periphery (see figure 2b,c). In case B2, all of the droplets are in the wake region of other preceding droplets, mimicking a droplet group inside a large droplet cloud (see figure 2a). The volume fraction of ethanol and the sub-droplet size affect the degree of breakup (Shinjo et al. 2014). These are input parameters for the system, basically determined at the stage of fuel preparation. In this study, a realistic fuel spray is considered for which a moderate volume fraction is practical. Based on the results of Shinjo et al. (2014), the sub-droplet number and size are chosen to be moderate so that they characterise a representative puffing-dominant mode. For all of the cases, the diameters of the parent decane droplet and ethanol sub-droplets are D = 30 µm and D sub = 4.6 µm respectively. The number of ethanol sub-droplets is 19 in a parent decane droplet, and the volume fraction of ethanol is 7.4 %. The air velocity is U G = 10 m s −1 , the air temperature T G = 900 K, the air pressure p = 1 MPa and the droplet Reynolds number Re = 30. The number of nucleation in cases A1-A5 is varied to cover a range of the degree of breakup. Physical properties such as heat capacity c p , thermal conductivity λ, viscosity µ, latent heat of evaporation h l and surface tension σ are retrieved from the databases of NIST (National Institute of Standards and Technology) (NIST 2011) and Jasper (1974). For all of the physical properties, temperature dependence is considered.
The ambient pressure is set at p = 1 MPa. At this pressure, the boiling temperature of decane is 565 K, which is higher than that of ethanol (425 K). The superheat limit of ethanol is 477 K. The embedded ethanol sub-droplets can be superheated up to this limit, which is still below the boiling temperature of decane. As the liquid ethanol temperature approaches the superheat limit, nucleation (initial generation of vapour bubbles) occurs and triggers explosive boiling. Nucleation is strongly dependent on the local liquid temperature. In homogeneous nucleation theory, the nucleation rate J (m −3 s −1 ) is given by Avedisian & Andres (1978), Avedisian & Glassman (1981) and Avedisian (1985) as where the activation energy is ∆ = 16πσ 3 /3(p V − p L ) 2 , Γ is a pre-exponential factor, k f is a rate constant, N 0 is the molecular number density, k is the Boltzmann constant, σ is the surface tension coefficient, p V is the vapour pressure and p L is the liquid pressure. The nucleation rate rapidly increases above a certain temperature due to the exponential effects of temperature. It is impossible to deterministically predict the nucleation location and timing in the framework of the continuum dynamics; therefore, in this study nucleation is initiated by placing tiny vapour bubbles. Past findings show that the bubble nucleation is more likely to occur at a liquid-liquid interface rather than deep inside a sub-droplet ( FIGURE 3. (Colour online) Temperature field before vapour bubble nucleation. (a) Schematic of the liquid temperature inside the parent droplet predicted by the ECME model. (b) Temperature (K) field before nucleation. The liquid shape is represented by green isosurfaces and the temperature distribution is drawn on a two-dimensional cut plane. The solid black line shows the parent droplet surface on the cut plane. The convective air flow is from −y to +y.
1981; Avedisian 1985). To mimic this in the present simulation, bubbles are placed within the ethanol sub-droplets adjacent to the ethanol/decane interface. This has been also discussed in Shinjo et al. (2014). The droplet heating process prior to puffing is much slower than puffing. Inside an emulsion droplet, an internal circulation motion is developed due to convective heating, and therefore the inner-droplet temperature is not uniform. To properly approximate the heating outcome to save computational time, the model proposed in Shinjo et al. (2016) is utilised. In this previous study, the flow physics of convective heating of an emulsion droplet in a transient liquid Péclet number regime (100 < Pe L < 500) was elucidated. This regime corresponds to realistic flow conditions in a combustor at Re ∼ 30. The initial temperature distribution in the parent droplet is approximated by the ECME (effective conductivity with modified eccentricity) model, which calculates the inner-droplet liquid temperature by using an eccentric axisymmetric coordinate system to reproduce the internal circulation effect, as schematically shown in figure 3(a).
During the droplet heating process, coalescence of embedded ethanol sub-droplets is expected to occur gradually, which changes the number and size of the sub-droplets. This is due to thermocapillary migration induced by the non-uniform liquid temperature distribution inside a parent diesel droplet and reduced effectiveness of the biodiesel surfactant at high temperature (Kadota & Yamasaki 2002). However, this coalescence process is slow and only has a minor effect on a fuel spray, since the elapsed heating time is short for complete coalescence to occur in a combustor (Watanabe & Okazaki 2013;Shinjo et al. 2014). Even when puffing starts, the inner-droplet liquid temperature (except in the near-surface region) is still relatively low where the surfactant still maintains the repelling force between the sub-droplets. For the time scale of puffing studied here, sub-droplet coalescence does not occur even if the sub-droplets are clustered. Therefore, in this study, coalescence of sub-droplets is not simulated. Instead, the simulation cases start with finite-size ethanol sub-droplets, and the number of nucleation is varied to change the strength of puffing and thus the extent of droplet breakup. The numerical procedure for simulation initialisation in this study is as follows. First, the initial inner-droplet liquid temperature field is given by the ECME model. Then, the coupled gas/liquid heating simulation is started, initially with uniform gasphase temperature and velocity. This fully coupled heating simulation is continued until the gas-phase velocity and species concentration fields around the droplet reach a quasi-steady state, as shown in figure 3(b). Then, nucleation is triggered. Inside the parent droplet, the liquid temperature is high in the rear region (Shinjo et al. 2016). Therefore, for all of the cases, the nucleation starts from the rear region of a parent droplet.
To resolve the embedded sub-droplet structures, the grid spacing is equidistantly set as x = y = z = 0.19 µm, and the total number of grid points is 243 million for cases A1-A5 and 576 million for cases B1 and B2. A grid convergence study was conducted in Shinjo et al. (2014), and the current grid spacing is smaller than that used there. The simulations have been conducted on the UK national supercomputer ARCHER using 1728 and 4096 processor cores for the series-A and series-B cases respectively. The wall-clock time is 300-400 h for each case. Figure 4 shows the temporal sequence of puffing for case A2. In this figure, the parent decane droplet and embedded ethanol sub-droplets are shown by green isosurfaces, decane vapour by blue isosurfaces for Y decane = 0.1 and 0.2, and ethanol vapour by red isosurfaces for Y ethanol = 0.05 and 0.1. In this case, two distantly separated vapour bubbles are activated in the rear region of the emulsion droplet. As time passes, the bubbles grow to two substantially large ones before puffing occurs, because the two ethanol sub-droplets in which nucleation is triggered are rather deep inside the parent droplet. As the bubbles grow, the emulsion droplet finally ruptures and ejects boiled ethanol vapour.

Dynamics of puffing
From the ejection on the upper z side, which is directed towards the +xz direction, a ligament is pushed out due to the recoiling motion of the liquid surface after ethanol vapour ejection (figure 4a,b). The mechanism is schematically shown in figure 5. The ligament ejection occurs when the recoiling motion is strong. This is likely to occur FIGURE 5. Schematic of the puffing dynamics and after-puffing ligament ejection. when the bubble size is large before the emulsion droplet ruptures, namely when the boiling sub-droplet is deep inside the parent droplet as in the present case. When a boiling ethanol sub-droplet is close to the parent decane droplet surface, the recoiling motion is relatively weak and will lead to oscillations of the remaining sub-droplet (Shinjo et al. 2014). From the ejection in the lower z side in figure 4(b), only gaseous ethanol vapour is ejected and no liquid fragmentation occurs. The ejection is mostly directed towards the wake direction, where the mixture of decane/air is already fuel-rich. The ejected ethanol vapour (a secondary fuel) makes the mixture richer locally, and extends the wake in the downstream, but the width of the wake is not much affected in the side transverse direction.
The ejection velocity is examined, since it is related to the trajectory of an ethanol vapour pocket discussed in § 4.3.1. Although the bubble growth inside the parent droplet is quasi-steady, the ejection process is rather rapid and unsteady after a hole is generated. Figure 6 shows the sequence of the velocity field for the ejection in the lower z side. At t = 10.56 µs (figure 6a), the ejection hole opens and ethanol vapour starts to come out. At t = 11.22 µs (figure 6b), the ethanol vapour jet is accelerating. As schematically shown in figure 5, the periphery of the hole is shrinking due to surface tension and the hole is becoming larger quickly afterwards. At t = 11.88 µs (figure 6c), the magnitude of the ejection velocity is at maximum and the bubble size 456 J. Shinjo, J. Xia, L. C. Ganippa and A. Megaritis Ethanol vapour Vapour ejection After-puffing ligament ejection FIGURE 7. (Colour online) Temporal evolution of puffing. The frame rate is 4000 f.p.s.. The blend composition is decane 60 % and bioethanol 40 % in volume fraction. In (c), the ejecting ethanol vapour is outlined by the dashed lines. Puffing is repeated during the droplet lifetime.
has become smaller. Finally, at t = 13.20 µs (figure 6d), the jet stops as most of the ethanol vapour has ejected out. The jet velocity magnitude can be estimated from the vapour pressure p V of the bubble by using the quasi-steady Rayleigh-Plesset equation. The Rayleigh-Plesset equation predicts the growth of a spherical bubble and can be written as (Rayleigh 1917;Plesset & Zwick 1954;Mikic Rohsenow & Griffith 1970) where R is the bubble radius,Ṙ = dR/dt,R = d 2 R/dt 2 , p 0 is the ambient pressure, ρ L is the liquid density and σ is the surface tension coefficient. It should be stressed that the prediction of the Rayleigh-Plesset equation for the bubble growth is only used for comparison with the DNS results. In the DNS, the bubble growth is directly resolved. Although the actual bubble shape is not perfectly spherical, the result below indicates that the effect of non-sphericity is minor in the estimation. Assuming a quasi-steady balance between the pressure difference and the surface tension, the vapour pressure is p V ∼ p 0 + 2σ /R. The estimated pressure difference is p V − p 0 ∼ 2σ /R = 4.8 kPa, with σ ∼ 9 × 10 −3 N m −1 and R ∼ 0.25D, which is in good agreement with the observed pressure difference p V − p 0 ∼ 4.9 kPa in the present DNS just before the droplet rupture. Using ρ V U 2 jet /2 = p V − p 0 , the jet velocity is U jet = √ 2(p V − p 0 )/ρ V ∼ √ 4σ /ρ V R. Using the observed representative values in the DNS, p V − p 0 ∼ 4.9 kPa, σ ∼ 9 × 10 −3 N m −1 , ρ V ∼ 13.9 kg m −3 , R ∼ 0.25D, the jet velocity is estimated to be U jet = 26 m s −1 , which is close to the observed magnitude of U jet ∼ 20 m s −1 in figure 6. When puffing is accompanied by ligaments or liquid fragments, as in the ejection in the upper z region in figure 4, the estimation of the jet velocity is more complicated, but the puffing dynamics is the same.
The puffing dynamics has been also visualised as direct images in our experiment (Avulapati et al. 2016), as shown in figure 7, for a decane/bioethanol droplet. Similarly to figure 1, a droplet of D = 2 mm is suspended on a fibre under combustion conditions. This figure is only used for qualitative comparison. In figure 7(a), some part of the embedded liquid ethanol has already become gaseous ethanol vapour by boiling. In figure 7(c), the ethanol vapour ejects towards the lower-right direction. Since it is difficult to see the ethanol vapour in this photo, the dashed lines are added to outline the vapour pocket. Ligament ejection is observed after the ethanol vapour ejection as a result of the recoiling motion, as shown in figure 7(e). The puffing and after-puffing dynamics is similar to the simulation result shown in figure 4. This indicates that the present numerical results reproduce the dynamics of puffing properly, although it should be pointed out that the fibre used to anchor the droplet has an influence on nucleation, and the location and direction of boiled ethanol vapour ejection will also be affected. Therefore, the experimental result is only qualitatively referenced here.
Puffing-enhanced mixing of emulsion droplet  FIGURE 9. Increase of the liquid surface area due to puffing. Boiling is initiated at t = 5.5 µs; S 0 is the initial surface area of the droplet.
As the number of boiling sub-droplets increases, the degree of breakup increases. Figure 8 shows the ethanol vapour ejection and droplet breakup dynamics for cases A3 and A4, where more ethanol sub-droplets boil compared with case A2 in figure 4. The boiling behaviour of each sub-droplet is similar to that seen in figure 4. If neighbouring vapour bubbles merge during the bubble growth, the vapour ejection becomes more intense ( figure 8a,b). As a result, the recoiling motion becomes intense and ligament ejection occurs, as shown in figure 8(b). Subsequently, secondary-droplet pinch-off occurs due to the surface tension effect. In case A4, where more ethanol sub-droplets boil than in case A3, the additional vapour bubbles near the side periphery of the parent droplet do not merge with other ethanol vapour bubbles, and thus the effect of the additional vapour ejection is found to be minor, as indicated by the comparison between figures 8(b) and 8(c). This is similar for case A5. Therefore, in § § 4.2.3 and 4.2.4, A2 and A3 are chosen as two representative multi-nucleation cases for further analysis of vapour mixing due to puffing.
The profiles of the liquid surface area for cases A1-A5 are plotted in figure 9. For all of the cases, boiling is initiated at t = 5.5 µs. As the boiling progresses, the droplet breakup increases the liquid surface area. Generally, breakup of the parent emulsion droplet becomes more intense as the number of boiling sub-droplets increases. A large increase in the surface area of cases A2 and A3 indicates that the bubble merging motion has induced large breakup, as shown in figure 8. A relatively small difference in the surface area between A3 and A4 (and also A5) also confirms the limited effect of additional ethanol vapour ejection in cases A4 and A5, as stated above.

4.2.
Effects of puffing on fuel/air mixing 4.2.1. Mixture fraction Z and scalar dissipation rate χ The mixture fraction Z and the dissipation rate of the mixture fraction, or the scalar dissipation rate (SDR), χ quantify the mixing and micromixing of the reactants in the gas phase respectively, and are two key variables in a variety of non-premixed flame models. In the present study, two mixture fractions Z 1 and Z 2 are used to describe the mixing of the vapour of the primary fuel decane and the secondary fuel ethanol with air. They are defined as where Y i is the mass fraction of species i. Setting Y i,0 = 0 in the far field and Y i,1 = 1 at the droplet centre, the mixture fraction Z i reduces to the fuel mass fraction Y F,i (Zoby et al. 2011a), namely Z 1 = Y decane and Z 2 = Y ethanol in the gas phase. Without chemical reactions, the transport equation for each mixture fraction in the physical space where S Z i is the evaporation source term. Conversion into the Z space, following the framework of the flamelet model, yields (Peters 2000 where Le i = λ i /ρc p,i D dif is the Lewis number of species i and S Y i is the evaporation source term. A similar equation can be derived for the temperature T (Hasse & Peters 2005). In these equations, the SDRs characterise the diffusion effect. They are defined as Assuming that the mixing is dominant in the compressive strain direction (Hasse & Peters 2005), the SDR for the primary fuel decane χ 1 is governed by rate makes the local gradients steeper, which contributes positively to the rate of change of the SDR. Three diffusion terms exist since there are two fuels. The effect of the evaporation source term is confined in the vicinity of an evaporating/boiling liquid surface. An analogous equation can be formulated for χ 2 . It is clear from (4.4) and (4.6) that the additional terms due to the dissipation rate χ 2 and cross-SDR χ 12 quantify the direct puffing effects on the mixing and micromixing of the vapour of the primary fuel decane with air.

4.2.2.
Decane-vapour/air mixing field before puffing First, the decane-vapour distribution before puffing is examined (Z 2 = χ 2 = χ 12 = 0). The data samples are taken from a line in the transverse direction at two downstream positions of l 1 = 2.3D and l 2 = 1.3D from the droplet centre, as shown in figure 10(a). It is known that the fuel-vapour field in the Z space in the wake of a droplet takes a similar form to the gaseous fuel field of a laminar diffusion flame in the Z space (Zoby et al. 2011a). The conditional mean χ 1 | Z 1 , or simplified as χ | Z in this subsection, downstream of an evaporating fuel droplet has been modelled using an error function profile as (Peters 2000;Hasse & Peters 2005; where χ ref is a reference value, erfc −1 is the inverse complementary error function and Z max is the maximum mixture fraction value. Another proposed model is (Zoby et where α = U/4D dif l, A =ṁ d (Z d − Z ∞ )/4πρD dif l,ṁ d is the mass evaporation rate, Z d and Z ∞ are the mixture fractions inside and far from the droplet respectively, and l is the streamwise distance of the considered location from the droplet centre. This solution is based on the assumption where r tr is the transverse distance from the droplet centre. Figure 10(b) shows the comparison of the observed Z-χ correlation for decane from the DNS results and the fitted F(Z) and S(Z) at l 1 = 2.3D in case A2. Both . (Colour online) Doubly conditional mean SDRs χ 1 | Z 1 , Z 2 , χ 2 | Z 1 , Z 2 and χ 12 | Z 1 , Z 2 (s −1 ) for case A2. It should be noted that the contour legends are different. Data points only exist in Z 1 + Z 2 1 by definition: (a-c) t = 11.88 µs; (d-f ) t = 13.20 µs.
the F(Z) and S(Z) curves are in good agreement with the DNS results, and F(Z) shows slightly better agreement with regard to its curve shape. At a closer distance of l 2 = 1.3D, the observed Z-χ correlation and the F(Z) curve again show good agreement. This indicates that the two models, especially F(Z), can properly predict the fuel/air mixing before puffing occurs. This information is used as a reference for the following discussion.

Weak puffing effects on mixing
Puffing of boiled vapour of the secondary fuel ethanol is expected to enhance the gas-phase mixing in both space and time. Here, effects of weak puffing, with no ejection of liquid fragments, on mixing are first investigated. Since the primary fuel is decane, puffing effects on χ 1 and Z 1 are of particular interest. The data samples taken for analysis in this section are from the rectangular region in figure 4(b) in case A2. For this case, most of the ethanol vapour is ejected into the decane-rich wake region, and the shape of the stoichiometric mixture fraction (Z st ) isosurface is not strongly affected during the puffing. However, the ethanol vapour strongly affects the production of decane vapour from the parent droplet surface. It is therefore important and interesting to investigate the ethanol vapour puffing effects on the local fuel/air mixing. Figure 11 shows the SDRs doubly conditioned on Z 1 and Z 2 , χ 1 | Z 1 , Z 2 and χ 2 | Z 1 , Z 2 , and the cross-SDR, χ 12 | Z 1 , Z 2 , at t = 11.88 µs when the puffing is most intense (figure 4b), and at t = 13.20 µs when puffing-enhanced mixing has progressed further (figure 4c). It should be noted that the upper limit of the mixture fractions in the figure is bounded by Z 1 + Z 2 = 1 by definition. Immediately after the ejection starts, χ 1 and χ 2 increase, as shown in figure 11(a-c). Cross-mixing of the decane and ethanol fuel vapour is also evident, as χ 12 is negative. The SDRs in the region of large values of Z 2 are due to the ethanol vapour ejection and those in the region of large values of Z 1 are due to the evaporation of decane from the parent droplet. The distributions of χ 1 and χ 2 are qualitatively symmetric in this region where weak puffing occurs. As the mixing progresses further, the magnitude of the SDRs decreases, as shown in figure 11(d-f ) at t = 13.20 µs.
To focus on puffing effects on the mixing of the primary fuel decane with air, the mean SDR conditioned on Z 1 , χ 1 | Z 1 , is shown in figure 12(a). By comparing figure 12(a) with figures 10(b) and 10(c), in which χ 1 | Z 1 is presented before puffing occurs, two conclusions are evident. First, χ 1 | Z 1 during the puffing cannot be approximated by F(Z) or S(Z). Second, the magnitude of χ 1 | Z 1 has increased by one to two orders of magnitude after puffing, demonstrating that the micromixing has been considerably enhanced by puffing. At t = 11.88 µs, enhanced mixing due to the ejected ethanol vapour is evident in Z 1 = 0.1-0.25. The other peak in the region of large values of Z 1 is due to the evaporation of decane from the parent droplet. At t = 13.20 µs, the puffing jet has already stopped (figure 6d) and χ 1 is much smaller, which indicates that puffing-induced mixing has weakened. Therefore, puffing enhances mixing in a short period of the ethanol vapour ejection.
To better understand puffing effects on mixing in quantified detail, budget analysis of the conditionally averaged (4.6) is performed. The contributions of each conditionally averaged right-hand side term to the rate of change ∂ χ 1 | Z 1 /∂t are shown in figure 12(b). Such an analysis will shed crucial light on how puffing would affect the modelling of χ 1 | Z 1 in the framework of the conditional moment closure (CMC) approach. These budget terms are conditionally averaged in the indicated region in figure 4(b), excluding the region within 5 x of the liquid surface of the parent droplet to remove the local peak evaporation source terms for clarity. The production term 2aχ 1 in (4.6) is approximated as 2Sχ 1 , where S is the norm of the strain rate tensor S ij = (1/2)(∂u i /∂x j + ∂u j /∂x i ).
As can be seen in figure 12(b), the dissipation term is negative by definition, and its magnitude is small. One leading positive contribution is from the production term due to the strong strain rate of the ejecting ethanol vapour jet (figure 6). It is interesting to note that the cross-diffusion due to χ 12 is another strong source that is comparable to . (Colour online) Doubly conditional mean SDRs χ 1 | Z 1 , Z 2 , χ 2 | Z 1 , Z 2 and χ 12 | Z 1 , Z 2 (s −1 ) for case A3. Data points only exist in Z 1 + Z 2 1 by definition: (a-c) t = 9.24 µs; (d-f ) t = 10.56 µs.
the production term. This indicates that puffing-induced cross-diffusion of the primary and secondary fuels works to increase the magnitude of χ 1 in the full spectrum of Z 1 . The other two diffusion terms related to χ 1 and χ 2 are negative, which represents a normal diffusion process. At t = 11.88 µs, the puffing motion is strong with mixing in progress. The net contribution of all of the right-hand side terms in figure 12(b) is negative, which decreases χ 1 . Later, at t = 13.20 µs, the magnitudes of these terms are smaller, which means that the puffing jet has stopped and the fuel/air mixing has progressed further.

Strong puffing effects on mixing accompanying ejection of liquid-fuel ligaments
In case A3, the ejecting vapour and ligaments interact with each other in the region shown by the rectangle in figure 8(a). This ejection is oriented towards the side transverse direction and strongly influences the shape of the stoichiometric mixture fraction isosurface (see § 4.3). Figure 13 shows the doubly conditional mean SDRs at t = 9.24 µs and t = 10.56 µs, corresponding to the initial and middle stages of the puffing. The distributions of the SDRs are similar to those in case A2, for example between figures 13(d-f ) and 11(a-c). The magnitude of the SDRs is larger, since the puffing jet ejection motion is stronger than in case A2. One obvious difference between figures 13 and 11 is that the SDR values are more strongly disturbed in figure 13 due to multiple ejections and ejected liquid fragments.
The conditional mean χ 1 | Z 1 is shown in figure 14(a), and the budget analysis using (4.6) is presented in figure 14(b). Similar trends to those in figure 12 are found. One important difference in figure 14(b) compared with figure 12(b) is that the cross-scalar diffusion due to χ 12 now takes the dominant role in enhancing χ 1 | Z 1 , whereas in figure 12(b) the magnitudes of the production and cross-scalar diffusion are comparable.
Through the detailed analysis, it is clear that ejection of boiled ethanol vapour initially makes a strong stratification in the reactive fuel gas mixture of decane and air, followed by mixing that homogenises the stratification. This indicates that puffing, which is an intermittent ethanol vapour ejection motion, will change the local mixing characteristics around the emulsion droplet temporarily and spatially. In the next subsection, puffing effects on droplet grouping will be investigated.

Effects of puffing on droplet grouping 4.3.1. Region of influence of ejected ethanol vapour pockets
A flame around an evaporating fuel droplet is basically a non-premixed (diffusion) flame. The fuel vapour is supplied from the droplet surface and the oxidiser from the ambient air. The stoichiometric mixture fraction isosurface is determined by the mixing balance between the fuel and the oxidiser. The stoichiometric condition well represents the flame structure for non-premixed combustion. Therefore, in the present study it is indicative to investigate how puffing affects the stoichiometric mixture fraction isosurface. A variable transformation is useful to define the overall mixture fraction Z and a mixing parameter Y as (Hasse & Peters 2005) (4.10b) For the combustion of two fuels of decane and ethanol, the global reaction can be considered as reactants (C 10 H 22 , C 2 H 5 OH, O 2 ) → products (CO 2 , H 2 O).  where one degree of freedom (Y, which represents the ratio of the two fuels) remains. Figure 15 shows the stoichiometric mixture fraction isosurface and streamlines for case A3 at three time instants. Initially, when the effect of puffing is limited, the stoichiometric surface surrounds the droplet almost axisymmetrically (figure 15a). The fuel-rich area extends towards the wake direction due to the convective flow. After the ethanol ejection starts in the transverse direction, it is clear that the fuel-rich region is rapidly expanding in the transverse direction ( figure 15b,c). Although the ejection direction is not fully perpendicular to the air flow direction and is slightly inclined towards the downstream direction, the effect of puffing on expanding the stoichiometric surface is evident. The ejected ethanol vapour pocket gradually becomes affected by the momentum of the air flow, and the trajectory of the pocket becomes bent further downstream. In figure 15(c), the ethanol pocket has reached a few diameters away from the parent droplet surface in the transverse direction.
The trajectory of an ejected ethanol vapour pocket is determined by the convective effect. The convective length scale is l c = U jet t and the diffusion length scale is l dif = D dif t during the time scale t. With the estimation of t = 1 µs, D dif = 1.0 × 10 −5 m 2 s −1 and U jet = 20 m s −1 , the scale ratio is l c / l dif = 6.3. Therefore, the convective effect is dominant compared with diffusion. This implies that the trajectory of an ethanol vapour pocket has a similarity to that of a jet in cross-flow (Keffer & Baines 1963;Kamotani & Greber 1972;Broadwell & Breidenthal 1984;Smith & Mungal 1998;Chochua et al. 2000;Su et al. 2000;Muppidi & Mahesh 2005). A correlation for the jet-in-cross-flow trajectory obtained by several experiments and simulations is known to scale with a power law as (4.13) where x d and x tr are the downstream distance (in the air flow direction) and the transverse distance (in the flow-normal direction) respectively, and d is the jet diameter. Here, r is the effective velocity ratio r = (ρ j u 2 j /ρ mf u 2 mf ) 1/2 , (4.14) Puffing-enhanced mixing of emulsion droplet  Muppidi & Mahesh (2005). For the two simulation cases from Muppidi & Mahesh (2005), the difference lies in the upstream boundary layer (BL) thickness.
where the subscript j denotes the jet and mf the main free-stream air flow. For a wide range of r, it is known that the constant α j is α j = 1/3-1/2 (Keffer & Baines 1963;Kamotani & Greber 1972;Broadwell & Breidenthal 1984;Smith & Mungal 1998;Chochua et al. 2000;Su et al. 2000;Muppidi & Mahesh 2005). Figure 16 shows the trajectory profiles of jets in cross-flow, which are replotted from figures 1 and 8 in Muppidi & Mahesh (2005), and those of several ethanol vapour pockets in two cases in the present study. In the present cases, the pocket position is defined as that of the centre of gravity (CG) of the pocket. Most experiments in the literature (Keffer & Baines 1963;Kamotani & Greber 1972;Broadwell & Breidenthal 1984;Smith & Mungal 1998;Chochua et al. 2000;Su et al. 2000) are in turbulent conditions, and the simulation cases of Muppidi & Mahesh (2005) are in laminar conditions. The trajectories are affected to some extent by the effective velocity ratio r, the boundary layer thickness of the main free-stream flow and the jet velocity profile (Muppidi & Mahesh 2005). However, all of the trajectories scale well with (4.13) within a limited range of scatter. This indicates that the trajectories of the ethanol vapour pockets ejected by puffing in the present study are also similar to the trajectory of a jet in cross-flow. In puffing, d cannot be strictly defined and d ∼ D/3 is used. It should also be noted that puffing is not continuous ejection but finite-duration ejection, and that the initial ejection direction is slightly inclined towards the downstream (x d ) direction, resulting in a relatively smaller x tr /rd compared with the jet-in-cross-flow results. Still, it can be said that the trajectory of an ethanol vapour pocket is similar to that of a jet in cross-flow. For a wide range of r, the penetration length in the transverse (x tr ) direction is approximately x tr /rd < 2. Considering the ejected ethanol vapour and air densities, the effective velocity ratio becomes r ∼ 1.8-3.3 in the present puffing cases. Assuming d ∼ D/3, this relation becomes x tr /D < 2.2. Therefore, the ejected ethanol vapour may reach a few diameters away from the parent decane droplet in the transverse direction under the current convective conditions (see figure 15). For an isolated droplet that is far away from other droplets, such as those in the downstream dilute spray region, this length scale may not be significant in terms of droplet grouping because flames surrounding droplets may remain isolated. When other droplets exist in the vicinity, such as those in the dense and transient spray regions, however, the ejected ethanol vapour may reach the region of influence of neighbouring droplets. Therefore, puffing may induce mutual interactions between droplets and change the group combustion configuration. In § 4.3.2, the droplet grouping under puffing will be investigated in the multiple-droplet configuration.
The primary-fuel SDR at the stoichiometric condition χ 1,st is strongly related to the local non-premixed flame structure. In a turbulent combusting flow, it is known that the probability density function (p.d.f.) of χ 1,st can be approximated by a lognormal distribution (Pitsch & Steiner 2000). In the present study, without puffing, the flow field around a droplet is laminar and axisymmetric. Therefore, at a fixed position downstream of the droplet in the streamwise direction, the p.d.f. of χ 1,st has a sharp peak at a specific value. Figure 17 shows the p.d.f. profiles of χ 1,st for case A3 at the downstream position of l 1 = 2.3D from the droplet centre. It should be noted that this position is the same as that shown in figure 10(a) for case A2. The corresponding flow field can be seen in figure 15. At t = 10.56 µs, the p.d.f. has a peak around χ 1,st ∼ 190 s −1 , which is in agreement with the peak SDR shown in figure 10(b). After puffing (t = 13.20 µs and t = 14.52 µs), the decane vapour field gradually becomes asymmetric and the p.d.f. profiles are shifted towards a larger χ 1,st . This indicates that the motion of ejected ethanol vapour pockets causes an increase of χ 1,st . Puffing occurs from the parent decane droplet surface, which is in a decane-rich region. Affected by the puffing motion, the stoichiometric surface is expanded and χ 1,st increases. This indicates that when a flame is existent, both the flame topology and the reaction rate will be affected by puffing.
Puffing-enhanced mixing of emulsion droplet

Droplet grouping under puffing
The above discussion in § 4.3.1 indicates that mutual droplet interaction becomes significant when the inter-droplet distance l is on the scale of D in a multiple-droplet configuration with puffing. In this subsection, puffing effects on droplet grouping characteristics are investigated using cases B1 and B2. The present setting mimics a droplet group in the dense spray region. Case B1 is a baseline case and has four droplets in the front and four in the back. Case B2 has the same geometrical droplet locations, but the streamwise boundary conditions are periodic. Thus, all of the droplets are in the wake region in this case. The droplet positions (normalised by D) are (x, y, z) = (0, 0, 0), (+1.47, +0.13, +1.28) and (+1.47, +2.05, +0.96) for droplets D1, D2 and D3 respectively (see figure 18a). Droplets D2 and D3 are placed in tandem in the streamwise direction and named as a droplet pair P2. The inter-droplet distance between D1 and P2 in the transverse direction (+xz direction) is l = 1.95D. D1 and P2 are not grouped at the current droplet Reynolds number of Re = 30 without puffing. The interaction between D1 and P2 under puffing and its effects on the grouping characteristics will be the main interest in this subsection. The For a droplet group, the interaction occurs in both the transverse and streamwise directions. For the effect of puffing on droplet grouping under convective conditions, the transverse direction is more important, as can be expected from § 4.3.1. Puffing itself is, once initiated, strong, fast and mostly independent of the free-stream flow. Even if there are neighbouring droplets, the dynamics of puffing and vapour ejection is similar to that for an isolated single droplet. Therefore, whether there are other droplets in the region of influence of a puffing droplet becomes the sole factor determining puffing effects on droplet grouping. Figures 18(a) and 18(b) show the temporal evolution of puffing in case B1. The yellow isosurface represents the stoichiometric surface and the pink isosurface represents Y ethanol = 0.05. Puffing occurs for droplet D1 located in the centre and upstream. Under the droplet Reynolds number Re = 30, the stoichiometric surface does not spread much in the transverse direction. In the pair P2, D3 is in the wake of D2. Although the transverse distance l is not large, the wakes of droplet D1 and P2 do not affect each other before puffing, as indicated in figure 18(a). As already stated, puffing progresses almost independently of neighbouring droplets and expands the local stoichiometric surface surrounding D1 in the transverse direction. The ejected ethanol vapour pockets travel towards the droplet pair P2. Finally, the fuel vapour from D1 starts to interact with the stoichiometric surface surrounding P2, as shown by the dashed arrow in figure 18(b).
To further quantify puffing effects on droplet grouping, the profiles of Z 1 , Z 2 and Z st along the free-stream flow direction (+y) at the middle transverse distance between D1 and P2 are shown in figure 18(c). The position of the line of interest is indicated by the dashed line in figure 18(a). The stoichiometric mixture fraction (the thin solid line in figure 18c) represents the criterion whether D1 and P2 interact or not, namely if Z = Z 1 + Z 2 < Z st , D1 and P2 are separate, and if Z 1 + Z 2 Z st , they are in the same group. At t = 9.24 µs, D1 and P2 are separate at y/D < 2.8, as indicated by the dashed red arrow in figure 18(c), which is also clear in figure 18(a). At t = 13.20 µs, the ejected ethanol vapour leads to the increased Z 2 in y/D = 1.8-3.2, as shown in figure 18(c). Accordingly, this region becomes fuel-rich, and D1 and P2 now belong to the same group at y/D > 1.8, as shown by the solid red arrow in figure 18(c). It should be noted that the apparent rise in Z 1 is a combined result of its magnitude increase and the stoichiometric surface displacement due to puffing of the stratified decane vapour, as shown in figures 15 and 18.
The DNS results of case B1 indicate that puffing may change the droplet grouping characteristics, even if the geometrical configuration is not changed during the puffing. The effect is significant especially for puffing in the transverse direction. The discussion in § 4.3.1 (figure 16) has indicated that a critically influenced region exists within a distance of a few droplet diameters away from the puffing droplet in the transverse direction. This conclusion gives important information for modelling puffing effects on droplet grouping and group combustion.
In case B2, the boundary conditions in the free-stream flow direction (+y) are changed to periodic conditions, mimicking a droplet group inside a large droplet cloud. In this case, all of the droplets are already in the wake of their preceding droplets and surrounded by fuel-rich gas mixture of decane and air. Figure 19(a) shows the stoichiometric surface distribution at t = 13.20 µs. It is clear that the ethanol puffing dynamics of D1 is similar to that shown in figure 18(b). However, the mixing of the ejected ethanol vapour with the decane/air reactive gas mixture does not change the droplet grouping characteristics. As shown in figure 19(b), a similar puffing effect to that in figure 18(c) can be seen, as shown by the increased Z 2 (and Z 1 ) in y/D = 1.8-2.6. In this case, however, the investigated region already has a fuel-rich mixture in the upstream due to the wake of preceding droplets, and D1 and P2 have already been grouped before the puffing. Therefore, the puffing effect is not significant in terms of the droplet grouping. As schematically shown in figure 2(a), actual droplet groups in a dense spray can be categorised into the above two cases, namely droplet groups on the periphery of the spray or droplet groups inside the spray. Recent studies on the droplet interaction in a multi-dimensional droplet array also reflect the conclusion that the interaction is different between the most upstream droplets and inner droplets (Sirignano 2014). For emulsion droplets with puffing, it is more significant to consider the configuration of case B1 which includes front droplets on the periphery of a spray, since the puffing probability is strongly related to the droplet heating conditions and is therefore different for the two kinds of droplet groups. It is known that downstream droplets in the wake region undergo less heat and mass transfer compared with front droplets (Tal, Lee & Sirignano 1984;Tsai & Sterling 1991;Chiang & Kleinstreuer 1992;Umemura 1994). This suggests that emulsion fuel droplets on the periphery of a spray may be more prone to puffing than those in the wake region.
The local Nusselt number Nu, a normalised heat flux, is defined as where n denotes the droplet surface-normal direction, the subscript s indicates surface and 0 indicates free-stream air. In order to identify which droplet is more prone to heating prior to puffing, figure 20 shows the instantaneous local Nu along the surface of the droplets D1-D3 in case B1 at t = 9.24 µs when puffing has just started. The flow asymmetry due to puffing is not strong at this stage and, for simplicity, each plot shows the distribution for 0-180 • from the front stagnation point along a droplet surface line that includes the top position on the droplet surface, as shown in figure 20 (see also figure 18). It is clear that for the front droplets (D1 and D2) the Nusselt numbers are higher compared with that for the droplet (D3) in the wake, which is similar to the results in the literature. Therefore, the puffing probability for the front droplets is typically higher due to the higher heat flux at the droplet surface. The different puffing probabilities due to different heating conditions for droplets in different regions of a spray should also be incorporated, in addition to the change of droplet grouping due to puffing discussed above, in modelling of puffing in a biofuel-blended diesel spray.

Concluding remarks
Puffing dynamics and its effects on fuel/air mixing and micromixing for a single droplet and multiple droplets of a decane/bioethanol emulsion mixture under convective heating have been investigated by DNS, which solves all interface dynamics. Under convective heating, the inner-droplet temperature becomes high in the rear region of an emulsion droplet, where initiation of boiling occurs. As the ethanol vapour is ejected into the decane/air gaseous mixture, intense mixing occurs locally and transiently among the vapour of the primary fuel decane, the vapour of the secondary fuel ethanol and air. The mixing characteristics have been investigated using SDRs, which quantify the diffusion processes among the three scalars. The contribution of cross-fuel scalar mixing to the primary-fuel SDR is identified to be predominant due to ethanol vapour puffing. Since the ejected ethanol vapour increases the local mixture fraction, the stoichiometric mixture fraction isosurface is also affected by puffing. When the ethanol ejection is in the wake direction, the wake is elongated but the shape of the stoichiometric surface is not much affected. When the ejection is in the transverse region, the stoichiometric surface shape is considerably changed. The ejection dynamics in the transverse direction is similar to that of a jet in cross-flow. The ethanol fuel-vapour trajectory scales well with the correlation obtained for jets in cross-flow. The transient travelling distance of ethanol fuel-vapour pockets may reach a few diameters away from the parent decane droplet in the transverse direction. When other droplets are present in the neighbourhood of the puffing droplet, this ejection may change the droplet grouping transiently, even if the geometrical locations of the droplets are not changed. This effect is especially 207.241.231.81, on significant for droplets in the upstream of a droplet group on the periphery of a spray due to their direct exposure to the hot air and thus higher heat transfer, which can be quantified by the local Nusselt number. On the other hand, droplets in the wake of other droplets receive less heat and thus puffing is less likely. Even if puffing occurs for these downstream droplets, the ejected ethanol vapour may not change the droplet grouping significantly, since these droplets are already surrounded by fuel-rich gaseous mixture. In modelling puffing effects on an emulsion fuel spray, such a difference in puffing possibility also needs to be considered, in addition to puffing effects on droplet grouping characteristics.