A three-dimensional model of streamer discharges in unsteady airflow

Streamers are thin ionized channels that initiate large-scale electric discharges such as sparks and lightning, and a deeper understanding of their nature will shed light on the fascinating dynamics of long sparks. Streamers have been extensively studied both by numerical and experimental means [1–5]. Streamer dynamics are highly nonlinear and small changes in the discharge conditions can lead to large variations in streamer dynamics. The fast propagation and small dimensions of streamers do not allow experiments to study the microscopic characteristics of streamers with high resolution. Thus, numerical investigation of streamer discharges have been used to understand streamer dynamics, and numerous models of varying complexity have been developed. We can for instance mention particle, hybrid particle/fluid and fluid models [3, 6–9].

In this paper, we use a fluid description for cold atmospheric plasma. We are interested in unsteady fluctuations in air density and also the influence of the air velocity. A few streamer models have tried to study the streamer propagation in non-uniform air [10, 11]. There are two main physical mechanisms from which unsteady airflow might affect the electric discharges. Firstly, the breakdown electric field is modified due to variations in air pressure, which consequently affects the streamer dynamics [10, 12]. Secondly, the airflow affects the motion of charged particles by transferring momentum. The role of this mechanism in electric discharges can be significant, yet only a few research papers have investigated it [13, 14]. Therefore, we have included Navier–Stokes and turbulence equations to model the airflow.

The model can study these effects in an accurate and realistic manner by coupling electric discharge set of equations with Navier Stokes and turbulence equations. It is fully three dimensional, allowing calculation of realistic flows around complex structures, and of streamers branching and interactions with each other [15]. The model also allows the incorporation of structures with complicated geometries and is built on an unstructured mesh to capture the curves and complexities of the structure. Unstructured meshes have previously been used in streamer simulations [16].

The aim of the model is to capture the effects of unsteady airflow on electric discharges. Part of these effects involve the motion of ions, and therefore occur on longer timescales than the streamer timescales. Our model integrates implicitly over time, removing the Courant–Friedrichs–Lewy (CFL) restriction on the time step size, which allows us to model long-duration electric discharges of hundreds of microseconds.

The model is built on the commercial software Ansys Fluent. Ansys Fluent is based on the finite-volume method and allows the implementation and modification of equations and numerical schemes associated to electrical discharges. Furthermore, it supports adaptive refinement on unstructured mesh which is considered necessary for our model.

In the following sections, we describe the physics of the model, its implementation, and the numerical schemes used. Later, we verify the model accuracy by comparing to a test case from [17]. In the last section, we present a simulation of long-duration discharge which shows the ability of the code to couple the airflow and electric discharge dynamics, and explains the influence of the wind on the discharge.

Charged particles dynamics are modeled by drift–diffusion equations [18–20].

$$\begin{align}\hfill \frac{\partial {n}_{\text{e}}}{\partial t}& =\nabla \cdot \left({D}_{\text{e}}\enspace \nabla {n}_{\text{e}}\right)+\nabla \cdot \left({\boldsymbol{v}}_{\text{air}}+{\mu }_{\text{e}}\boldsymbol{E}{n}_{\text{e}}\right)\hfill \\ \hfill & \quad +{n}_{\text{e}}\left({\nu }_{\text{i}}-{\nu }_{\text{a}}\right)-{S}_{{r}_{\text{e}-\text{p}}},\hfill \end{align} \tag{ 1.1 }$$

$$\begin{align}\hfill \frac{\partial {n}_{\text{p}}}{\partial t}& =\nabla \cdot \left({D}_{\text{p}}\enspace \nabla {n}_{\text{p}}\right)+\nabla \cdot \left({\boldsymbol{v}}_{\text{air}}+{\mu }_{\text{p}}\boldsymbol{E}{n}_{\text{p}}\right)\hfill \\ \hfill & \quad +{n}_{\text{e}}{\nu }_{\text{i}}-{S}_{{r}_{\text{e}-\text{p}}}-{S}_{{r}_{\text{n}-\text{p}}},\hfill \end{align} \tag{ 1.2 }$$

$$\begin{equation}\frac{\partial {n}_{\text{n}}}{\partial t}=\nabla \cdot \left({D}_{\text{n}}\enspace \nabla {n}_{\text{n}}\right)+\nabla \cdot \left({\boldsymbol{v}}_{\text{air}}+{\mu }_{\text{n}}\boldsymbol{E}{n}_{\text{n}}\right)+{n}_{\text{e}}{\nu }_{\text{a}}-{S}_{{r}_{\text{n}-\text{p}}},\end{equation} \tag{ 1.3 }$$

where n_e, n_p, and n_n are the densities of electrons, positive ions and negative ions and E is the electric field vector. The loss terms ${S}_{{r}_{\text{e}-\text{p}}}$ and ${S}_{{r}_{\text{n}-\text{p}}}$ describe the charged particle losses due to electron–ion and ion–ion recombination processes [21]. The scalar parameters D and μ are the diffusion and mobility coefficients, and ν_i, and ν_a are the ionization and attachment rate coefficients. These swarm parameters are approximated based on the local field and they are calculated using the Bolsig + software or determined from experiments [22–28].

We assume that the superposition law holds and the velocities of particles can be approximated by adding the air velocity v _air to the drift velocity. The drift–diffusion assumption for the velocity is valid for electrons as the time constant of momentum transfer is short compared to timescales of electric field variations in the streamer head [7, 29]. Although the ions have a higher inertia than electrons, we included the drift–diffusion model for ions in order to provide accurate results at microsecond timescales. The model is justified since the ion velocities are smaller than that of electrons and the ion displacement in streamer timescales is negligible. The effect of eventual delay to relax to the drift–diffusion velocity does not have significance in long timescales.

Poisson's equation for electrostatics (1.4) is used for the calculation of the electric potential:

$$\begin{equation}{\nabla }^{2}\phi =-\frac{{\rho }_{\text{c}}}{{{\epsilon}}_{0}},\end{equation} \tag{ 1.4 }$$

where ϕ is the electric potential, ₀ is the vacuum permittivity, and ρ_c is the charge density. The electric field is the gradient of the electrical potential and is calculated through:

$$\begin{equation}\boldsymbol{E}=-\nabla \phi .\end{equation} \tag{ 1.5 }$$

Photoionization effects are modeled based on the approximation proposed by Luque et al [30] of the integral model of photoionization from Zhelezniak et al [31]:

$$\begin{equation}{S}_{\text{ph}}\left({\boldsymbol{r}}_{\boldsymbol{1}}\right)={\int }_{{V}_{1}}\mathrm{d}{V}_{1}{\int }_{{V}_{2}}\frac{1}{4\pi }\frac{I\left({\boldsymbol{r}}_{\boldsymbol{2}}\right)f\left(\vert {\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\vert \right)}{\vert {\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}{\vert }^{2}}\enspace \mathrm{d}{V}_{2},\end{equation} \tag{ 1.6 }$$

where I( r ₂ ) is the emission of photons from differential volume at r ₂ , and S_ph( r ₁ ) is the photoionization source term at r ₁ . The method proposed by Luque et al [30] approximates the absorption function, f(| r ₁ − r ₂ |), by a series of exponential functions:

$$\begin{equation}\frac{f\left(\vert {\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\vert \right)}{{p}_{{\text{O}}_{2}}}={p}_{{\text{O}}_{2}}\vert {\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\vert \sum\limits _{j=1}^{3}{A}_{j}{e}^{-{\lambda }_{j}{p}_{{\text{O}}_{2}}\vert {\boldsymbol{r}}_{\boldsymbol{1}}-{\boldsymbol{r}}_{\boldsymbol{2}}\vert },\end{equation} \tag{ 1.7 }$$

where ${p}_{{\text{O}}_{2}}$ is the partial pressure of molecular oxygen. The integral equation, (1.6), can then be calculated by solving a set of Helmholtz differential equations. The fitting parameters are taken from [32]. The set of Helmholtz equations are [32]:

$$\begin{equation}{\nabla }^{2}{S}_{{\text{ph}}_{j}}\left(\boldsymbol{r}\right)-{\left({\lambda }_{j}{p}_{{\text{O}}_{2}}\right)}^{2}{S}_{{\text{ph}}_{j}}\left(\boldsymbol{r}\right)=-{A}_{j}\enspace {p}_{{\text{O}}_{2}}^{2}I\left(\boldsymbol{r}\right)\quad j=\left\{1,2,3\right\},\end{equation} \tag{ 1.8 }$$

with I( r ) defined as:

$$\begin{equation}I\left(\boldsymbol{r}\right)=\frac{{p}_{\text{q}}}{p+{p}_{\text{q}}}\xi \frac{{\nu }_{u}}{{\nu }_{\text{i}}}{S}_{\text{i}}\left(\boldsymbol{r}\right),\end{equation} \tag{ 1.9 }$$

where p_q/(p + p_q) is a quenching factor, ξ is the photionization efficiency, ν_u is the electron impact excitation frequency of level u, and S_i( r ) is the electron impact ionization source term [32]. Finally, the photoionization source term is:

$$\begin{equation}{S}_{\text{ph}}=\sum\limits _{j=1}^{3}{S}_{{\text{ph}}_{j}}.\end{equation} \tag{ 1.10 }$$

The general framework for implementing equations comes from Ansys Fluent which allows us to implement user-defined scalar (UDS) transport equations [33]. We used this capability to implement the drift–diffusion equations for charged particles, Poisson's equation for the electric potential and the Helmholtz equations for photoionization. The general form of the transport equation provided by Ansys Fluent is [33]:

$$\begin{equation}\frac{\partial \phi }{\partial t}+\frac{\partial }{\partial {x}_{\text{i}}}\left(v\phi -D\frac{\partial \phi }{\partial {x}_{\text{i}}}\right)={S}_{\phi },\end{equation} \tag{ 2.1 }$$

where, v is the velocity of the transport of the variable ϕ, D is the diffusion coefficient, and S_ϕ is the source term.

The drift–diffusion equations of charged particles have been implemented using the general format of the transport equation (2.1).

Poisson's equation and the Helmholtz equations have been implemented using the transport equation (2.1). We removed the unsteady and advection terms from (2.1) and replaced the source term with the source terms from (1.4) and (1.8). The diffusion term in (2.1) is discretized by a central difference scheme.

2.1. Spatial discretization

Ansys Fluent allows the use of several different schemes to calculate the advection flux. In our simulations, we used a second order upwind scheme and a slope limiter to bound the face value to the values of the cells straddling the face.

$$\begin{equation}{n}_{f}={n}_{{\mathrm{c}}_{\mathrm{0}}}+\zeta \enspace \nabla {n}_{{\mathrm{c}}_{\mathrm{0}}}\cdot \boldsymbol{\delta }{\boldsymbol{r}}_{{\mathrm{c}}_{\mathrm{0}}-\boldsymbol{f}},\end{equation} \tag{ 2.2 }$$

where ${n}_{{\mathrm{c}}_{\mathrm{0}}}$ and n_f are respectively the value of the variable at the upwind cell center and face center, and vector $\boldsymbol{\delta }{\boldsymbol{r}}_{{\mathrm{c}}_{\mathrm{0}}-\boldsymbol{f}}$, connects the cell center to the face center. ζ is a scalar between 0 and 1 and is calculated as follows:

$$\begin{equation}\zeta =\text{Max}\left[\text{Min}\left(1,\frac{{n}_{c1}-{n}_{c0}}{\nabla {n}_{{\mathrm{c}}_{\mathrm{0}}}\cdot \boldsymbol{\delta }{\boldsymbol{r}}_{{\mathrm{c}}_{\mathrm{0}}-\boldsymbol{f}}}\right),0\right],\end{equation} \tag{ 2.3 }$$

where ζ is chosen such that the value of the n_f would always be bounded between n_c0 and n_c1. When ζ = 1 we have a 2nd order upwind scheme, and when ζ = 0 we have a 1st order upwind scheme.

The discretization of the diffusion flux is done using the central difference, a second order scheme. This scheme assumes a linear profile for the variable between the center of the two cells straddling the face.

The gradient calculation is carried out using the least square cell based method [34]. This method provides comparable accuracy with lower computation time than Green Gauss node based method [35, 36]. It preserves a 1st order spatial accuracy in an irregular mesh and a 2nd order spatial accuracy in a regular mesh [35].

2.2. Time discretization

Ansys Fluent solves the UDS transport equations in an implicit way, and allows the implementation of user defined time integration schemes. In this paper, we present results from two different integration schemes. The first is the classical 2nd order backward difference scheme that provides accurate results in small time steps. The second is an implicit ELP scheme designed to allow the use of long time steps.

The 2nd order backward difference scheme is not accurate in high electric field regions especially if the time step size is large. To remedy this shortcoming, we implemented an implicit ELP integration scheme.

We have implemented an integration scheme for electrons and ions following [37], which uses an exact integration for the linear part describing charged particles growth and decay (ELP). The motivation behind this is the error in high electric field regions around streamer tip [1, 4, 38]. If the time step is not small enough, as compared to ionization timescale, the second order scheme introduces a relatively large error. The implicit ELP integration scheme integrates the source term due to direct electron impact ionization in an exact manner, and therefore does not suffer the inaccuracies in high electric field regions at large time steps.

To derive the ELP scheme, we write the equation describing the electron dynamics:

$$\begin{equation}\dot {{n}_{\text{e}}}=\left({\nu }_{\text{i}}-{\nu }_{\text{a}}\right){n}_{\text{e}}+{N}_{\text{e}},\end{equation} \tag{ 2.4 }$$

where n_e is the electron density and N represents the terms due to photoionization, recombination and transport fluxes. We integrate the linear growth and decay term, (ν_i − ν_a)n_e, in an exact manner and the other term, N, is integrated using an implicit 2nd order scheme. The derivation of the ELP integration scheme for electron density, n_e, is similar to the one described in [37], which we modified to incorporate adaptive time stepping. The derivation of the scheme and the terms that are required to be passed to Fluent are summarized in appendix A.

ELP schemes are not directly applicable to ions. Ion growth and decay dynamics include terms that are not linear w.r.t. ion densities, but linear to electron density.

$$\begin{equation}\dot {{n}_{\text{i}}}={\nu }_{\text{i}}{n}_{\text{e}}+{N}_{\text{i}}.\end{equation} \tag{ 2.5 }$$

The ion dynamics has a large term due to electron impact ionization, and is solved by taking into account the ionization from the linear part of (2.4). The linear part of the growth and decay dynamics is integrated in an exact manner, and the other parts of the dynamics are integrated implicitly with second order accuracy. A detailed derivation of the integration scheme for ions is also provided in appendix A. The advantages of using the implicit ELP scheme over the second order backward difference scheme are presented in section 3.

The spatial and temporal discretization presented above lead to a pseudo-linear system of equations:

$$\begin{equation}\boldsymbol{A}\left({\boldsymbol{x}}_{\boldsymbol{n}+\mathrm{1}}\right){\boldsymbol{x}}_{\boldsymbol{n}+\mathrm{1}}=\boldsymbol{b}\left({\boldsymbol{x}}_{\boldsymbol{n}+\mathrm{1}},{\boldsymbol{x}}_{\boldsymbol{n}}\right),\end{equation} \tag{ 2.6 }$$

where the matrix A includes the coefficients of variables and vector b includes the source terms and the explicit contributions from previous time steps. Both A and b are dependant on the solution variables at the current time step. Fluent solves this system of equations in an iterative way, allowing A and b to update after each iteration. The algebraic multigrid method (AMG) has been used to solve the set of linearized equations at each iteration [39].

The solution convergency at the end of each time step is checked based on the scaled residuals defined by Ansys Fluent [33].

2.3. Adaptive mesh refinement

Due to the multiscale nature of streamer discharges [40], it is crucial for 3D streamer simulations to use adaptive mesh refinement. A typical length of charge separation in a streamer head is 4 μm which can be several orders of magnitude smaller than the size of the domain, and therefore without adaptive mesh refinement the computation time would be extremely high. In this work, we use the adaptive mesh refinement from Ansys Fluent that is based on the hanging node method [33].

An important aspect of adaptive mesh refinement is the criterion to choose the refinement region. Our criterion for refinement is based on the parameter:

$$\begin{equation}{R}_{\text{p}}=\frac{{S}_{{\text{ph}}_{2}}\enspace \vert \boldsymbol{E}\vert }{\mathrm{max}\left({S}_{{\text{ph}}_{2}}\enspace \vert \boldsymbol{E}\vert \right)},\end{equation} \tag{ 2.7 }$$

where ${S}_{{\text{ph}}_{2}}$ is the 2nd term of the photoionization source (1.10) and | E | is the electric field magnitude. The contour plots of R_p is shown in figure 1 (left). Including both photoionization and electric field has two advantages. Firstly, the refiner does not refine the streamer body due to low electric field inside the streamer body. Secondly, it refines the high electric field region around the electrode only when there is high electron density near the electrode in the form of a patch or a streamer, otherwise it skips the regions around the electrode. Consequently, the high density mesh is focused in the streamer head and the region in front of it.

Zoom In Zoom Out Reset image size

The algorithm works as follows:

• We evaluate the value of max$\left({S}_{{\text{ph}}_{2}}\enspace \vert \boldsymbol{E}\vert \right)$ in the domain, and then we calculate the value of R_p from (2.7) for every cell.
• If the value of R_p in a cell is larger than some pre-specified values $\left\{{R}_{{\text{p}}_{1}},\enspace {R}_{{\text{p}}_{2}},\enspace {R}_{{\text{p}}_{3}}\right\}$, then the cell's volume should be smaller than a pre-specified value $\left\{{V}_{{\text{p}}_{1}},\enspace {V}_{{\text{p}}_{2}},\enspace {V}_{{\text{p}}_{3}}\right\}$, otherwise the cell is considered large, and is marked for refinement. The mesh refinement layers in front of the streamer are shown by {p₁, p₂, p₃}, and $\left\{{V}_{{\text{p}}_{1}},\enspace {V}_{{\text{p}}_{2}},\enspace {V}_{{\text{p}}_{3}}\right\}$ are the corresponding maximum cell volume allowed in that refinement level.
• The volumes $\left\{{V}_{{\text{p}}_{1}},\enspace {V}_{{\text{p}}_{2}},\enspace {V}_{{\text{p}}_{3}}\right\}$, are chosen in advance such that the maximum dimension of the cell is smaller than the ionization length (1/α), where α is the townsend coefficient of ionization.

The above algorithm leads to different refinement levels covering the streamer tip, resulting in a mesh density that gradually increases toward the streamer tip. Figure 2 shows the mesh around the streamer tip and the different mesh density levels.

Zoom In Zoom Out Reset image size

Figure 2 (left) shows that the refinement region fully covers the streamer tip. However, the refinement region also spreads to the side of the streamer tip which is considered unnecessary, as the streamer is not propagating in that direction. To improve the algorithm, we have added a second criterion that prevents refinement of cells in these regions, which is based on the optimization parameter ${R}_{{\text{p}}_{\text{opt}}}$ defined as:

$$\begin{equation}{R}_{{\text{p}}_{\text{opt}}}=\frac{\nabla {S}_{{\text{ph}}_{2}}\cdot {\boldsymbol{v}}_{\text{e}}}{\vert \nabla {S}_{{\text{ph}}_{2}}\vert \enspace \vert {\boldsymbol{v}}_{\text{e}}\vert }\end{equation} \tag{ 2.8 }$$

and is in fact the cosine of the angle between the gradient of the photoionization source term and electrons velocity vector. Figure 1 (right) shows the contours of ${R}_{{\text{p}}_{\text{opt}}}$ around the streamer, and that the unnecessary mesh refinement on streamer sides can be avoided if we prevent refinement of cells where $\vert {R}_{{\text{p}}_{\text{opt}}}\vert {< }0.75$. Figure 2 compares the mesh around the streamer tip before (left) and after (right) applying this improvement. We note that the above criteria do not depend on the direction of the electron velocity, the refinement algorithm is therefore applicable for both positive and negative streamers.

We also have an extra criterion that refines the streamer body up to 10 µm. The criterion is based on the positive ion density such that only cells having a density above a given threshold are chosen for refinement.

Finally, a criterion has been added to moderately refine the mesh at the vicinity of the electrode; it is simply based on the distance between the cell and the electrode.

The accuracy of the model is verified by comparing the model results to other streamer codes. Bagheri et al [17] have performed a comparison of six streamer codes from different research groups. They compared the accuracy and performance of the codes with respect to three test cases. Here, we have simulated the third test case, which is an axisymmetric simulation of a streamer discharge in air. The configuration of the simulation domain is shown in figure 3. In order to enhance the electric field, an spherical Gaussian patch of positive ions with a max density of 5 × 10¹⁸ m⁻³ and an e-folding radius of 0.4 mm is placed in the domain. The top plane electrode is at potential ϕ = 18.75 kV and the lower electrode is grounded. The transport and rate coefficients, and photoionization parameters are calculated in the same way as in [17]. It is to be noted, that we performed the simulation using a 3D model instead of an axisymmetric one. Furthermore, we adopt the same notation as Bagheri et al [17] to refer to different streamer codes either by the name of the institute or the country of origin.

Zoom In Zoom Out Reset image size

We performed three simulations with different time step and minimum mesh sizes. The first one uses a time step size of 2.5 ps and a minimum mesh size of 2.5 µm; it is tagged 'DTU–2 ps–2.5 μm' in figure 4. In the second simulation ('DTU–2 ps–5 μm'), we used a time step size of 2 ps, and a minimum mesh size of 5 μm. In the third simulation ('DTU–10 ps–5 μm'), the time step and the minimum mesh size were respectively 10 ps and 5 μm. Figure 3 (2nd and 3rd panels) shows the contour plots of electron density and electric field at 15 ns from 'DTU–2 ps–5 μm' simulation.

Zoom In Zoom Out Reset image size

In figure 4, we compare some important characteristics of streamer simulations with the study by Bagheri et al [17]. Figure 4 (top left) compares the maximum electric field at the tip of the streamer at different lengths. The maximum electric field from our model is in good agreement with the results from CWI and DE. The results from 'DTU–2 ps–2.5 μm' is almost identical to the results from the DE model. Interestingly, the result from 'DTU–10 ps–5 μm" simulation is closer to the results from 'DTU–2 ps–2.5 μm' than the result from 'DTU–2 ps–5 μm' simulation. The reason is that the 10 ps time step size corresponds to a CFL value in the streamer tip region that is close to 1, which leads to lower numerical diffusion and dispersion, and hence increases the accuracy of the solution.

Figure 4 (top right) shows the variation of L(t) − vt where L(t) is the distance between the maximum electric field and the upper electrode and v is an estimate of the propagation velocity set to 0.06 cm ns⁻¹. The results from CWI, DE and 'DTU–2 ps–2.5 μm' are almost identical. In figure 4 (bottom left) the total number of electrons in the domain is shown vs time. The results from DTU model are in agreement with the results from the CWI and DE models. And finally, figure 4 (bottom right) compares the conservation of charge in different codes by plotting the total charge Q_tot in the domain. The DTU model conserves the electric charge even when very large time steps are used. This attribute is mostly due to the use of the implicit ELP integration scheme in time, where the scheme uses exact integration for ionization and attachment source terms.

The comparison of the models showed that the results from our model are in agreement with other streamer codes, specifically with the codes from CWI and DE. Our code did exceptionally well in conserving the charges when using long time steps. This comparison demonstrate the accuracy of the model to simulate streamer simulation in steady air. In the next section, we simulate a long-duration electric discharge in unsteady air.

The settings and the run times of the DTU model and of the models from [17] are reported in tables 1 and 2. For DTU model, we performed the same simulation in an axisymmetric configuration with both fixed mesh and adaptive mesh refinement. It should be mentioned that implicit schemes require more computational time per time step; however, they are not limited by the CFL requirement or the dielectric relaxation time as opposed to explicit models, and can use longer time steps [41]. Specifically, as we will see in the next section, the implicit models can be very effective in pulsed streamers simulations by using long time steps between the streamer pulses.

In order to compare the performance of the time integration schemes, we repeated the previous simulations using the 2nd order backward difference scheme.

Figure 5 compares the results of the implicit integration scheme to the results obtained using the implicit ELP scheme along with the results from CWI and DE models. Figure 5 (left) compares the streamer tip electric field vs streamer length. Figure 5 (right) shows the variation of L(t) − vt in time. The implicit and implicit ELP schemes perform similarly at the small time step size of 2 ps. However, increasing the time step size to 10 ps leads to an inaccurate solution for the 2nd order implicit integration scheme, while the implicit ELP scheme remains accurate.

Zoom In Zoom Out Reset image size

To conclude the comparison, both schemes are stable when taking large time steps. However, inaccuracies develop when using the 2nd order implicit integration scheme. The implicit ELP integration scheme remains stable and accurate while using long time steps, which is essential to model long-duration discharges. Furthermore, the implicit ELP scheme allows the solution to converge faster in each time step, leading to a lower computation cost.

We modeled the behavior of positive streamers in the presence of lateral wind, based on an experiment performed in [14]. The experiment shows that the positive streamer corona when subjected to a transverse wind, tilts in the direction of the wind. Here, we have devised a simulation to capture the effect of the airflow, and investigate the main mechanisms that affect positive streamers in low speed airflow. This simulation demonstrates the ability of our model to simulate long discharges coupled with unsteady airflow, and also helps us understand the physical mechanism explaining the tilting.

First, we need to point out some peculiarities of this simulation:

• (a)  
  The effect of the airflow on electric discharges reveals itself at the timescales corresponding to heavy particles, and therefore the duration of this simulation needs to be long (around 0.1 ms), which renders the simulation very challenging. However, the implicit ELP integration scheme allows us to take very large time steps and makes this simulation feasible.
• (b)  
  We use a simplified chemical kinetic model to lower the computation time knowing that the duration of our simulation is in the order of milliseconds, and that intermediate ionic species (such as ${\mathrm{N}}_{2}^{+}$, ${\mathrm{O}}_{2}^{+}$, ${\mathrm{N}}_{4}^{+}$, ${\mathrm{O}}_{2}{\mathrm{N}}_{2}^{+}$, ...) are converted to ${\mathrm{O}}_{4}^{+}$ on a timescale of few nanoseconds [21, 42–44]. Therefore, the charged particles considered in the model (1.1) and (1.2) are electrons, positive ${\mathrm{O}}_{4}^{+}$ ions and negative ${\mathrm{O}}_{2}^{-}$ ions, for which we calculate the ion mobility, diffusion coefficient and recombination rate. The effects of detachment processes and water content are neglected in this simulation in order to focus on understanding the effects of airflow.
• (c)  
  In this simulation, we have streamers that reach the stagnation point where the electric field in their tip fades away, the propagation stops and the streamer vanishes. Previously, Pancheshnyi [43] and Starikovskiy [44], using a first order fluid model and local field and density approximation, simulated positive streamers that eventually came into stagnation. In their simulations, the streamer head became narrower while the electric field at the streamer tip grew indefinitely and finally the simulation could not continue [43, 44]. We performed an axisymmetric simulation of a positive streamer using the local field and density approximation for ionization and experienced the same behavior, see figure 6. The electric field in the streamer tip grows unboundedly and eventually the simulation stops, similar to the behavior observed in [43, 44]. However, the local field and density approximation of the ionization source term is not accurate at high fields, as a large contribution to ionization comes from fast, high energy electrons which will spread the ionization over a significant length. Therefore, we smooth the ionization source term in high electric field regions to mimic the non-local aspect of electron ionization. We solved the following equation to smooth the ionization source term:

  $$\begin{equation}{L}^{2}{\nabla }^{2}{S}_{\text{ni}}-{S}_{\text{ni}}=-{n}_{\text{e}}{\nu }_{\text{i}},\end{equation} \tag{ 4.1 }$$

  where L is the smoothing length scale, and S_ni is the smoothed ionization source term, used in the drift–diffusion equations for electrons and positive ions (1.1) and (1.2). We have set the boundary condition at electrode surface for the smoothing equation to Neumann zero (∂S_ni/∂n = 0). We can note that the total number of ionization events in the domain is not altered by this smoothing process. By integrating both sides of (4.1) over the entire domain, we have:

  $$\begin{equation}{\int }_{V}{L}^{2}{\nabla }^{2}{S}_{\text{ni}}\enspace \mathrm{d}V-{\int }_{V}{S}_{\text{ni}}\enspace \mathrm{d}V=-{\int }_{V}{n}_{\text{e}}{\nu }_{\text{i}}\enspace \mathrm{d}V\end{equation} \tag{ 4.2 }$$

  and,

  $$\begin{equation}{\int }_{V}{S}_{\text{ni}}\enspace \mathrm{d}V={\int }_{V}{n}_{\text{e}}{\nu }_{\text{i}}\enspace \mathrm{d}V,\end{equation} \tag{ 4.3 }$$

  since the integral of the diffusion term over the volume cancels. Figure 7 compares the ionization source term S_i using the local field approximation in the streamer tip with the smoothed ionization source term S_ni. The smoothing length scale, L, is limited to 4 μm in high field regions, and as shown in figure 7 (top), the smoothing process slightly changes the ionization source term near the tip, and as we get away from the streamer tip, the smoothing effects disappear. The smoothed ionization source term solves the issue of the unbounded electric field for a non-propagating streamer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 8, we show the simulation domain in a point-plane electrode configuration with the point anode at 8.5 kV electric potential. In order to experiment with wind speed magnitude coherent with the ones found near the tip of wind turbine blades, we set a lateral wind of velocity v_air = 50 m s⁻¹, blowing transversely. To initiate the first streamer, we place a neutral Gaussian plasma patch of electrons and positive ions on the tip of the electrode to initiate the streamer. The maximum density of the patch is 10¹⁴ m⁻³ and the e-folding radius is 0.14 mm. The charged particles boundary condition at the point electrode surface is set to absorb all the incoming charge fluxes, and no out flux.

Zoom In Zoom Out Reset image size

The air streamlines around the electrode are shown in figure 9. Compressible Navier Stokes equations for an ideal gas along with shear-stress transport k–ω turbulent model [45] is used to calculate the airflow around the electrode.

Zoom In Zoom Out Reset image size

The 8.5 kV electric potential applied to the point electrode leads to a maximum electric field of 110 kV cm⁻¹ on the electrode tip. The streamer emerges from the plasma patch and propagates in the inhomogeneous field created by the anode while a uniform lateral wind is blowing to the streamer. During the initial phase after the streamer initiation, the streamer tip electric field rises to around 190 kV cm⁻¹ and the streamer velocity V_s reaches 0.45 mm ns⁻¹ (see figure 10). As the streamer length increases, the streamer thickness shrinks and the electric potential in the streamer tip decreases, such that there is a steady decline in the electric field and streamer velocity. The streamer propagation slows down till it reaches stagnation, at this point as seen in figure 10 (right) the electric field in the streamer tip fades out.

Zoom In Zoom Out Reset image size

The plots of the electron density in figure 11, show the propagation of the streamer over time. The last panel in figure 11 shows the electron density at 50 μs, which is long after the streamer has stopped propagating and clearly shows the effect of the airflow and electron decay.

Zoom In Zoom Out Reset image size

The streamer reaches a maximum channel length of 3.15 cm. In order to compare the channel length to the one measured in experiments, we try to choose an experiment with a single streamer to avoid the streamer interaction effects. However, it seems most experiments in atmospheric pressure lead to a collection of streamers rather than single streamers [5]. Here, we have chosen an experiment from [46] that produces relatively few streamers, and hence the streamer interactions are low. In one of their experiments, they captured positive streamer discharges in a gas composition of 20% oxygen and 80% nitrogen at a pressure of 200 mbar. In the experiment, the electrode with an 11 kV potential initiated a few streamers that crossed a 16 cm gap. If we scale down the streamer length using the Townsend similarity law to account for the low pressure, we will have an equivalent streamer length of 3.2 cm in atmospheric pressure. This result in an electrode potential to streamer length ratio of 3.4 kV cm⁻¹, the corresponding ratio in our simulation is around 2.7 kV cm⁻¹, which is in good agreement. Our results are also in good agreement with the results from [44], where they report an average electric field inside the streamer of 3.25 kV cm⁻¹.

After the streamer stops propagating and the electric field in the streamer tip vanished, the simulation enters the quiescent phase. In this phase, the density of the charged particles created mostly by the 1st streamer are continuously decreasing. The main mechanisms removing charged particles are electron absorption by the anode, attachment, and ion–ion and electron–ion recombination. The total number of charged particles in the discharge gap is plotted in figure 12 (left), and we can notice a large decrease in the total number of particles after the 1st streamer. The electron density inside the domain is reduced to a greater extent than the positive ions, mostly due to absorption by the anode. A similar behavior has been observed in [18]. The positive charge density remainder from the 1st streamer dampens the field around the electrode, inhibiting the emergence of a new streamer. Figure 12 (right), shows that the electric field near the electrode rises over time as the positive ions density is reduced and pushed further away from the anode. Note that the implicit ELP integration scheme was essential in performing this simulation, as we had to use very long time steps.

Zoom In Zoom Out Reset image size

During the quiescent phase, we observe a very interesting effect. The plasma channel remaining from the initial streamer, remains attached to the electrode surface while being moved by the air, as shown in figure 13. The streamer remains attached to the surface, notwithstanding the separation of the airflow from the electrode surface. The reason for this behavior is that the electrons inside the streamer channel are continuously moving toward the electrode, and they ionize the air in their path and maintain the channel (see figure 13). Therefore, the streamer channel remains constantly connected to the electrode surface. It should be noted that the electric field inside the streamer channel near the electrode, due to charge removal, is high enough to provide sufficient ionization, which is shown in the ionization source plots in figure 13 (bottom row). The significance of this sustained plasma channel connection will be understood in the next section, where we discuss the initiation of a new streamer from the remaining charges in this plasma channel.

Zoom In Zoom Out Reset image size

During the quiescent phase the electric field keeps rising as shown in figure 12 (right), and eventually it will reach the Laplacian field in the absence of screening charges allowing the possible initiation of a 2nd streamer. However, the rise in the electric field was slow and in order to speed up the simulation, after 50 μs, we increased artificially the electrode potential to 10.5 kV with a rise time of 100 ns. Then the charge density remaining from the 1st streamer and transported by the wind over the surface of the electrode, developed into a new streamer. The plots of electron density in figure 14 clearly show the initiation of a streamer from the residual charges of the 1st streamer.

Zoom In Zoom Out Reset image size

Initially, the streamer starts to follow the pre-ionized path remaining from the 1st streamer, as positive streamers tend the follow the region with a higher free electron density [47]. Later, as the streamer move away from the electrode, the electron density in the pre-ionized path decreases and the streamer starts to change direction and follows the field lines instead. In figure 15, the plots of electron density and electric field clearly demonstrate the curved path of the 2nd streamer. In the final 1 cm of the 2nd streamer, it shifts slightly toward left due to the field created by the residual charges of the 1st streamer. The new streamer propagates for 3.75 cm, which is longer than the 1st streamer as the potential of the electrode has been increased from 8.5 kV to 10.5 kV. Figure 16, shows the evolution of the electric field and streamer velocity for the 2nd streamer. The maximum length and the propagation time of the 2nd streamer is longer than the 1st one due to a higher electrode potential. Similarly to the 1st streamer, the 2nd streamer stops propagating, and another quiescent phase starts with a sharp reduction of the charged particle densities and a steady increase in electric field near the electrode. We modeled this phase until 104 μs, which show an even more pronounced tilting of the 2nd streamer plasma channel in the direction of the wind. The channel also remains attached to the electrode.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A 3D model is presented to couple the dynamics of streamer discharges and unsteady airflow. The model allows us to investigate the effects of unsteady airflow on electric discharge characteristics at different timescales. The effects of unsteady airflow on streamer discharges through charge redistribution manifest in longer timescales than the typical timescales of a streamer discharge. We therefore implemented a model that allows large time step sizes by using an implicit ELP time integration scheme. The model also uses an unstructured mesh allowing the incorporation of solid bodies with complex geometries. Comparing the results of the model with a test case from literature, shows a good agreement. The ability of the model to couple streamer discharges and unsteady airflow was demonstrated by a simulation of the initiation and termination of two successive positive streamers in transverse airflow.

The results of the simulation can be summarized by several points:

• During its propagation, the 1st positive streamer had only a small tilting in the direction of the wind due to its short duration. It shows that the effect of the low speed airflow on the electric discharges manifests itself on longer timescales than the ones of the streamer propagation driven by electron timescales.
• After the stagnation of the 1st streamer, the streamer channel is moved by the air. However, it remains attached to the surface of the electrode, even though the air streamlines suggest that the streamer channel should separate from the electrode surface. This is an intrinsic characteristic of positive streamers and is due to the flow of electrons from the streamer channel toward the electrode and subsequent ionization inside the plasma channel near the electrode surface. It is uncertain that negative streamers remain similarly attached since electrons move away from the electrode.
• The subsequent streamer was initiated from the residual charges of the 1st streamer close to electrode surface. These charges have been swept across the electrode surface by the wind. The results suggest that positive streamers in lateral wind have a tendency to emerge from the trailing edge of the electrode.
• The 2nd streamer, emerging from the trailing edge of the electrode, shows a clear tilting in the direction of the wind. The tilting later intensifies, as the streamer channel keeps moving in the direction of the wind.

This simulation shows that the effects of the low speed airflow on positive streamers mainly result from the charged particle redistribution in the domain which will affect the subsequent streamer discharges.

This project has received funding from the European Unions Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement 722337. The computational resources were provided by DTU Computing Center (DCC).

The data that support the findings of this study are available upon reasonable request from the authors.

This appendix gives a derivation of the integration scheme for charged particles and describe how it is implemented in Ansys Fluent. The derivation of the integration scheme for electrons is similar to [37]. In the following, we show the modifications we brought to the scheme to allow a varying time step size and its implementation in Ansys Fluent. The derivation of the integration scheme for ions also follows [37] with a higher level of complexity coming from the fact that the linear term in ion dynamics is proportional to electron density and not the ion density.

A.1. 2nd order implicit ELP integration scheme for electrons

Starting from the Taylor series expansion of the electron density ${n}_{{\text{e}}_{n+1}}$ at the time step n + 1 given by:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}=\sum\limits _{k=0}^{\infty }{n}_{{\text{e}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{k}}{k!},\end{equation} \tag{ A.1 }$$

where the jth time derivative of the electron density is given by:

$$\begin{equation}{n}_{{\text{e}}_{n}}^{\left(1\right)}={\lambda }_{n+\frac{1}{2}}{n}_{{\text{e}}_{n}}+{N}_{{\text{e}}_{n}},\end{equation} \tag{ A.2 }$$

$$\begin{equation}{n}_{{\text{e}}_{n}}^{\left(2\right)}={\lambda }_{n+\frac{1}{2}}^{2}{n}_{{\text{e}}_{n}}+{\lambda }_{n+\frac{1}{2}}{N}_{{\text{e}}_{n}}+{N}_{{\text{e}}_{n}}^{\left(1\right)},\end{equation} \tag{ A.3 }$$

$$\begin{align}{\vdots}\hfill \\ \hfill {n}_{{\text{e}}_{n}}^{\left(j\right)}& ={\lambda }_{n+\frac{1}{2}}^{j}{n}_{{\text{e}}_{n}}+\sum\limits _{m=0}^{j-1}{\lambda }_{n+\frac{1}{2}}^{j-m-1}{N}_{{\text{e}}_{n}}^{\left(m\right)},\hfill \end{align} \tag{ A.4 }$$

where ${\lambda }_{n+\frac{1}{2}}=\left({\nu }_{{\text{i}}_{n}}+{\nu }_{{\text{i}}_{n+1}}-{\nu }_{{\text{a}}_{n}}-{\nu }_{{\text{a}}_{n+1}}\right)/2$ is the effective ionization frequency averaged between steps n and n + 1, and ${N}_{{\text{e}}_{n}}$ represents the other source terms, which include the variations due to fluxes, photoionization and recombination. Replacing (A.4) inside the Taylor series expression (A.1), we have:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}=\sum\limits _{k=0}^{\infty }{\lambda }_{n+\frac{1}{2}}^{k}{n}_{{\text{e}}_{n}}\frac{{\Delta}{t}_{n}^{k}}{k!}+\sum\limits _{k=1}^{\infty }\sum\limits _{m=0}^{k-1}{\lambda }_{n+\frac{1}{2}}^{k-m-1}{N}_{{\text{e}}_{n}}^{\left(m\right)}\frac{{\Delta}{t}_{n}^{k}}{k!}.\end{equation} \tag{ A.5 }$$

The first term is the familiar exponential growth term and the second term can be simplified by changing the order of summation such that we have:

$$\begin{align}\hfill & \sum\limits _{m=0}^{\infty }{N}_{{\text{e}}_{n}}^{\left(m\right)}\sum\limits _{k=m+1}^{\infty }{\lambda }_{n+\frac{1}{2}}^{k-m-1}\frac{{\Delta}{t}_{n}^{k}}{k!}\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{m=0}^{\infty }\frac{{N}_{{\text{e}}_{n}}^{\left(m\right)}}{{\lambda }_{n+\frac{1}{2}}^{m+1}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-\sum\limits _{k=0}^{m}\frac{{\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)}^{k}}{k!}\right)\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{m=0}^{\infty }{N}_{{\text{e}}_{n}}^{\left(m\right)}{\Delta}{t}_{n}^{m+1}{Q}_{m+1}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right),\hfill \end{align} \tag{ A.6 }$$

with Q_j defined as [37] by:

$$\begin{equation}{Q}_{j}\left(x\right)=\frac{{e}^{x}-{E}_{j}\left(x\right)}{{x}^{j}},\hspace{25.0pt}{E}_{j}\left(x\right)=\sum\limits _{k=0}^{j-1}\frac{{x}^{k}}{k!}\end{equation} \tag{ A.7 }$$

In order to have a 2nd order implicit ELP integration scheme the integration scheme needs to equalize the Taylor series up to the 2nd order terms:

Assuming the following form for the integration scheme:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}={n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}+\gamma {\Delta}{t}_{n}{N}_{{\text{e}}_{n+1}}+\beta {\Delta}{t}_{n}{N}_{{\text{e}}_{n}}\end{equation} \tag{ A.8 }$$

and expanding ${N}_{{\text{e}}_{n+1}}$ in Taylor series, we have:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}={n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}+\left(\gamma +\beta \right){\Delta}{t}_{n}{N}_{{\text{e}}_{n}}^{\left(0\right)}+\sum\limits _{k=1}^{\infty }\gamma {\Delta}{t}_{n}^{k+1}{N}_{{\text{e}}_{n}}^{\left(k\right)}\end{equation} \tag{ A.9 }$$

(A.9) and (A.6) are then equal up to the 2nd order term if

$$\begin{equation}\gamma ={Q}_{2}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad \beta ={Q}_{1}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)-\gamma \end{equation} \tag{ A.10 }$$

The integration scheme to calculate ${n}_{{\text{e}}_{n+1}}$ now depends on N_n , which includes the flux terms and the source terms from the time step n. In order to have an implicit solver, we derive a similar expression for ${n}_{{\text{e}}_{n-1}}$ and combine the two expressions to remove the dependency on N_n . We proceed as follows:

The Taylor series expansion of the electron density at previous time step, ${n}_{{\text{e}}_{n-1}}$:

$$\begin{equation}{n}_{{\text{e}}_{n-1}}=\sum\limits _{k=0}^{\infty }{n}_{{\text{e}}_{n}}^{\left(k\right)}{\left(-1\right)}^{k}\frac{{\Delta}{t}_{n-1}^{k}}{k!},\end{equation} \tag{ A.11 }$$

after following a similar derivation path as that of ${n}_{{\text{e}}_{n+1}}$, we have:

$$\begin{equation}{n}_{{\text{e}}_{n-1}}={n}_{{\text{e}}_{n}}{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}+{\gamma }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n+1}}+{\beta }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n}},\end{equation} \tag{ A.12 }$$

where ${\lambda }_{n-\frac{1}{2}}=\left({\nu }_{{\text{i}}_{n}}+{\nu }_{{\text{i}}_{n-1}}-{\nu }_{{\text{a}}_{n}}-{\nu }_{{\text{a}}_{n-1}}\right)/2$ is the effective ionization frequency averaged between steps n and n − 1, and

$$\begin{equation*}{\gamma }^{\prime }=\frac{{\Delta}{t}_{n-1}}{{\Delta}{t}_{n}}{Q}_{2}^{\prime }\left({\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}\right),\end{equation*}$$

and

$$\begin{equation}{\beta }^{\prime }={Q}_{1}^{\prime }\left({\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}\right)-{\gamma }^{\prime },\end{equation} \tag{ A.13 }$$

with Q_j ' defined as:

$$\begin{equation}{Q}_{j}^{\prime }\left(x\right)=\frac{{e}^{-x}-{E}_{j}^{\prime }\left(x\right)}{{x}^{j}},\hspace{25.0pt}{E}_{j}^{\prime }\left(x\right)=\sum\limits _{k=0}^{j-1}{\left(-1\right)}^{k}\frac{{x}^{k}}{k!}.\end{equation} \tag{ A.14 }$$

We combine the approximations of ${n}_{{\text{e}}_{n+1}}$ and ${n}_{{\text{e}}_{n-1}}$ in such a way to remove the explicit dependence on ${N}_{{\text{e}}_{n}}$. The expressions of the ${n}_{{\text{e}}_{n+1}}$ and ${n}_{{\text{e}}_{n-1}}$ are:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}={n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}+\gamma {\Delta}{t}_{n}{N}_{{\text{e}}_{n+1}}+\beta {\Delta}{t}_{n}{N}_{{\text{e}}_{n}}.\end{equation} \tag{ A.15 }$$

$$\begin{equation}{n}_{{\text{e}}_{n-1}}={n}_{{\text{e}}_{n}}{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}+{\gamma }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n+1}}+{\beta }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n}}.\end{equation} \tag{ A.16 }$$

We multiply the first equation by scalar A and subtract it from the second equation, and then we choose the scalar A such that the coefficient of ${N}_{{\text{e}}_{n}}$ in the final equation would be zero. The final equation solved in terms of ${N}_{{\text{e}}_{n+1}}$ would be:

$$\begin{equation}{N}_{{\text{e}}_{n+1}}=\frac{A{n}_{{\text{e}}_{n+1}}-{n}_{{\text{e}}_{n}}\left(A{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}\right)-{n}_{{\text{e}}_{n-1}}}{\left(A\gamma {\Delta}{t}_{n}-{\gamma }^{\prime }{\Delta}{t}_{n-1}\right)},\end{equation} \tag{ A.17 }$$

where,

$$\begin{equation}A=\frac{{\beta }^{\prime }{\Delta}{t}_{n-1}}{\beta {\Delta}{t}_{n}}.\end{equation} \tag{ A.18 }$$

The drift–diffusion equation for electrons (1.1) can now be written as:

$$\begin{equation}{N}_{{\text{e}}_{n+1}}=\nabla \cdot \left({D}_{\text{e}}\enspace \nabla {n}_{\text{e}}\right)+\nabla \cdot \left({\boldsymbol{v}}_{\text{air}}+{\mu }_{\text{e}}\boldsymbol{E}{n}_{\text{e}}\right)-{S}_{{r}_{e-p}},\end{equation} \tag{ A.19 }$$

where the ionization and attachment source term are removed from the equation, and the unsteady term, $\frac{\partial {n}_{\text{e}}}{\partial t}$, is approximated by ${N}_{{\text{e}}_{n+1}}$.

At the initial time step we do not have the value of ${n}_{{\text{e}}_{n-1}}$. Therefore, we have to integrate the terms represented by ${N}_{{\text{e}}_{n+1}}$ using a first order scheme. The scheme to calculate the ${n}_{{\text{e}}_{n+1}}$ would be:

$$\begin{equation}{n}_{{\text{e}}_{n+1}}={n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}+{Q}_{1}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right){\Delta}{t}_{n}{N}_{{\text{e}}_{n+1}}.\end{equation} \tag{ A.20 }$$

Solving for ${N}_{{\text{e}}_{n+1}}$, we have:

$$\begin{equation}{N}_{{\text{e}}_{n+1}}=\frac{{n}_{{\text{e}}_{n+1}}-{n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}}{{Q}_{1}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right){\Delta}{t}_{n}}.\end{equation} \tag{ A.21 }$$

Ansys Fluent expects the unsteady term in the UDS transport equations to be decomposed into implicit and explicit terms, they are respectively called 'apu' and 'su' in the Fluent manual [48], and are given by:

$$\begin{equation}\mathrm{a}\mathrm{p}\mathrm{u}=\frac{-A}{\left(A\gamma {\Delta}{t}_{n}-{\gamma }^{\prime }{\Delta}{t}_{n-1}\right)}.\end{equation} \tag{ A.22 }$$

$$\begin{equation}\mathrm{s}\mathrm{u}=-\frac{-{n}_{{\text{e}}_{n}}\left(A{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}\right)-{n}_{{\text{e}}_{n-1}}}{\left(A\gamma {\Delta}{t}_{n}-{\gamma }^{\prime }{\Delta}{t}_{n-1}\right)}.\end{equation} \tag{ A.23 }$$

We have limited the use of the implicit ELP integration scheme to the regions of the domain where the electric field magnitude is strictly above the breakdown field, i.e. where the electron density grows exponentially. In the other regions, the integration is done by a classical 2nd order backward difference scheme. The two schemes will provide the same level of accuracy near the breakdown field, therefore there is no need for a special treatment of the interface between the cells using the exponential and the 2nd order backward difference schemes.

A.2. 2nd order implicit ELP integration scheme for ions

We now describe the derivation of the integration scheme for the ions. Starting from the Taylor series expansion of the ions density at the next time step, ${n}_{{\text{i}}_{n+1}}$:

$$\begin{equation}{n}_{{\text{i}}_{n+1}}=\sum\limits _{k=0}^{\infty }{n}_{{\text{i}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{k}}{k!}\end{equation} \tag{ A.24 }$$

where the jth time derivative of the ions density is given by:

$$\begin{equation}{n}_{{\text{i}}_{n}}^{\left(1\right)}={\nu }_{{\text{i}}_{n+\frac{1}{2}}}{n}_{{\text{e}}_{n}}+{N}_{{\text{i}}_{n}}^{\left(0\right)},\end{equation} \tag{ A.25 }$$

$$\begin{equation}{n}_{{\text{i}}_{n}}^{\left(2\right)}={\nu }_{{\text{i}}_{n+\frac{1}{2}}}{n}_{{\text{e}}_{n}}^{\left(1\right)}+{N}_{{\text{i}}_{n}}^{\left(1\right)},\end{equation} \tag{ A.26 }$$

$$\begin{align}\hfill & {\vdots}\hfill \\ \hfill {n}_{{\text{i}}_{n}}^{\left(j\right)}& ={\nu }_{{\text{i}}_{n+\frac{1}{2}}}{n}_{{\text{e}}_{n}}^{\left(j-1\right)}+{N}_{{\text{i}}_{n}}^{\left(j-1\right)}\hfill \end{align} \tag{ A.27 }$$

where ${\nu }_{{\text{i}}_{n+\frac{1}{2}}}$ is the averaged ionization frequency between steps n and n + 1, and ${N}_{{\text{i}}_{n}}$ represents the other terms which includes fluxes, photoionization and recombination. As it can be seen in (A.27), we also need to have the general expression for jth time derivative of the electron density:

$$\begin{equation}{n}_{{\text{e}}_{n}}^{\left(j\right)}={\lambda }_{n+\frac{1}{2}}^{j}{n}_{{\text{e}}_{n}}+\sum\limits _{m=0}^{j-1}{\lambda }_{n+\frac{1}{2}}^{j-m-1}{N}_{{\text{e}}_{n}}^{\left(m\right)}.\end{equation} \tag{ A.28 }$$

Replacing (A.27) in the Taylor series expansion (A.24), we have:

$$\begin{align}\hfill {n}_{{\text{i}}_{n+1}}& ={n}_{{\text{i}}_{n}}+{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}\sum\limits _{k=1}^{\infty }{\lambda }_{n+\frac{1}{2}}^{k-1}\frac{{\Delta}{t}_{n}^{k}}{k!}+{\nu }_{{\text{i}}_{n+\frac{1}{2}}}\sum\limits _{k=2}^{\infty }\hfill \\ \hfill & \quad {\times}\sum\limits _{j=0}^{k-2}{\lambda }_{n+\frac{1}{2}}^{k-j-2}{N}_{{\text{e}}_{n}}^{\left(j\right)}\frac{{\Delta}{t}_{n}^{k}}{k!}+\sum\limits _{k=1}^{\infty }{N}_{{\text{i}}_{n}}^{\left(k-1\right)}\frac{{\Delta}{t}_{n}^{k}}{k!},\hfill \end{align} \tag{ A.29 }$$

$$\begin{align}\hfill & ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)+{\nu }_{{\text{i}}_{n+\frac{1}{2}}}\sum\limits _{k=0}^{\infty }\qquad \qquad \hfill \\ \hfill & \quad {\times}\sum\limits _{j=0}^{k}{\lambda }_{n+\frac{1}{2}}^{k-j}{N}_{{\text{e}}_{n}}^{\left(j\right)}\frac{{\Delta}{t}_{n}^{k+2}}{\left(k+2\right)!}+\sum\limits _{k=0}^{\infty }{N}_{{\text{i}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{\left(k+1\right)}}{\left(k+1\right)!}.\hfill \end{align} \tag{ A.30 }$$

There are four different terms in the rhs The first term is the ion density in the current time step. The second term is a closed form analytical expression and needs no further simplification. The fourth term also needs no further simplification. The third term is simplified by changing the order of summation:

$$\begin{align}\hfill & \sum\limits _{j=0}^{\infty }{N}_{{\text{e}}_{n}}^{\left(j\right)}{\left({\Delta}{t}_{n}\right)}^{j+2}\enspace \frac{1}{{\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)}^{j+2}}\sum\limits _{k=j+2}^{\infty }\frac{{\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)}^{k}}{\left(k\right)!}\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{j=0}^{\infty }{N}_{{\text{e}}_{n}}^{\left(j\right)}{\left({\Delta}{t}_{n}\right)}^{j+2}{Q}_{j+2}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right),\hfill \end{align} \tag{ A.31 }$$

with Q_j (x) defined as [37] by:

$$\begin{equation}{Q}_{j}\left(x\right)=\frac{{e}^{x}-{E}_{j}\left(x\right)}{{x}^{j}},\hspace{25.0pt}{E}_{j}\left(x\right)=\sum\limits _{k=0}^{j-1}\frac{{x}^{k}}{k!}\end{equation} \tag{ A.32 }$$

Then ${n}_{{\text{i}}_{n+1}}$ reads:

$$\begin{align}\hfill {n}_{{\text{i}}_{n+1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\hfill \\ \hfill & \quad +{\nu }_{{\text{i}}_{n+\frac{1}{2}}}\sum\limits _{j=0}^{\infty }{N}_{{\text{e}}_{n}}^{\left(j\right)}{\left({\Delta}{t}_{n}\right)}^{j+2}{Q}_{j+2}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)\hfill \\ \hfill & \quad +\sum\limits _{k=0}^{\infty }{N}_{{\text{i}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{\left(k+1\right)}}{\left(k+1\right)!}.\hfill \end{align} \tag{ A.33 }$$

In order to have a 2nd order implicit ELP integration scheme the integration scheme needs to equalize the Taylor series up to the 2nd order terms. Assuming the following form of the integration scheme:

$$\begin{align}\hfill {n}_{{\text{i}}_{n+1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n}^{2}\left(\zeta {N}_{{\text{e}}_{n+1}}+\eta {N}_{{\text{e}}_{n}}\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n}\left({\gamma }_{\text{i}}{N}_{{\text{i}}_{n+1}}+{\beta }_{\text{i}}{N}_{{\text{i}}_{n}}\right).\hfill \end{align} \tag{ A.34 }$$

Expanding the ${N}_{{\text{i}}_{n+1}}$ and ${N}_{{\text{e}}_{n+1}}$ in Taylor series we have:

$$\begin{align}\hfill {n}_{{\text{i}}_{n+1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\qquad \qquad \hfill \\ \hfill & \quad +\left(\zeta +\eta \right){N}_{{\text{e}}_{n}}^{\left(0\right)}{\Delta}{t}_{n}^{2}+\zeta \sum\limits _{k=1}^{\infty }{N}_{{\text{e}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{k+2}}{k!}\hfill \\ \hfill & \quad +\left({\gamma }_{\text{i}}+{\beta }_{\text{i}}\right){N}_{{\text{i}}_{n}}^{\left(0\right)}{\Delta}{t}_{n}+{\gamma }_{\text{i}}\sum\limits _{k=1}^{\infty }{N}_{{\text{i}}_{n}}^{\left(k\right)}\frac{{\Delta}{t}_{n}^{k+1}}{k!}.\hfill \end{align} \tag{ A.35 }$$

And by equating (A.35) and (A.33) up to the 2nd order term, we have:

$$\begin{equation}\zeta ={\nu }_{{\text{i}}_{n+\frac{1}{2}}}{Q}_{3}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right),\end{equation} \tag{ A.36 }$$

$$\begin{equation}\eta ={\nu }_{{\text{i}}_{n+\frac{1}{2}}}\left({Q}_{2}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)-{Q}_{3}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)\right),\end{equation} \tag{ A.37 }$$

$$\begin{equation}{\gamma }_{\text{i}}=\frac{1}{2}\end{equation} \tag{ A.38 }$$

$$\begin{equation}{\beta }_{\text{i}}=\frac{1}{2}.\end{equation} \tag{ A.39 }$$

The integration scheme to calculate ${n}_{{\text{i}}_{n+1}}$ depends on ${N}_{{\text{e}}_{n+1}}$, ${N}_{{\text{e}}_{n}}$ and ${N}_{{\text{i}}_{n}}$. In order to have an implicit solver, we derive a similar expression for ${n}_{{\text{i}}_{n-1}}$ and combine this expression with the ones for ${n}_{{\text{i}}_{n+1}}$, ${n}_{{\text{e}}_{n+1}}$, and ${n}_{{\text{e}}_{n-1}}$ to remove the dependency on ${N}_{{\text{e}}_{n+1}}$, ${N}_{{\text{e}}_{n}}$ and ${N}_{{\text{i}}_{n}}$. We proceed as follow:

The Taylor series expansion of the ions density at the previous time step, ${n}_{{\text{i}}_{n-1}}$:

$$\begin{equation}{n}_{{\text{i}}_{n-1}}=\sum\limits _{k=0}^{\infty }{n}_{{\text{i}}_{n}}^{\left(k\right)}{\left(-1\right)}^{k}\frac{{\Delta}{t}_{n-1}^{k}}{k!},\end{equation} \tag{ A.40 }$$

after following a similar derivation path as that of ${n}_{{\text{i}}_{n+1}}$, we have:

$$\begin{align}\hfill {n}_{{\text{i}}_{n-1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n-\frac{1}{2}}}}{{\lambda }_{n-\frac{1}{2}}}\left({e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}-1\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n-1}^{2}\left({\zeta }^{\prime }{N}_{{\text{e}}_{n+1}}+{\eta }^{\prime }{N}_{{\text{e}}_{n}}\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n-1}\left({\gamma }_{\text{i}}^{\prime }{N}_{{\text{i}}_{n+1}}+{\beta }_{\text{i}}^{\prime }{N}_{{\text{i}}_{n}}\right)\hfill \end{align} \tag{ A.41 }$$

where ${\nu }_{{\text{i}}_{n-\frac{1}{2}}}$ is the averaged ionization frequency between steps n and n − 1, and:

$$\begin{equation}{\zeta }^{\prime }={\nu }_{{\text{i}}_{n-\frac{1}{2}}}\frac{{\Delta}{t}_{n-1}}{{\Delta}{t}_{n}}{Q}_{3}^{\prime }\left({\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}\right)\end{equation} \tag{ A.42 }$$

$$\begin{equation}{\eta }^{\prime }={\nu }_{{\text{i}}_{n-\frac{1}{2}}}{Q}_{2}^{\prime }\left({\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}\right)-{\zeta }^{\prime },\end{equation} \tag{ A.43 }$$

$$\begin{equation}{\gamma }_{\text{i}}^{\prime }=\frac{{\Delta}{t}_{n-1}}{2{\Delta}{t}_{n}},\end{equation} \tag{ A.44 }$$

$$\begin{equation}{\beta }_{\text{i}}^{\prime }=-1-{\gamma }_{\text{i}}^{\prime },\end{equation} \tag{ A.45 }$$

with Q_j ' defined as:

$$\begin{equation}{Q}_{j}^{\prime }\left(x\right)=\frac{{e}^{-x}-{E}_{j}^{\prime }\left(x\right)}{{x}^{j}},\hspace{25.0pt}{E}_{j}^{\prime }\left(x\right)=\sum\limits _{k=0}^{j-1}{\left(-1\right)}^{k}\frac{{x}^{k}}{k!}.\end{equation} \tag{ A.46 }$$

We combine the approximations of ${n}_{{\text{i}}_{n+1}}$, ${n}_{{\text{i}}_{n-1}}$, ${n}_{{\text{e}}_{n+1}}$, and ${n}_{{\text{e}}_{n-1}}$ to remove the explicit dependence on ${N}_{{\text{e}}_{n+1}}$, ${N}_{{\text{e}}_{n}}$ and ${N}_{{\text{i}}_{n}}$. The expressions of ${n}_{{\text{i}}_{n+1}}$, ${n}_{{\text{i}}_{n-1}}$, ${n}_{{\text{e}}_{n+1}}$, and ${n}_{{\text{e}}_{n-1}}$ are:

$$\begin{align}\hfill {n}_{{\text{i}}_{n+1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n}^{2}\left(\zeta {N}_{{\text{e}}_{n+1}}+\eta {N}_{{\text{e}}_{n}}\right)+{\Delta}{t}_{n}\left({\gamma }_{\text{i}}{N}_{{\text{i}}_{n+1}}+{\beta }_{\text{i}}{N}_{{\text{i}}_{n}}\right),\hfill \end{align} \tag{ A.47 }$$

$$\begin{align}\hfill {n}_{{\text{i}}_{n-1}}& ={n}_{{\text{i}}_{n}}+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n-\frac{1}{2}}}}{{\lambda }_{n-\frac{1}{2}}}\left({e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}-1\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n-1}^{2}\left({\zeta }^{\prime }{N}_{{\text{e}}_{n+1}}+{\eta }^{\prime }{N}_{{\text{e}}_{n}}\right)\hfill \\ \hfill & \quad +{\Delta}{t}_{n-1}\left({\gamma }_{\text{i}}^{\prime }{N}_{{\text{i}}_{n+1}}+{\beta }_{\text{i}}^{\prime }{N}_{{\text{i}}_{n}}\right),\hfill \end{align} \tag{ A.48 }$$

$$\begin{equation}{n}_{{\text{e}}_{n+1}}={n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}+\gamma {\Delta}{t}_{n}{N}_{{\text{e}}_{n+1}}+\beta {\Delta}{t}_{n}{N}_{{\text{e}}_{n}}\end{equation} \tag{ A.49 }$$

$$\begin{equation}{n}_{{\text{e}}_{n-1}}={n}_{{\text{e}}_{n}}{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}+{\gamma }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n+1}}+{\beta }^{\prime }{\Delta}{t}_{n-1}{N}_{{\text{e}}_{n}}.\end{equation} \tag{ A.50 }$$

We multiply the above expressions with scalar coefficients, A₁, A₂ and A₃, chosen such that there will not be any explicit dependence on ${N}_{{\text{i}}_{n}},\enspace {N}_{{\text{e}}_{n+1}}$ and ${N}_{{\text{e}}_{n}}$ in the final formulation. We have:

$$\begin{align}\hfill {N}_{{\text{i}}_{n+1}}& =\left[{n}_{{\text{i}}_{n+1}}+{A}_{1}{n}_{{\text{i}}_{n-1}}+{A}_{2}{n}_{{\text{e}}_{n+1}}+{A}_{3}{n}_{{\text{e}}_{n-1}}\right.\hfill \\ \hfill & \quad -\left.\left(1+{A}_{1}\right){n}_{{\text{i}}_{n}}-{A}_{2}{n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}\right.\hfill \\ \hfill & \quad -\frac{{A}_{1}{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n-\frac{1}{2}}}}{{\lambda }_{n-\frac{1}{2}}}\left({e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}-1\right)\hfill \\ \hfill & \quad +\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\hfill \\ \hfill & \quad -\left.{A}_{3}{n}_{{\text{e}}_{n}}{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}\right]{\left({\gamma }_{\text{i}}{\Delta}{t}_{n}+{A}_{1}{\gamma }_{\text{i}}^{\prime }{\Delta}{t}_{n-1}\right)}^{-1}\hfill \end{align} \tag{ A.51 }$$

with:

$$\begin{equation}{A}_{1}=-\frac{{\Delta}{t}_{n}{\beta }_{\text{i}}}{{\Delta}{t}_{n-1}{\beta }_{\text{i}}^{\prime }},\end{equation} \tag{ A.52 }$$

$$\begin{equation}{A}_{2}=\frac{-{\Delta}{t}_{n}^{2}\zeta -{A}_{1}{\Delta}{t}_{n-1}^{2}{\zeta }^{\prime }-{A}_{3}{\gamma }^{\prime }{\Delta}{t}_{n-1}}{\gamma {\Delta}{t}_{n}},\end{equation} \tag{ A.53 }$$

$$\begin{align}\hfill {A}_{3}& =\left(-{\Delta}{t}_{n}^{2}\eta -{A}_{1}{\Delta}{t}_{n-1}^{2}{\eta }^{\prime }+\beta \frac{{\Delta}{t}_{n}^{2}\zeta +{A}_{1}{\Delta}{t}_{n-1}^{2}{\zeta }^{\prime }}{\gamma }\right)\hfill \\ \hfill & \quad {\times}{\left[{\Delta}{t}_{n-1}\left({\beta }_{n}-\frac{{\gamma }^{\prime }}{\gamma }\beta \right)\right]}^{-1}.\hfill \end{align} \tag{ A.54 }$$

The drift–diffusion equation for positive ions (1.2) can now be written as:

$$\begin{equation}{N}_{{\text{i}}_{n+1}}=\nabla \cdot \left({D}_{\text{p}}\enspace \nabla {n}_{\text{p}}\right)+\nabla \cdot \left({\boldsymbol{v}}_{\text{air}}+{\mu }_{\text{p}}\boldsymbol{E}{n}_{\text{p}}\right)-{S}_{{r}_{\text{e}-\text{p}}}-{S}_{{r}_{\text{n}-\text{p}}},\end{equation} \tag{ A.55 }$$

where the ionization source term has been removed from the equation, and the unsteady term, ∂n_p/∂t, is approximated by ${N}_{{\text{i}}_{n+1}}$ (A.51).

In the initial time step there is no information about the previous time step so we have to use the first order method:

$$\begin{align}\hfill {N}_{{\text{i}}_{n+1}}& =\left[{n}_{{\text{i}}_{n+1}}+A{n}_{{\text{e}}_{n+1}}-{n}_{{\text{i}}_{n}}-A{n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}\right.\hfill \\ \hfill & \quad -\left.\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\right]\enspace {\Delta}{t}_{n}^{-1},\hfill \end{align} \tag{ A.56 }$$

where,

$$\begin{equation}A=-\frac{{\Delta}{t}_{n}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}{Q}_{2}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)}{{Q}_{1}\left({\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}\right)}.\end{equation} \tag{ A.57 }$$

The implicit and explicit parts needed by Ansys Fluent are respectively given by:

$$\begin{equation}\mathrm{a}\mathrm{p}\mathrm{u}=-{\left({\gamma }_{\text{i}}{\Delta}{t}_{n}+{A}_{1}{\gamma }_{\text{i}}^{\prime }{\Delta}{t}_{n-1}\right)}^{-1},\end{equation} \tag{ A.58 }$$

$$\begin{align}\hfill \mathrm{s}\mathrm{u}& =-\left[{A}_{1}{n}_{{\text{i}}_{n-1}}+{A}_{2}{n}_{{\text{e}}_{n+1}}+{A}_{3}{n}_{{\text{e}}_{n-1}}-\left(1+{A}_{1}\right){n}_{{\text{i}}_{n}}\right.\hfill \\ \hfill & \quad -\left.{A}_{2}{n}_{{\text{e}}_{n}}{e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}\right.-\frac{{A}_{1}{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n-\frac{1}{2}}}}{{\lambda }_{n-\frac{1}{2}}}\hfill \\ \hfill & \quad {\times}\left({e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}-1\right)+\frac{{n}_{{\text{e}}_{n}}{\nu }_{{\text{i}}_{n+\frac{1}{2}}}}{{\lambda }_{n+\frac{1}{2}}}\left({e}^{{\lambda }_{n+\frac{1}{2}}{\Delta}{t}_{n}}-1\right)\hfill \\ \hfill & \quad -\left.{A}_{3}{n}_{{\text{e}}_{n}}{e}^{-{\lambda }_{n-\frac{1}{2}}{\Delta}{t}_{n-1}}\right]{\left({\gamma }_{\text{i}}{\Delta}{t}_{n}+{A}_{1}{\gamma }_{\text{i}}^{\prime }{\Delta}{t}_{n-1}\right)}^{-1}.\hfill \end{align} \tag{ A.59 }$$

Similar to electrons, the implicit ELP scheme for ions is applied to the regions of the domain where the electric field magnitude is strictly above the breakdown field. In the other regions, the integration is done by a classical 2nd order backward difference scheme.