Response of spray number density and evaporation rate to velocity oscillations

Abstract

A theoretical investigation of the effect of gas velocity oscillations on droplet number density and evaporation rate is presented. Oscillations in gas velocity cause a number density wave, i.e. an inhomogeneous, unsteady variation of droplet concentration. The number density wave, as it propagates downstream at the mean flow speed, causes modulation of the local evaporation rate, creating a vapour wave with corresponding oscillations in equivalence ratio. The present work devises an analytical formulation of these processes. Firstly, the response of a population of droplets to oscillations in the gas velocity is modelled in terms of a number density wave. Secondly, the formulation is extended to incorporate droplet evaporation, such that an analytical expression for the evaporation rate modulation is obtained. Subsequently, the droplet 1D convection-diffusion transport equation with the calculated evaporation source term is solved using an appropriate Green’s function to determine the resulting equivalence ratio perturbations. The dynamic response of equivalence ratio fluctuations to velocity oscillations is finally characterized in terms of a frequency-dependent transfer function. The aforementioned analytical approach relies on a number of simplifying approximations, nevertheless it was validated with good agreement against 1D Euler-Lagrange CFD simulations.

Introduction

The economic and environmental requirements for aero-engines continue to push for combustors operating in lean combustion regimes to reduce $N{O}_X$ emission. Due to energy density requirements, liquid fuels are used for aero-engines. Operating liquid fuelled engines in lean conditions may give rise to combustion instabilities.^1–6 Combustion instabilities occur due to the coupling of pressure and velocity oscillations with unsteady heat release rate. In liquid fuelled systems affected by combustion instabilities, the gas velocity oscillations modulate the evaporation rate and thereby, the equivalence ratio. Equivalence ratio fluctuations cause heat release fluctuations,^7,8 which make the overall system prone to combustion instabilities. Investigating spray combustion involves the modelling of several interacting phenomena – such as atomization, evaporation and transport of droplets – that occur over a wide range of spatial and temporal scales. Among the various sub-processes of spray formation and transport, droplet evaporation is particularly important, because it is often rate controlling. Additionally, the unsteady evaporation rate may lead to various burning regimes – from premixed to non-premixed – thus impacting the flame structure. ⁹ As a result, the dynamic response of evaporation is a key constituent of spray combustion dynamics that needs to be modelled accurately.

As droplet evaporation plays an important role in spray combustion, previous works have assessed the impact of acoustics on a single evaporating droplet. Tong et al. ¹⁰ numerically studied the effect of travelling acoustic waves on droplets that are injected at regular intervals into a combustor. It was observed that the overall droplet evaporation rate is enhanced and that the gain of the droplet vaporization response function is sufficiently large to sustain instability in the combustor. In a similar study, Duvvur et al. ¹¹ numerically investigated the effect of acoustic standing waves on the droplet vaporization process with an improved droplet evaporation model. It was seen that only in certain frequency ranges the droplet vaporization process was able to drive combustion instabilities in the combustor and the amplitude of the oscillation had minimal impact on the vaporization rate compared to changes in the fuel volatility. Sujith et al. ¹² examined the effect of an axial acoustic field on the behaviour and evaporation of the droplet. The conclusion from this study was that the droplet evaporation rate increases with the increase of acoustic velocity frequency. Prud’Homme et al. ¹³ conducted a theoretical investigation of the droplet dynamic vaporization response to acoustic oscillations with a droplet evaporation model accounting for finite thermal conductivity inside the liquid phase of the droplet. Enhanced evaporation rate due to the application of acoustic waves was observed. It was concluded that the internal thermal exchange inside the droplet is an important factor in modelling the droplet evaporation rate in the context of combustion instabilities.

Adequate studies have been carried out on the response of a single droplet to acoustic oscillations. However, in realistic situations we are more concerned with response of a population of droplets. The following discussion reviews some works concerned with response of a spray to acoustic oscillations. The effect of acoustic field on an ethanol spray flame was experimentally studied by Dubey et al. ¹⁴ . The Sauter mean and arithmetic mean diameters of the spray were seen to decrease in the presence of an acoustic field by 15% and 20%, respectively, due to enhanced evaporation rate of the spray. Gurubaran et al. ¹⁵ observed droplet clustering and a variation of droplet mean diameter downstream of the injector and concluded that the acoustic pressure amplitude has a significant effect on the particle size distribution. Gajan et al. ¹⁶ experimentally investigated the behaviour of a spray downstream of an aero-engine injector submitted to acoustic excitation and noticed the formation of a spray droplet density wave. With numerical analysis, the droplet density wave was shown to have two origins: due to the atomization process and due to the effect of oscillating flow on the transport of droplets. Apeloig et al. ¹⁷ experimentally studied the unsteady interaction of a kerosene spray downstream from a multi-point injector with the flame, where the liquid jet was injected into the air as a crossflow. During an instability cycle, different droplet concentration were seen, leading to heat release rate oscillations. In a recent, similar experimental study by Bodoc et al. ¹⁸ the formation of an alternating dense and diluted zones of two-phase flow was observed downstream of the multipoint jet-in-crossflow injector.

On the one hand, experimental and numerical studies^19–22,2 of the response of a single droplet or a population of droplets to an acoustic field have been carried out. Nevertheless, general conclusions from these studies cannot be drawn unless an elaborate parametric analysis is performed. Such parametric analysis in experimental campaigns and numerical studies might be expensive. On the other hand, theoretical study on the response of a population of droplets to an acoustic field may lead to a simple mathematical formulation and facilitate a rich and inexpensive parametric analysis of the number density wave and subsequent modulated evaporation rate. In the present work, we aim to devise an analytical formulation of the response of a population of droplets with and without evaporation to an acoustic field. The goal of formulating an analytical description of the response of population of droplets to acoustic oscillation comprises of three objectives: (1) analytical formulation of number density wave without evaporation, (2) analytical formulation of the evaporation spray response, and finally (3) description of the evaporation spray dynamics using a transfer function.

Previous studies have observed droplet clustering ¹⁵ and the formation of a droplet density wave^23,16 downstream of the injector through experiments and simulations. Achury and Polifke ²⁴ proposed a theoretical approach to study the response of a single droplet to acoustic oscillations. The current work extends this work from a single droplet to a population of droplets by proposing an analytical formulation of the Number Density (ND) wave. Such a ND wave results from modulating the gas velocity with continuous injection of mono-disperse droplets.

Prior studies^11,10,25 have dealt with the response of droplet(s) to large acoustic pressure oscillations leading to an oscillatory response of the evaporation. However, in those studies the implicit role of velocity fluctuations was not taken into account. Velocity oscillations cause an ND wave, where droplets exhibit a heterogeneous distribution, which leads to a spatio-temporal oscillation of evaporation rate and the formation of a fuel vapour wave. The vapour wave is convected downstream at the mean flow speed causing equivalence ratio fluctuations. A theoretical development, based on the concept of solving the droplet motion equation by perturbation theory, was used to incorporate droplet evaporation to the spray number density wave by Katoshevski et al. ²⁶ . In the present work, instead, we extend the analytical formulation of the ND wave by incorporating the droplet evaporation. Furthermore, we model the equivalence ratio fluctuations caused by the propagation of the vapour wave by analytically solving the 1D convection-diffusion equation by means of a Green’s function. The oscillatory evaporation response is then characterized by an analytical transfer function of the evaporation. The results of the analytical approach are validated against 1D Euler-Lagrange CFD simulations.

Modelling the spray as a collection of sub-systems comprising of a given droplet ND wave, evaporation and fuel vapour transport facilitates the development of a so-called network model ²⁷ for the response of a flame to acoustic wave. A Network model is a low-order representation of all the involved sub-systems characterised by blocks with associated inputs and outputs, which subsequently can be interconnected forming a global system (investigated spray) that can be described analytically by a transfer function. Such a method of studying acoustic response via network modelling may be useful when studying realistic sprays with increased complexity.

In the next section, an analytical formulation of the number density wave is deduced from the droplet motion equation. This is followed by the Evaporative Spray Response section, which shows the formation of evaporation rate fluctuations caused by the ND wave. Finally, the Equivalence Ratio Fluctuations section describes the formation of equivalence ratio fluctuations and subsequent characterisation of the evaporation response in terms of an analytical transfer function.

Droplet Dynamics

In the present study, we aim at analytically determining the response of a population of droplets to velocity oscillations that typically occur in liquid-fuelled combustors. In pressure-swirled atomizers used in gas turbine combustors, due to the hydrodynamic instabilities involved in the breakup process, primary atomization occurs. In this region, the spray is quite dense and its velocity and temperature do not substantially change from their injection values. ²⁸ Further break-up of the droplets occurs downstream due to the interaction with the turbulent gaseous flow in the secondary atomization process. Droplet heating and vaporization occurs only downstream from the primary atomization region, once the droplet dispersion causes the spray to become dilute to allow for the droplets to receive thermal energy from the surrounding gas. As a consequence, atomization, spray formation, vaporization and combustion can be studied independently. In the following discussion, we focus on the dynamics of the droplets’ motion and simultaneous vaporization behaviour due to fluctuating gas velocity occurring in the dilute downstream region where the liquid volume fraction is small. Due to small liquid volume fraction, the droplet-droplet interaction can be neglected and the modelling approach can be carried out with a “one-way” coupling approach as shown in the work of Elghobashi. ²⁹ As seen in the previous studies,^16,21 due to the oscillation of the gaseous field, droplets downstream of the injector concentrate inhomogeneously in space and time, thus giving rise to a droplet number density wave. Such a droplet number density can be analysed by studying the response of a population of droplets to acoustic forcing. In CFD, as well as theoretical approaches, the Lagrangian or mass-point perspective is well-suited to model dilute spray problems, as it is computationally inexpensive and easier to model than approaches that resolve droplet/fluid interfaces. In the Lagrangian framework, the droplet acceleration is balanced by the forces acting on it, which need to be modelled for accurate droplet dynamics. The Lagrangian equation of motion for a spherical droplet with mass $m_d$ and velocity ${\mathbf{u}}_d$ is given by:

$$m_d\frac{\mathrm{d}{\mathbf{u}}_d}{\mathrm{d}t}={\mathbf{F}}_B+{\mathbf{F}}_S,$$

(1)

where ${\mathbf{F}}_B$ refers to body forces, such as gravity acting on the droplet, and ${\mathbf{F}}_S$ refers to surface forces, such as pressure gradient and shear stress. Maxey and Riley ³⁰ presented an equation of motion for a small rigid particle of diameter $D$ immersed in a non-uniform flow field that resolves all the forces acting on the particle. In their hypothesis, forces can be split broadly into two categories: undisturbed and disturbed flow forces. The former represent the surface forces required to accelerate the fluid that would occupy the volume of the particle if the particle were absent. ³¹ Maxey and Riley further classified undisturbed flow forces into pressure gradient force and shear stress. The latter ${\mathbf{F}}_{\mathbf{u}\mathbf{f}}$ can be written as:

$${\mathbf{F}}_{\mathbf{u}\mathbf{f}}=V_d(\underset{\text{,press. grad. }}{\underbrace{-\mathrm{\nabla }p}}+\underset{\text{shear stress }}{\underbrace{{\nu }_c{\mathrm{\nabla }}^2{\mathbf{u}}_c}})=\frac{1}{\gamma }(\frac{D{\mathbf{u}}_c}{Dt}-\mathbf{g}),$$

(2)

where $\mathbf{g}$ corresponds to buoyancy force, $\gamma ={\rho }_d/{\rho }_c$ represents the density ratio of the liquid to carrier fluid, ${\nu }_c$ represents the kinematic viscosity of the carrier fluid and $V_d$ is the velocity of the droplet. The disturbed flow forces ${\mathbf{F}}_{\mathbf{d}\mathbf{f}}$ correspond to the forces exerted due to the presence of the droplet. The disturbed flow forces can be sub-divided into steady-state drag, virtual mass force and Basset or history integral forces. The complete Maxey-Riley equation including all the surface forces is given as:

$$\begin{array}{cc}\frac{\mathrm{d}{\mathbf{u}}_d}{\mathrm{d}t}=& \underset{\text{body, : , gravity }}{\underbrace{\mathbf{g}}}+\underset{\text{undisturbed flow }}{\underbrace{\frac{1}{\gamma }(\frac{D{\mathbf{u}}_c}{Dt}-\mathbf{g})}}\\ & +\underset{\text{steady state drag }}{\underbrace{\frac{3C_D}{4\gamma D}|{\mathbf{u}}_c-{\mathbf{u}}_d|({\mathbf{u}}_c-{\mathbf{u}}_d)}}\\ & +\underset{\text{virtual mass }}{\underbrace{\frac{{\mathrm{\Delta }}_A}{2\gamma }(\frac{D{\mathbf{u}}_c}{Dt}-\frac{\mathrm{d}{\mathbf{u}}_d}{\mathrm{d}t})}}+{\mathbf{F}}_{\text{Basset }}\end{array}$$

(3)

where ${\mathrm{\Delta }}_A$ is an empirical constant used to extend the approach to high Reynolds number flows. Added or virtual mass effects become important only when the particle/droplet density is comparable to the fluid density, for example when the ‘particle’ is a gas bubble in a liquid, and therefore can be neglected in this investigation. Similarly, Basset forces, which represent history acceleration effects, become insignificant when the density of the particle is very large compared to the gaseous phase density. ³² In the current investigation, the particle density (water) is higher than the density of the gaseous medium (air) and therefore Basset forces can be neglected. The drag force depends on the coefficient of drag, $C_D=24f_1/R{e}_d$ , where $R{e}_d$ is the droplet Reynolds number defined as $R{e}_d=|{\mathbf{u}}_d-{\mathbf{u}}_c|D/{\nu }_c$ . A non-linearity is introduced by the drag force ³³ if the flow is outside the Stokes regime, then $f_1$ becomes a function of $R{e}_d$ as shown below:

$$f_1=\{\begin{array}{cc}1& \to \text{ Stokes ;flow},R{e}_d<1\\ 1+0.15R{e}_d^{0.687}& \to \text{ Schiller \& Naumann},R{e}_d<800\end{array}$$

(4)

In a 1D flow, with ${\ddot{x}}_d={\dot{u}}_d$ and ${\dot{x}}_d=u_d$ and particle relaxation time ${\tau }_d=\gamma D^2/18{\nu }_c$ , the equation of motion for particle position $x_d$ in the Stokes flow regime ( $f_1=1$ ) after neglecting the gravity, Basset and virtual mass forces reduces to:

$${\ddot{x}}_d+\frac{f_1{\dot{x}}_d}{{\tau }_d}=\frac{f_1{u}_c}{{\tau }_d}+\frac{1}{\gamma }\frac{\mathrm{\partial }{u}_c}{\mathrm{\partial }t}$$

(5)

Droplet Population Response

In this section, an analytical description of the response of a population of droplets to velocity oscillations is introduced. Consider droplets with velocity $u_d$ being continuously injected into the one dimensional domain with carrier gas velocity $u_c$ (Eq. (6)) resembling droplets being sheared from an atomizer.

$$u_c={\overline{u}}_c+{\hat{u}}_c\sin (\omega t+\varphi ),$$

(6)

where ${\overline{u}}_c$ is the mean velocity, ${\hat{u}}_c$ is the fluctuating component of the gas velocity, $\omega =2\pi f$ is the angular frequency and $\varphi$ is the phase angle of the flow oscillation, which is included to generalize the solution. In this work, droplets are injected into a pulsating flow (oscillation superposed on mean flow) of an incompressible medium at the mean gas velocity, corresponding to zero slip velocity with respect to mean gas velocity ( $u_d={\overline{u}}_c$ ). Due to the application of pulsating flow, the domain under investigation is assumed to be acoustically compact, i.e. the droplet-acoustic interactions occur over a length considerably smaller than the acoustic wavelength. Therefore, acoustic velocity and displacement are homogeneous in the domain.

Due to the acoustic compactness of the problem, a non-dimensional equation can be derived from Eq. (5) by normalizing time with the oscillation time period, $\tilde{t}=t/{\tau }_d$ and velocity by the oscillation amplitude, ${\tilde{u}}_d=u_d/{\hat{u}}_c$ . Under the dilute flow regime the governing equation can be written with relative amplitude of oscillation, $\epsilon ={\hat{u}}_c/{\overline{u}}_c$ , non-dimensional frequency $\tilde{\omega }=\omega \tilde{t}$ along with other dimensionless parameters which control the response to velocity oscillation ¹² :

$$\frac{\mathrm{d}{\tilde{u}}_d}{\mathrm{d}\tilde{t}}=-\frac{f_1}{2\pi }({\tilde{u}}_d-\sin (\tilde{\omega }\tilde{t})-\frac{1}{\epsilon })+\frac{\tilde{\omega }}{\gamma }\cos (\tilde{\omega }\tilde{t})$$

(7)

An analytical solution of Eq (7) is possible as the droplets are injected at the fluid mean flow velocity ( $R{e}_d\sim 0$ , $f_1=1$ ) and due to the application of oscillating flow (infinite acoustic wavelength). If the Schiller-Naumann extension for the drag law is accounted for, an analytical solution is not possible and is resolved numerically as discussed in Achury et al. ²⁴ .

In an acoustically compact domain, the number density wave, as reported previously in the works of Gurubaran et al. ¹⁵ and Gajan et al. ¹⁶ , results from the interaction of the pulsating flow with the droplets at the injection time. When droplets are continuously injected into a pulsating flow, droplets at an asymptotic state (after a long time) convect at mean flow speed ${\overline{u}}_c$ and oscillate around that for a given frequency, amplitude of oscillation, droplet size and fluid and gas properties. When droplets are injected into the flow domain, each droplet experiences a different gas velocity corresponding to the relative phase of the fluctuating gas velocity. For example, in Figure 1 a droplet injected at instance $t_2$ (blue) sees a higher gas velocity, whereas a droplet injected at $t_1$ (red) sees a lower gas velocity. As droplets experience different gas velocity when injected, the trajectory of each droplet is different as shown by the dotted lines in Figure 2. Differences in droplet trajectories due to different injection times gives rise to differences in spacing of the droplets in the domain at any given observation time. This inhomogeneous distribution of droplets in space, as shown by black crosses in Figure 2, determines the Number Density wave.

Figure 1.

Droplet injection times $t_i$ (i=0,1,…6) marked by vertical coloured solid lines on the time axis

Figure 2.

Individual trajectories of six droplets injected at $t_i$ (c.f. Figure 1) according to Eq. (8) . Black crosses show droplet positions at time $\tilde{t}$ = 280. Dimensional parameters for this example: Droplet size 35  $\mu$ m, (Stokes number = 0.8), mean flow velocity ${\overline{u}}_c=5$   m/s, gas phase oscillation amplitude ${\hat{u}}_c=2.5$  m/s and frequency $f=250$  Hz (T = 1/f = 0.004 s)

In this work, we seek to obtain a closed form solution for the proposed number density wave. This is achieved by first assessing an individual droplet trajectory and then analysing the evolution of an ensemble of droplet positions with observation time to determine the ND wave function. In 1-D, to assess the response of a droplet to a given fluid excitation, one may solve Eq. (3) with flow excitation $u_c(x,t)$ . Using the ansatz for $u_c(x,t)$ in the droplet motion equation, the analytical solution of the droplet positions $x_d$ in the Stokes flow regime ( $f_1=1$ ) can be obtained by integrating over the lifetime $t$ of the particle. To identify relevant parameters that govern the droplet trajectories, the analytical solution is represented in a non-dimensional form with the droplet relaxation time ${\tau }_d$ and mean velocity ${\overline{u}}_c$ as characteristic time and velocity scales.

$$\begin{array}{rl}& {\tilde{x}}_d(\tilde{t},{\tilde{t}}_i)={\tilde{C}}_1+{\tilde{C}}_2{e}^{-(\tilde{t}-\widetilde{t_i})}+\tilde{t}\\ & +{\tilde{C}}_3[\frac{\gamma -1}{\gamma }\sin (\tilde{\omega }\tilde{t}+\varphi )-(\frac{1}{\tilde{\omega }}+\frac{\tilde{\omega }}{\gamma })\cos (\tilde{\omega }\tilde{t}+\varphi )]\end{array}$$

(8)

The coefficients ${\tilde{C}}_1$ and ${\tilde{C}}_2$ calculated using the initial condition $x_d(0)=0$ and ${\dot{x}}_d(0)=u_{d0}={\overline{u}}_c$ are also shown in non-dimensional form as:

$$\frac{{\tilde{C}}_3={\hat{u}}_c/{\overline{u}}_c}{1+{\tilde{\omega }}^2=\frac{\epsilon }{1+{\tilde{\omega }}^2}}$$

(9)

$$\begin{array}{rl}{\tilde{C}}_2=& \frac{{\overline{u}}_c-u_{d0}}{{\overline{u}}_c}\\ & +{\tilde{C}}_3[\frac{1-\gamma }{\gamma }\cos (\tilde{\omega }{\tilde{t}}_i-\varphi )-(\frac{1}{\tilde{\omega }}+\frac{\tilde{\omega }}{\gamma })\sin (\tilde{\omega }\widetilde{t_i}-\varphi )]\end{array}$$

(10)

$$\begin{array}{rl}{\tilde{C}}_1=& -{\tilde{C}}_2{e}^{-\widetilde{t_i}}\\ & +{\tilde{C}}_3[(\frac{1}{\tilde{\omega }}+\frac{\tilde{\omega }}{\gamma })\cos (\tilde{\omega }\widetilde{t_i}-\varphi )+\frac{1-\gamma }{\gamma }\sin (\tilde{\omega }\widetilde{t_i}-\varphi )]\end{array}$$

(11)

The first term of the non-dimensional droplet position equation (Eq. (8)) is a constant, the second term describes an exponential decay where the droplets “loose” the influence of the initial condition as time progresses, the third term represents the mean convection of droplets which increases linearly with time and the fourth term describes the modulation of the droplet positions due to flow oscillation. The trajectory of the individual droplets injected at different injection times: $t_i$ = $0,T/5,2T/5,3T/5,4T/5,T$ of the $250Hz$ gas oscillation signal $u_c$ , and for a given phase angle $\varphi =0$ , is shown in Figure 2 by dotted lines at an observation non-dimensional time of 280. It is to be noted that, a high excitation amplitude (50%) and a low convection velocity, which is not typical of values seen in gas turbines, are used to distinctly visualize individual droplet trajectories in Figure 2 which otherwise would not visible. However, higher velocity can be considered in the current analytical framework. Differences in the droplet trajectories due to different injection times can also be seen in the zoomed in view (inset) in Figure 8. Since the gas phase excitation is periodic, a periodic behaviour of the droplet trajectories is also reached asymptotically where they convect at mean speed and fluctuate at the oscillation frequency. This asymptotic behaviour can also be seen by zeroing the exponential term in Eq. (8), which results in an expression with terms involving convection at ${\overline{u}}_c$ and oscillation ( ${\tilde{C}}_1$ and ${\tilde{C}}_3$ without ${\tilde{C}}_2$ ), as given in Eq. (12):

$$\begin{array}{cc}{\tilde{x}}_d(\tilde{t},{\tilde{t}}_i)=& \tilde{t}+{\tilde{C}}_3[\frac{1-\gamma }{\gamma }(\sin (B)-\sin (A))\\ & -(\frac{1}{\tilde{\omega }}+\frac{\tilde{\omega }}{\gamma })(\cos (B)-\cos (A))],\end{array}$$

(12)

where $A=\tilde{\omega }\tilde{t}+\varphi$ and $B=-\tilde{\omega }\widetilde{t_i}+\varphi$ .

Definition of the Number Density Wave

Figure 2 shows that individual droplet trajectories differ from each other leading to different inter-droplet spacing at any given observation time, as shown by black crosses. Such a heterogeneous distribution of droplets in space results in a number density wave. The evolution of the number density wave can be obtained by tracking the change in spacing of droplets as time progresses. The change in spacing of droplets is inversely dependent on injection time $t_i$ . Thus, the number density wave can be mathematically written as the inverse of time derivative $(d{x}_d/d{t}_i{)}^{-1}$ . In order to obtain the inverse time derivative, we first differentiate the dimensional form of Eq. (8) to get $d{x}_d/d{t}_i$ (Eq. (13)) and then take the inverse to get the expression for number density wave ( $\varrho (x_d(t_i),\varphi )$ ):

$$\begin{array}{rl}& {\displaystyle \frac{d{\tilde{x}}_d}{d{\tilde{t}}_i}=1-\left({\displaystyle \frac{{\overline{u}}_c-u_{d0}}{{\overline{u}}_c}}\right)e^{-{\tilde{t}}_i}+\epsilon \left(1-e^{-\tilde{t}}\right)}\\ & \left({\displaystyle \frac{\tilde{\omega }}{\gamma }\cos \left(-\tilde{\omega }{\tilde{t}}_i+\varphi \right)+\sin \left(-\tilde{\omega }{\tilde{t}}_i+\varphi \right)}\right)\end{array}$$

(13)

$$\begin{array}{rl}\tilde{\varrho }({\tilde{x}}_d,\varphi )& =\dot{n}{\tau }_d[1-\left(\frac{{\overline{u}}_c-u_{d0}}{{\overline{u}}_c}\right)e^{-{\tilde{t}}_i}\\ & +\epsilon (\frac{\tilde{\omega }}{\gamma }\cos (\tilde{\omega }{\tilde{t}}_i-\varphi )-\sin (\tilde{\omega }{\tilde{t}}_i-\varphi ))(1-e^{-{\tilde{t}}_i}){]}^{-1}\end{array}$$

(14)

where, $\dot{n}$ is the rate of droplet injection and $\varphi$ is the phase angle of the gas oscillation corresponding to the time when the droplet is injected into the domain. The number density wave fundamentally depends on $x_d$ and time of injection $t_i$ . It also depends on time $t$ through the inverse function $t_i$ and $x_d$ . Equation. (14) is not yet useful as $t_i$ is not known explicitly. By inspection of Eq. (8), it can be recognized that $x_d(t_i)$ is composed by the superposition of the trajectories due to the relaxation of the droplets to the mean flow and the oscillation induced by ${\hat{u}}_c$ when ${\hat{u}}_c>0$ . An approximation for the inverse function to obtain $t_i$ from Eq. (8) assuming small oscillation amplitudes is given by:

$${\tilde{t}}_i\left(x_d\right)={\tilde{x}}_d+\frac{{\overline{u}}_c-u_{d0}}{{\overline{u}}_c}+W(z),$$

(15)

where $W(z)$ is a Lambert W function, which can be expressed as a power series according to the Lagrangian inverse theorem in the following way:

$$\begin{array}{cc}& W(z)={\displaystyle \sum _{n=1}^{\mathrm{\infty }}\frac{(-n{)}^{n-1}}{n!}{z}^n,\text{ with }}\\ & z\left({\tilde{x}}_d\right)=e^{-\frac{1}{{\overline{u}}_c}({\overline{u}}_c-u_{d0}+{\tilde{x}}_d)}\left(\frac{u_{d0-{\overline{u}}_c}}{{\overline{u}}_c}\right).\end{array}$$

(16)

As $t_i(x_d)$ is obtained, a closed form expression for the Number Density is obtained by combining Eqs. (14), (15) and (16). Before proceeding further, the analytical formulation of the obtained ND wave expression is validated against a 1D CFD simulation carried out in OpenFOAM.

Validation of the ND wave model

In this section, the number density wave is validated against a 1D CFD simulation with Lagrangian particle tracking scheme to obtain droplet statistics. Droplet position and liquid mass flow data from the simulation are post-processed to validate the number density wave and evaporation rate modulation, which will be discussed in the next section.

The CFD setup consists of a 1D duct of length 0.3 m in the $x$ direction and is discretized with 500 grid points. Water is used as the dispersed phase and air is used as the carrier phase. As the focus of the work is to demonstrate an analytical framework to capture the evaporation response of droplets to velocity oscillations, water is chosen as the material for dispersed phase, because compared to other fuels there is little ambiguity concerning the available thermophysical properties. The current framework can be used with typical fuels such as heptane, kerosene by considering respective thermophysical properties without altering the analytical framework. Mono-disperse droplets are injected continuously at the inlet of the domain ( $x=0$  m) as shown in the Figure 3, which are tracked using Lagrangian particle tracking method. To assess the response of a spray to velocity oscillations, the gas velocity at the inlet is modulated at a given frequency and amplitude. Other boundary conditions used for the simulation and analytical procedure are given in Tab. 1.

Figure 3.

Schematic of the CFD domain with droplet injection

Table 1.

Boundary conditions for the CFD simulations

Diameter, $\mu$ m	${\overline{u}}_c$ (m/s)	${\hat{u}}_c$ (m/s)	$u_{d0}$ (m/s)	$T_{gas}$ (K)	$T_{d0}$ (K)
30	5	2.5	5	700	300

The transient simulation with time step of $\mathrm{\Delta }t=5\times {10}^{-4}$ of the laminar flow is carried out for 0.1 s. The Lagrangian data in terms of velocity, position in the domain at any given observation time, and number of droplets is recorded at each time step. The number density is obtained by using the histogram and the position information of the droplets. The normalized ND wave ( $\overline{\varrho }=\varrho /\varrho ({\hat{u}}_c=0)$ ) calculated from the closed expression (Eq. 14) is in good agreement with the CFD simulation as shown in Figure 4.

Figure 4.

Comparison of analytical (—) and CFD (- - -) Normalized Number Density profile

Evaporative Spray Response

This section focuses on the characterisation of the evaporative spray response to a pulsating flow, where the previously developed theory regarding the formation of a number density wave for non-evaporating droplets is recalled and droplet evaporation is incorporated. This can be simply seen as the progressive reduction of the droplet size, along with the droplet convection. To include the effect of evaporation, we consider the $D^2$ model for evaporation of droplets injected with zero slip velocity with respect to mean velocity. We show later that for moderate amplitudes of excitation, the $D^2$ law can recover quite accurately the evaporation response as obtained with an advanced evaporation model employed in CFD. For a spherical droplet the rate of evaporation from $D^2$ law is expressed as:

$${\dot{m}}_d=2\pi D{\rho }_g{\mathcal{D}}_{\mathcal{A}\mathcal{B}}\ln (1+B_M)$$

(17)

which yields the following equation for the droplet diameter evolution:

$$\frac{d{D}^2}{\mathrm{d}t}=-\frac{8{\rho }_g{\mathcal{D}}_{\mathcal{A}\mathcal{B}}}{{\rho }_d}\ln (1+B_M)$$

(18)

where ${\rho }_g$ is the average gas density evaluated at film temperature and composition obtained using the $1/3^{rd}$ rule, ³⁴ ${\rho }_l$ is the density of the liquid, ${\mathcal{D}}_{\mathcal{A}\mathcal{B}}$ is the binary diffusion coefficient and $B_M$ is the Spalding mass transfer number. The droplet lifetime can then be calculated as:

$$t_d=\frac{D_0^2}{K}$$

(19)

where $D_0$ is the initial droplet size and $K$ is the evaporation constant. As a first approach, the simple aforementioned $D^2$ law can be used in the number density wave to incorporate size variation. In this case the droplet size evolution can be resolved for initial injection size $D_0$ and constant evaporation coefficient $K$ (Eq. (20)) which does not depend on the droplet Reynolds number. Please note that a constant evaporation coefficient can be justified in this case as the droplet is injected close to the wet-bulb temperature of water for the given gas phase temperature and pressure. Thus droplet heating time is negligible and the evaporation coefficient can be considered constant in time and can be deduced analytically. However, future work will consider extension of the theoretical framework to include non-negligible droplet heat up time analytically based on the liquid injection temperature, liquid material and gas phase temperature and pressure.

$$K=\frac{8{\rho }_g{\mathcal{D}}_{\mathcal{A}\mathcal{B}}}{{\rho }_d}\ln (1+B_M)$$

(20)

The droplet size evolution can be further used in the Number Density equation (Eq. (14)) via the droplet relaxation time parameter ${\tau }_d$ . The function $x_d(t,t_i)$ gives the position of the droplets at a given observation time $t$ , that were injected at injection time $t_i$ . In order to obtain the droplet size evolution for all the droplets injected in the domain, the $D^2$ law for size evolution can be written as:

$$D(t_i)=\sqrt{D_0^2-K{t}_i}$$

(21)

After having obtained the droplet size evolution equation for all the droplets, the droplet size evolution function is introduced in the ND wave Eq. (14) via parameter ${\tau }_d(D(t_i))$ which recursively depends on $t_i$ .

The number density wave equation (Eq. (14)) contains an implicit function in ${\tau }_d(t_i)$ , which is solved using a predictor-corrector scheme. In the first step, the initial droplet size $D_0$ is used to determine the ${\tau }_d(D_0)$ as done in the non-evaporative case. ${\tau }_d(D_0)$ is then used to obtain the $t_i$ values. In the next step, the injection time $t_i$ estimated in the first step is substituted in Eq. (21) to give the droplet size evolution.

Eq. (14) for the number density wave was derived for a constant droplet relaxation time ${\tau }_d$ . Inclusion of a variable ${\tau }_d(D(t_i))$ in the derivation of the number density equation (14) results in a problem that appears analytically intractable. Fortunately, given that the time scales of droplet relaxation and oscillation are smaller than the droplet life time $t_d$ , a quasi-steady approximation for the oscillatory evaporation appears justified, which is obtained by simply inserting the time-varying ${\tau }_d$ in Eq. (14). The fact that results from the analytical approach are in good agreement with CFD (see below) provides a posteriori validation of the quasi-steady ansatz. The limits of the validity of this assumption will be explored in future work.

Finally, the parameter ${\tau }_d(D(t_i))$ is updated and then used in the evaluation of the ND wave (Eq. (14)). The profile of the evaporation rate (Eq. (17)) is extended with ND wave Eq. (14) to give the evaporative spray response:

$$\dot{m}(x_d(t_i))=\underset{\text{individual droplet }}{\underbrace{\frac{\pi }{4}{\rho }_d\stackrel{D}{\overbrace{\sqrt{D_0^2-K{t}_i\left(x_d\right)}}}\stackrel{d{D}^2/dt}{\overbrace{K}}}}\underset{\text{number density }}{\underbrace{\varrho (x_d(t_i),\varphi )}}$$

(22)

The constant evaporation coefficient $K$ for the analytical expression in Eq. (22) is obtained by the $D^2$ law for $R{e}_d=0$ . In the Lagrangian CFD solver, the film theory of Abramzon and Sirignano ³⁵ is employed for the calculation of the evaporation of droplets. The evaporation rate in the CFD case is estimated as the local sum of the evaporated liquid which is added to the gas-phase system in the given time interval. This evaporation rate source term is extracted using a user-defined routine in OpenFOAM. The oscillating component of the gas velocity ${\hat{u}}_c>0$ , which elicits the formation of a ND wave, produces a local fluctuation of the droplet rate of evaporation. Droplet grouping caused by the ND wave in space and coupled with simultaneous evaporation of droplets gives rise to the spatio-temporal modulation of the evaporation rate. The evaporation rate in non-dimensional form ( ${\tilde{m}}_d={\dot{m}}_d/{\rho }_d{A}_{domain}{\overline{u}}_c$ ) with cross-sectional area $A_{domain}$ , is shown in Figure 5.

Figure 5.

Comparison of analytical (—) and CFD (- - -) evaporation rate of droplets, D = 30  $\mu \mathrm{m}$ , f = 250  Hz, K = 7.7e-8  ${\mathrm{m}}^2/\mathrm{s}$ , $T_{air}=700$   K, $T_{p0}=300$   K.

Figure. 5 compares the evaporation rate profile according to Eq. (22) with the CFD simulation. A good agreement is seen, except for the fact that the CFD case exhibits a gradual increase in oscillations at the start of the domain $\sim x$ = 0, which is not reproduced by the analytical result. This is plausible due to the absence of droplet heating process in the analytical formulation. Another important takeaway message from Figure 5 is that the $D^2$ law used for the analytical prediction recovers the group evaporation response obtained with the more realistic Abramzon and Sirignano evaporation model. In summary, the timing of the injection of droplets and fluctuation velocity lead to the formation of a number density wave. As a result of the ND wave and simultaneous evaporation of droplets, the evaporation rate $\dot{m}$ is also modulated.

Equivalence Ratio Fluctuations

In this section, we seek an analytical solution for the evaporation response of a population of droplets. We have seen from the previous section that due to the imposition of acoustic excitation, an inhomogeneous distribution of droplets occurs downstream of the injector to form a ND wave. As a result of the ND wave the rate of evaporation oscillates giving rise to the fuel vapour wave, which propagates downstream at mean convection speed.^36,37 Oscillations in fuel vapour production lead to oscillations in the equivalence ratio. The fact that air velocity oscillations lead to oscillations in the equivalence ratio allows to construct a transfer function characterizing the evaporation response to velocity oscillations. Such a transfer function block for the evaporation process coupled with a transfer function for other processes, such as atomization and heat release rate dynamics, completes the low-order model of a spray flame for the study of combustion instabilities.

The evaporated mass flow feeds the gaseous phase with fuel vapour. Thus, the expression Eq. (22) for the evaporated mass flow rate obtained in the previous section can be utilized as a source term in a transport equation for gaseous fuel to determine the mass fraction modulation of fuel vapour:

$$\rho \frac{\mathrm{\partial }Y}{\mathrm{\partial }t}+\rho {\overline{u}}_c\frac{\mathrm{\partial }Y}{\mathrm{\partial }x}-\rho \mathcal{D}\frac{{\mathrm{\partial }}^{2}Y}{\mathrm{\partial }{x}^2}=\dot{m}.$$

(23)

Essentially, Eq. (23) captures the effect of acoustic modulation via the source term, while the transport of mass fraction ignores the fluctuations in the convective-diffusive transport. This is justified as the evaporated vapour convects at mean flow speed which allows to evaluate Eq. (23) using only the mean flow speed without oscillation. It is to be noted that the density $\rho$ is out of the differential operators due to the combination of the conservative form of the transport equation and continuity equation, along with the assumption that density fluctuations are negligible due to compact acoustics and small vapour loading. This assumption is justified a posteriori, as the results from the analytical approach matches well with the CFD, where density variations are taken into account. The limit of the validity of such an assumption will be investigated in the future work.

We solve the 1D transport equation via the Green’s function method for convection-diffusion problem with a source term.

$$Y(x,t)={\int }_0^t{\int }_0^{\xi }G(\xi ,\tau ;x,t)\dot{m}(\xi ,\tau )\mathrm{d}\xi \mathrm{d}\tau$$

(24)

Now, the problem turns to find a Green’s function form that can solve Eq. (23). For the given 1D transport equation, where the source term $\dot{m}$ is set to act in the domain $x$ $\in$ $0-\xi$ and during time $0-\tau$ , the Green’s function ³⁸ with mean flow can be derived as:

$$G(\xi ,\tau ;x,t)=\frac{H(t-\tau )}{\sqrt{4\pi \mathcal{D}(t-\tau )}}exp[-\frac{[\xi -x+{\overline{u}}_c(t-\tau ){]}^2}{4\mathcal{D}(t-\tau )}],$$

(25)

where $\mathcal{D}$ is the diffusion coefficient and $H$ is the Heaviside function. It is to be noted here that, first, the solution $Y$ is evaluated at time $t$ , which means the source (evaporation of droplets) is active in $0-\xi$ until time $\tau =t$ . The diffusion coefficient $\chi$ is set equal to the value used in the OpenFOAM CFD simulation to enable direct comparison of the results.

The mass fraction result obtained from solving Eq. (23) is validated against the mass fraction result obtained from CFD by post-processing the mass fraction in the flow domain. Figure 6 shows the mass fraction obtained from CFD and the analytical Green’s function method with only mean flow at a given time instance. The result of the analytical solution in terms of the oscillation amplitude and the diffusion of fuel vapour is in good agreement with the CFD results. The discrepancies between the analytical solution and CFD result could be due to the absence of an oscillating component in the Green’s function solution of the transport equation Eq. (25).

Figure 6.

Comparison of analytical (—) and CFD (- - -) mass fraction profiles for given boundary conditions: $D_0$ = 30  $\mu \mathrm{m}$ , f = 250 Hz, ${\overline{u}}_c=5$  m/s, $u^{\prime }=2.5$   m/s, $T_{air}=700$   K, $T_{p0}=300$   K

The dynamic response of the equivalence ratio to air velocity oscillations can be studied by a evaporation transfer function as shown below:

$$F(\omega )=\frac{Y_{H_{2}O}^{\prime }/{\overline{Y}}_{H_{2}O}}{u_c^{\mathrm{\prime }}/{\overline{u}}_c}=G(\omega )e^{i\phi (\omega )}$$

(26)

where $G(\omega )$ and $\phi (\omega )$ represent the gain and phase delay of the evaporation response. Such a transfer function enables one to characterize the response of the evaporation process across the relevant frequency spectrum.

In this study, the gain of the transfer function for both analytical approach and CFD simulations is calculated by taking the ratio of the normalized mass fraction fluctuation and air velocity fluctuation time series signals imposed at the reference inlet position. The phase is calculated by taking the phase difference between mass fraction fluctuation and air velocity fluctuation time series signals at the reference location.

The gain and phase are calculated for the boundary conditions given in Tab. 1 at different Strouhal number (Sr = $fL/{\overline{u}}_c$ ) where $L$ is the axial length of the domain as shown in Figure 8. The evaporation process describes a low-pass behaviour meaning that the droplet evaporation process becomes unresponsive to rapid changes in the gas velocity similar to the observations made by Chaussonnet et al. ³⁹ for the atomization process. Such an observation needs to be rigorously studied for the realistic poly-disperse spray, which can exhibit different liquid evaporation time scales that will affect their frequency response. An estimate for the time delay of the transport process of vapour can be given by $\tilde{t}=\tilde{x}=11.2$ , which is almost equal to the value obtained by measuring the peak to peak signal from Figure 7. Furthermore, the evaporation process introduces an additional time delay due to the evaporation and transport processes, as seen in Figure 7, which can be a critical factor during combustion instability process. The time delay from only the convection process is plotted in Figure 8 phase plot as a solid red line. While the phase calculated from convection time delay agrees well for lower frequencies, there is a mismatch for higher frequencies. Therefore, an accurate characterisation of the evaporation response to acoustic excitation is necessary for robust thermoacoustic stability analysis of realistic spray flames.

Figure 7.

Time series of the air velocity (top) and fuel mass fraction measured at $\tilde{x}=11.2$ for $D_0=30$   $\mu \mathrm{m}$ , f = 100  Hz, ${\overline{u}}_c=5$   m/s, $u^{\prime }=2.5$   m/s

Figure 8.

Calculated Analytical (—) and CFD (- - -) evaporation transfer function. The convective phase decay is shown in red.

Summary and Outlook

This work proposes an analytical solution of the response of a mono-disperse population of droplets to acoustic excitation in terms of the number density. Due to relative timing of the injection of droplets with gas velocity, a ND wave is formed. The analytical expression for a particle population response is extended by incorporating the evaporation of the droplets. The resulting analytical formulation describes the oscillatory evaporation rate for the linear drag regime. It is shown that the evaporation rate profile inherits the characteristics of the ND wave without any phase lag. The oscillatory evaporation rate gives rise to a fuel vapour wave, which convects downstream with the mean flow speed.

The propagation of vapour wave manifests itself in the form of equivalence ratio fluctuations. The latter are determined by solving a 1D convection-diffusion equation for gaseous phase using an appropriate Green’s function. Results from the analytical approach are in good agreement with 1D Euler-Lagrange CFD simulations. The evaporation response of monodisperse evaporative spray to acoustic oscillations has been calculated. The resulting transfer function exhibits a low-pass behaviour. Additionally, the evaporation process introduces a characteristic time delay in the equivalence ratio fluctuation. An accurate determination of such a time delay is critical for thermoacoustic stability analysis.

Future work will explore the limits of the validity of the assumptions made in this work. Subsequent work will also consider the inclusion of transient droplet heating, which can be seamlessly integrated to the current formalism. Droplet heating may introduce a non-negligible time delay in the evaporation process. As such, a simple estimate of the global time delay associated with convection process ( $t=x/V$ as done above) will not be accurate. Furthermore, the present work can be easily extended to account for polydisperse flows. For that purpose, the current evaporation response has to be evaluated over the size distribution of the droplets. Finally, it will be important to perform the evaporation response study for turbulent flows, where increased mixing and turbulent diffusion might play a key role in determining the amplitude of equivalence ratio oscillations.