A spectral deferred correction strategy for low Mach number reacting flows subject to electric fields

TL;DR

This paper introduces a spectral deferred correction algorithm for low Mach number reacting flows with electric fields, enabling stable, second-order accurate simulations of complex electrochemical interactions.

Contribution

It extends the MISDC framework to include electric field effects and charged species transport with a Jacobian-Free Newton Krylov approach for implicit treatment.

Findings

Demonstrates second-order convergence in space and time.

Shows stability in steady and unsteady flow problems.

Designed for potential multidimensional applications.

Abstract

We propose an algorithm for low Mach number reacting flows subjected to electric field that includes the chemical production and transport of charged species. This work is an extension of a multi-implicit spectral deferred correction (MISDC) algorithm designed to advance the conservation equations in time at scales associated with advective transport. The fast and nontrivial interactions of electrons with the electric field are treated implicitly using a Jacobian-Free Newton Krylov approach for which a preconditioning strategy is developed. Within the MISDC framework, this enables a close and stable coupling of diffusion, reactions and dielectric relaxation terms with advective transport and is shown to exhibit second-order convergence in space and time. The algorithm is then applied to a series of steady and unsteady problems to demonstrate its capability and stability. Although developed in a one-dimensional case, the algorithmic ingredients are carefully designed to be amenable to multidimensional applications.

Tables

Table 1. Table 1: List of ions included in the chemical mechanism along with their molecular weight.

Cation	H3O+	HCO+	C2H3O+	CH5O+
$W_m$ [g/mol]	19.02	29.02	43.05	33.05
Anions	OH-	O-	O${}_{}{}^{-}{}_2{}^{}$	CO${}_{}{}^{-}{}_3{}^{}$	HCO${}_{}{}^{-}{}_2{}^{}$	HCO${}_{}{}^{-}{}_3{}^{}$
$W_m$ [g/mol]	17.01	16.00	32.00	60.01	44.01	61.02

Table 2. Table 2: Characteristic of the 1D laminar premixed flame

Operating conditions
Tin [K]	$U_{adv,in}$ [m/s]	Pressure [Pa]	Yfuel,in	YO2,in	YN2,in
350.0	0.371	101325.0	0.055	0.220	0.725
Numerical parameters
L [m]	$n_x$	${\Delta }_x$ [$\mu$m]	$k_{max}$	$\gamma$	${\gamma }_{MG}$
0.01	[128,2048]	[156,9.77]	4	1.0e-4	1.0e-4

Table 3. Table 3: Estimated relaxation time of electrons and H 3 O + for increasing external forcing intensity Δ ​ V Δ 𝑉 \Delta V .

Relaxation time [s]
$\Delta V$ = 100V		$\Delta V$ = 1000V		$\Delta V$ = 2500V
${\tau }_{e^{-}}$	${\tau }_{H_3{O}^{+}}$	${\tau }_{e^{-}}$	${\tau }_{H_3{O}^{+}}$	${\tau }_{e^{-}}$	${\tau }_{H_3{O}^{+}}$
2.5e-6	2.5e-4	2.5e-7	2.5e-5	1.0e-7	1.0e-5

Table 4. Table 4: Considered forcing frequencies and half periods for the AC cases along with ratios to the characteristic relaxation time of electrons and H 3 O + at A A ​ C = 1000 subscript 𝐴 𝐴 𝐶 1000 A_{AC}=1000 V.

Case	$f_{AC}$ [Hz]	${\tau }_{AC}/2$ [s]	${\tau }_{AC}/2{\tau }_{e^{-}}$	${\tau }_{AC}/2{\tau }_{H_3{O}^{+}}$
AC1k	1000	5e-4	2000	20
AC2.5k	2500	2e-4	800	8
AC5k	5000	1e-4	400	4
AC10k	10000	5e-5	200	2
AC25k	25000	2e-5	80	0.8
AC50k	50000	1e-5	40	0.4
AC100k	100000	5e-6	20	0.2

Equations

$$p(\mathbf{x},t)=p_0+\pi (\mathbf{x},t)$$

$$\frac{\partial (\rho Y_m)}{\partial t}+\nabla \cdot ({\text{emU}}_{adv}\rho Y_m)$$

$$\frac{\partial (n_e)}{\partial t}+\nabla \cdot ({\text{emU}}_{adv}{n}_e)$$

$$\frac{\partial (\rho h)}{\partial t}+\nabla \cdot ({\text{emU}}_{adv}\rho h)$$

$$\frac{\partial (\rho {\text{emU}}_{adv})}{\partial t}+\nabla \cdot (\rho {\text{emU}}_{adv}{\text{emU}}_{adv})$$

$$p_0=\rho \frac{R_u}{\overline{W}}T=\rho R_u\sum _m\frac{Y_m}{W_m}$$

$$\displaystyle\begin{split}\boldsymbol{\nabla}\cdot\boldsymbol{U}_{adv}&=\frac{1}{\rho c_{p}T}\bigg{(}\boldsymbol{\nabla}\cdot\lambda\boldsymbol{\nabla T}+\sum_{m}\mathbf{\Gamma}_{m}\cdot\boldsymbol{\nabla h_{m}}+\sum_{m+e}z_{m}Y_{m}\mathbf{\Gamma}_{m}\cdot\boldsymbol{E}\bigg{)}\\ &+\frac{1}{\rho}\sum_{m}\frac{\overline{W}}{W_{m}}\boldsymbol{\nabla}\cdot\mathbf{\Gamma}_{m}+\frac{1}{\rho}\sum_{m}\bigg{(}\frac{\overline{W}}{W_{m}}-\frac{h_{m}}{c_{p}T}\bigg{)}\dot{\omega}_{m}\equiv S\end{split}$$

$$\boldsymbol{\nabla}\cdot\boldsymbol{U}_{adv}=S+\underbrace{\frac{f}{p_{eos}}\bigg{(}\frac{p_{eos}-p_{0}}{\Delta t}+\boldsymbol{U}_{adv}\cdot\boldsymbol{\nabla p_{eos}}\bigg{)}}_{\delta\chi}\equiv\hat{S}$$

$$\boldsymbol{\tau}=\mu\bigg{[}\boldsymbol{\nabla}\textbf{{\em U}}_{adv}+(\boldsymbol{\nabla}\textbf{{\em U}}_{adv})^{T}-\frac{2}{3}\mathcal{I}(\boldsymbol{\nabla}\cdot\textbf{{\em U}}_{adv})\bigg{]}$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho {\text{emU}}_{adv})=0$$

$$\mathbf{\Gamma}_{m}=\rho Y_{m}\boldsymbol{{\cal V}}_{m}\,,\hskip 20.00003pt\boldsymbol{{\cal V}}_{m}=-\varUpsilon_{m}\,\textbf{{\em d}}\hskip 0.50003pt_{m}=-\varUpsilon_{m}\,\big{(}\textbf{{\em d}}\hskip 0.50003pt_{m,d}+\textbf{{\em d}}\hskip 0.50003pt_{m,\textit{ef}}\big{)}=\boldsymbol{{\cal V}}_{m,d}+\boldsymbol{{\cal V}}_{m,\textit{ef}}$$

$$\{\begin{array}{rl}\text{emd}{}_{m,d}& ={\nabla \mathbf{X}}_m\\ \text{emd}{}_{m,\text{ef}}& =\frac{{\nu }_m{q}_e{N}_{\text{scalebox0.75emA}}}{R_{u}T}\frac{\overline{W}}{W_m}{Y}_m\nabla \mathbf{\varphi }=z_m{Y}_m\frac{\overline{W}}{R_{u}T}\nabla \mathbf{\varphi }\end{array}$$

$$\nabla \cdot \mathbf{E}=\frac{q_t}{{ϵ}_0{ϵ}_r}$$

$$\mathbf{E}=-\nabla \mathbf{\varphi }$$

$$-{ϵ}_0{ϵ}_r{\nabla }^2\varphi =q_t$$

$${\kappa }_m={\mathcal{D}}_m\frac{{\nu }_m{q}_e{N}_{\text{scalebox0.75emA}}}{R_{u}T}$$

$${\nabla \mathbf{X}}_m=\frac{\overline{W}}{W_m}{\nabla \mathbf{Y}}_m+\frac{Y_m}{W_m}\nabla \overline{\mathbf{W}}$$

$$\rho Y_{m}\boldsymbol{{\cal V}}_{m,d}=-\rho\frac{W_{m}}{\overline{W}}{\cal D}_{m}\Big{(}\frac{\overline{W}}{W_{m}}\boldsymbol{\nabla Y}_{m}+\frac{Y_{m}}{W_{m}}\boldsymbol{\nabla{\overline{W}}}\Big{)}=-\rho{\cal D}_{m}\boldsymbol{\nabla Y}_{m}-\rho{\cal D}_{m}\frac{Y_{m}}{\overline{W}}\boldsymbol{\nabla{\overline{W}}}$$

$${\tilde{\mathbf{\Gamma }}}_m\equiv -\rho {\mathcal{D}}_m{\nabla \mathbf{Y}}_m-\rho {\mathcal{D}}_m\frac{Y_m}{\overline{W}}\nabla \overline{\mathbf{W}}-\rho Y_m{\kappa }_m\nabla \mathbf{\varphi }$$

$$c=\frac{T-T_{in}}{T_{max}-T_{in}}$$

$${\mathbf{\phi }}_t=\mathbf{F}(t,\mathbf{\phi }(t)),t\in [t^n,t^{n+1}];$$

$$\mathbf{\phi }(t^n)={\mathbf{\phi }}^n$$

$$\mathbf{\phi }(t)={\mathbf{\phi }}^n+{\int }_{t^n}^t\mathbf{F}(\tau ,\mathbf{\phi }(\tau ))\text{d}\tau$$

$${\mathbf{\phi }}^{(k+1)}(t)={\mathbf{\phi }}^n+{\int }_{t^n}^t\left[\mathbf{F}({\mathbf{\phi }}^{(k+1)})-\mathbf{F}({\mathbf{\phi }}^{(k)})\right]\text{d}\tau +{\int }_{t^n}^t\mathbf{F}({\mathbf{\phi }}^{(k)})\text{d}\tau$$

$$\mathbf{F}(t,\mathbf{\phi }(t))\equiv \mathbf{A}(t,\mathbf{\phi }(t))+\mathbf{D}(t,\mathbf{\phi }(t))+\mathbf{R}(t,\mathbf{\phi }(t))$$

$${\mathbf{U}}_{ef,m}^{(k)}={\mathbf{U}}_{adv}^{(k)}-{\nu }_m{\kappa }_m{\nabla \mathbf{\varphi }}^{n+1,(k)}$$

$${\widetilde{\rho Y}}_m^{n+1,k+1}=\rho Y^n+\Delta t\left[{\mathbf{A}}_m^{n+1/2,(k)}+\frac{1}{2}\left({\mathbf{D}}_m^n-{\mathbf{D}}_m^{n+1,(k)}\right)+{\mathbf{D}}_{m,AD}^{n+1,(k+1)}+{\mathbf{I}}_{R,m}^{(k)}\right]$$

$$\frac{\partial (n_e)}{\partial t}=-\nabla \cdot n_e(\text{emU}-{\kappa }_e\nabla \mathbf{\varphi })+\nabla \cdot D_e\nabla {\mathbf{n}}_{\mathbf{e}}+I_{R,e}$$

$${\epsilon }_0{\epsilon }_r{\nabla }^2\varphi =-\sum _m{z}_m{\widetilde{\rho Y}}_m+e{n}_e,$$

$$-n_e^{n+1}-\Delta t\nabla \cdot n_e^{n+1}(\text{emU}-{\kappa }_e{\nabla \mathbf{\varphi }}^{n+1})-\Delta t\nabla \cdot D_e{\nabla {\mathbf{n}}_{\mathbf{e}}}^{n+1}+f_e=0$$

$$-n_e^{n+1}+\frac{{\epsilon }_0{\epsilon }_r}{e}{\nabla }^2{\varphi }^{n+1}+f_{\varphi }=0,$$

$${\mathcal{J}}^{(l)}\delta {\mathbf{X}}^{(l)}=-F({\mathbf{X}}^{(l)})$$

$$(S_F^{-1}{\mathcal{J}}^{(l)}{S}_X)(S_X^{-1}\delta {\mathbf{X}}^{(l)})=-S_F^{-1}F({\mathbf{X}}^{(l)})$$

Full text

1 Introduction

Experiments have shown that applying electric fields to flames can provide an effective control of the combustion process by enhancing flame propagation speed, improving flame stabilization and reducing pollutant emissions [1, 2]. However, the development of such technology has proven difficult without a clear understanding of the interaction mechanisms between the flame and the electric field, and the use of electric fields is currently limited to flame detection sensors [3].

The chemical decomposition of hydrocarbons proceeds mainly through reactions involving neutral intermediate radicals. However, some reactions, called chemi-ionization reactions, also produce small quantities of charged chemical species and electrons [4, 5, 6]. These particles undergo a force when subjected to an electric field and their interactions with the surrounding gas can result in a global flame response to the electric field. Three major effects have been advanced in the literature [1]: 1) the collision of charged particles with neutral ones induces a bulk convective transport in the gas called the ionic wind effect; 2) the transport of highly reactive charged particles from the reactive layer of the flame to the low temperature zone enhances the fuel oxidation rate; and 3) for strong electric fields ohmic heating increases the flame temperature, resulting in a higher flame speed. These processes were found to have an effect on flame speed [7, 8, 9, 10, 11], flame stabilization [12] and NOx and soot formation [13, 14, 15]. The extent to which each process is important depends on the applied potential difference, the polarity, the distance between the electrodes and the flame, and the operating conditions, making it difficult to compare results from different experiments and to provide clear design guidelines for engineers.

Over the last decade, several groups have developed numerical methods to analyze the interactions of an electric field with charged particles in a flame. In most applications, the flame can be considered weakly ionized, i.e., the number density of electrons is much smaller than that of neutrals. However, the presence of charged particles, especially light electrons, results in challenging numerical issues associated with the wide scale separation between the electron dielectric relaxation scale and the comparatively slow hydrodynamic scale. Consequently, early studies focused mainly on steady-state one-dimensional flames [16, 17, 18] without an external applied electric field and identified the main chemical pathways associated with ions as well as the role of the ambipolar diffusion in the charged species spatial distribution. More recently, these steady-state numerical studies have been used to provide a more complete characterization of the flame response to external forcing (also called the $i-V$ curve, relating the current drawn from the flame to the applied voltage difference) [19, 20, 21, 22]. In agreement with experimental evidence, the effect of the external electric field is found to strongly depend on its polarity. The current is found to increase linearly with the potential difference before reaching a saturation current for high (positive) voltage. These studies highlight the dependence of the numerical results on the choices of the chemical mechanism and, to a lesser extent, on the modeling of electron and ion transport properties [23, 24]. Steady-state multi-dimensional simulations have also been reported [25, 26], showing that the simulations are able to capture qualitatively the change in flame shape and position resulting from the ionic wind. Due to the aforementioned multi-scale nature of the problem, fewer unsteady simulations are reported in the literature [27, 28, 29, 21, 30]. These simulations capture the effect of the electric field on the flame base position and investigate both direct current (DC) and alternative current (AC) conditions. To partially alleviate the fast electron drift velocity constraint on the stability of the numerical method, Belhi et al. [28, 29] employed a small value of the electron mobility $\kappa_{e}$ and a linearized approximation of the charged species transport equation. The effects of these assumptions on the flame response was not evaluated and this approach cannot be extended to more realistic values of $\kappa_{e}$ and higher intensity external electric fields without significant reduction of the simulation time step. In the plasma community, semi-implicit methods have been developed to overcome the electron time scale constraint [31]. However these approaches allow at best a couple orders of magnitude increase of the time step ($\sim 10^{-11}-10^{-13}$ s depending on the intensity of the electric field), which remains several order of magnitudes smaller than the hydrodynamic time scale in typical turbulent combustion applications ($\sim 10^{-7}-10^{-8}$ s).

In this paper, we propose a strategy based on multi-implicit spectral deferred correction (MISDC) method [32] to include the coupling between charged species and an electric field in a low Mach number combustion framework. The MISDC approach allows tight coupling between the different physical processes in a multi-scale simulation by including the effect of each process in their separate integration (in contrast to Strang splitting methods that consider each process sequentially and independently [32]). To alleviate the electron dielectric relaxation time scale constraint, the non-linear system formed by the coupled electron conservation equation and electrostatic potential equation is solved implicitly using a preconditioned Jacobian-free Newton Krylov (JFNK) method.

The paper is organized as follows. In Section 2 we introduce the low Mach number conservation equations including the electrostatic potential equation as well as the chemical and transport models. In Section 3 we discuss the changes implemented in the MISDC algorithm and details of the solution of the implicit non-linear system. We then provide a skeletal description of the time advance procedure. In Section 4 we present results for premixed flames in 1D under DC and AC conditions. Finally, the paper finishes with the main take-away of our approach and discusses future work.

2 Low Mach number equation set

2.1 Low Mach number equation set

This paper builds on the low Mach number equations set reported in previous work [33, 32], with the addition of an electrical drift contribution in the momentum, species and enthalpy equations, a separate conservation equation for the electron number density and a Poisson equation for the electrostatic potential to obtain an electric field consistent with the charged species distribution.

In the low Mach number regime, the characteristic velocity of the fluid $\textbf{{\em U}}_{adv}$ is much smaller than the speed of sound $\mathbf{a}$ (typically $|\textbf{{\em U}}_{adv}|/a=M\sim 0.1$ or even smaller), so the effect of acoustic wave propagation can be neglected since it does not affect the dynamics of the system. In numerical simulations, this effect is mathematically removed from the equations of motion and the system evolves subject to a time step based on the advective CFL condition. In low Mach number conditions, the total pressure can be decomposed into a spatially uniform (thermodynamic) component $p_{0}$, and a perturbational term, $\pi$, that drives the flow:

[TABLE]

Although the formulation supports a time varying $p_{0}$ (arising, for example, in closed chamber applications [43]), we assume an open domain here to simplify the exposition.

The set of equations describing species, electrons, enthalpy and momentum conservation in the low Mach number limit [33] are given by:

[TABLE]

where $N$ is the total number of species (excluding the electrons), $\rho$ is the density, $Y_{m}$ is the mass fraction of species $m$, $n_{e}$ is the electron number density, $\textbf{{\em U}}_{adv}$ is the fluid advective velocity, $\mathbf{\Gamma}_{m}$ (resp. $\mathbf{\Gamma}_{e}$) is the diffusion mass flux of species $m$ (electrons), $h=\sum_{m}(Y_{m}h_{m})$ is the mixture total (sensible and chemical) enthalpy with $h_{m}(T)$ the enthalpy of species $m$, $\dot{\omega}_{m}$ (resp. $\dot{\omega}_{e}$) is the production rate of species $m$ (electrons) due to chemical reactions, $\lambda$ is the thermal conductivity, $z_{m}$ is the electric charge per unit mass of species $m$, $\boldsymbol{E}$ is the electric field, and $\pi$ is the perturbational pressure arising from the low Mach number approximation. The evolution equations are closed by an equation of state, $p_{0}=p(\rho,T,Y)$ for the thermodynamic pressure. Note that the low Mach assumption requires that the flow evolve subject to a constant $p_{0}$. This DAE system can be solved by differentiating the equation of state in the frame of the fluid and requiring that the evolution be constrained to satisfy constant pressure in this frame [34]. Here we assume a mix of ideal gases:

[TABLE]

where $p_{0}$ is the ambient pressure, $\overline{W}$ is the mean molecular weight of the mixture, $W_{m}$ is the molecular weight of species $m$ and $R_{u}$ is the universal gas constant. Expanding in partial derivatives and using the conservation equations, the constant $p_{0}$ condition can be recast as a constraint on the velocity [33]:

[TABLE]

where $c_{p}$ is the specific heat at constant pressure for the mixture. Since this constraint is a linearization of the equation of state, the thermodynamic variables will not remain consistent with $p_{0}$ numerically; in order to prevent this thermodynamic drift, a correction term $\delta\chi$ has been added to to the constraint equation 7:

[TABLE]

where $0<f<1$ is a damping factor (see Day et al. [33] for details on the iterative implementation of this equation). Compared to classical low Mach number reactive flows, two additional source terms appear in the conservation equations: 1) the Lorentz volumetric forces (last term in Eq. (5) and 2) the ohmic heating, corresponding to the work of the Lorentz forces (last term in Eq. (4)).

The stress tensor in the momentum Eq. 5 is defined as:

[TABLE]

where $\mu(Y_{m},T)$ is the dynamic viscosity and $\mathcal{I}$ is the identity tensor (we ignore the bulk viscosity here). Since neither species diffusion nor chemistry redistribute total mass, we have $\sum_{m}{\mathbf{\Gamma}_{m}}=0$ and $\sum_{m}{\dot{\omega}_{m}}=0$. Noting that $\sum_{m}Y_{m}=1$ (ignoring the mass of electrons), the continuity equation can be derived summing up the species continuity equations:

[TABLE]

The diffusion flux of species $m$ can be expressed as:

[TABLE]

where $\boldsymbol{{\cal V}}_{m}$ defines the diffusion velocity of species $m$ in terms of EGLIB’s “flux diffusion vector” $\varUpsilon_{m}=\frac{W_{m}}{\overline{W}}{\cal D}_{m}$ [35] and the driving forces $\textbf{{\em d}}\hskip 0.50003pt_{m}$; ${\cal D}_{m}$ is the mixture-averaged diffusion coefficient of species $m$. Ignoring Dufour, Soret and barodiffusion terms, the diffusion and electric driving forces are, respectively:

[TABLE]

where $X_{m}$ is the mole fraction of species $m$, $\phi$ is the electric potential, $\nu_{m}$ is the valence ($e$ charges per molecule) of species $m$, $q_{e}$ is the elementary electron charge, $N_{\textbf{\scalebox{0.75}{{\em A}}}}$ is Avogadro’s number, and $z_{m}=\nu_{m}q_{e}N_{\textbf{\scalebox{0.75}{{\em A}}}}/W_{m}$ is the charge per unit mass of species $m$.

Under the electrostatic assumption, the local electric field $\boldsymbol{E}$ is obtained from Gauss’ law:

[TABLE]

where $q_{t}=\sum_{m}z_{m}\rho Y_{m}+q_{e}n_{e}$ is the local total charge number density of the mixture and $\epsilon_{0}$ and $\epsilon_{r}=1$ are the vacuum permittivity and the relative permittivity of the gaseous medium, respectively. The electric field is the negative gradient of the electrostatic potential $\phi$, i.e:

[TABLE]

Inserting Eq. (13) in Eq. (14) we obtain the electrostatic potential equation:

[TABLE]

The drift velocity $\boldsymbol{{\cal V}}_{m,ef}$ can also be written as $\boldsymbol{{\cal V}}_{m,ef}=\kappa_{m}\boldsymbol{E}$, where $\kappa_{m}$ is the mobility of species $m$ in the mixture. Thus, consistent with the Einstein relation [36], the mobility is defined as:

[TABLE]

The right-hand side of the diffusive driving force in Eq. 12 can be rewritten as:

[TABLE]

and so the diffusion fluxes can be rewritten in terms of mass fractions gradients plus $\overline{W}$ corrections:

[TABLE]

We will use this form of the transport equation to build an iterative time-implicit update scheme based on lagging the corrections and sweeping through the species with decoupled linear solves for the Crank-Nicolson update. The resulting form of the diffusive species flux will be:

[TABLE]

However, since we employed mixture averaged diffusion coefficients, equation 19 will not in general satisfy that $\sum_{m}\tilde{\boldsymbol{\Gamma}}_{m}=0$; to conserve mass, we introduce a correction velocity [33] that guarantees that these fluxes sum to zero. Since we use an implicit approach to compute the diffusion fluxes, we first solve the implicit system to evaluate the original fluxes $\tilde{\boldsymbol{\Gamma}}_{m}$, then we conservatively correct $\tilde{\boldsymbol{\Gamma}}_{m}$ so that they sum to zero on each cell face (we will denote the corrected fluxes as $\mathbf{\Gamma}_{m}$), and finally we modify the time-advanced values of the mass fractions $Y_{m}$ to be consistent with the corrected fluxes.

2.2 Chemical mechanism and species transport properties

The chemical mechanism employed in this work combines the GRI3.0 [37] for the oxidation of methane with the reaction mechanism for charged species reported in Belhi et al. [38]. The combined mechanism contains 61 species (not including electrons) and 386 reactions, and includes 10 ions (4 cations and 6 anions) as listed in Table 1. In the remainder of the paper, charged species refers to the ions whereas charged particles also includes the electron. Several studies have showed that anions are only present in very small quantities in freely evolving flames; electrons account for most of the negative charges. However, Belhi et al. [30] recently showed that including the anions (especially large anions such as CO${}_{3}^{-}$ and HCO${}_{3}^{-}$) is essential to reproduce the ionic wind motion observed experimentally.

The thermodynamic data for the charged species from the Burcat [39] database were used. The computation of the transport properties for the charged particles listed in Table 1 uses the EGLIB library. Specific treatment of the ion/neutral or ion/ion collision is not investigated in this work, the use of (n,6,4) and Coulomb [40] interaction potentials for ions/neutrals and ions/ions collisions as described in Han et al. [24] will be studied in future work.

The electron transport coefficients require a more detailed treatment. For low values of the reduced electric field $|\boldsymbol{E}|/\mathcal{N}$, where $\mathcal{N}$ is the background gas number density, the electrons are in thermal equilibrium with the mixture. In these conditions, the electron temperature is equal to that of the mixture: electrons are accelerated by the electric field $\boldsymbol{E}$, but the collision frequency with neutral species (represented by $\mathcal{N}$) is high enough to prevent the electrons from reaching high kinetic energy conditions. For higher values of $|\boldsymbol{E}|/\mathcal{N}$, the electrons gain sufficient kinetic energy that their energy (temperature) is higher than the remainder of the mixture. For this case the electrons are said to be non-thermal. Under these conditions, the evaluation of the electron transport coefficients require the computation of the evolution of the electron energy distribution function (EEDF) by solving the Boltzmann equation [41]. Additionally, the chemi-ionization reaction $\text{CH}+\text{O}\rightarrow\text{HCO}^{+}+e^{-}$ is no longer the only chemical pathway producing electrons since impact ionization rates become important [41]. However, this last effect is not included in our framework at the present time and its relevance will be the subject of future studies.

Previous studies employed a constant value of the electron mobility $\kappa_{e}=0.4$ mJ*-1s-1* [21, 26] since this constant value was found to provide a good agreement with simulations obtained from more detailed thermal electron transport calculations [23]. The framework developed in this work aims at simulating realistic engineering applications characterized by relatively high external voltages, conditions at which electrons can no longer be assumed thermal. In order to include a non-thermal electron transport coefficient without explicitly computing the evolution of the EEDF, the mixture composition and temperature are extracted from the simulation as function of the progress variable $c$:

[TABLE]

Note that other definitions of the progress variable could be used for fuels exhibiting more complex behaviors. This information is then used in the BOLSIG+ [42] code to estimate the EEDF and the corresponding value of the electron mobility and diffusion coefficient at different values of $|\boldsymbol{E}|/\mathcal{N}$. The resulting two-dimensional tables are shown in Fig. 1 and electron transport properties are extracted from these tables during the simulation using $c$ and $|\boldsymbol{E}|/\mathcal{N}$. Note that

3 MISDC strategy

3.1 MISDC strategy

The present strategy builds upon the MISDC methodology developed in Nonaka et al. [32, 43]. As a brief reminder, the spectral deferred correction (SDC) method [44] solves a system of ordinary differential equation:

[TABLE]

using the integral form:

[TABLE]

The SDC method generates successive approximation $\boldsymbol{\varphi}^{(k)}(t)$ of $\boldsymbol{\varphi}(t)$ using the update equation:

[TABLE]

where the explicit dependence of $\boldsymbol{F}$ and $\boldsymbol{\varphi}$ on $t$ in the integrals has been dropped for simplicity. By using a low-order approximation of the first integral and a more accurate quadrature rule for the second integral, the SDC method effectively constructs an arbitrary order (the order of the quadrature rule) solution by successive low-order corrections of the approximation $\boldsymbol{\varphi}^{(k)}(t)$. In MISDC [45, 46], $\boldsymbol{F}$ is decomposed into distinct processes, that can be treated separately in their own time scales:

[TABLE]

with $\boldsymbol{A}$, $\boldsymbol{D}$ and $\boldsymbol{R}$ referring here to the advection, diffusion and reaction processes, respectively. Here, following [32], $\boldsymbol{A}(t,\boldsymbol{\varphi}(t))$ and $\boldsymbol{D}(t,\boldsymbol{\varphi}(t))$ are piece-wise constant over each time step. The former is evaluated using a second-order Godunov method while the latter is evaluated using a midpoint rule. To accommodate the stiffness of hydrocarbon chemical reactions, the update equation for the reaction is formulated as an ODE and integrated using a stiff ODE package such as CVODE. The effects of advection and diffusion are taken into account as temporally constant forcing terms in the chemical ODE integration (see [32] for more details on the integration procedure).

The main steps of the integration algorithm are summarized in Algorithm 1. The set of transported thermodynamic scalars is written as $\boldsymbol{\varphi}=(\rho,\rho h,T,n_{e},\rho Y_{m})^{T}$. The superscript $n$ indicates the timestep and $(k)$ is the SDC iteration index. The diffusion operator for scalar $\varphi$ at time $t^{n}$ is written $\boldsymbol{D}_{\varphi}^{n}$, the $k$-th approximation of this operator at time $t^{n+1}$ is written $\boldsymbol{D}_{\varphi}^{n+1,(k)}$ and the $k$-th approximate of the advection operator obtained with the Godunov procedure is $\boldsymbol{A}^{n+1/2,(k)}_{\varphi}$. The charged species drift flux appearing in Eq. 19 is non-symmetric and introduces numerical instabilities when discretized with the species diffusion flux using a second order centered scheme. To overcome this difficulty, the drift flux is treated along with the advective flux in the second-order Godunov procedure by constructing an effective velocity for each species, $m$:

[TABLE]

The resolution of the coupled electron/electrostatic potential non-linear system requires provisional charged species mass fraction $\widetilde{\rho Y}^{n+1,(k+1)}_{m}$:

[TABLE]

where $\boldsymbol{I}^{(k)}_{R,m}$ is the integrated representation of the reaction term for species $m$ from the previous SDC iteration.

3.2 Non-linear implicit solution

At each SDC iteration, we solve the non-linear system formed by the electron conservation equation (Eq. (3)) and the electrostatic potential equation (Eq. (15)):

[TABLE]

where $\widetilde{\rho Y}^{n+1,k+1}_{m}$ is the provisional charged species mass fraction at the current SDC iteration and $I_{R,e}$ is the last evaluation of the electron chemical source term.

Using a first order backward Euler time discretization, the implicit non-linear system can be written as:

[TABLE]

where $f_{e}=-I_{R,e}^{n}+n_{e}^{n}$ and $f_{\phi}=\sum_{m}z_{m}\widetilde{\rho Y}^{n+1,k+1}_{m}/e$. Introducing $\boldsymbol{X}=(\boldsymbol{n_{e}},\boldsymbol{\phi})$, Eq. (29) can be written as $F(\boldsymbol{X})=0$, where $F(\boldsymbol{X})$ is the non-linear residual. This system is solved using a Jacobian-free Newton-Krylov (JFNK) method [47].

The basis of JFNK is the iterative non-linear Newton solution, where at each iteration $l$, a linear system of the form:

[TABLE]

is solved. Here, $\delta\boldsymbol{X}^{(l)}=\boldsymbol{X}^{(l+1)}-\boldsymbol{X}^{(l)}$ is the Newton update and $\mathcal{J}^{(l)}=\mathcal{J}(\boldsymbol{X}^{(l)})=\partial F(\boldsymbol{X}^{(l)})/\partial\boldsymbol{X}$ is the system Jacobian matrix. In practice, the components of $\boldsymbol{X}^{(l)}$ and $F(\boldsymbol{X}^{(l)})$ can have entries than span a large range of values which can affect the solution of the linear system (30) and destroy the convergence properties of Newton’s method. To address this issue, Eq. (30) is scaled by two diagonal matrices $S_{F}$ and $S_{X}$:

[TABLE]

where $S_{F}$ contains typical values of $F(n_{e})$ and $F(\phi)$ respectively, and $S_{X}$ contains typical values of $n_{e}$ and $\phi$. The typical values are evaluated at the beginning of the non-linear iterations since the apropriate values may evolve with the solution. The Newton iterations are stopped when the norm of scaled residual is reduced by $\epsilon_{F}$ orders of magnitude or the scaled magnitude of the Newton step drops below a certain value $\epsilon_{X}$:

[TABLE]

These tolerances must be chosen to ensure that the non-linear solution residual remains smaller than the truncation error of the numerical schemes. A backtracking linesearch algorithm is employed for globalization of the Newton method [48].

For large non-linear systems encountered in multi-dimensional simulations, the computational cost and memory requirements of solving a linear system with a direct solver at each Newton iteration are prohibitive. Our implementation thus uses the GMRES Krylov method [49] to solve the scaled linear system (31). For clarity, the left and right scaling matrices will not be carried in the following description and the outer (Newton) iteration index $l$ is dropped (the scaling is implemented in the code; however). The GMRES starts with an initial guess $\delta\boldsymbol{X}_{0}$, and the corresponding residual $\boldsymbol{r_{0}}=-F(\boldsymbol{X}_{0})-\mathcal{J}\delta\boldsymbol{X}_{0}$. In the context of a Newton-Krylov method, $\delta\boldsymbol{X}_{0}=0$ is used since the Newton step tends toward zero as we go through the Newton iterations. At the $p^{th}$ iteration of the GMRES method, we construct an approximation $\delta\boldsymbol{X}_{p}$ of the solution by solving a minimization problem in the Krylov subspace $\mathcal{K}_{p}$ of $\mathcal{J}$:

[TABLE]

It can be seen that the GMRES method only needs the action of the Jacobian matrix on a vector. For large linear systems, the construction and storage of matrix $\mathcal{J}$ can hinder the performance and the scalability of the algorithm. In the JFNK context, the explicit construction of $\mathcal{J}$ is dropped in favor of a finite difference approximation of the the matrix/vector product $\mathcal{J}\boldsymbol{v}$:

[TABLE]

where $\varepsilon_{FD}$ is a small number. The quality of the approximation of $\mathcal{J}.\boldsymbol{v}$ depends on the choice of $\varepsilon_{FD}$. Here we use the method employed in the Trilinos package [50]:

[TABLE]

where $\lambda_{FD}=\varepsilon_{mach}^{1/3}$ is a small parameter related to the machine precision $\varepsilon_{mach}$. The linear solver is iterated until:

[TABLE]

A constant value of $\gamma$ is kept throughout the simulation and the effect of the choice of $\gamma$ on the non-linear system solution will be assessed in Section 4.2.

The performance of the JFNK depends strongly on the number of GMRES iterations required to solve (31). If $\mathcal{J}$ has a large condition number, the Krylov method requires a large number of iterations to converge. In this case, it is necessary to apply a preconditioner to the linear system:

[TABLE]

where $P$ is an approximation of $\mathcal{J}$, such that $P^{-1}\mathcal{J}\sim\mathcal{I}$. The main objective of the preconditioner is to cluster the eigenvalues of the resulting $P^{-1}\mathcal{J}$ matrix, allowing the GMRES method to find a good $\delta\boldsymbol{X}_{p}$ in a small Krylov space (i.e. small number of iterations). To construct the preconditioner, we start by linearizing Eq. (29):

[TABLE]

This allows us to write the block matrix form of the Jacobian $\mathcal{J}$ resulting from the spatio-temporal discretization of Eq. (29):

[TABLE]

where the block matrices $\boldsymbol{\mathcal{D}_{f}}$, $\boldsymbol{\mathcal{D}_{r}}$ and $\boldsymbol{\mathcal{L}}$ are the spatial operators underlined in Eq. 39. Note that $\boldsymbol{\mathcal{I}_{e}}$ actually differs from the identity matrix because of the scaling applied to the linear system (31). Schur factorization of the inverse of the 2 $\times$ 2 block Jacobian is written as:

[TABLE]

where $\boldsymbol{\mathcal{S}}=\boldsymbol{\mathcal{L}}-\boldsymbol{\mathcal{I}_{e}}(\Delta t\boldsymbol{\mathcal{D}_{f}}-\boldsymbol{\mathcal{I}})^{-1}\Delta t\boldsymbol{\mathcal{D}_{r}}$ is the Schur complement of $\mathcal{J}$. Here, $P^{-1}$ is the exact inverse of the Jacobian matrix and it only requires $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\boldsymbol{\mathcal{I}})^{-1}$ and $\boldsymbol{\mathcal{S}}^{-1}$, both of which are easier to invert than $\mathcal{J}$. However, computing $\boldsymbol{\mathcal{S}}^{-1}$ is still difficult since the construction of $\boldsymbol{\mathcal{S}}$ requires the solution of $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\boldsymbol{\mathcal{I}})^{-1}$. Instead, we use an approximation $\widetilde{\boldsymbol{\mathcal{S}}}=\boldsymbol{\mathcal{L}}+\boldsymbol{\mathcal{I}_{e}}\Delta t\boldsymbol{\mathcal{D}_{r}}$ of $\boldsymbol{\mathcal{S}}$ that is easier to solve. It can be seen that for small time steps, $\widetilde{\boldsymbol{\mathcal{S}}}$ is a good approximation of $\boldsymbol{\mathcal{S}}$. Both $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\boldsymbol{\mathcal{I}})$ and $\widetilde{\boldsymbol{\mathcal{S}}}$ are then diagonally dominant and can be solved effectively using a multi-grid (MG) approach. The present implementation uses a standard V-cycle approach with red-black Gauss-Siedel relaxation to solve both linear systems to a tolerance $\gamma_{MG}$. The effect of $\gamma_{MG}$ on the performance of the JFNK is evaluated in Section 4.2. Applying $P^{-1}$ to any vector $\boldsymbol{v}$ requires the application of the successive matrices of Eq. (41) to $\boldsymbol{v}$. In its classical Schur factorization form (41), this entails four MG solves (three solves of $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\boldsymbol{\mathcal{I}})^{-1}$ and one of $\widetilde{\boldsymbol{\mathcal{S}}}^{-1}$). To save one MG solve, the block factorization (41) is rewritten in the following form:

[TABLE]

The solution of the implicit non-linear system is summarized in Algorithm 2. Superscript $l$ corresponds to the Newton iteration index while subscript $p$ is the GMRES iteration index.

4 Numerical experiments

We first evaluate the robustness and performance of the proposed algorithm in order to optimize the numerical parameters and tolerances employed in the JFNK. Then, simulations of steady one-dimensional premixed methane/air flames subject to DC electric fields are performed in order to estimate the accuracy of the complete algorithm and provide comparisons with experimental data. Finally, the behavior of flames subjected to AC electric fields at various frequencies is analyzed.

4.1 Numerical set-up

Throughout this section, we consider an unstrained one-dimensional burner-stabilized premixed methane/air flame. The operating conditions correspond to the experimental study of Speelman et al. [20]: the inlet velocity is set to the flame speed of a stoichiometric methane/air flame at T = 300 K while the inlet temperature is set to T = 350 K, such that the flame is stabilized on the left boundary of the domain. Simulations are initialized from a resolved CANTERA [51] solution ($\sim$ 4000 unequally-spaced grid points), that does not include the effect of the electric field. The CANTERA solution is interpolated onto a set of uniform grids with varying resolution, and simulations are evolved initially without external electric forcing for 5 ms in order to eliminate any spurious artifacts introduced by the initialization. Subsequently, the external electric field is activated and set to the desired values. The main characteristics of the simulations are listed in Table 2. Unless otherwise specified, the numerical parameters ( $k_{max}$, $\gamma$, …) listed in Table 2 are employed.

The interactions of the electric field with the charged particles in the flame introduces additional time scales compared to classical reactive flow simulations. The following is an overview of the relevant characteristic time scales and summarizes the specific treatment used here in the numerical strategy:

• •

bulk advective time scale:

[TABLE]

• •

species/electron diffusive time scale:

[TABLE]

where $d$ is the number of dimensions.

• •

chemical reaction time scale:

[TABLE]

where $\dot{\omega}_{m}$ is chemical production rate per volume of species $m$.

• •

charged species/electron effective convective time scale:

[TABLE]

[TABLE]

where the drift velocity of the species is considered. Note that for large values of the external electric field, the drift velocity can oppose the convective velocity and its magnitude can exceed it. Additionally, the large mobility of the electrons results in a more stringent time step constraint, compared to ions.

• •

electron dielectric relaxation time scale characterizes the response of the electric field to a change in the electron distribution:

[TABLE]

The first three time scales are common in reactive flow simulations. In most combustion simulations using detailed chemical kinetics, the chemical time-step constraint is alleviated in the numerical implementation by using a stiff ODE integrator. Additionally, we use a semi-implicit Crank-Nicholson method for conduction and species diffusion which enables time steps larger than the fast diffusive time scales of light species. The advection of charged species is treated time explicitly so that the advective time scale constrains the overall simulation time-step. Although this often results in time steps smaller than $\tau_{bulk}$, the charged species time scales are still several orders of magnitude larger than that of the electrons. Typical values of the various time scales are plotted against the external electric forcing $\Delta V$ in Fig. 2. The data is based on a $n_{x}=512$ grid points simulation, corresponding to $\Delta_{x}=19.5\mu$m. The bulk advective time scale is only shown as a reference for the classical low Mach number time constraint. Both $\tau_{conv,e}$ and $\tau_{conv,m}$ decrease with increasing values of $\Delta V$; $\tau_{conv,e}$ is approximately four orders of magnitude smaller than $\tau_{bulk}$. Both advective time scales exhibit a plateau at around $\Delta V=250$ V.cm*-1*, corresponding to the saturation voltage. At the same location, the dielectric time scale jumps to exceed $\tau_{bulk}$. This is behavior is related to the drop in peak electron number density as the external voltage exceeds the saturation value. Across the range of $\Delta V$ considered, the time scales associated with electrons are several orders of magnitude more stringent than the others, thus highlighting the need for an implicit treatment of the electrons.

4.2 Iterative solvers performance

Solution of the implicit non-linear electron/electrostatic potential system with JFNK involves several tolerances, which can have a significant impact on both the robustness and performance of the proposed methodology. A series of tests are performed in order to evaluate the optimal settings.

At the lower level of the JFNK algorithm is the application of the inverse of the preconditioner $\widetilde{P}$ on the GMRES basis vectors (see Eq. (42)). This requires three MG solutions (two solutions of $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\mathcal{I})^{-1}$ and one solution of $\widetilde{\boldsymbol{\mathcal{S}}}^{-1}$) using a standard V-cycle approach [52]. Two relaxation operations are applied going down and up each level of the V-cycle based on a Red-Black Gauss Siedel. Figure 3(a) shows the total number of GMRES iterations per SDC iteration, as function of the V-cycle tolerance $\gamma_{MG}$. Figure 3(b) shows the total number of V-cycle as function of $\gamma_{MG}$. To separate the effect of each block on the performance of the preconditioner, the MG tolerance is tested for one block while the other is solved exactly (using a tri-diagonal solver in the present one-dimensional case). The number of GMRES iterations is only marginally affected by the multigrid tolerance on $\widetilde{\boldsymbol{\mathcal{S}}}^{-1}$, while loose tolerance on $(\Delta t\boldsymbol{\mathcal{D}_{f}}-\mathcal{I})^{-1}$ results in a large increase of the number of iterations. However, the number of V-cycles directly increases the CPU cost of the algorithm and, in the present cases, the trade-off between the MG tolerance and the total number of V-cycles shows V-cycles minimized around $\gamma_{MG}\sim 1.0e^{-4}$ as can be observed in Figure 3(b).

Given these settings, the efficiency of the preconditioner can be directly evaluated by comparing the convergence of the GMRES solver with and without preconditioning. Figure 4 shows the GMRES residual as function of the GMRES iteration count for different values of $\Delta V$ both with and without the preconditioner. The preconditioned systems converge 20 to 50 times faster, regardless of the operating conditions. Note that simulations are performed at a constant CFLconv,m, so that the time step is reduced as $\Delta V$ increases, resulting in a more efficient preconditioning ($\widetilde{\boldsymbol{\mathcal{S}}}^{-1}$ tends towards $\boldsymbol{\mathcal{S}}^{-1}$ as $\Delta_{t}$ decreases). Additionally, the preconditioned system convergence is only marginally affected by the size of the system whereas the unpreconditioned system convergence rate decreases with $n_{x}$ (not shown).

4.3 Method convergence

In order to evaluate the convergence properties of the complete algorithm, simulations are performed halving the inlet velocity and evolving the system to a fixed time with increasing resolution, decreasing $\Delta x$ by a factor two with each refinement. The simulations are performed at a constant CFLm,conv. The error is obtained by comparing the results at resolution $\Delta x$ with those computed with twice the resolution $\Delta x/2$. The $L^{2}$ norm of the error for a simulation with $n_{x}$ cells is:

[TABLE]

where $\varphi_{i}^{c-f}$ is the average of the fine results onto the coarser grid. Figure 5 shows the $L^{2}$ norm of the error at four grid resolutions for 6 scalars: $\rho$, $\rho h$, $Y_{CH_{4}}$, $Y_{H_{2}}$, $n_{e}$ and $Y_{H_{3}O^{+}}$. The slope of the error shows that second order is reached for all variables across the range of external forcing considered. The error on neutral species and mixture averaged quantities is not affected by the external forcing whereas the error on $n_{e}$ and $Y_{H_{3}O^{+}}$ decreases for high $\Delta V$ values (but remains a second order convergence rate). This decrease in errors indicates that the applied voltage is higher than the saturation voltage at which the charged species are drawn away from the reaction zone by the electric drift as fast as they are produced by chemical reactions.

4.4 Steady premixed flame under DC

Burner-stabilized, steady-state premixed methane/air flames subjected to DC electric fields have been studied using the PREMIX program in previous studies [20, 38]. Fig. 6 shows the temperature as well as oxygen, CH, electrons, H3O+ and C2H3O+ profiles across the flame in the absence of an external electric field. Oxygen and CH are the key neutral species controlling the production rate of electrons, i.e. the number of charged particles in the flame and consequently the maximum current that can be drawn from the flame [20, 22]. Accordingly, the peak electron density in Fig. 6 is located near the corresponding maximum of CH. Note that the number density of charged species is about 5 orders of magnitude smaller than that of an intermediate radical such as O. In the absence of an external electric field, the sum of number densities of the two major cations (H3O+ and C2H3O+) equals that of electrons, as ambipolar diffusion tends to balance charge separation resulting in a near electro-neutral gas.

The peak value of electron and H3O+ is higher than that reported in a previous study [38] where the neutral chemical mechanism was optimized to better reproduce the CH distribution. This study showed that the GRI3.0 mechanism over-predicts the CH mass fraction, resulting in higher chemi-ionization rate and electron maximum number density.

Figure 7 shows comparisons between experimental $i$-V curves [20] and the present simulations. The current $i$ is evaluated by computing the charge flux carried by the charged species $m$:

[TABLE]

and summing over positive and negative species:

[TABLE]

Finally, $i=\boldsymbol{J}S_{expe}$, where $S_{expe}=7.04$ cm*-2* is the experimental cross section of the burner [20]. Note that from the species conservation equations, in steady-state conditions $\boldsymbol{J}_{m}=\frac{z_{m}}{q_{e}}\dot{\omega}_{m}$, showing that the current drawn from the flame is directly related to the production rate of charged particles. The simulation results are consistent with the experimental data: the current increases for positive voltage until it reaches a plateau as the applied voltage exceeds a saturation value. In contrast, higher negative voltage is required to reach saturation conditions. This effect of the polarity, known as diodic effect, results from the large difference in distance between the flame and each electrode [20]. The over-prediction of the saturation current is consistent with the fact that the mass fraction of CH is over-predicted by the GRI3.0 mechanism.

The profiles of charge particles at different values of the external voltage $\Delta V$ is presented in more detail in Figures 8, 9 and 10, which also show the steady-state profiles of electrostatic potential and electric field, for both positive and negative $\Delta V$. For sub-saturation voltages, the electric field profiles show the existence of a ’dead zone’, where the electric field is close to zero and the particles are not affected by electric forces. As the external voltage intensity increases, the electrode sheath develops, eventually penetrating into the reaction zone of the flame. As saturation conditions are reached, the peak number densities of charged particles drop since they are convected away from the reaction zone as fast as they are produced through chemi-ionization. This drop is responsible for the jump of $\Delta t_{diel}$ in Fig. 2 and the drop in $L^{2}$-norm of the error on charged particles in Fig. 5.

The charged particles profiles in Fig. 10 show that the proposed algorithm is able to accommodate very sharp profiles in the charged species distribution, without introducing numerical noise that would eventually lead to unstable numerical oscillations due to the strong non-linear coupling between the electrons and the electrostatic potential.

4.5 Unsteady premixed flame

In order to demonstrate the potential of the method to tackle unsteady simulations, the behavior of the burner-stabilized flame subjected to AC conditions is studied. Three forcing amplitudes $A_{AC}$ are considered: 100V, 1000V and 2500V respectively corresponding to conditions below both positive and negative saturation voltages, above positive and below negative saturation voltages and above both positive and negative saturation voltage. In order to estimate a characteristic relaxation time of the electrical structure of the flame, the time required for the electron and H3O+ to travel across the computational domain is evaluated using an averaged effective velocity $\overline{\boldsymbol{U}}_{ef}$ based on a representative mobility value for each particle and the external forcing intensity $\Delta V$. Table 3 summarizes these relaxation times for the three forcing amplitudes considered, showing that electrons relax within a few micro-seconds, whereas it takes about a 100 times longer for the H3O+ to relax.

These relaxation time scales can be compared to the half period of an AC forcing $\tau_{AC}/2$ to distinguish several regimes (for a fixed value of $\Delta V$): 1) for low forcing frequency, $f_{AC}$, both electrons and ions remain in quasi-equilibrium with the instantaneous potential difference and the charged particles profiles are close to the corresponding steady-states; 2) for higher $f_{AC}$, the electrons are close to equilibrium, but the slower ions do not reach steady-state, changing the current drawn from the flame and possibly inducing an asymmetric ionic wind due to the diodic effect; and 3) for very high $f_{AC}$ the ions are too slow to respond to the change in external electric potential and the ionic wind (mainly due to the motion of ions) becomes negligible. In practice, only the first two regimes are of interest to study the effect of ionic wind on the flame behavior. Additionally, in the first regime, the flame structures are expected to remain close to the ones described in Section 4.4 so that we will focus on the second regime by considering $f_{AC}$ listed in Table 4.

The temporal evolution of the H3O+ profile during a statistically steady period of the AC forcing is shown in Fig. 11 for the seven values of $f_{AC}$ considered at a constant forcing amplitude of 1000 V. Each plot shows the evolution of the one-dimensional H3O+ number density profile (horizontal direction) as function of the normalized time $t^{*}=t/\tau_{AC}$ (vertical direction, from top to bottom). A few periods are necessary before reaching statistically steady oscillations. Note that, these plots confirm that the proposed algorithm is able to smoothly capture the fast motion of steep charged species fronts.

Figure 11 shows that for the initially positive polarity, the development of the cathode sheath is qualitatively similar to the steady states depicted in Fig. 10(a). For low forcing frequency, positive saturation conditions are reached for most of the cycle first half-period. As the forcing frequency increases, the cathode sheath is no longer able to fully develop and H3O+ depletion near the right boundary of the computational domain remains during part the second half of the cycle, even though the polarity is reversed. Additionally, the peak value of H3O+ is found to decrease with increasing frequency, indicating that the charged particle profiles are no longer able to relax to the forcing free profiles while the instantaneous voltage is close to zero.

To analyze the effect of the forcing frequency and amplitude on the ionic wind effect, the integral of the Lorentz forces appearing the momentum equation (5) across the computational domain is computed:

[TABLE]

corresponding to the force per unit area. Additionally, to evaluate the effect of this force on the flame, the average work of the Lorentz forces $W_{Lorentz}$ over a period is evaluated from the Ohmic heating term appearing in the energy equation (4). Figures 12(a-c) show the temporal evolution of $F_{Lorentz}$ during one AC forcing period for the three values of the forcing amplitude while Fig. 12(d) shows the evolution of $W_{Lorentz}$ as function of the frequency for different values of $A_{AC}$. For small forcing amplitude (below both positive and negative saturation values) the proximity of the flame to the anode, also responsible for the diodic effect, results in an overall negative Lorentz force, the work of which increases with increasing frequency. As the forcing amplitude increases, the positive Lorentz force during the second half of forcing cycle becomes more important. In these conditions, increasing the forcing frequency results in an increase of the averaged work generated by the Lorentz force, up to a critical frequency above which the electric field is no longer able to penetrate into the flame and the work begin to decrease. These results indicate that, as in the DC cases, the effect of the AC electric field not only depends upon the forcing frequency and amplitude, but also the flame position compared to the electrodes.

5 Conclusion

This work proposes a new numerical strategy to include the motion of charged particles in simulations of low Mach number reactive flows in the presence of electric fields. We have found that to overcome the stringent time-step constraint imposed by fast electrons and their coupling with the electrostatic potential equation, a non-linear implicit solution of the system of equations governing these two quantities is necessary. Keeping in mind the need for an efficient methodology in large scale (multi-dimensional) computations, we have developed an algorithm that introduces a JFNK solver within the SDC iterations developed for classical reactive flow simulation. To obtain good performance, we constructed a preconditioner based on the Schur decomposition of the Jacobian matrix for the electrons/electrostatic potential system. An approximation of the Schur complement of the Jacobian matrix is proposed enabling use of multi-grid method to approximate the inverse of the preconditioner in the iterative linear solve.

We demonstrated on one-dimensional burner-stabilized premixed flame configurations, that second-order accuracy is reached for all the transported variables and for a wide range of external electric forcing. The numerical results compare well with experimental data regarding the current-voltage characterization of the flame (given the uncertainty on the chemical mechanism) and detailed analysis of the charged particles profiles are consistent with previous studies using steady-state one-dimensional solvers.

The proposed strategy is currently being implemented in the low Mach number reactive flow solver PeleLM, which is based on the block-structured adaptive mesh refinement library AMReX. The resulting unique numerical tool will allow us to investigate realistic engineering applications of electric field controlled flames that have so far not been possible.

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program (under Award Number DE-SC0008271 and under contract No. DE-AC02-05CH11231).