2 Vortex sheet formation of flickering buoyant diffusion flames

As stated above, the essence of the flickering flame lies in the dynamics of the toroidal vortices. The first question to ask is how the toroidal vortex forms. In vortex dynamics, the formation and evolution of a toroidal vortex, formally known as the ‘vortex ring’, have been studied extensively (Maxworthy Reference Maxworthy1972, Reference Maxworthy1977; Saffman Reference Saffman1978; Didden Reference Didden1979; Glezer Reference Glezer1988; Shariff & Leonard Reference Shariff and Leonard1992; Gharib, Rambod & Shariff Reference Gharib, Rambod and Shariff1998). Physically, the appearance of a vortex ring may be considered as the outcome of the growing and rolling-up of a cylindrical-shaped vortex sheet. A well-known example is the starting vortex jet (Didden Reference Didden1979; Nitsche & Krasny Reference Nitsche and Krasny1994; Gharib et al. Reference Gharib, Rambod and Shariff1998; Mohseni & Gharib Reference Mohseni and Gharib1998), where the vortex sheet is continuously supplied by pushing a fluid slug out of a circular jet nozzle to form a cylindrical shear layer. From this perspective, the toroidal vortex of a diffusion flame should not be fundamentally different and its formation must involve a growing vortex sheet.

A schematic of vortex sheet in a laminar diffusion flame is shown in figure 2. The vorticity growth inside the vortex sheet is governed by the vorticity transport equation,

(2.1) $$\begin{eqnarray}\frac{\text{D}\unicode[STIX]{x1D74E}}{\text{D}t}=(\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-\unicode[STIX]{x1D74E}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})+\frac{1}{\unicode[STIX]{x1D70C}^{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p)+\frac{\unicode[STIX]{x1D70C}_{A}}{\unicode[STIX]{x1D70C}^{2}}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \boldsymbol{g})+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D74E},\end{eqnarray}$$

where $\boldsymbol{u}$ and $\unicode[STIX]{x1D74E}$ are the velocity and vorticity vectors, $\unicode[STIX]{x1D70C}$ the local density and $p$ the gauge pressure, $\boldsymbol{g}$ the gravitational acceleration vector, $\unicode[STIX]{x1D708}$ the kinematic viscosity and $\unicode[STIX]{x1D70C}_{A}$ the gas density of the ambient environment. On the right-hand side of (2.1), the first vortex tilting/stretching term vanishes for either two-dimensional flows or axisymmetric flows without swirling, the second dilatation term vanishes for incompressible flows and the fifth diffusion term describes the redistribution of vorticity and thus is not a source of vorticity production. Vorticity generation are attributed to the third baroclinic term and the fourth gravitational term, both of which entail the presence of variable density. In this case, the growth of the vortex sheet is caused by the flame-induced vorticity addition mechanisms, which differ from the vorticity flux supplied by the inflow of a jet.

Figure 2. Schematic of the vortex sheet for the growth of a toroidal vortex in an axisymmetric laminar diffusion flame. Note that the zero thickness of flame and vortex sheets is exaggerated for illustration.

To evaluate the formation of a flame-induced vortex sheet (or viscous shear layer) and its roll-up into a toroidal vortex, we focus on the circulation ( $\unicode[STIX]{x1D6E4}$ ) of a control mass, defined by

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\oint _{\unicode[STIX]{x2202}A}\boldsymbol{u}\boldsymbol{\cdot }\text{d}\boldsymbol{s}.\end{eqnarray}$$

As illustrated in figure 2, the dashed box $\unicode[STIX]{x2202}A$ is a material contour around the vortex sheet segment between $s_{v1}$ and $s_{v2}$ , and $A$ is the area encircled by the contour so $A$ is a control mass; $\text{d}\boldsymbol{s}$ represents a material line element along $\unicode[STIX]{x2202}A$ . The vorticity vector has only one azimuthal component as $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D714}\hat{\unicode[STIX]{x1D73D}}$ under the assumption of axisymmetric flow without swirling (the axial vorticity is important in swirling flames as studied by Klimenko & Williams (Reference Klimenko and Williams2013) and Yu & Zhang (Reference Yu and Zhang2017a ,Reference Yu and Zhang b )). Applying the divergence theorem to (2.2) yields $\unicode[STIX]{x1D6E4}=\iint \unicode[STIX]{x1D714}\,\text{d}A$ , meaning $\unicode[STIX]{x1D6E4}$ is a measure for the total vorticity inside  $A$ .

The rate of total change of $\unicode[STIX]{x1D6E4}$ can be derived from (2.2) as

(2.3) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}}{\text{d}t}=\oint _{\unicode[STIX]{x2202}A}\frac{\text{D}\boldsymbol{u}}{\text{D}t}\boldsymbol{\cdot }\text{d}\boldsymbol{s}+\oint _{\unicode[STIX]{x2202}A}\boldsymbol{u}\boldsymbol{\cdot }\frac{\text{D}(\text{d}\boldsymbol{s})}{\text{D}t}.\end{eqnarray}$$

According to the identity $\text{D}(\text{d}\boldsymbol{s})/\text{D}t=\text{d}\boldsymbol{u}$ (Wu, Ma & Zhou Reference Wu, Ma and Zhou2007), the second integral term of (2.3) vanishes as $\oint _{\unicode[STIX]{x2202}A}\boldsymbol{u}\boldsymbol{\cdot }\text{d}\boldsymbol{u}=\oint _{\unicode[STIX]{x2202}A}\,\text{d}(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u})/2=0$ ; the integrand $\text{D}\boldsymbol{u}/\text{D}t$ in the first integral term is given by the Navier–Stokes equation,

(2.4) $$\begin{eqnarray}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-\frac{1}{\unicode[STIX]{x1D70C}}[\unicode[STIX]{x1D735}p+(\unicode[STIX]{x1D70C}_{A}-\unicode[STIX]{x1D70C})\boldsymbol{g}]+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\end{eqnarray}$$

where the gauge pressure $p$ is related to the absolute pressure $p_{ab}$ as $\unicode[STIX]{x1D735}p_{ab}=\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D70C}_{A}\boldsymbol{g}$ . Here, $\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}$ is the shear stress (per unit mass) and approaches zero outside the viscous shear layer, where the flow is effectively inviscid. This implies $\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}=\mathbf{0}$ on $\unicode[STIX]{x2202}A$ . The term $|\unicode[STIX]{x1D735}p|/|(\unicode[STIX]{x1D70C}_{A}-\unicode[STIX]{x1D70C})\boldsymbol{g}|$ physically represents a characteristic Froude number. For buoyancy-dominated flames, we adopt the small Froude number assumption ( $\mathit{Fr}\ll 1$ ) so the pressure gradient is negligible compared with the gravity term in (2.4). With $\text{d}\boldsymbol{s}$ expressed as $\hat{\boldsymbol{s}}\,\text{d}s$ , where $\hat{\boldsymbol{s}}$ is the unit tangential vector along the contour $\unicode[STIX]{x2202}A$ , equation (2.3) can be rewritten as

(2.5) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}}{\text{d}t}=\oint _{\unicode[STIX]{x2202}A}-\frac{\unicode[STIX]{x1D70C}_{A}-\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}}\boldsymbol{g}\boldsymbol{\cdot }\hat{\boldsymbol{s}}\,\text{d}s=\int _{s_{v1}}^{s_{v2}}\unicode[STIX]{x1D70C}_{A}\left(\frac{1}{\unicode[STIX]{x1D70C}_{1}}-\frac{1}{\unicode[STIX]{x1D70C}_{2}}\right)(\hat{\boldsymbol{s}}_{v}\boldsymbol{\cdot }\boldsymbol{g})\,\text{d}s,\end{eqnarray}$$

where $\hat{\boldsymbol{s}}_{v}$ is the unit tangential vector of the vortex sheet as shown in figure 2. Assuming the gas density at either side of the vortex sheet is constant, which means $\unicode[STIX]{x1D70C}_{1}(s)=\unicode[STIX]{x1D70C}_{f}$ where $\unicode[STIX]{x1D70C}_{f}$ is the density of the gas at the flame sheet and $\unicode[STIX]{x1D70C}_{2}(s)=\unicode[STIX]{x1D70C}_{A}$ , we have

(2.6) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}}{\text{d}t}=\unicode[STIX]{x1D70C}_{A}\left(\frac{1}{\unicode[STIX]{x1D70C}_{1}}-\frac{1}{\unicode[STIX]{x1D70C}_{2}}\right)\left(\int _{s_{v1}}^{s_{v2}}\,\text{d}\hat{\boldsymbol{s}}_{v}\boldsymbol{\cdot }\boldsymbol{g}\right)=-g\unicode[STIX]{x0394}z(r^{\ast }-1),\end{eqnarray}$$

where $\unicode[STIX]{x0394}z$ is the vertical length of the vortex sheet associated with the control mass  $A$ , as shown in figure 2. The density ratio, $r^{\ast }=\unicode[STIX]{x1D70C}_{A}/\unicode[STIX]{x1D70C}_{f}$ , is a measurable quantity for a given flame. For a rough estimation of  $r^{\ast }$ , assuming the flame is isobaric and following the ideal gas law, we have $r^{\ast }=\unicode[STIX]{x1D70C}_{A}/\unicode[STIX]{x1D70C}_{f}=T_{f}/T_{A}$ , where $T_{A}$ and $T_{f}$ are the temperatures of the ambient air and the flame sheet, respectively. For common fuels burned in air, the flame temperature at normal atmospheric conditions (1 atm, $T_{A}\approx 300$  K) varies in a wide range roughly between 1200 K and 2400 K, corresponding to the $r^{\ast }$ range of 4–8.

Equation (2.6) indicates an important feature of the buoyancy-induced vortex sheet that the generation rate of total vorticity (circulation) is independent of the geometric shape of the sheet, and only dictated by the vertical length of the sheet.

3 Formation of periodic toroidal vortices

It is evident from the flow visualization of Chen et al. (Reference Chen, Seaba, Roquemore and Goss1989) and numerical simulation of Katta et al. (Reference Katta, Goss and Roquemore1994), among others, that the periodicity of the flame structures is strongly correlated with the toroidal vortices, which are rolled up from the buoyancy-induced vortex sheet. The study of unsteady diffusion flames by Cantwell, Lewis & Chen (Reference Cantwell, Lewis and Chen1989) further illustrated the sequential evolution of an individual toroidal vortex and its coupling with the dynamics of the flame. Following these studies, we present the formation process of a toroidal vortex in figure 3 to illustrate the relation between the flickering flame and the periodic toroidal vortices. At $t=0$ , a previous toroidal vortex (marked by grey-dashed line) matures and detaches, leading to the generation of a new toroidal vortex near the base of the flame. This new vortex core together with its trailing vortex sheet (marked by blue-solid line) then grow and advect downstream under buoyancy. At $t=\unicode[STIX]{x1D70F}$ , where $\unicode[STIX]{x1D70F}=1/f$ is the period of flame flickering, the new toroidal vortex becomes fully developed and its circulation is defined as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{TV}=\unicode[STIX]{x1D6E4}_{B}(\unicode[STIX]{x1D70F})=\int _{0}^{\unicode[STIX]{x1D70F}}\frac{\text{d}\unicode[STIX]{x1D6E4}_{B}}{\text{d}t}\,\text{d}t,\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{B}$ is the total circulation of the moving control volume $B$ (marked by the red-dashed box), which encloses the vortex core and its trailing vortex sheet.

Figure 3. Schematic of the periodic formation process of a toroidal vortex. The vortex sheets associated with the formation of toroidal vortex are tracked by the blue-dashed lines.

The rate of change of $\unicode[STIX]{x1D6E4}_{B}$ is given by

(3.2) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}_{B}}{\text{d}t}=\frac{\text{d}}{\text{d}t}\int _{B}\,\text{d}\unicode[STIX]{x1D6E4}=\frac{\text{d}}{\text{d}t}\int _{0}^{h(t)}\unicode[STIX]{x1D6FE}_{z}(z,t)\,\text{d}z,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{z}(z,t)=\text{d}\unicode[STIX]{x1D6E4}/\text{d}t$ . $h(t)$ is the height of the upper boundary of  $B$ , which rises following the convection of the toroidal vortex. Applying the Leibniz integral rule to (3.2), we have

(3.3) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}_{B}}{\text{d}t}=\unicode[STIX]{x1D6FE}_{z}(h(t),t)\frac{\text{d}h(t)}{\text{d}t}+\int _{0}^{h(t)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}_{z}(z,t)}{\unicode[STIX]{x2202}t}\,\text{d}z.\end{eqnarray}$$

Considering the control mass ${\mathcal{V}}_{B}$ that instantaneously overlaps with the control volume $B$ at time  $t$ , and applying (2.6) to ${\mathcal{V}}_{B}$ , we can obtain the rate of change of the circulation in  ${\mathcal{V}}_{B}$ :

(3.4) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}_{{\mathcal{V}}_{B}}}{\text{d}t}=-gh(t)(r^{\ast }-1).\end{eqnarray}$$

Further noting that $\unicode[STIX]{x1D6E4}_{{\mathcal{V}}_{B}}=\int _{{\mathcal{V}}_{B}}\unicode[STIX]{x1D6FE}_{z}(z,t)\,\text{d}z$ , $\unicode[STIX]{x1D6E4}_{{\mathcal{V}}_{B}}$ can be also derived using the Reynolds transport theorem (Batchelor Reference Batchelor1967) as

(3.5) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D6E4}_{{\mathcal{V}}_{B}}}{\text{d}t} & = & \displaystyle \unicode[STIX]{x1D6FE}_{z}(h(t),t)\left.\frac{\text{d}z}{\text{d}t}\right|_{z=h(t)}-\unicode[STIX]{x1D6FE}_{z}(0,t)\left.\frac{\text{d}z}{\text{d}t}\right|_{z=0}+\int _{0}^{h(t)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}_{z}(z,t)}{\unicode[STIX]{x2202}t}\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6FE}_{z}(h(t),t)\frac{\text{d}h(t)}{\text{d}t}-\unicode[STIX]{x1D6FE}_{z}(0,t)\left.\frac{\text{d}z}{\text{d}t}\right|_{z=0}+\int _{0}^{h(t)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}_{z}(z,t)}{\unicode[STIX]{x2202}t}\,\text{d}z.\end{eqnarray}$$

Comparing (3.3) and (3.5) and using (3.4), we obtain

(3.6) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6E4}_{B}}{\text{d}t}=-gh(t)(r^{\ast }-1)+\dot{\unicode[STIX]{x1D6E4}}_{i}(t),\end{eqnarray}$$

where $\dot{\unicode[STIX]{x1D6E4}}_{i}(t)=(\text{d}\unicode[STIX]{x1D6E4}/\text{d}t)|_{z=0}$ is the rate of circulation addition by the inflow, which enters the system through the lower boundary of  $B$ .

In the current study, $\dot{\unicode[STIX]{x1D6E4}}_{i}(t)$ reflects the possible presence of an initial vortex sheet associated with the fuel stream at the fuel inlet, either from a jet or induced by evaporation of a pool flame. The fuel stream may be considered as a jet flow, the rate of circulation addition by which was given by Didden (Reference Didden1979) as

(3.7) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D6E4}}_{i}(t)=-C_{j}U_{0}^{2},\end{eqnarray}$$

where $U_{0}$ is the initial jet velocity and $C_{j}$ is a constant relating to the configuration and boundary condition of the jet exit. According to Krieg & Mohseni (Reference Krieg and Mohseni2013), $C_{j}$ is 0.5 for an ideal parallel-nozzle jet, but it could be inaccurate if the jet entrains a radial flow and causes additional circulation generation. For example, the effective $C_{j}$ for a converging jet can be under-predicted by a factor of 3 (Krieg & Mohseni Reference Krieg and Mohseni2013). Consequently, we shall compare the results of $C_{j}=0.5$ and $C_{j}=1.5$ in the theoretical predictions.

Applying (3.6) and (3.7) to (3.1), we obtain

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{TV}=-gH\unicode[STIX]{x1D70F}(r^{\ast }-1)-C_{j}U_{0}^{2}\unicode[STIX]{x1D70F},\end{eqnarray}$$

where $H=\unicode[STIX]{x1D70F}^{-1}\int _{0}^{\unicode[STIX]{x1D70F}}h(t)\,\text{d}t$ is the time-averaged height of the growing toroidal vortex. The integration of (3.8) is based on the assumption that the vortex sheet grows in a quasi-steady manner, given that the flame flickering period is so short that the densities of the gases across the vortex sheet remain approximately unchanged meanwhile.

In the current study, the two common types of diffusion flames are pool flames, where the fuel vapour enters the system from the evaporation of liquid pool, and jet flames, where the fuel stream is supplied via a gaseous jet. The initial fuel velocity $U_{0}$ is a given boundary condition for jet flames, whereas for pool flames it is an eigenvalue of the thermochemical-coupled equations (Yu & Zhang Reference Yu and Zhang2017a ,Reference Yu and Zhang b ). The effect of fuel velocity on the flickering frequency of jet flames has been investigated by Hamins et al. (Reference Hamins, Yang and Kashiwagi1992), Sato, Amagai & Arai (Reference Sato, Amagai and Arai2000) and Fang et al. (Reference Fang, Wang, Guan, Zhang and Wang2016), among others, which indicates the existence of two distinct regimes, namely momentum driven and buoyancy driven. Since buoyancy-driven flames have a characteristic velocity of $U=\sqrt{gD}$ , the motion of the upper boundary of the toroidal vortex may be roughly estimated as $h(t)=C_{h}Ut$ with $C_{h}$ being a constant prefactor; so $H=C_{h}U\unicode[STIX]{x1D70F}/2$ . With $\unicode[STIX]{x1D6E4}_{TV}$ scaled by $-UD$ , equation (3.8) can be written in the non-dimensional form as

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{TV}^{\ast }=\frac{C_{h}}{2}\mathit{Ri}\unicode[STIX]{x1D70F}^{\ast 2}+C_{j}\sqrt{\mathit{Fr}}\unicode[STIX]{x1D70F}^{\ast },\end{eqnarray}$$

where $\unicode[STIX]{x1D70F}^{\ast }=U_{0}\unicode[STIX]{x1D70F}/D$ . The Richardson number and the Froude number are defined as $\mathit{Ri}=(r^{\ast }-1)gD/U_{0}^{2}$ and $\mathit{Fr}=U_{0}^{2}/gD$ , respectively. It is noted that, in theory, equation (3.9) is only valid for buoyancy-driven flames that are characterized by $U_{0}\ll U$ , corresponding to $\mathit{Fr}\ll 1$ or $\mathit{Ri}\gg 1$ .

In this study, we propose that flame flickering is caused by the alternate formation and detachment of the toroidal vortices, the fundamental mechanism of which does not conflict with the Kelvin–Helmholtz instability proposed in the literature (Chen & Roquemore Reference Chen and Roquemore1986; Buckmaster & Peters Reference Buckmaster and Peters1988; Chen et al. Reference Chen, Seaba, Roquemore and Goss1989). In fact, these mechanisms share two key similarities. First, both mechanisms are characterized by vortex roll-up. Second, they are caused by velocity shear, and are closely related to strength of the shear layer. The difference between them can be understood as spatial or temporal development of the instability. Specifically, a typical Kelvin–Helmholtz instability can be thought as a spatial one where multiple vortices roll up simultaneously along a uniform shear layer in space; whereas the periodic formation and detachment of vortices can be considered as a temporal instability where local growth of vortex sheet causes the vortices to develop in sequence. Thus, a series of fully developed vortices together render the feature of the Kelvin–Helmholtz instability.

4 Theoretical predictions

The above analysis hitherto has addressed the formation and growth of a toroidal vortex, which leaves us another question: how does the toroidal vortex shed? The answer to a similar problem on vortex ring formation has been given by Gharib et al. (Reference Gharib, Rambod and Shariff1998) and Mohseni & Gharib (Reference Mohseni and Gharib1998), who pointed to the existence of a universal non-dimensional formation number, above which a vortex ring would be too strong to maintain growth and consequently detach from its vorticity-feeding shear layer. In other words, the formation number can be considered as a dimensionless measure of the upper limit of the total circulation of an attached vortex ring. This mechanism can be interpreted as a unified condition for the detachment of a rolled-up vortex ring. We hypothesize that the shedding of the current toroidal vortex is controlled by the same mechanism, because the strength of the vortex sheet is constantly growing under the buoyancy-induced vorticity generation, and will roll up into a vortex once the vortex sheet becomes sufficiently strong. However, this should be distinguished from the case of constant laminar flow through a nozzle, where the strength of the vortex sheet remains unchanged. As the result, the initial vortex sheet is too weak to roll up unless otherwise being disturbed or instability develops in the downstream. This explains why a continuous laminar jet generally does not have alternate formation and detachment of vortex rings.

From the above discussion, the continuous growing shear layer causes gradual accumulation of the circulation inside the toroidal vortex until it reaches a threshold denoted by  $C$ . The typical value of $C=4$ for an ideal starting vortex jet (Gharib et al. Reference Gharib, Rambod and Shariff1998) is adopted here, although it could change notably under various conditions (Dabiri & Gharib Reference Dabiri and Gharib2005; Krueger, Dabiri & Gharib Reference Krueger, Dabiri and Gharib2006; Lawson & Dawson Reference Lawson and Dawson2013; Xia & Mohseni Reference Xia and Mohseni2015). Applying $\unicode[STIX]{x1D6E4}_{TV}^{\ast }=C$ to (3.9) and solving for $\unicode[STIX]{x1D70F}$ yields

(4.1) $$\begin{eqnarray}f=\frac{1}{\unicode[STIX]{x1D70F}}=\frac{1}{2C}\sqrt{\frac{g}{D}}\left(C_{j}\mathit{Fr}+\sqrt{C_{j}^{2}\mathit{Fr }^{2}+2CC_{h}(r^{\ast }-1)}\right).\end{eqnarray}$$

This completes the derivation of the frequency relation for buoyancy-driven diffusion flames. Equation (4.1) is similar to the scaling formula given by (1) of Cetegen & Ahmed (Reference Cetegen and Ahmed1993), with different power on the $\mathit{Fr}$ term. An extended discussion on the differences between the two formulas will be provided at the end of this section. Now, let $\mathit{Fr}\rightarrow 0$ , equation (4.1) can be further simplified to

(4.2) $$\begin{eqnarray}f=\sqrt{\frac{C_{h}(r^{\ast }-1)}{2C}\cdot \frac{g}{D}},\end{eqnarray}$$

which recovers the prominent scaling law, $f\propto (g/D)^{1/2}$ , obtained by Byram & Nelson (Reference Byram and Nelson1970), among others.

Figure 4. (a) Flickering frequencies ( $f$ ) of pool and jet flames as a function of $\sqrt{g/D}$ , with data collected from Byram & Nelson (Reference Byram and Nelson1970), Sibulkin & Hansen (Reference Sibulkin and Hansen1975), Schönbucher et al. (Reference Schönbucher, Arnold, Banhardt, Bieller, Kasper, Kaufmann, Lucas and Schiess1988), Baum & McCaffrey (Reference Baum and McCaffrey1989), Weckman & Sobiesiak (Reference Weckman and Sobiesiak1989), Hamins et al. (Reference Hamins, Yang and Kashiwagi1992), Cetegen & Ahmed (Reference Cetegen and Ahmed1993), Durox et al. (Reference Durox, Yuan, Baillot and Most1995), Ghoniem et al. (Reference Ghoniem, Lakkis and Soteriou1996), Yoshihara, Ito & Torikai (Reference Yoshihara, Ito and Torikai2013), Fang et al. (Reference Fang, Wang, Guan, Zhang and Wang2016). The gravitational acceleration was adjusted in the range of $1.5g{-}6.0g$ for Durox et al.’s (Reference Durox, Yuan, Baillot and Most1995) jet flames, and $0.5g{-}1.0g$ for Yoshihara et al.’s (Reference Yoshihara, Ito and Torikai2013) pool flames. ‘PF’ and ‘JF’ in the legend represent pool flame and jet flame, respectively. (b) Theoretical predictions of $\mathit{St}$ versus $1/\mathit{Fr}$ for all jet flames in (a).

Next, we compare the theoretical formula with experiment for validation. Figure 4(a) shows the relation of $f$ against $\sqrt{g/D}$ for various pool flames and jet flames from existing literature, which have diverse fuel types, fire source dimensions and gravities. The experimental data are in good agreement with (4.2), and confirms the unification of the theory for buoyant diffusion flames. Similar comparisons of data were reported by Hamins et al. (Reference Hamins, Yang and Kashiwagi1992) and Cetegen & Ahmed (Reference Cetegen and Ahmed1993) to demonstrate the relation $f\sim \sqrt{D}$ . However, we further verified the relation of $f\sim \sqrt{g/D}$ with pool flames of varying gravity between $0.5g$ and $1.0g$ (Yoshihara et al. Reference Yoshihara, Ito and Torikai2013) and jet flames of varying gravity between $1.5g$ and $6.0g$ (Durox et al. Reference Durox, Yuan, Baillot and Most1995). The trend line obtained from figure 4(a), $f=0.48\sqrt{g/D}$ with $g$ being $9.8~\text{m}~\text{s}^{-2}$ , is equivalent to the scaling relation, $f=1.5\sqrt{1/D}$ , given by Cetegen & Ahmed (Reference Cetegen and Ahmed1993). This suggests that the relation,

(4.3) $$\begin{eqnarray}\sqrt{\frac{C_{h}(r^{\ast }-1)}{2C}}=0.48,\end{eqnarray}$$

applies to all buoyant diffusion flames. Consequently, we can use (4.3) to eliminate the undetermined prefactor $C_{h}$ and rewrite (4.1) as

(4.4) $$\begin{eqnarray}\mathit{St}=\frac{C_{j}}{2C}\sqrt{\mathit{Fr}}+\sqrt{\left(\frac{C_{j}}{2C}\right)^{2}\mathit{Fr}+\frac{1}{4.34\mathit{Fr}}},\end{eqnarray}$$

where the Strouhal number is defined as $\mathit{St}=fD/U_{0}$ . We have arrived at the non-dimensional frequency formula of flame flickering, where the only undetermined scaling factor is  $C_{j}$ , related to the boundary condition of the fuel inlet.

Figure 4(b) plots the experimental data from all jet flames in figure 4(a), showing generally promising agreement with the predictions by (4.4). It is seen that the theoretical predictions with $C_{j}=0.5$ and $C_{j}=1.5$ are almost identical for $\mathit{Fr}<0.1$ , and $C_{j}=0.5$ shows a better agreement with the experimental data for $\mathit{Fr}>0.1$ . This is beyond the range of $\mathit{Fr}\ll 1$ , based on which the present theory is derived. The result does suggest that the applicable range of this theory can be extended to $\mathit{Fr}\leqslant 1$ . Equation (4.4) should also be valid for pool flames; however, since $U_{0}$ values for those pool flames in figure 4(a) were not measured or reported, a direct validation is not available here. As a reference, we also plot the scaling relation of Hamins et al. (Reference Hamins, Yang and Kashiwagi1992), which was obtained based on fitting jet flame data in the range of $10^{-6}<\mathit{Fr}<10^{8}$ . It is seen that their fitting line, with an exponent of $-0.57$ , does not yield the best match with data in the $\mathit{Fr}<1$ region, whereas the current theory (4.4) predicts that in the limit $\mathit{Fr}\rightarrow 0$ the exponent of the scaling law should be exactly  $-0.5$ .

Finally, we contrast our theory (4.4) with (1) of Cetegen & Ahmed (Reference Cetegen and Ahmed1993). Their equation was derived by applying the Bernoulli equation to calculate the convective velocity of the flame puff, which was then integrated to obtain the period as a function of the flame puff height, $H_{f}$ , in the form of ((A12) of Cetegen & Ahmed Reference Cetegen and Ahmed1993)

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D70F}^{\ast }=\frac{1}{\mathit{St}}=\frac{1}{C_{f}\mathit{Ri}}\left(\sqrt{\frac{2\mathit{Ri }H_{f}}{D}+1}-1\right),\end{eqnarray}$$

where $C_{f}$ is a constant prefactor. It should be noted that Cetegen and Ahmed made two key assumptions. The first is a linear correlation between the inner flame velocity and the convective velocity of the flame puff, which was loosely implied from their previous work (Zukoski et al. Reference Zukoski, Cetegen and Kubota1985). The second is an implicit assumption that $H_{f}=D/2$ , which removes possible dependence of the height of the flame puff on the dynamics of the toroidal vortex. So far, this assumption has not been supported by the existing literature. The current theory does not require these assumptions. For quantitative comparison, we rewrite (1) of Cetegen & Ahmed (Reference Cetegen and Ahmed1993) in the $\mathit{St}$ – $\mathit{Fr}$ form,

(4.6) $$\begin{eqnarray}\mathit{St}=K\left(\frac{1}{\sqrt{r^{\ast }-1}}+\sqrt{\frac{1}{r^{\ast }-1}+\frac{1}{\mathit{Fr}}}\right),\end{eqnarray}$$

where the fitting parameter $K$ was suggested by Cetegen and Ahmed to be 0.5 based on their experimental data. As seen in figure 5, both theories recover the same scaling law of $\mathit{St}\sim \mathit{Fr}^{-0.5}$ in the limit $\mathit{Fr}\rightarrow 0$ . It is noted that (4.6) is explicitly dependent on $r^{\ast }$ , which varies in a wide range for diffusion flames. This might explain why Cetegen & Ahmed’s (Reference Cetegen and Ahmed1993) formula predicted different $K$ values for different jet flames as shown in their figure 14. In contrast, the present theory (4.4) is a unified one for buoyancy-driven jet flames, and relies on two physical constants $C$ and  $C_{j}$ , the former accounting for the toroidal vortex detachment and the latter for the initial circulation of inflow.

Figure 5. Comparison between the current theory (4.4) with $C_{j}=0.5$ and Cetegen & Ahmed’s (Reference Cetegen and Ahmed1993) model (4.6).  $r^{\ast }$ for a common diffusion flame varies between 4 and 8.