Methane/Air Flame Control in Non-Premixed Bluff Body Burners Using Ring-Type Plasma Actuators

Abstract

Enhancing the combustion efficiency and flame stability in conventional systems is essential for reducing carbon emissions and advancing sustainable energy solutions. In this context, electrohydrodynamic plasma actuators offer a promising active control method for modifying and regulating flame characteristics. This study presents a numerical investigation into the effects of a ring-type plasma actuator positioned on the co-flow air side of a non-premixed turbulent methane/air combustion system—an approach not previously reported in the literature. The ring-type plasma actuator was designed by placing electrodes along the perimeter of the small diameter wall of the air duct. The impact of the plasma actuator on the reacting flow field within the burner was analyzed, with a focus on its influence on the flow dynamics and flame structure. The results, visualized through velocity and temperature contours, as well as flow streamlines, provide insight into the actuator’s effect on flame behavior. Two operating modes of the plasma actuators were evaluated: co-flow mode, where the aerodynamic effect of the plasma actuators was directed downstream; and counter-flow mode, where the effects were directed upstream. The findings indicate that the co-flow actuation positively reduces the flame height and enhances the flame anchoring at the root, whereas counter-flow actuation slightly weakens the flame root. Numerical simulations further revealed that co-flow actuation marginally increases the energy release by approximately 0.13%, while counter-flow actuation reduces the energy release by around 7.8%.

1. Introduction

The urgent need to combat global warming and climate change has driven global efforts towards significantly reducing atmospheric emissions. A major contributor to these emissions is the combustion of hydrocarbon fossil fuels, which releases substantial quantities of carbon dioxide and other greenhouse gases into the atmosphere. Emissions from combustion systems not only exacerbate environmental problems but also pose challenges to energy efficiency and sustainability. To mitigate these negative effects, the development of more efficient combustion systems and the exploration of alternative green fuels are critical.

One of the key challenges in improving combustion efficiency lies in addressing the inherent inefficiencies associated with the turbulent mixing of fuel and oxidant, particularly in non-premixed combustion systems. Poor mixing can lead to incomplete combustion, greater emissions of pollutants, and reduced energy output. To overcome these challenges, various flow control techniques have been explored, broadly categorized into passive and active methods.

Active flow control techniques have gained considerable attention due to their ability to dynamically adapt and enhance combustion processes in real time. Among these, surface dielectric barrier discharge (DBD) plasma actuators have emerged as a promising tool for flame and combustion control [1,2]. These actuators are advantageous due to their lack of mechanical components, rapid response times, low energy consumption rates, and ability to be applied locally where needed [3].

DBD plasma actuators consist of two electrodes—one exposed and one grounded or embedded—separated by a dielectric material. These electrodes are placed asymmetrically on either side of the dielectric layer, with a high-voltage waveform typically applied to the exposed electrode while the grounded electrode remains at electrical ground level. When a sufficiently high voltage with the appropriate frequency is applied, the air adjacent to the actuator surface becomes ionized. This ionized region experiences a force due to the electric field, which accelerates the ionized air, generating what is known as an “ionic wind” [4].

The application of plasma-induced ionic wind for flow control has been explored in a variety of contexts, from aerospace to combustion systems. A comprehensive review of the application, physics, and modelling of the DBD plasma actuators for flow control was provided by Corke et al. [5]. Chen and Liao [2] experimentally investigated flow interactions downstream of a bluff body burner equipped with an annular plasma actuator. Qingya et al. [6] conducted an experimental study on the influence of a DBD plasma actuator on diffusion flame combustion characteristics. Wang and Roy [7,8] performed numerical simulations of the combustion stabilization and manipulation of the recirculation zones inside a gas turbine combustor using serpentine plasma actuators. The geometry and operating parameters of the plasma actuator, such as the electrode configuration and applied voltage, can be tailored to optimize the flow control effects for specific applications. In combustion systems, particularly those involving non-premixed flames, the ability to manipulate the flow field and enhance the mixing using plasma actuators holds significant promise for improving the combustion efficiency and reducing emissions. Li et al. [9,10] developed a plasma swirl injector with helical electrodes and compared its performance to a design using a straight electrode configuration. Their results demonstrated that the new design can adjust both the axial and azimuthal velocities and can be used to mitigate the flame blow-off or flashback. Furthermore, Li and Jiang [11] investigated the effects of electrical parameters on plasma swirler characteristics using both experimental and numerical methods. They observed that the swirl number increases almost linearly with the electrode length. While the literature highlights the potential of plasma actuators to influence flame and combustion control, the challenge remains to optimize their operating conditions and configurations to maximize their effectiveness.

Various models with differing levels of complexity have been proposed for modelling the influence of DBD plasma actuators [12]. Among these, the plasma fluid model and particle-in-cell model are highly sophisticated, offering detailed insights into physics of DBD plasma actuators. Shang and Huang [13] employed the drift–diffusion model to study plasma flow control and analyzed the effects of Joule heating, periodic electrostatic force, and Lorentz acceleration. The plasma fluid models generally require significant computational resources due to the large difference between the spatial and temporal scales of the flow and the plasma. Consequently, their application to large-scale problems is limited. In such cases, phenomenological or semi-empirical models of the plasma actuators provide viable alternatives. For example, Shyy et al. [14] proposed a simplified linearized body force model to represent the effects of a plasma actuator. Orlov et al. [15] developed a lumped-element circuit model to estimate the plasma body force, accounting for the spatial and temporal variations of the discharge.

Despite these advancements, there has been limited numerical research on the influence of DBD plasma actuators on combustion and flame control. Specifically, no numerical studies have analyzed the effects of a ring plasma actuator positioned in the co-flow stream of a non-premixed bluff body burner. This gap presents an opportunity for the further exploration of DBD plasma actuators in flame and combustion control. In this study, the potential of a ring-type plasma actuator to induce controlled perturbations in the flow field was investigated, focusing on its ability to alter the recirculation zones and flame shape. The geometry of the ring plasma actuator is straightforward, consisting of electrodes placed along the periphery of the air duct in a circular or ring configuration. This simple design enables two significant operating modes, co-flow and contour flow actuation, achieved by adjusting the positions of the exposed and embedded electrodes. The objective of the present study is to assess the efficacy of this configuration in achieving more stable flame and efficient combustion.

2. Methodology

2.1. Governing Equations of Fluid Dynamics

The flow was assumed to be steady-state, incompressible, viscous, and single-phase, with gravity effects neglected. Additionally, only the aerodynamic effects of the surface DBD plasma actuator were considered. Under these assumptions, the governing equations for continuity, momentum, and energy are as follows [16,17]:

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}\left(\rho u_i\right)=0$$

(1)

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left(\rho u_i{u}_j\right)=-{\displaystyle \frac{\mathsf{\partial }p}{\mathsf{\partial }{x}_i}}+{\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left[\mu \left({\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_j}}+{\displaystyle \frac{\mathsf{\partial }{u}_j}{\mathsf{\partial }{x}_i}}-{\displaystyle \frac{2}{3}}{\delta }_{ij}{\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_i}}\right)\right]+{\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left(-\rho u_i^{\prime }{u}_j^{\prime }\right)+f_{b,i}$$

(2)

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}\left[u_i\left(\rho E+p\right)\right]={\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left(K_{eff}{\displaystyle \frac{\mathsf{\partial }T}{\mathsf{\partial }{x}_j}}-{\displaystyle \sum }_k{h}_k{J}_{k,j}\right)+S_h+S_r$$

(3)

In Equation (2), $f_{b,i}$ represents the volumetric force caused by the plasma actuator. Here, only the plasma body force is considered, as the production of the active or radical species and the heat dissipation by the plasma actuator is assumed to be negligible due to the small size of the plasma discharge region. The term $J_{k,j}$ in Equation (3) represents the diffusion flux of species k, expressed as:

$$J_{k,j}=-\left(\rho D_{m,k}+{\displaystyle \frac{{\mu }_t}{S{c}_t}}\right) {\displaystyle \frac{\mathsf{\partial }{Y}_k}{\mathsf{\partial }{x}_j}}-D_{T,k}{\displaystyle \frac{1}{T}}{\displaystyle \frac{\mathsf{\partial }T}{\mathsf{\partial }{x}_j}}$$

(4)

where $S{c}_t$ is the turbulent Schmidt number, ${\mu }_t$ is the turbulent viscosity, $D_{m,k}$ is the molecular diffusivity, and $D_{T,k}$ is the thermal diffusivity. The term $S_h$ in Equation (3) represents the energy source due to chemical reactions and is defined as:

$$S_h=-{\displaystyle \sum }_k{\displaystyle \frac{h_k^0}{M_{w,k}}}{R}_k$$

(5)

where $h_k^0$ is the enthalpy of the formation, $M_{w,k}$ is the molecular weight of species k, and $R_k$ is the reaction rate of species k. The term $S_r$ represents the heat losses due to thermal radiation. The discrete ordinates (DO) radiation model is used to estimate these losses. The radiative transfer equation (RTE) for this model is expressed as [18]:

$$\nabla ·\left[I\left(\overrightarrow{r},\overrightarrow{s}\right)\overrightarrow{s}\right]+\left(\alpha +{\sigma }_s\right)I\left(\overrightarrow{r},\overrightarrow{s}\right)=\alpha n^2{\displaystyle \frac{\sigma T^4}{\pi }}+{\displaystyle \frac{{\sigma }_s}{4\pi }}{{\displaystyle \int }}_0^{4\pi }I\left(\overrightarrow{r},\overrightarrow{s}\right)d{\Omega }^{\prime }$$

(6)

In the above equation, $I$ is the radiation intensity, $\overrightarrow{r}$ is the position vector, $\overrightarrow{s}$ is the direction vector, $\alpha$ is the absorption coefficient, ${\sigma }_s$ is the scattering coefficient, $\sigma$ is the Stefan–Boltzmann constant ($5.669\times {10}^{-8} \mathrm{W}/({\mathrm{m}}^2·{\mathrm{K}}^4$)), $n$ is the refractive index, and ${\Omega }^{\prime }$ is the solid angle. In this work, the absorption coefficient is calculated using the weighted sum of gray gases model (WSGGM). Using the DO radiation model, the radiative heat source $S_r$ is estimated as:

$$S_R=\underset{A}{\int }\left(\int I\overrightarrow{s}·\overrightarrow{n} d{\Omega }^{\prime }\right)dA$$

(7)

In the above equation, $\overrightarrow{n}$ is the unit normal vector to the area $A$. In Equation (2), the Reynolds averaging method for simulating turbulent flows requires modelling of the Reynolds stresses $-\rho u_i^{\prime }{u}_j^{\prime }$. These stresses are commonly expressed in terms of time-averaged velocity gradients using the Boussinesq hypothesis [19]:

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left({\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_j}}+{\displaystyle \frac{\mathsf{\partial }{u}_j}{\mathsf{\partial }{x}_i}}\right)-{\displaystyle \frac{2}{3}}\left(\rho k+{\mu }_t{\displaystyle \frac{\mathsf{\partial }{u}_k}{\mathsf{\partial }{x}_k}}\right){\delta }_{ij}$$

(8)

Similarly, in Equation (3), the turbulent heat flux term $-\rho {u^{\prime }}_i{T}^{\prime }$ can be modelled analogously to the Reynolds stresses. The heat flux due to turbulence is expressed using the turbulent thermal conductivity coefficient $K_t$ caused by the turbulent mixing by establishing an analogy similar to Fourier’s law as follows:

$$-\rho {u^{\prime }}_i{T}^{\prime }=-K_t{\displaystyle \frac{\mathsf{\partial }T}{\mathsf{\partial }{x}_j}}$$

(9)

Using the above relationship, the total heat flux in the energy equation is defined as a combination of the heat flux caused by the mixing and molecular interactions using the effective thermal conductivity coefficient $K_{eff}=K+K_t$. A common way to calculate $K_t$ is to use the turbulent Prandtl number:

$$K_t={\displaystyle \frac{c_p{\mu }_t}{P{r}_t}}$$

(10)

In this study, the standard $k$-$\epsilon$ model was used. After examining alternative RANS turbulence models, the standard $k$-$\epsilon$ model was chosen for its balance between accuracy and computational efficiency. In this model, two transport equations (Equations (11) and (12)) are solved for the kinetic energy of the turbulence $k$ and dissipation rate of the turbulent kinetic energy $\epsilon$ [19]:

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}\left(\rho k{u}_i\right)={\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right){\displaystyle \frac{\mathsf{\partial }k}{\mathsf{\partial }{x}_j}}-{\displaystyle \frac{2}{3}}{\delta }_{ij}{\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_i}}\right]+G_k-\rho \epsilon$$

(11)

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}\left(\rho \epsilon u_i\right)={\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{\epsilon }}}\right){\displaystyle \frac{\mathsf{\partial }\epsilon }{\mathsf{\partial }{x}_j}}-{\displaystyle \frac{2}{3}}{\delta }_{ij}{\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_i}}\right]+C_{1\epsilon }{\displaystyle \frac{\epsilon }{k}}{G}_k-C_{2\epsilon }\rho {\displaystyle \frac{{\epsilon }^2}{k}}$$

(12)

In the above equations, $G_k$ represents the production of turbulent kinetic energy and is directly linked to the Reynolds stresses and the mean velocity gradients as:

$$G_k=-\overline{u_i^{\prime }{u}_j^{\prime }}{\displaystyle \frac{\mathsf{\partial }{u}_j}{\mathsf{\partial }{x}_i}}$$

(13)

By substituting the Boussinesq assumption for the Reynolds stresses, the production of turbulent kinetic energy can be expressed as:

$$G_k={\mu }_t{\tilde{S}}^2$$

(14)

$$\tilde{S}=\sqrt{2S_{ij}{S}_{ij}} with S_{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathsf{\partial }{u}_i}{\mathsf{\partial }{x}_j}}+{\displaystyle \frac{\mathsf{\partial }{u}_j}{\mathsf{\partial }{x}_i}}\right)$$

(15)

According to the turbulence model, the following relationship holds for ${\mu }_t$:

$${\mu }_t=\rho C_{\mu }{\displaystyle \frac{k^2}{\epsilon }}$$

(16)

The constants of the standard $k$-$\epsilon$ model based on the work of Dally et al. [20] and Hossain and Malalasekera [21] are as follows:

$$C_{\mu }=0.09 C_{1\epsilon }=1.6 C_{2\epsilon }=1.92 {\sigma }_k=1.0 {\sigma }_{\epsilon }=1.3$$

(17)

2.2. Combustion Modeling

The combustion process of methane and air is modelled using a simplified two-step chemical reaction:

$$\begin{gathered} Step 1\text{:} C{H}_4+1.5O_2\to CO+2H_{2}O \\ Step 2\text{:} CO+0.5O_2\to C{O}_2 \end{gathered}$$

(18)

These reactions are expressed in the general form:

$${\displaystyle \sum }_{k=1}^N{\upsilon }_{k,r}^{\prime }{\chi }_k\to {\displaystyle \sum }_{k=1}^N{\upsilon }_{k,r}^{″}{\chi }_k$$

(19)

where ${\upsilon }_{k,r}^{\prime }$ and ${\upsilon }_{k,r}^{″}$ are the stoichiometric coefficients for reactants and products, respectively, and ${\chi }_k$ represents the chemical species. The transport of species in this model is governed by the following equation:

$${\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}\left(\rho u_i{Y}_k\right)=-{\displaystyle \frac{\mathsf{\partial }}{\mathsf{\partial }{x}_i}}{J}_{k,i}+R_k,\hspace{1em}\hspace{1em}k=C{H}_4, O_2, CO, C{O}_2, H_{2}O$$

(20)

where $Y_k$ is the mass fraction of the k-th species, $J_{k,i}$ is the diffusion flux, and $R_k$ is the net rate of production or consumption. It is assumed that the chemical kinetics is significantly faster than turbulent mixing, so the reaction rate is controlled by turbulent mixing. This scenario is common in non-premixed flames with fast-burning fuel, where turbulence dictates the mixing of fuel and oxidizers into the reaction zones. The turbulence–chemistry interaction is modelled using the eddy dissipation model, where the net rate of production or consumption of $k$-th species in the $r$-th reaction is [22]:

$${\mathrm{R}}_{k,r}=\min \left({\upsilon }_{k,r}^{\prime }{M}_{w,k}A\rho {\displaystyle \frac{\epsilon }{k}}mi{n}_R\left({\displaystyle \frac{Y_R}{{\upsilon }_{R,r}^{\prime }{M}_{w,R}}}\right), {\upsilon }_{k,r}^{\prime }{M}_{w,k}AB\rho {\displaystyle \frac{\epsilon }{k}}\left({\displaystyle \frac{{{\displaystyle \sum }}_P{Y}_R}{{{\displaystyle \sum }}_j^N{\upsilon }_{j,r}^{\prime }{M}_{w,j}}}\right)\right)$$

(21)

In the above equation, $A$ and $B$ are empirical coefficients with values of 4 and 0.5, respectively. The latter relation assumes that the chemical reaction rate is a function of the large eddy mixing time scale value ${\displaystyle \frac{k}{\epsilon }}$.

2.3. Phenomenological Model of the Plasma Actuator

To model the flow control effect of the plasma actuator, a body force term is introduced into the momentum equation. This body force is calculated using the approach proposed by Suzen et al. [23]. In this model, two equations are solved for the electric potential ($\varphi$) and the charge density (${\rho }_c$), as follows:

$$\nabla · \left({ϵ}_r\nabla \varphi \right)=0$$

(22)

$$\nabla · \left({ϵ}_r\nabla {\rho }_c\right)={\displaystyle \frac{{\rho }_c}{{\lambda }_d^2}}$$

(23)

In the above equation, ${ϵ}_r$ represents the relative dielectric permittivity and ${\lambda }_d$ denotes the Debye length. These equations are solved using proper boundary conditions for the electric potential and charge density. As shown in Figure 1, the alternating voltage is applied to the exposed electrode and the embedded electrode is assumed to be grounded. At the interface of the dielectric layer and the fluid zone above the embedded electrode, a time-varying spatial distribution is considered for the charge density. The applied voltage at the exposed electrode and the charge density distribution at the dielectric surface are expressed as:

$${\varphi }_{electrode}\left(t\right)={\varphi }_{max}f\left(t\right)$$

(24)

$${\rho }_{c,w}\left(x,t\right)={\rho }_{c,max}G\left(x\right)f\left(t\right)$$

(25)

In these equations, ${\rho }_{c,max}$ and ${\varphi }_{max}$ are the maximum values of the charge density and electric potential, respectively. The functions $f\left(t\right)$ and $G\left(x\right)$ describe the temporal and spatial distributions of the charge density, respectively. Suzen and Huang [23] proposed using a half-Gaussian profile for the spatial distribution of the charge density as follows:

$$G\left(x\right)=\exp \left[-{\displaystyle \frac{{\left(x-\xi \right)}^2}{2{\gamma }^2}}\right]$$

(26)

In the above function, $\gamma$ controls the rate of decay of the charge density, while $\xi$ determines the location of its maximum at the dielectric interface. The temporal variation of the charge density ($f\left(t\right)$) is assumed to match the wave form of the applied voltage. For a sinusoidal waveform:

$$f\left(t\right)=\sin \left(2\pi ft\right)$$

(27)

In the above, $f$ represents the frequency of the applied voltage. The body force term (force per unit volume) can then be calculated as:

$$\stackrel{\rightharpoonup }{f_b}={\rho }_c\times \stackrel{\rightharpoonup }{E}$$

(28)

where $\stackrel{\rightharpoonup }{E}$ is the electric field vector calculated as $\stackrel{\rightharpoonup }{E}=-\nabla \varphi$. Equations (22) and (23) are normalized using the maximum values of the applied voltage and charge density and along with the temporal profile $f\left(t\right)$ as follows:

$${\varphi }^{\ast }={\displaystyle \frac{\varphi }{{\varphi }_{max}f\left(t\right)}}, {\rho }_c^{\ast }={\displaystyle \frac{{\rho }_c}{{\rho }_{c,max }f\left(t\right)}}$$

(29)

Using the normalized definitions, the boundary conditions described in Equations (24) and (25) simplify to:

$${{\varphi }^{\ast }}_{electrode}=1, {\rho }_{c,w}^{\ast }\left(x\right)=G\left(x\right)$$

(30)

This formulation enables solving Equations (22) and (23) in a steady-state manner. The boundary conditions given by Equation (30) and normalized forms of the Equations (22) and (23) are shown schematically in Figure 1.

Examining the body force definition reveals that its magnitude primarily depends on ${\varphi }_{max}$ and ${\rho }_{c,max }$. These parameters can be adjusted to control the strength of the plasma body force. Experimental studies have shown that the plasma body force or the induced velocity exhibits a non-linear dependence on the applied voltage magnitude. To account for this non-linearity, the correct scaling of the charge density must be incorporated.

Following the work of Abdollahzadeh et al. [24,25], the scaling the charge density can be determined through analytical or semi-empirical estimations of the plasma-induced thrust. Assuming that the plasma is generated within a region characterized by a length $l_p$, height $h_p$, and width $L$, the body force magnitude can be approximated as:

$$\left|f_b\right|~{\displaystyle \frac{Thrust/L}{l_p{h}_p}}$$

(31)

Furthermore, since the plasma body force is dependent on ${\varphi }_{max}$ (or equivalently the maximum electric field $E_{max}$) and ${\rho }_{c,max }$, the charge density scaling is given by:

$${\rho }_{c,max }~{\displaystyle \frac{Thrust/L}{{\varphi }_{max}{l}_p{h}_p}}=\chi {\displaystyle \frac{Thrust/L}{E_{max}L{l}_p{h}_p}}$$

(32)

Here, $\chi$ is a proportionality constant that relates the thrust and electric field to the maximum charge density. The estimation of the ${\rho }_{c,max }$ requires knowledge of the thrust $T$, the maximum electric field $E_{max}$, and the constant $\chi$. Experimental evidence indicates that plasma is generated only when the applied voltage significantly exceeds the plasma ignition voltage $V_0$. This implies that the thrust and maximum charge density are functions of the difference (${\varphi }_{max}-V_0$). In other words, charged particle generation and accumulation occur only when the ${\varphi }_{max}>V_0$.

The plasma actuator can be modelled using an equivalent capacitance and the thrust it generates can be linked to the energy consumption of the actuator when plasma is formed ${\varphi }_{max}>V_0$. According to Yoon and Han [26], the plasma actuator can be represented by an equivalent circuit comprising two capacitors, $C_1$ and $C_2$, connected in series. The equivalent capacitance is defined as:

$$C_{eq}={\displaystyle \frac{C_1{C}_2}{C_1+C_2}}$$

(33)

$${\displaystyle \frac{C_1}{L}}={\displaystyle \frac{\pi {ϵ}_{r,d}{ϵ}_0}{\ln \frac{2t_d+\frac{t_e}{2}}{\frac{t_e}{2}}}}$$

(34)

$$\begin{gathered} {\displaystyle \frac{C_2}{L}}={\displaystyle \frac{{ϵ}_{0}A}{t_d}} \\ A={\displaystyle \frac{k_2{\left(\frac{f}{1000}\right)}^{\frac{1}{4}}\left(\frac{{\varphi }_{max}}{V_0}-1\right)}{\ln \frac{2t_d+\frac{t_e}{2}}{\frac{t_e}{2}}}}\left(\ln {\displaystyle \frac{\frac{t_e}{2}}{r_{pl}}}-\ln {\displaystyle \frac{\frac{t_e}{2}}{\sqrt{r_{pl}^2+{\left(2t_d+t_e\right)}^2}}}\right) \end{gathered}$$

(35)

In the above equation, $t_e$ and $t_d$ represent the thickness of the electrode and the dielectric layer, respectively; ${ϵ}_{r,d}$ also stands for the relative dielectric permittivity of the dielectric layer and ${ϵ}_0$ is the vacuum permittivity (${ϵ}_0=8.8541\times {10}^{-12} \mathrm{F}·{\mathrm{m}}^{-1}$); $k_2$ depends on the wave form. For an AC sinusoidal voltage, it will be equal to unity. For a saw tooth, it will have a value equal to 0.66, and for rectangular wave form, it will have a value equal to 2. For the equivalent capacitance of $C_{eq}$, Yoon and Han assumed that the thrust produced by the plasma actuator is proportional to its power consumption. They developed the following semi-empirical relation for the thrust per unit length produced by the actuator:

$${\displaystyle \frac{Thrust}{L}}=k_1{\left(2\pi f\right)}^{{\alpha }_1}{C}_{eff}{\left(1-{\displaystyle \frac{{\varphi }_{max}}{V}}\right)}^2{\varphi }_{max}^2$$

(36)

This equation provides a reliable estimate of the net thrust generated by the plasma actuator, which includes the combined effects of the plasma body force and shear stress acting on the fluid [27]. To account for the shear stress acting on the fluid near the surface of the actuator, an analogy can be considered between the induced plasma wall jet and a plane wall jet [28]. Previous studies in the literature have shown the validity of this analogy. For a steady, incompressible flow, the boundary layer solution for the plane wall jet provides the following expressions for the wall jet momentum flux $K$ (interpreted as thrust) and the wall shear stress [29,30]:

$$K=0.884{\left({\displaystyle \frac{\rho {\mathcal{F}}^3}{\upsilon x}}\right)}^{{\displaystyle \frac{1}{4}}}$$

(37)

$${\tau }_w=0.221\rho {\left({\displaystyle \frac{{\mathcal{F}}^3}{{\rho }^3\upsilon x^5}}\right)}^{{\displaystyle \frac{1}{4}}}$$

(38)

In the above, $\mathcal{F}$ is the wall jet constant, which is also called the “flux of momentum flux” for wall jets, which is a conserved quantity [31]. The net thrust can be assumed to be equal to K and the wall shear force per unit of length is estimated as:

$${\displaystyle \frac{F_{friction}}{L}}\cong 0.25K$$

(39)

Thus, the average plasma force per unit length can be approximated as:

$${\displaystyle \frac{F_{plasma}}{L}}={\displaystyle \frac{Thrust}{L}}+{\displaystyle \frac{F_{friction}}{L}}\cong 1.25{\displaystyle \frac{Thrust}{L}}$$

(40)

Using the above relation, the average body force can be scaled as:

$$\left|f_b\right|~{\displaystyle \frac{1.25k_1{\left(2\pi f\right)}^{{\alpha }_1}{C}_{eff}{\left(1-\frac{V_0}{{\varphi }_{max}}\right)}^2{\varphi }_{max}^2}{l_p{h}_p}}~{\rho }_{c,max }{\varphi }_{max}{\chi }_1{\chi }_2$$

(41)

In the above, ${\chi }_1$ and ${\chi }_2$ are related to the distribution of the charge density and electric field, respectively. Assuming that the maximum electric field depends on the voltage amplitude, ${\chi }_2$ is obtained automatically as a result of solution of Equation (22), while ${\chi }_1$ is derived from the solution of the charge density equation under the influence of the Gaussian profile for the normalized surface charge density. These parameters can be estimated as:

$${\chi }_1{\chi }_2~{\displaystyle \frac{\int {\rho }_{c}EdV}{l_p{h}_p}}$$

(42)

The maximum charge density ${\rho }_{c,max }$ is then given by:

$${\rho }_{c,max }~{\displaystyle \frac{1.25k_1{\left(2\pi f\right)}^{{\alpha }_1}{C}_{eff}{\left(1-\frac{V_0}{{\varphi }_{max}}\right)}^2{\varphi }_{max}}{l_p{h}_p{\chi }_1{\chi }_2}}$$

(43)

The approximate length and height of the plasma region can be estimated using the following relations proposed by Yoon et al.:

$$l_p~0.005{\left({\displaystyle \frac{k_1{\left(2\pi f\right)}^{{\alpha }_1}{C}_{eff}{\left({\varphi }_{max}-V_0\right)}^2}{0.015}}\right)}^{1/3}$$

(44)

$$h_p~0.0013{\left({\displaystyle \frac{k_1{\left(2\pi f\right)}^{{\alpha }_1}{C}_{eff}{\left({\varphi }_{max}-V_0\right)}^2}{0.015}}\right)}^{1/3}$$

(45)

3. Results and Discussion

3.1. Validation of Model

To validate the proposed model for the linear plasma actuator, we compared our numerical results with the experimental studies conducted by Thomas et al. [32] and Durscher and Roy [33]. The geometrical and operational parameters of these validation cases are summarized in

Table 1. Geometrical and operational parameters of the validation cases for the plasma actuator model.

Parameter	Value		Dimension
	Thomas et al. [32]	Durscher and Roy [33]
Dielectric layer thickness	6.35	3.0	mm
Electrode thickness	0.04	0.07	mm
Relative permittivity of the dielectric	2	3	-
Width of the exposed electrode	12.7	5.0	mm
Width of the embedded electrode	25.4	20.0	mm
Asymmetric gap between electrodes	12	0	mm
Operating voltage	40	20	KVpp
Operating Frequency	2000	7000, 14,000	Hz
Voltage shape	Sawtooth	Sinusoidal	-

.

Figure 2 illustrates the variation of the net thrust generated by the plasma actuator as a function of the applied voltage, computed using Equation (36). The numerical results are compared with the experimental data from Thomas et al. [32] and Durscher and Roy [33]. The comparison encompasses various dielectric materials, applied frequencies, dielectric thicknesses, and voltage waveforms, demonstrating the general validity of Equation (36) for estimating the net thrust produced by the plasma actuator.

In Figure 3, the velocity profiles induced by the plasma actuator are numerically computed and compared against the experimental results of the studies by Thomas et al. and Durscher and Roy. The proposed model exhibits acceptable accuracy in predicting the induced velocities. It should be mentioned that for higher operating voltages, as in the case of the Thomas et al.’s study, the plasma-induced velocities are formed in such a way as to produce a turbulent boundary layer [34]. The same behaviour is also reported by Maden et al. [35]. Consequently, simulations corresponding to Durscher and Roy’s experiments were conducted under laminar flow conditions, whereas those replicating Thomas et al.’s setup assumed a turbulent flow state.

For validation of the base flow case in the burner under the non-reacting condition, the experimental results of Dally et al. [20], Tong et al. [36], and Caetano and da Silva [37] were considered. The validation cases involved comparing the axial and radial velocity profiles at various points within the burner. These cases span a wide range of fuel and oxidant velocities, as well as bluff body geometries, providing a robust framework for model validation. Figure 4 illustrates the comparisons of the axial and radial velocity profiles obtained from the numerical simulations with the corresponding experimental data. Consequently, the agreement observed in the velocity profiles presented in Figure 4 shows the general validity of the current numerical solution.

For the validation of the reacting flow regime inside the burner, we also considered the experimental results of Dally et al. [20]. For this case, the temperature and velocity profiles were compared with the experimental results, as shown in Figure 5. The simulated temperature profiles show excellent agreement with the experimental measurements, capturing both the peak flame temperatures and gradients. The axial and radial velocity profiles further validate the numerical model’s ability to replicate the flow dynamics in a reacting environment. The predicted flow acceleration and deceleration patterns near the flame and the recirculation regions closely match the experimental data.

3.2. Influence of Ring Plasma Actuator

To investigate the influence of the ring dielectric barrier discharge (DBD) plasma actuator on the reacting flow and flame structure, a non-premixed bluff body burner with the geometry schematically depicted in Figure 6 was considered. This geometry was adapted from the work of Manrique et al. [38]. The plasma actuator is placed in the co-flow air stream in the vicinity of the bluff body, as shown in Figure 6. The plasma actuator consists of an exposed electrode and an embedded electrode, each with a width of 15 cm. These electrodes are asymmetrically positioned with no gap between them, and they are separated by a dielectric layer with a thickness of 3 mm made of Kapton. The plasma actuator operates with a high-voltage sinusoidal waveform of 20 kV peak-to-peak (kVpp) at 20 kHz. Additionally, the position of the embedded electrode can be adjusted to generate an ionic wind that opposes the air co-flow. In Figure 6, the diameter of the bluff body $D_B$, the diameter of the co-flow air tube $D_{co\text{-}flow}$, the diameter of the central jet duct $D_j$, and the distance of the edge of the exposed electrode from the tip of the bluff body $\Delta x$ are 60 mm, 200 mm, 7.1 mm, and 10 mm, respectively.

At the inflow boundaries, fuel and air enter the computational domain with constant uniform velocities. A fuel velocity of 8 m/s and air velocity of 4 m/s were considered, corresponding to a fuel-to-air velocity ratio of 2. The mass fraction of the methane at the fuel inlet was set to unity, and the mass fraction of oxygen at the air inlet was set to 0.23. The turbulent intensity at the fuel inlet was 5.8%, while the turbulent intensity at the air inlet was 4.4%.

Due to the axisymmetric geometry of the burner equipped with the ring plasma actuator, a 2D axisymmetric computational domain was used. It is important to note that the configuration of the plasma actuator considered here will only induce axial and radial body force components. If the actuator were installed at any angle relative to the axis, a tangential force component would be induced and a 3D computational domain would need to be considered. A mesh was generated for the 2D axisymmetric computational domain using structured cells. The mesh was refined near the surface of the bluff body to ensure that the $y^{+}$ value remained below 1, ($y^{+}<1$). In addition, the mesh near the surface of the DBD plasma actuator was also refined such that the cell wall distance was smaller than the Debye length. A schematic of the computational grid is shown in Figure 7. A grid sensitivity analysis was performed by monitoring the variation of the maximum temperature in the computational domain for various grid sizes. The results are presented in Figure 8. As shown in this figure, the differences between predictions on different grids tend to diminish as the grid resolution increases, and a grid with about 100K cells was sufficient to provide mesh-independent results. The grid convergence index (GCI) was also computed for grid 1, grid 2, and grid 3, leading to a GCI₂₁ value of 1.44%. This indicates that the uncertainty in the predicted temperature field is approximately 1.44%.

The influence of the ring DBD plasma actuator on the reacting flow and flame structure inside the burner is shown in Figure 9 and Figure 10. Figure 9 presents the velocity contours superimposed with flow streamlines for three cases—without the plasma actuator, with the plasma actuator in co-flow mode, and with the plasma actuator in counter-flow mode. When the DBD plasma actuator is active in the co-flow mode (Figure 9b), the internal recirculation is slightly enhanced. When the DBD plasma actuator operates in the counter-flow mode (Figure 9c), a clockwise-rotating separation bubble forms upstream of the bluff body within the air duct, deflecting the airflow direction. This separation bubble acts as a virtual modification to the shape of the bluff body. It can be expected that with the changes in the location and the strength of the plasma actuator in the counter-flow mode, various modifications to the size of the bluff body can be achieved. In Figure 9c, the internal and external recirculation zones are weakened due to the effect of the DBD actuator operating in counter-flow mode.

Figure 10 qualitatively shows the flame structures for the cases with and without plasma actuators using temperature contours. In this figure, regions with high temperatures (e.g., >1000 K) can be assumed to represent the flame structure. When the plasma actuator operates in the co-flow mode, the flame root is slightly strengthened, suggesting enhanced combustion stability (Figure 10b). However, for the case that the plasma actuator is active in counter-flow mode (Figure 10c), the plasma actuator’s influence on the flame root is slightly negative, as the flame appears less anchored near the bluff body.

To further illustrate the temperature differences in the flames between various cases, Figure 11 presents a plot of the differences in temperature between the co-flow and no actuator cases and a plot of the differences in temperature between the counter-flow and no actuator cases. In Figure 11a, the temperature near the bluff body is notably enhanced for the co-flow case, indicating improved combustion efficiency in this configuration. In contract, Figure 11b shows that in the counter-flow mode, the temperature increase occurs only in the very narrow small zone near the surface of the bluff body, and moving downstream the temperature is reduced, consistent with the observed weakening of the flame root in this configuration.

Figure 12 presents the radial temperature profiles at x = 10 mm, 80 mm, and 250 mm (equivalent to x/D_b = 0.167, 1.34, and 4.167, respectively) downstream of the bluff body surface. The results show that when the plasma actuator is active in the co-flow configuration, the maximum temperature of the flame near the bluff body surface (x = 10 mm) increases. This position can be assumed as the flame root. When the plasma actuator is activated in the counter-flow mode, the maximum flame temperature is reduced. This observation is directly related to the evolution of the internal and external recirculation zones downstream of the bluff body. When the plasma actuator is in counter-flow mode, the internal recirculation is significantly weakened and the rate and extent to which the fuel and air are mixed is influenced. This result indicates that among the two plasma cases, the co-flow actuation is acting favourably. This behaviour can also be observed in Figure 12b,c, further downstream of the bluff body. In Figure 12a, the radial distance where the temperature is above 1000 K can be considered as the flame root thickness. The results in the Figure 12a show that the thickness of the flame is reduced when the plasma actuator is activated in the counter-flow mode. Overall, the results in Figure 12 indicate that the co-flow actuation enhances the flame structure, improving both the flame stability and root anchoring.

The primary objective of using plasma actuators is to alter the flame structure to achieve greater stability or enhance the combustion efficiency by increasing the energy released during combustion. For a bluff body burner, flame stabilization can be observed through a stronger, closer-to-root flame. The results in Figure 10, Figure 11 and Figure 12 indicate that this goal was partially achieved, particularly in the co-flow mode. To assess the energy released, the integral of the heat of the reaction source (Equation (5)) was calculated for all three cases—no actuator, co-flow actuation, and counter-flow actuation. The simulation results confirm that co-flow actuation slightly increases the energy release by approximately 0.13%, while the counter-flow mode reduces the energy release by around 7.8%.

4. Conclusions

In this paper, a numerical study of the influence of a ring DBD plasma actuator positioned on the annular air co-flow side of a non-premixed bluff body burner is presented. First, a novel modified phenomenological model for the plasma actuator was developed and validated. The plasma actuator was then analyzed in two operational modes—co-flow and counter-flow. The results showed that under the tested inflow conditions and plasma operating conditions (20 kVpp and 20 kHz), the activation of the plasma actuator in the co-flow mode strengthens both the internal and outer recirculation zones downstream of the bluff body, leading to an increase in the flame root thickness. However, when the plasma actuator is activated in counter-flow mode, the recirculation zones are modified in a slightly negative mode under the conditions tested here. In this mode, a separation bubble is formed within the annular air duct, virtually modifying the shape of the bluff body. For this case, the flame root shown by the temperature profiles is weakened. Analyzing the energy release through the heat of the reaction for various cases showed that the co-flow mode actuation slightly increased the energy release by 0.13%, whereas the contour-flow mode actuation reduced the energy release by 7.8%. The practical implications of these findings are significant for the design and optimization of advanced combustion systems. The ability to control the flame stability and energy release through strategic placement and the operational modes of plasma actuators can inform the development of more efficient and stable combustion devices.

It is important to note that the conclusions drawn in this study are based on the specific operational conditions tested. Nevertheless, better effects could potentially be achieved at higher applied voltages, as the plasma-induced thrust depends on the voltage and frequency. Future research studies should focus on exploring a broader range of operating conditions, such as varying the inflow conditions, applied voltage, and frequency, as well as optimizing the position of the plasma actuator to maximize its aerodynamic effects and for improved flame stabilization.