3.1. Theoretical model

To understand the spatial evolutions of the mean ISL and OSL in a swirling flame, figure 3(a) presents a diagram of the shear layers, where a circular control volume with infinitesimal height $\Delta z$ is defined at the axial location $z$ of the ISL. To preserve generality, we next develop the theoretical framework for a general axisymmetric shear layer, with its inner and outer boundaries described by $\delta _1(z)$ and $\delta _2(z)$ , respectively.

Figure 3.

Schematics of (a) the control volume (CV) for an ISL segment and (b) its radial discretization.

We start with the Reynolds transport theorem (Batchelor Reference Batchelor2000) in the form

(3.1) \begin{equation} \frac {\mathrm {d}}{\mathrm {d}t}\int _{V}^{}\psi\, \mathrm {d}V=\int _{V^*}^{}\left [\dfrac {\mathrm {d}\psi }{\mathrm {d}t}+\psi\, \boldsymbol {\nabla \cdot u}\right ] \mathrm {d}V=\int _{V^*}^{}\frac {\partial \psi }{\partial t}\,{\rm d}V+\int _{A^*}^{}\psi\, \boldsymbol {u \cdot n}\,\mathrm {d}A, \end{equation}

where $\psi$ represents any intensive quantity, $V$ the material volume, and $V^*$ and $A^*$ the corresponding stationary control volume and control surface, respectively. To account for the mean flow, we apply the decompositions $\psi =\bar {\psi }+\psi ^\prime$ and $\boldsymbol {u}=\bar {\boldsymbol {u}}+\boldsymbol {u}^\prime$ to express (3.1) in terms of mean parts and fluctuating parts:

(3.2) \begin{align} \begin{split} \int _{V^*}^{}\left [\frac {\partial \bar {\psi }}{\partial t}+\bar {\boldsymbol {u}}\boldsymbol {\cdot } \boldsymbol {\nabla }\bar {\psi }+\bar {\psi }\,\boldsymbol{\nabla} \boldsymbol{\cdot} \bar{ \boldsymbol{u}}+\overline {\boldsymbol {u}^\prime \boldsymbol {\cdot} \boldsymbol{\nabla}\psi ^\prime }+\overline {\psi ^\prime\, \boldsymbol {\nabla} \boldsymbol{\cdot} \boldsymbol{u}^\prime }\right ] \mathrm {d}V \\ =\int _{V^*}^{}\frac {\partial \bar {\psi }}{\partial t}\,{\rm d}V+\int _{A^*}^{}\left [\bar {\psi }\,\bar{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{n}+{\overline {{\psi ^\prime }\, \boldsymbol{u}^\prime\boldsymbol{\cdot}\boldsymbol{n}}}\right ]\mathrm {d}A. \end{split} \end{align}

The fluctuation terms can be cancelled out to yield the mean equation

(3.3) \begin{equation} \int _{V^*}^{}\left [\dfrac {\mathrm {d}\psi }{\mathrm {d}t}+\psi\, \boldsymbol {\nabla \cdot u}\right ] \mathrm {d}V=\int _{A^*}^{}\psi\, \boldsymbol {u \cdot n}\,\mathrm {d}A, \end{equation}

where the $\bar {\ }$ symbol denoting the mean quantities has been dropped hereafter for simplicity. Note that although $\mathrm {d}\psi /\mathrm {d}t$ here only has the convection term, we retain the total derivative for its physical meaning. For an axisymmetric shear layer with swirl, we divide the circular control volume radially into $n$ ring-like segments, each with infinitesimal width $\Delta r$ , as shown in figure 3(b). For the $i$ th segment with average radius $r_i$ , substituting $\psi$ in (3.3) with the azimuthal vorticity $\omega _\theta$ yields

(3.4) \begin{align} \begin{split} 2\pi r_i\, \Delta z\, \Delta r \left [\frac {\mathrm {d} \omega _\theta }{\mathrm {d}t}+ \omega _\theta \left (\boldsymbol {\nabla \cdot u}\right )\right ] = 2\pi r_i\,\Delta r\, u_z \omega _\theta \left|_{z}^{z+\Delta z}+2\pi r\,\Delta z\, u_r \omega _\theta \right|_{r_i-{\Delta r}/{2}.}^{r_i+{\Delta r}/{2},} \end{split} \end{align}

Dividing (3.4) by $2\pi r_i$ and then integrating it for $r$ from $\delta _1$ to $\delta _2$ with the boundary condition $\omega _\theta (\delta _1)=\omega _\theta (\delta _2)=0$ , the last term in (3.4) can be derived as

(3.5) \begin{align} \begin{split} \sum _{i=1}^{n-1}\frac {1}{2}\left(\frac {1}{r_i}+\frac {1}{r_{i+1}}\right)u_r\omega _\theta\, \Delta z\, \Delta r =\Delta z\int _{\delta _1}^{\delta _2}\frac {u_r\omega _\theta }{r}\,\mathrm {d}r. \end{split} \end{align}

So the integral of (3.4) becomes

(3.6) \begin{align} \begin{split} \Delta z\int _{\delta _1}^{\delta _2} \left [\frac {\mathrm {d} \omega _\theta }{\mathrm {d}t}+ \omega _\theta \left (\boldsymbol {\nabla \cdot u}\right )\right ]\mathrm {d}r = \left.\int _{\delta _1}^{\delta _2} u_z \omega _\theta \,\mathrm {d}r \right|_{z}^{z+\Delta z}+\Delta z\int _{\delta _1}^{\delta _2}\frac {u_r\omega _\theta }{r}\,\mathrm {d}r. \end{split} \end{align}

The total change of vorticity $\mathrm {d}\omega _\theta /\mathrm {d}t$ in (3.6) can be further calculated based on the Reynolds-averaged vorticity transport equation (for axisymmetric swirling flow),

(3.7) \begin{align} \frac {\mathrm {d}\omega _\theta }{\mathrm {d}t} &=\underbrace {u_r\frac {\partial \omega _\theta }{\partial r}+u_z\frac {\partial \omega _\theta }{\partial z}+\frac {\omega _r u_\theta }{r}}_{\mathrm {convection}} =\underbrace {\omega _r\frac {\partial u_\theta }{\partial r}+\omega _z\frac {\partial u_\theta }{\partial z}}_{\mathrm {tilting}}+\underbrace {\frac {u_r\omega _\theta }{r}}_{\mathrm {stretching}}\underbrace {{}-\omega _\theta \left [\frac {1}{r}\frac {\partial \left (ru_r\right )}{\partial r} +\frac {\partial u_z}{\partial z}\right ]}_{\mathrm {thermal \: dilatation}} \nonumber\\&\quad+\underbrace {\frac {1}{\rho ^2}\left (\frac {\partial \rho }{\partial z}\frac {\partial p}{\partial r}-\frac {\partial \rho }{\partial r}\frac {\partial p}{\partial z}\right )}_{\mathrm {baroclinic \: torque}}\underbrace {{}-\frac {\rho _A}{\rho ^2}\left (\frac {\partial \rho }{\partial r}g_z\right )}_{\mathrm {buoyancy}}\nonumber\\ &\quad +\underbrace {\left (\nu +\nu _t\right )\left (\frac {\partial ^2 \omega _\theta }{\partial r^2}+\frac {1}{r}\frac {\partial \omega _\theta }{\partial r}+\frac {\partial ^2 \omega _\theta }{\partial z^2}-\frac {\omega _\theta }{r^2}\right )}_{\mathrm {diffusion}}, \end{align}

which is obtained by taking the curl of the Reynolds-averaged momentum equation (Pope Reference Pope2000); the detailed derivation is skipped for brevity. Equation (3.7) contains the key terms, including vortex stretching and tilting, baroclinic torque, thermal dilatation, buoyancy and vorticity diffusion, consistent with the conventional vorticity equation (Danaila, Kaplanski & Sazhin Reference Danaila, Kaplanski and Sazhin2021). However, a major distinction here is the emergence of the turbulent eddy viscosity $\nu _t$ , which accounts for the unsteady/turbulent effect on vorticity variation in the mean flow. Considering the difficulty in taking the curl of the Reynolds stress term, introducing the eddy viscosity concept allows us to interpret the vorticity change by Reynolds stress as a turbulent diffusion effect, which is of physical meaning; the estimation of $\nu _t$ will be discussed later.

At this point, it is interesting to note that when combining (3.7) and (3.6), the terms associated with $\omega _\theta (\boldsymbol {\nabla \cdot u} )$ and $u_r\omega _\theta /r$ in (3.6) exactly cancel with the thermal dilatation term and vortex stretching terms in $\mathrm {d}\omega _\theta /\mathrm {d}t$ . As a result, in the limit $\Delta z \rightarrow 0$ , (3.6) may take the simplified form

(3.8) \begin{align} \begin{split} \int _{\delta _1}^{\delta _2} \frac {\mathrm {d} \omega _\theta ^*}{\mathrm {d}t}\,\mathrm {d}r =\frac {\mathrm {d}}{\mathrm {d}z}\int _{\delta _1}^{\delta _2} u_z \omega _\theta \,\mathrm {d}r. \end{split} \end{align}

where ${\mathrm {d}\omega _\theta ^*}/{\mathrm {d}t}$ denotes (3.7) void of the thermal dilatation and vortex stretching terms. Physically, $\textstyle\int _{\delta _1}^{\delta _2} u_z \omega _\theta \,\mathrm {d}r$ in (3.8) can be interpreted as a radial vortex force applied on the total vorticity across the shear layer. Then (3.8) implies that the vorticity transport in the mean shear layer gives rise to the axial variation of the shear layer’s radial vortex force.

Recall that the goal of this study is to understand the axial variation of the shear layer intensity quantified by ${\mathrm {d}\gamma }/{\mathrm {d}z}$ . We next introduce a characteristic axial velocity associated with the above vortex force:

(3.9) \begin{equation} u_z^*=\frac {\int _{\delta _1}^{\delta _2}u_z\omega _\theta \,\mathrm {d}r}{\int _{\delta _1}^{\delta _2}\omega _\theta \,\mathrm {d}r}. \end{equation}

So (3.8) can be derived as

(3.10) \begin{align} \begin{split} \int _{\delta _1}^{\delta _2} \frac {\mathrm {d} \omega _\theta ^*}{\mathrm {d}t}\,\mathrm {d}r =\frac {\mathrm {d}(u^*_z\gamma )}{\mathrm {d}z}=\gamma\, \frac {\mathrm {d}u^*_z}{\mathrm {d}z}+u^*_z\,\frac {\mathrm {d\gamma }}{\mathrm {d}z}, \end{split} \end{align}

further yielding the calculation of ${\mathrm {d}\gamma }/{\mathrm {d}z}$ in the form

(3.11) \begin{equation} \frac {\mathrm {d}\gamma }{\mathrm {d}z} = \frac {1}{ u_z^*}\int _{\delta _1}^{\delta _2} \left (\frac {\mathrm {d}\omega _\theta ^*}{\mathrm {d}t} - \omega _\theta \frac {\mathrm {d}u^*_z}{\mathrm {d}z} \right ) \mathrm {d}r. \end{equation}

This establishes a theoretical model for interpreting the streamwise variation of the shear layer intensity. Equation (3.11) confirms that $\mathrm {d}\gamma /\mathrm {d}z$ is closely related to vorticity generation in the mean shear layer. Note that the second term in the integral of (3.11) is the vorticity generation owing to the axial gradient of the shear layer’s characteristic axial velocity, which can be understood physically as a thermal dilatation effect associated with the shear layer.

Before proceeding to model validation, we note that most previous works have considered $\mathrm {d}\omega _\theta /\mathrm {d}t$ to study the Lagrangian vorticity transport in an axisymmetric shear layer. However, in the Eulerian analysis of the shear layer’s intensity variation, our work reveals that the original thermal dilatation term should be replaced by a shear-layer thermal dilatation, and the vortex stretching effect also diminishes.

3.2. Estimation of key parameters

In practice, the model implementation of (3.11) depends on the calculation of ${\mathrm {d}\omega _\theta ^*}/{\mathrm {d}t}$ , which involves the estimations of pressure gradient $\boldsymbol{\nabla} {p}$ , $\nu _t$ and $\rho$ . Here, $\boldsymbol{\nabla} {p}$ is estimated based on the time-averaged momentum equations (Durst & Arnold Reference Durst and Arnold2008; Vanderkindere et al. Reference Vanderkindere, Laskari, Ganapathisubramani and Dekat2019) as

(3.12) \begin{equation} \boldsymbol{\nabla} {p} = -\rho [-(\bar{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla})\bar{\boldsymbol{u}}+\boldsymbol{\nabla} \boldsymbol{\cdot} \overline {\boldsymbol{u}^\prime \boldsymbol {u}^\prime }-\nu \,\nabla ^2\bar{\boldsymbol{u}}]. \end{equation}

For $\nu _t$ , we adopt the Boussinesq relation (Pope Reference Pope2000) to express the Reynolds stress as

(3.13) \begin{equation} \overline {\boldsymbol{u}^\prime\boldsymbol{u}^\prime} = \frac {2}{3}\bar {k}\,\boldsymbol {I}-2\nu _t\bar {\boldsymbol {S}}, \end{equation}

where $\bar {\boldsymbol {S}}$ is the mean rate-of-strain tensor, $\bar {k}$ is the TKE, and $\boldsymbol {I}$ is the identity tensor. Experimentally, the six different components of the symmetric Reynolds stress tensor $\overline {\boldsymbol{u}^\prime\boldsymbol{u}^\prime}$ can be computed based on the velocity fluctuation data. By taking the Frobenius product ( $\boldsymbol {:}$ ) between (3.13) and $\bar {\boldsymbol {S}}$ (Tammisola & Juniper Reference Tammisola and Juniper2016), we can obtain the eddy viscosity as

(3.14) \begin{equation} \nu _t = \frac {1}{2}\left [\frac {\left (\frac {2}{3}\bar {k}\boldsymbol {I}-\overline {\boldsymbol{u}^\prime\boldsymbol{u}^\prime}\right )\boldsymbol {:}\bar {\boldsymbol {S}}}{\bar {\boldsymbol {S}}\boldsymbol {:}\bar {\boldsymbol {S}}}\right ]. \end{equation}

Figure 4(a) plots the $\nu _t$ field for the same case of $Re=3600$ and $\phi=0.6$ , showing that $\nu _t$ is relatively low in the near field, but is significantly enhanced in the far field. This result is consistent with the flow unsteadiness evident from the transient measurements – the flow and flame exhibit quasi-steady features in the near field, and gradually evolve into a periodically unsteady state in downstream locations near the OSL. Figure 4(b) plots the radial variation of $\nu _t$ at $z=30$ mm together with the error bars, with maximum uncertainty $4.4\, \%$ . The error is estimated from the propagation of uncertainties of $\bar {k}$ , $\bar {\boldsymbol {S}}$ and $\overline {\boldsymbol{u}^\prime\boldsymbol{u}^\prime}$ in (3.14), based on the measured fluctuating velocities. Figure 4(c) further plots the corresponding distribution of the ratio between the eddy viscosity $\nu _t$ and the molecular viscosity $\nu$ . The relatively small $\nu _t$ confirms that the shear layers are essentially laminar in the near field, while the vortex shedding induces a weak turbulence downstream. In this study, $\nu$ is approximated using the Sutherland’s law for air (White & Majdalani Reference White and Majdalani2006); the error of this approximation is estimated to be smaller than $2.53\, \%$ , through the comparison with the component-weighted viscosity of the burnt gas mixture at 1700 K.

For the present weak-turbulence mild-swirling flow, (3.14) provides an overall reasonable and robust estimation of the eddy viscosity. However, its application to strong turbulence or high-swirling flow could be limited. For the former, the measured fluctuating velocity might not adapt to the rapid change of local velocity, resulting in an underestimated $\nu _t$ . For the latter, flow rotation could break the alignment assumption between the mean shear rate and the turbulent stress tensor in the Boussinesq relation (Pope Reference Pope2000).

Figure 4.

(a) Distribution of the estimated $\nu _t$ . The dashed red lines indicate the boundaries of the ISL and OSL. (b) Radial variation of $\nu _t$ at $z=30$ mm (marked by the yellow dashed line in (a)) and error bars. (c) Distribution of $\nu _t/\nu$ .

Figure 5.

(a) Time-averaged particle image. The red line represents the flame front extracted from the mean OH* field. (b) Contour of the radial gradient of particle density. (c) Distribution of the estimated gas density. The red dashed lines in (b) and green dashed lines in (c) represent the $1\, \%$ peak contours of the particle-density gradient. (d) Radial variations of the non-dimensional particle density (NPD), temperature $T$ , and density $\rho$ , with error bars at $z=20$ mm (marked by the yellow dashed line in (c)).

To estimate $\rho$ , we take inspiration from the approaches of Kazbekov, Kumashiro & Steinberg (Reference Kazbekov, Kumashiro and Steinberg2019) and Douglas et al. (Reference Douglas, Emerson, Hemchandra and Lieuwen2021). As the first step, the entire field is divided into three regions – burned, unburned and transition – where the densities are estimated separately based on the different temperatures using the ideal gas law. The temperature of the burned region is calculated using our previously proposed heat loss model (Wang et al. Reference Wang, Ren, Gu, Yu, Feng, Qi and Xia2023) for unconfined swirling flames, whereas the unburned region simply has the ambient temperature. The transition region essentially corresponds to a narrow preheat zone next to the flame front, where the temperature distribution is approximated using cubic spline interpolation. To this end, it is crucial to identify the transition region. From the mean particle density distribution in figure 5(a), the boundary line corresponding to the low-to-high density transition clearly overlaps with the flame front, which can be explained by the thermal diffusion effect that renders the flame front a barrier to the tracer particles. This provides the rationale to correlate the preheat zone to the region of dramatic particle density variation, based on the $1\, \%$ peak contour of the particle density gradient in figure 5(b). Accordingly, the resultant density distribution is plotted in figure 5(c). Figure 5(d) further presents the radial variations of the non-dimensional particle density (NPD), temperature, and density at the transverse plane $z=20\ \rm mm$ , with the error bars indicating relatively small uncertainties (less than $4.9\, \%$ ) in $\rho$ estimation. It should be noted that the correlation between the gradients of temperature and particle density becomes questionable in the far field, as the vortex-induced mixing dilutes the particle density near the flame front, which is evident from figure 5(a). Therefore, the result is plotted only up to $z=30\ \rm mm$ in figure 5(c).

3.3. Model validation and discussion

For model validation, figures 6(a,b) compare the predicted ${\rm d}\gamma /{\rm d}z$ of the ISL and OSL from (3.11) with that calculated directly from the experimental data in figure 2(e). Note that $\omega _\theta$ in the cylindrical coordinates has opposite sign to that in the Cartesian coordinates, so $\gamma$ is negative for the ISL and positive for the OSL. The results show promising agreement between the experimental calculation and the model prediction. In particular, the inflection point ${\rm d}\gamma /{\rm d}z=0$ for the ISL is predicted accurately, well capturing the observed non-monotonic variation trend of $\gamma$ . Similarly, the model predicts ${\rm d}\gamma /{\rm d}z\approx 0$ in the near field of the OSL, which also agrees with the observation of the near-constant $\gamma$ when $z\lt15$ mm. The $\gamma$ of the OSL starts to decrease further downstream, the rate at which is also effectively captured by the model. These verify the correctness of the theoretical derivation as well as the rationality of the parameter estimations.

Figure 6.

Comparison between ${\rm d}\gamma /{\rm d}z$ calculated directly from figure 2(e) and the prediction from (3.11) for (a) the ISL and (b) the OSL.

Figure 7.

Mean fields of (a) the total rate of change of vorticity and the four contributing terms, including (b) vortex tilting, (c) baroclinic torque, (d) vorticity diffusion and (e) shear-layer dilatation. The dashed red lines are the $\omega _\theta=0$ contour representing the ISL and OSL boundaries. Variations of the radially integrated contributing terms (3.15) in (f) the ISL and (g) the OSL along the $z$ axis. Here, $\dot {\gamma }_t$ , $\dot {\gamma }_b$ , $\dot {\gamma }_d$ and $\dot {\gamma }_s$ represent the vortex tilting, baroclinic torque, vorticity diffusion terms and shear-layer dilatation, respectively.

Next, we access the contributions of the different terms of (3.11) to the ISL and OSL intensities. Figures 7(a–e) plot the spatial distributions of the total rate of change of vorticity together with the four generation terms, including vortex tilting, baroclinic torque, vorticity diffusion and shear-layer dilatation, where the buoyancy term is neglected considering that the Froude number ( ${\sim}15$ ) is much larger than unity. In general, the vortex tilting is important primarily in the near field of the ISL where the swirl intensity is the strongest. The baroclinic torque contributes only to increasing the ISL intensity as it is related to the sharp density gradient across the flame surface. And the vorticity diffusion always serves as a decay factor. The roles of vortex tilting, baroclinic torque and vorticity diffusion in tuning the intensities of the shear layers, observed from figures 7(b–d), are largely consistent with the Konle & Sattelmayer (Reference Konle and Sattelmayer2009) numerical studies and Kazbekov & Steinberg (Reference Kazbekov and Steinberg2021) experimental studies on swirling flames. However, our analysis shows that the vortex stretching has no effect on the shear-layer intensity, and the thermal dilatation effect can be recast into the shear-layer dilatation. The latter tends to cause the ISL intensity to increase slightly in the near field and then to decrease in the far field, as shown in figure 7(e), whereas its effect on the OSL intensity is exactly the opposite. This distinction likely results from the flow reversal in the CRZ, which causes $u_z^*$ to first decrease and then increase along the ISL.

The four generation terms can be further integrated across the shear layer as

(3.15) \begin{align} \begin{split} \begin{cases} &\dot {\gamma }_t=\displaystyle\int _{\delta _1}^{\delta _2}\left ({\omega _r\,\dfrac {\partial u_\theta }{\partial r}+\omega _z\,\dfrac {\partial u_\theta }{\partial z}}\right )\mathrm {d}r,\\[3pt] &\dot {\gamma }_b=\displaystyle\int _{\delta _1}^{\delta _2}{\dfrac {1}{\rho ^2}\left (\dfrac {\partial \rho }{\partial z}\,\dfrac {\partial p}{\partial r}-\dfrac {\partial \rho }{\partial r}\,\dfrac {\partial p}{\partial z}\right )}\mathrm {d}r,\\[3pt] &\dot {\gamma }_d=\displaystyle\int _{\delta _1}^{\delta _2}{\left (\nu +\nu _t\right )\left (\dfrac {\partial ^2 \omega _\theta }{\partial r^2}+\dfrac {1}{r}\,\dfrac {\partial \omega _\theta }{\partial r}+\dfrac {\partial ^2 \omega _\theta }{\partial z^2}-\dfrac {\omega _\theta }{r^2}\right )}\mathrm {d}r,\\[10pt] &\dot {\gamma }_s=\displaystyle\int _{\delta _1}^{\delta _2}{\left (- \omega _\theta \,\dfrac {\mathrm {d}u^*_z}{\mathrm {d}z}\right )}\mathrm {d}r. \end{cases} \end{split} \end{align}

So the total rate of change of the shear layer intensity $\dot {\gamma }$ can be estimated as $\dot {\gamma }=\dot {\gamma }_t+\dot {\gamma }_b+\dot {\gamma }_d+\dot {\gamma }_s$ . The axial variations of $\dot {\gamma }$ , $\dot {\gamma }_t$ , $\dot {\gamma }_b$ , $\dot {\gamma }_d$ and $\dot {\gamma }_s$ are plotted in figures 7(f,g) for the ISL and OSL, respectively. Figure 7(f) shows that $\dot {\gamma }$ is negative in the ISL, which contributes to enhanced shear-layer intensity in the near field ( $z\lt15$ mm). This enhancement is mainly attributed to the combined effects of vortex tilting and baroclinic torque, which outweighs that of vorticity diffusion. In the far field, the ISL intensity begins to decrease ( $\dot {\gamma } \gt 0$ ) owing to the dominance of the vorticity diffusion and shear-layer dilatation effects. As for the OSL, figure 7(g) shows that $\dot {\gamma }$ is negative in the far field, indicating an intensity decay that results from the vorticity diffusion prevailing over the shear-layer dilatation. In contrast, the vorticity generation is weak, and the OSL intensity barely changes in the near field. This can be understood as the OSL being away from the flame root so the gradients of velocity and density are relatively small.