# An implicit time integration approach for simulation of corona discharges

Nguyen Tuan Dung, Christophe Besse, François Rogier

► To cite this version:

Nguyen Tuan Dung, Christophe Besse, François Rogier. An implicit time integration approach for simulation of corona discharges. Computer Physics Communications, 2023, 294, pp.108906. 10.1016/j.cpc.2023.108906. hal-04102798

HAL Id: hal-04102798

https://hal.science/hal-04102798

Submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An implicit time integration approach for simulation of corona discharges

Christophe Besse<sup>a</sup>, Nguyen Tuan Dung<sup>a,b</sup>, François Rogier<sup>b</sup>

<sup>a</sup>Institut de Mathématiques de Toulouse, UMR5219, Université de Toulouse 3, CNRS, Toulouse, 31400, France

<sup>b</sup>ONERA-The French Aerospace Lab, Toulouse, 31400, France

## Abstract

The simulation of corona discharges usually requires an excessive amount of time. The objective of this work is therefore to reduce the CPU time of corona discharges. The profound disparity between physical time scales in a gas discharge renders explicit time integration impractical in large time scale applications, for example in flow controls. We propose an implicit time strategy in the scope of the Local Field Approximation plasma model. We investigate some algorithms to solve the discretized plasma model, which are adapted to take into consideration a constraint on the electron density. The constraint is introduced to compensate for unaccounted sources of electrons; it is expressed by a minimum/floor density imposed on the electrons. Therefore, simplified plasma models can be employed for discharge simulations, for instance a four-specie, four-reaction model is used in this article. The investigation is done on a one-dimensional wire-to-wire corona discharge. After having identified the most suitable algorithm, we carry on to two-dimensional wire-to-wire discharges and compare the numerical results with experiment data for different floor densities.

**Keywords:** implicit method, corona discharge, LFA model, differential inclusion

## 1. Introduction

In the past two decades plasma actuators have been extensively studied for their applications in airflow control [1, 2], anti-icing [3] and flight sustaining [4], among many others. These devices are robust, light, low energy
consuming and easy to install, thus they have a potential interest for in-
dustrial use. A typical actuator usually consists of two or more electrodes
and/or a dielectric barrier that are capable to sustain a stable discharge. The
newborn ions from the ionization process are drifted along the field and col-
liding with neutral particles in air, creating an electro-hydrodynamic (EHD)
force. A lot of experimental efforts have been made to characterize the EHD
force [5, 6, 1, 7, 8, 9, 10, 11, 12, 13], to optimize the actuator geometry
[14, 15, 16, 17, 18] or to measure the dielectric surface charge [19, 20, 21].
Also during these years, a large numerical research has been conducted to
get more insights from a plasma discharge such as distribution density of
particles, how they interact with each other and with the potential field, how
they accumulate on the dielectric surface or how they dictate the onset of
different discharge phases [22, 23, 24, 25, 26, 27, 28, 29, 30].

For simulation of corona discharges, many fluid models have been pro-
posed throughout the years. Morrow [31] studied the discharge for concentric-
sphere electrodes with a set of evolution continuity equations for five species
(electrons, positive ions, negative ions, neutrals) coupled with a Poisson equa-
tion for the electrostatic potential. This approach took into account also the
electron source from photoionization but was not validated experimentally.
Adamiak & Atten [32] proposed a simplified, stationary model consisting of
the Poisson equation and a nonlinear equation for the space charge, at the
expense of capturing the discharge dynamics. This approached was validated
experimentally for different actuators such as point-to-plane [32], needle-to-
plane [33], wire-cylinder-plate [34] or needle-to-ring [35]. A similar model, but
evolutionary, was proposed by Guan et al. [36] but the boundary condition
for the charge equation depends on the experimental electric current. Chen
et al. [37] used a model similar to Morrow’s with three species (electrons,
positive ions, negative ions) to study negative needle-to-cylinder discharges
and validated experimentally. Photoionization was however neglected since
secondary emission from the cathode surface has more contribution for neg-
ative corona discharges [38].

Simulation of corona discharges proves to be very challenging in terms of
waiting time [31] since there is a disparity among characteristic time scales
of several phenomena within a discharge. For example, in flow control, the
time scale of air flow motion is on milliseconds, but the time scale of elec-
tron motion is only on picoseconds. Many numerical efforts have been made
to reduce the simulation time. Ventzek et al. [39] proposed a semi-implicit
scheme which allowed to compute the field implicitly and separately from the
specie densities. Hagelaar [40] proposed a scheme in the same category in the
presence of the electron energy equation. Seimandi et al. [24] studied some
asymptotic models to reduce the problem dimension in the ionization lay-
ers near the electrodes. Dufour & Rogier developed a sub-cycling algorithm
in their multiscale plasma solver COPAIER [41]. The idea was to repeat-
edly “freeze” the field and ions motion and update the electron density in
the meantime. There are other numerical methods developed for studying
streamer propagation, such as fully coupled implicit schemes [42, 43, 44] or
asynchronous time stepping [45], but their application on corona discharges
does not seem to be obvious.

The description of all physicochemical processes in a non-equilibrium plasma discharge is very complex. In air, the number of plasma-chemical reactions can be as high as around 300 [46]. The integration of all processes in the numerical model would be impractical for large time scale applications. Furthermore, some processes, such as photoionization, require a huge computation effort. The effects of photoionization are very important to discharge regimes such as streamer propagation because they are structurally non-uniform and extremely sensitive to simulation conditions [47, 48], but not for the glow regime in which we are interested [38] in this article. For these reasons, we aim to drop most of the plasma-chemical processes, including photoionization for the sake of simulation time. We enforce instead a *minimum/floor density* on the electron density to compensate for unaccounted sources of electrons. This approach was implemented in ONERA's in-house plasma solver COPAIER [41] with explicit time integration. Simulations with floor density have been contributing to the study of plasma actuators [49, 50, 51, 52]. Recently, simulation results from [53] with a floor density around $10^{11}-10^{12}m^{-3}$ showed good agreement with experiment data for a wire-to-wing corona discharge.

The structure of this work is presented as follows. In section 2, we clarify
the integration of the floor densities to the plasma model. The enforcement
of a minimum density $\psi$ on the charged particle density $n$ can be seen as
a constraint: $n \ge \psi$, which is not automatically satisfied by the solution of
the particle continuity equation, but necessitates the rewriting of the charged
specie conservation law under the form of a differential inclusion. Section 3
introduces the implicit time discretization of the plasma model and investi-
gates some algorithms to solve the differential inclusions numerically. Since
the problem is nonlinear, direct methods are not suitable to use and hence
splitting methods such as Lie operator splitting or iterative methods such
as Douglas-Rachford or Gauss-Seidel algorithms are employed. However, we shall see in section 4 that the Lie splitting is not a good candidate. In fact, it does not conserve the steady state solution of a DC corona discharge since the numerical electric current depends on the chosen time step. On the other hand, the studied iterative methods do conserve the stationary solution and therefore they should be preferred for simulations of corona discharges or even of the corona phase in DBD discharges. Finally, section 5 discusses the capability of the proposed model to reproduce the electric currents measured in real physics experiments [7]. The studied case is a corona discharge generated by a wire-to-wire actuator and the simulations are conducted in two dimensions with COPAIER.

## 2. Modeling of gas discharges with floor density

### 2.1. Simplified fluid model

We use the Local Field Approximation (LFA) model [54] to describe the plasma dynamics. It is composed of a set of conservation laws coupled with a Poisson equation. Let us assume that the discharge is contained in an open bounded domain $\Omega \subset \mathbb{R}^d$ ($d = 1, 2, 3$) (see fig. 1) and evolves on a time interval $(0, T)$. Let us denote as $\mathcal{L}$ the set of particle species that are considered in the discharge and as $n_i(t, \mathbf{x})$ (measured in $m^{-3}$) the particle density of the species $i \in \mathcal{L}$. The densities are grouped in a vector denoted as $\mathbf{U} = (n_i)_{i \in \mathcal{L}}$. The electrostatic field is denoted as $\mathbf{E}(t, \mathbf{x})$ ($Vm^{-1}$).

The generic form of the continuity equations for the species reads as follows,

$$ \partial_t \mathbf{U} + \nabla \cdot \mathbf{F}(\mathbf{E}, \mathbf{U}) = \mathbf{S}(E/N, \mathbf{U}), \quad (t, \mathbf{x}) \in (0, T) \times \Omega, \quad (1) $$

where $\mathbf{F} = (\mathbf{f}_i)_{i \in \mathcal{L}}$, $\mathbf{S} = (S_i)_{i \in \mathcal{L}}$, $E = |\mathbf{E}|$ and $N$ ($m^{-3}$) is the air density. More precisely, $\mathbf{f}_i(\mathbf{E}, \mathbf{U})$ is the particle flux density of species $i$ and $S_i(E/N, \mathbf{U})$ is the creation/decay rate per unit volume of $i$ via plasma-chemical processes. For $\mathbf{x} \in \Omega$, $\mathbf{f}_i = \mathbf{u}_i n_i - D_i \nabla n_i$ where $\mathbf{u}_i(\mathbf{E})$ ($ms^{-1}$) is the drift velocity of $i$ and $D_i(E/N)$ ($m^2s^{-1}$) is its diffusion rate. The charged species are subject to a mobility law which reads as

$$ \mathbf{u}_i = \mathrm{sign}(z_i) \mu_i \mathbf{E}, $$

where $z_i$ is the charge number of $i$ and $\mu_i(E/N)$ ($m^2V^{-1}s^{-1}$) is its mobility coefficient. The field is derived from the electric potential $\phi$ as

$$ \mathbf{E} = -\nabla \phi. $$
The electric potential $\phi(t, \boldsymbol{x})$ ($V$) is given by Gauss' law,

$$ -\nabla \cdot (\varepsilon \nabla \phi) = \rho, \qquad (2) $$

where $\varepsilon(\boldsymbol{x})$ ($CV^{-1}m^{-1}$) is the electric permittivity of the domain $\Omega$, $\rho = q \sum_{i \in \mathcal{L}} z_i n_i$ ($Cm^{-3}$) is the space charge density and $q \approx 1.6 \times 10^{-19}C$ is the elementary charge. We assume for simplicity that $n_i$, $\phi$, $\boldsymbol{E}$, $S_i$ are smooth functions defined on $(0, T) \times \Omega$.

<figure><img src="image_2.png" /><figcaption>Figure 1: A two-dimensional open domain $\Omega$ with boundary $\Gamma = \bar{\Omega} \setminus \Omega$. The coordinates of an arbitrary point is denoted as $\boldsymbol{x}$.</figcaption></figure>

The number of species and physicochemical reactions is generally very numerous [46, 55]. In order to limit the complexity of the system, we adopt a simplified model from J.-P. Bœuf et al.'s work [27] where the composition of the ionized gas is lumped into four species: electrons (denoted as $e$), positive ions ($p$), negative ions ($n$)$^1$ and neutral particles ($N$). The reactions consist only of Townsend ionization $e+N \rightarrow 2e+p$, electron attachment $e+N \rightarrow n$, electron-ion recombination $e+p \rightarrow N$ and ion-ion recombination $p+n \rightarrow N$. In the following, we suppose that the neutral particle density, still denoted as $N$, is constant and $N = 2.5 \times 10^{25}m^{-3}$ at atmospheric pressure. The plasma dynamics is thus decoupled from airflow dynamics, the discharge model simply considers three continuity equations for $n_e$, $n_p$, $n_n$; therefore, $\mathcal{L} = \{e, p, n\}$. The specie source terms are written as,

$$ \begin{aligned} S_e &= (\alpha - \eta) N n_e - k_{ep} n_e n_p, \\ S_p &= \alpha N n_e - k_{ep} n_e n_p - k_{np} n_n n_p, \\ S_n &= \eta N n_e - k_{np} n_n n_p, \end{aligned} $$

where $\alpha$, $\eta$, $k_{ep}$ and $k_{np}$ ($m^3s^{-1}$) are resp. the ionization, electron attachment, electron-ion recombination and ion-ion recombination coefficients, which are

<footnote>$^1z_e = -1, z_p = +1, z_n = -1$</footnote>
positive. In the LFA model, these coefficients depend only on the reduced field strength $E/N$. In the simulations, $\alpha$, $\eta$ are extracted from look-up tables provided by the BOLSIG+ application [56]; $k_{ep} = 2 \times 10^{-13}$ and $k_{np} = 1.7 \times 10^{-13}$.

As discussed in section 1, we introduce a positive and smooth function on $\Omega$,

$$ \psi(\boldsymbol{x}) \geq 0, $$

which we dub floor density, to represent the smallest electron population that exists because of unaccounted plasma-chemical processes. We enforce the following constraint on the electron density,

$$ n_e(t, \boldsymbol{x}) \geq \psi(\boldsymbol{x}), \quad \forall (t, \boldsymbol{x}) \in (0, T) \times \Omega. \qquad (3) $$

The problem as a whole is ill-posed since the solution of eq. (1) (for $i = e$) may not respect eq. (3). Therefore, a new reformulation of the problem is necessary.

Let us denote $\mathcal{K}$ as the closed convex set of smooth functions defined on $\Omega$ that satisfy eq. (3), i.e.

$$ \mathcal{K} = \{ v | v \geq \psi \}^2. $$

The characteristic function of $\mathcal{K}$ is defined as

$$ I_{\mathcal{K}}(v) = \begin{cases} 0 & \text{if } v \in \mathcal{K}, \\ +\infty & \text{if } v \notin \mathcal{K}. \end{cases} $$

The subdifferential [57, Proposition 3] of $I_{\mathcal{K}}$ at $v$ is a set of functions defined as

$$ \partial I_{\mathcal{K}}(v) = \begin{cases} \emptyset & \text{if } v \notin \mathcal{K}, \\ \{0\} & \text{if } v > \psi, \\ \{u | \forall w \in \mathcal{K}, (u, w - v) \leq 0\} & \text{otherwise,} \end{cases} \qquad (4) $$

where $(\cdot, \cdot)$ is a scalar product defined on a Hilbert space $H$ such that for every $t \in (0, T)$, $n_i(t, \cdot) \in H$, $i \in \mathcal{L}$. For example, $H = L^2(\Omega)$ and $(u, v) = \int_{\Omega} u(\boldsymbol{x})v(\boldsymbol{x})d\boldsymbol{x}$.

<footnote n="2">²without ambiguity, all inequalities are point-wise, i.e. $v(\boldsymbol{x}) \geq \psi(\boldsymbol{x}) \quad \forall \boldsymbol{x} \in \Omega$</footnote>
As we have discussed, the constraint eq. (3) may not be compatible with eq. (1), but if we assume for now the contrary, i.e. $n_e(t, \cdot) \in K \ \forall t \in (0, T)$, then we would have from eq. (1) that

$$S_e - \partial_t n_e - \nabla \cdot \mathbf{f}_e + n_e = n_e \ge \psi, \quad \forall (t, \mathbf{x}) \in (0, T) \times \Omega,$$

Having these relations suggests that the projection on $K$ of the left-hand side is $n_e$. On the other hand, we have the following functional analysis theorem [58].

**Theorem 1** (Projection on a closed convex). Let $K \subset H$ be a non empty closed convex. Then for all $v \in H$, there exists a unique $u \in K$ (its projection) such that,

$$|v - u| = \min_{w \in K} |v - w|,$$

where $|\cdot|$ is the norm induced from $(\cdot, \cdot)$. Furthermore, $u$ is characterized by the property

$$(v - u, w - u) \le 0, \quad \forall w \in K.$$

Thus, we have $(S_e - \partial_t n_e - \nabla \cdot \mathbf{f}_e + n_e - n_e, w - n_e) \le 0$ for all $w \in K$ and $t \in (0, T)$. This and eq. (4) imply that $S_e(t, \cdot) - \partial_t n_e(t, \cdot) - \nabla \cdot \mathbf{f}_e(t, \cdot)$ is an element of the set $\partial I_K(n_e(t, \cdot))$ for all $t \in (0, T)$. Let us denote $S_\psi(t) = -\partial I_K(n_e(t, \cdot))$. The electron conservation law which integrates the floor density $\psi$ now transforms into a differential inclusion [59] which reads as

$$\partial_t n_e(t, \cdot) + \nabla \cdot \mathbf{f}_e(t, \cdot) \in S_e(t, \cdot) + S_\psi(t)^3, \quad \forall t \in (0, T). \quad (5)$$

Finally, the simplified LFA model reads as follows,

$$\left\{
\begin{array}{l}
\partial_t \mathbf{U} + \nabla \cdot \mathbf{F} \in \mathbf{S} + \begin{pmatrix} S_\psi \\ 0 \\ 0 \end{pmatrix}, \\
\mathbf{U}(0, \cdot) = \mathbf{U}^0, \\
\nabla \cdot (\varepsilon \mathbf{E}) = \rho,
\end{array}
\right.
\qquad (6)$$

where $\mathbf{U}^0 = (n_e^0 \ n_p^0 \ n_n^0)^t$ and $n_i^0$ are smooth functions defined on $\Omega$.

We remark that the above system is not charge conservative, i.e. $\partial_t \rho + \nabla \cdot \mathbf{j} \ne 0$. However, since the floor density $\psi$ that we choose in the simulations

<footnote>$^3S_e(t, \cdot) + S_\psi(t) \equiv \{ S_e(t, \cdot) + v \mid v \in S_\psi(t) \}$</footnote>

<footnote>$^4j = q(\mathbf{f}_p - \mathbf{f}_n - \mathbf{f}_e)$ is the conduction current density</footnote>
is much less than the charge density in corona discharges (see section 5.1), the charge conservation issue does not have much influence on the studied discharges.

Another remark is that the sign of $S_e$ and $\partial_t n_e + \nabla \cdot \mathbf{f}_e$ in eq. (5) is not arbitrarily chosen (we could have otherwise $S_e \in \partial_t n_e + \nabla \cdot \mathbf{f}_e + S_\psi$). The existence and uniqueness of solution of the nonlinear problem eq. (5) are subject of classical studies based on the theory of monotone operators [60] and variational inequalities [61, 62, 63].

## 2.2. Boundary conditions

The boundary conditions for the conservation laws are imposed on the fluxes $\mathbf{f}_i$. Different boundary conditions exist and each is defined on a non-overlapping boundary portion $\Gamma_k$. The union of these portions is the domain boundary $\Gamma$: $\Gamma = \cup_k \Gamma_k$ (see fig. 10 for example).

1. For free-flow boundaries $\Gamma_f$, which are fictive interfaces between the discharge domain and the outer air space, we assume that they are far enough from the discharge so that the potential and density gradients are null, i.e.

$$ (\mathbf{E} \cdot \boldsymbol{\nu})|_{\Gamma_f} = 0, \qquad (7) \\ f_{i|\Gamma_f} \equiv (\mathbf{f}_i \cdot \boldsymbol{\nu})|_{\Gamma_f} = (n_i \mathbf{u}_i \cdot \boldsymbol{\nu})|_{\Gamma_f} = 0, \quad i \in \mathcal{L}, $$

where $\boldsymbol{\nu}(\mathbf{x})$ is the unit outward normal of $\Omega$ at $\mathbf{x} \in \Gamma$.

2. For wall boundaries $\Gamma_w$, which are physical interfaces between the discharge domain and the electrodes<fnref n="5" />, we determine the species fluxes using the two-stream approach [64]. We assume that there is no particle reflection on the walls and only electron secondary emission due to ion bombardment is involved. The boundary conditions read as

$$ \phi|_{\Gamma_w} = \phi_0, \\ f_{i|\Gamma_w} = \left( n_i \max(\mathbf{u}_i \cdot \boldsymbol{\nu}, 0) \right)|_{\Gamma_w}, \quad i \in \mathcal{L} \setminus \{e\}, \qquad (8) \\ f_e|_{\Gamma_w} = \frac{1}{2} \left( \frac{8k_B T_e}{\pi m_e} \right)^{\frac{1}{2}} n_e|_{\Gamma_w} - 2\gamma f_{p|\Gamma_w}. $$

<footnote n="5">⁵we do not consider interfaces with dielectric material: they do not appear in the discharges studied in this article</footnote>
where $\phi_0$ is the electrode voltage, $k_B$ is the Boltzmann constant, $T_e$ is the electron temperature, $m_e$ is the electron mass and $\gamma$ is the secondary emission coefficient. We set $T_e = 2eV$ and $\gamma = 10^{-4}$.

3. For symmetric boundaries $\Gamma_s$, which are fictive interfaces between symmetric discharge domains, the boundary conditions are

$$ (\mathbf{E} \cdot \boldsymbol{\nu})|_{\Gamma_s} = 0, \\ f_{i|\Gamma_s} = 0, \quad i \in \mathcal{L}. \qquad (9) $$

### 2.3. Natural sources of charged species in normal condition

We assume that all the discharges in this paper are simulated in open space air. The level of pre-ionization charged species is dictated by natural radioactive decays, especially those of radon [65]. The appearance rate of an electron-positive ion pair is up to $10^9 m^{-3} s^{-1}$ which yields an equilibrium density of $10^9 - 10^{10} m^{-3}$. The floor density $\psi$ can also be envisioned to include the naturally formed electrons.

## 3. Numerical methods

### 3.1. Key features

A complete description of the COPAIER solver is found in [41], so we only highlight here some useful elements.

Numerical grids are locally refined near the electrodes where the ionization process is most active. This proves to be vital to the simulations as pointed out in section 4, where the electric currents do not exhibit artifact pulses like those in [23] for instance. The discharges studied in this paper do not develop streamers so only fixed grids are sufficient. The grids used in section 5 are triangular and generated by the Gmsh application [66].

The cell-centered finite volume approach is employed to discretized the continuity equations/inclusions. The Scharfetter-Gummel scheme [67] as well as the SGCC scheme [68] are used to flux approximation, although the latter has not yet been implemented in COPAIER. The P1-Lagrange finite element method is used to solve the Poisson equation.

Explicit time integration is customary in COPAIER. The numerical time step suffers some upper limitations. For example, on an one-dimensional uniform grid with step size $\Delta x$ and at the time level $t^l$, the time step $\Delta t^l \equiv$
$t^{l+1} - t^l$ must satisfy<fnref n="6" />

$$ \Delta t^l \le \min_x \left( \frac{\Delta x}{|\mathbf{u}_e(\mathbf{E}(t^l, \mathbf{x}))|} \tanh \frac{|\mathbf{u}_e(\mathbf{E}(t^l, \mathbf{x}))|\Delta x}{2D_e(\mathbf{E}(t^l, \mathbf{x}))} \right) \equiv \Delta t_e^l, $$

$$ \Delta t^l \le \min_{i, \mathbf{x}} \left( \frac{\Delta x}{|\mathbf{u}_i|} \tanh \frac{|\mathbf{u}_i|\Delta x}{2D_i} \right) \equiv \Delta t_{\text{ions}}^l, \quad i \in \mathcal{L} \setminus \{e\}, $$

$$ \Delta t^l \le \min_x \left( \frac{\varepsilon(\mathbf{x})}{\sigma(t^l, \mathbf{x})} \right) \equiv \Delta t_{\phi}^l. $$

In corona discharges, the dielectric relaxation time $\Delta t_{\phi}^l$ is much larger than the CFL condition $\Delta t_e^l$ (a factor $10^4 - 10^7$). The latter is around $10^{-13} - 10^{-12}$s comparing to the target physical time $10^{-3}$s. This huge discrepancy results in unbearable waiting time which can amount to weeks or even months in two-dimensional simulations, even on multi MPI processes. Therefore we could treat the conservation laws eq. (5) implicitly and then update the potential by eq. (2) after the charge density is up-to-date.

## 3.2. Discrete plasma model

For $l \ge 0$, let $\mathbf{U}^l(\mathbf{x})$ be an approximation of $\mathbf{U}(t^l, \mathbf{x})$ at discrete time $t^l$ (and similarly for $\phi$, $\mathbf{u}_i$, $D_i$, $\alpha$, $\eta$, $k_{ep}$, $k_{np}$) and $\partial_t \mathbf{U}(t^{l+1})$ be approximated by $\frac{\mathbf{U}^{l+1} - \mathbf{U}^l}{\Delta t^l}$, we consider the following time discretization of eq. (6),

$$ \left\{ \begin{array}{l} \frac{\mathbf{U}^{l+1} - \mathbf{U}^l}{\Delta t^l} + \nabla \cdot \mathbf{F}(\mathbf{E}^l, \mathbf{U}^{l+1}) \in \mathbf{S}(E^l/N, \tilde{\mathbf{U}}) + \begin{pmatrix} S_{\psi}^{l+1} \\ 0 \\ 0 \end{pmatrix}, \quad (10) \\ \\ \nabla \cdot (\varepsilon \nabla \mathbf{E}^{l+1}) = \rho^{l+1}, \end{array} \right. $$

with

$$ \mathbf{S}(E^l/N, \tilde{\mathbf{U}}) = \begin{pmatrix} (\alpha^l - \eta^l) N n_e^{l+1} - k_{ep}^l n_e^{l+1} n_p^l + S_{\psi}^{l+1} \\ \alpha^l N n_e^{l+1} - k_{ep}^l n_e^{l+1} n_p^l - k_{np}^l n_n^l n_p^l \\ \eta^l N n_e^{l+1} - k_{np}^l n_n^l n_p^l \end{pmatrix} $$

<footnote n="6">⁶$\sigma(t, \mathbf{x}) \equiv q \sum_{i \in \mathcal{L}} |z_i| \mu_i(E(t, \mathbf{x})) n_i(t, \mathbf{x})$ is the electrical conductivity</footnote>
and $S_{\psi}^l \equiv S_{\psi}(t^l) = -\partial_K(n_e(t^l, \cdot))$. We remark that the coefficients $u_e$, $D_i$, $\alpha$, $\eta$, $k_{ep}$, $k_{np}$ are computed from the potential $\phi^l$, so the Poisson equation is decoupled from the specie conservation laws. Also, in the recombination terms, only the electron density is taken at $t^{l+1}$. This choice was made to reduce the non-linearity of the system and because the recombination rates are much smaller than the ionization rate in corona discharges. The proposed discrete model is quite similar to the one proposed in [69] except the non-linearity only stemming from $S_{\psi}$.

Let us denote as $\mathcal{I}$ the identity operator,

$$ \mathcal{A}^l : U \mapsto \begin{pmatrix} \nabla \cdot \mathbf{f}_e(\mathbf{E}^l, U) - (\alpha^l - \eta^l)Nn_e + k_{ep}^l n_e n_p^l \\ \nabla \cdot \mathbf{f}_p(\mathbf{E}^l, U) - \alpha^l Nn_e + k_{ep}^l n_e n_p^l \\ \nabla \cdot \mathbf{f}_n(\mathbf{E}^l, U) - \eta^l Nn_e + k_{np}^l n_n^l n_p^l \end{pmatrix} \equiv \begin{pmatrix} (\mathcal{A}^l U)_e \\ (\mathcal{A}^l U)_p \\ (\mathcal{A}^l U)_n \end{pmatrix}, $$

$\mathcal{B}_{\psi} : U \mapsto \begin{pmatrix} -S_{\psi} \\ 0 \\ 0 \end{pmatrix}$ and $\mathbf{R}^l = \begin{pmatrix} n_e^l \\ n_p^l - \Delta t^l k_{np}^l n_n^l n_p^l \\ n_n^l - \Delta t^l k_{np}^l n_n^l n_p^l \end{pmatrix} \equiv \begin{pmatrix} \mathbf{R}_e^l \\ \mathbf{R}_p^l \\ \mathbf{R}_n^l \end{pmatrix}$. Then at $t^l$, the semi-discrete (in time) conservation laws of eq. (10) can be rewritten as

$$ (\mathcal{I} + \Delta t^l \mathcal{A}^l + \Delta t^l \mathcal{B}_{\psi}) \mathbf{U}^{l+1} \ni \mathbf{R}^l. \qquad (11) $$

<figure><img src="image_3.png" /><figcaption>Figure 2: A triangulation $\mathcal{T}$ of a polygonal domain $\Omega$ in $\mathbb{R}^2$</figcaption></figure>

Furthermore, let us consider a polygonal domain $\Omega$ in $\mathbb{R}^2$ and a mesh $\mathcal{T}$ consisting of closed polygonal convex subsets of $\Omega$ (see fig. 2). For an element
$K \in \mathcal{T}$, we denote as $\mathcal{E}_K$ the set of its edges. Let us consider the electron conservation law of eq. (10),

$$ \frac{n_e^{l+1} - n_e^l}{\Delta t^l} + \nabla \cdot \mathbf{f}_e(\mathbf{E}^l, \mathbf{U}^{l+1}) - (\alpha^l - \eta^l) N n_e^{l+1} + k_{ep}^l n_e^{l+1} n_p^l = h, $$

where $h$ is an element of $S_{\psi}^{l+1}$. Now integrating this equation on $K$ and dividing it by $|K|$, we have

$$ \frac{\bar{n}_{e,K}^{l+1} - \bar{n}_{e,K}^l}{\Delta t^l} + \frac{1}{|K|} \sum_{\lambda \in \mathcal{E}_K} \int_{\lambda} \mathbf{f}_e \cdot \boldsymbol{\nu}_K \mathrm{d}\boldsymbol{l} - \left( \overline{(\alpha^l - \eta^l) N n_e^{l+1}} \right)_K + \left( \overline{k_{ep}^l n_e^{l+1} n_p^l} \right)_K = \bar{h}_K, $$

where $\bar{v}_K \equiv \frac{1}{|K|} \int_K v(\boldsymbol{x}) \mathrm{d}\boldsymbol{x}$ for any function $v$, $\mathrm{d}\boldsymbol{l}$ is the unit length element and $\boldsymbol{\nu}_K(\boldsymbol{x})$ is the unit outward normal of $K$ at $\boldsymbol{x} \in \lambda \in \mathcal{E}_K$. We approximate $\alpha^l(\boldsymbol{x})$ for $\boldsymbol{x} \in K$ (and similarly for $\eta$, $k_{ep}$) by $\alpha^l(\boldsymbol{x}_K) \equiv \alpha_K^l$ with $\boldsymbol{x}_K$ the center of $K$, $\int_{\lambda} \mathbf{f}_e \cdot \boldsymbol{\nu}_K \mathrm{d}\boldsymbol{l}$ by $|\lambda| \mathbf{f}_{e,K,\lambda} \cdot \boldsymbol{\nu}_{K,\lambda} \equiv |\lambda| \mathbf{f}_e(\boldsymbol{x}_{K,\lambda}) \cdot \boldsymbol{\nu}_K(\boldsymbol{x}_{K,\lambda})$ with $\boldsymbol{x}_{K,\lambda}$ the center of $\lambda \in \mathcal{E}_K$ and then $\bar{n}_e \bar{n}_p$ by $\bar{n}_e \bar{n}_p$. Finally we have

$$ \frac{\bar{n}_{e,K}^{l+1} - \bar{n}_{e,K}^l}{\Delta t^l} + \frac{1}{|K|} \sum_{\lambda \in \mathcal{E}_K} |\lambda| \mathbf{f}_{e,K,\lambda} \cdot \boldsymbol{\nu}_{K,\lambda} - (\alpha_K^l - \eta_K^l) N \bar{n}_{e,K}^{l+1} + k_{ep,K}^l \bar{n}_{e,K}^{l+1} \bar{n}_{p,K}^l \approx \bar{h}_K, \quad (12) $$

for every $K \in \mathcal{T}$. Now for any smooth function $v$ defined on $\Omega$, we define a piecewise constant function $\bar{v}$ on $\Omega$ such that

$$ \bar{v}(\boldsymbol{x}) = \bar{v}_K, \quad \forall \boldsymbol{x} \in K \setminus \cup_{\lambda \in \mathcal{E}_K} \lambda. $$

It is not evident how $\bar{h}$ with $h \in S_{\psi}^{l+1}$ relates to $\psi$ and $\bar{n}_e$. Therefore, we make an assumption that

$$ \bar{h} \in S_{\bar{\psi}}(\bar{n}_e^l) \equiv -\partial I_{\bar{K}}(\bar{n}_e^l), \quad \bar{K} = \{ v \mid v \ge \bar{\psi} \}. \quad (13) $$

Finally, using the Scharfetter-Gummel scheme<fnref n="7" /> [67] to approximate the fluxes $\mathbf{f}_{i,K,\lambda} \cdot \boldsymbol{\nu}_{K,\lambda}$, $i \in \mathcal{L}$ yields a linear relation of $(\bar{n}_{i,K})_{K \in \mathcal{T}}$. The left-hand

<footnote n="7">⁷see Appendix A</footnote>
side of eq. (12), combining with similar equations for other species, can be written in the same way of eq. (11) as follows,

$$ (\bar{\mathcal{I}} + \Delta t^l \bar{\mathcal{A}}^l + \Delta t^l \bar{\mathcal{B}}_{\bar{\psi}}) \bar{\mathcal{U}}^{l+1} \ni \bar{\mathcal{R}}^l, \quad (14) $$

where $\bar{\mathcal{U}}^l = (\bar{n}_{e,K}^l \ \bar{n}_{p,K}^l \ \bar{n}_{n,K}^l)_{K \in \mathcal{T}}^t \in \mathbb{R}^{|\mathcal{L}||\mathcal{T}|}$, $\bar{\mathcal{I}}$ is the identity matrix, $\bar{\mathcal{A}}^l$ is a $|\mathcal{L}||\mathcal{T}| \times |\mathcal{L}||\mathcal{T}|$ matrix which is the discrete counterpart of $\mathcal{A}^l$ and $\bar{\mathcal{R}}^l$ is the discrete counterpart of $\mathcal{R}^l$. We remark that writing $\bar{\mathcal{B}}_{\bar{\psi}}\bar{\mathcal{U}}^{l+1}$ is a notation abuse since $\bar{\mathcal{U}}^{l+1}$ is a discrete vector; eq. (14) should be understood in the sense of eq. (12).

We shall investigate in the next section some algorithms to solve the (semi-)discrete conservation laws eq. (11)/eq. (14).

### 3.3. Solving the conservation laws

The complex structure of the operator $\bar{\mathcal{I}} + \Delta t^l \bar{\mathcal{A}}^l + \Delta t^l \bar{\mathcal{B}}_{\bar{\psi}}$ makes the evaluation of $(\bar{\mathcal{I}} + \Delta tA + \Delta tB)^{-1} \bar{\mathcal{R}}^l$ hard to achieve. Therefore, we explore some splitting and iterative methods to solve this problem, namely the Lie operator splitting, the Douglas-Rachford method [70] and a Gauss-Seidel-inspired algorithm. The first one is a classical and very simple method to implement. The second one is also classical and well known in convex optimization [71]. The last method is derived from the Gauss-Seidel algorithm which is used for solving linear systems. We shall compare these methods further in section 4.

Let us remark first of all that the operator $\partial I_K$ is maximal monotone on the Hilbert $H$. We recall that a (multivalued) operator $A$ is monotone if

$$ (y - u, w - v) \ge 0, \quad \forall v, w \in \mathcal{D}(A) \subset H, u \in A(v), y \in A(w), $$

where $\mathcal{D}(A)$ is the domain of $A$; and is maximal monotone if its resolvent $J_A^{\Delta t} = (\bar{\mathcal{I}} + \Delta tA)^{-1}$ is a contraction defined on $H$ [70]. In our case, if $(\bar{\mathcal{I}} + \Delta t^l \partial I_K)w \ni y$ then from the definition eq. (4) we have,

$$ \frac{y-w}{\Delta t^l} \in \partial I_K(n) \iff \forall v \in \mathcal{K}, \quad (y-w, v-w) \le 0. $$

It is shown that $w$ is the projection of $y$ on $\mathcal{K}$, i.e. $w(\boldsymbol{x}) = \max(y(\boldsymbol{x}), \psi(\boldsymbol{x}))$, which is unique according to theorem 1 and easily shown to be a contraction. Therefore, if $\boldsymbol{Y} = (y_e \ y_p \ y_n)^t$ then

$$ J_{\bar{\psi}}^{\Delta t^l} \boldsymbol{Y} = (\bar{\mathcal{I}} + \Delta t^l \bar{\mathcal{B}}_{\bar{\psi}})^{-1} \boldsymbol{Y} = (\max(y_e, \psi) \ y_p \ y_n)^t. \quad (15) $$
Hence, from eq. (11) the electron density at $t^{l+1}$ is $n_e^{l+1} = \max(\bar{R}_e^l - \Delta t^l (\mathcal{A}^l \bar{U}^{l+1})_e, \psi)$ which implies that $\bar{n}_e^{l+1} = \max(\bar{R}_e^l - \Delta t^l (\mathcal{A}^l \bar{U}^{l+1})_e, \bar{\psi})$. By this point, making the assumption eq. (13) is in fact to approximate the right-hand side with $\max(\bar{R}_e^l - \Delta t^l (\mathcal{A}^l \bar{U}^{l+1})_e, \bar{\psi})$. Further approximations of flux and chemical coefficients lead finally to $\bar{n}_e^{l+1} = \max(\bar{R}_e^l - \Delta t^l (\bar{\mathcal{A}}^l \bar{U}^{l+1})_e, \bar{\psi})$, which is the electron density solution of eq. (14).

Another remark is that $J_{\bar{B}_{\psi}}^{\Delta t^l}$ is independent of $\Delta t^l$, which means that the electrons are created instantaneously to reach the floor density. Therefore, we can omit the time step and write $J_{\bar{B}_{\psi}}$ instead of $J_{\bar{B}_{\psi}}^{\Delta t^l}$.

### 3.3.1. Lie splitting

The simplest operator splitting method is the Lie algorithm which is first-order in time. At $t^l$, we compute the numerical solution, still denoted as $\bar{U}^{l+1}$, by

$$ \bar{U}^{l+1} = J_{\bar{B}_{\psi}} J_{\bar{\mathcal{A}}^l}^{\Delta t^l} \bar{R}^l. \qquad (16) $$

The evaluation of $J_{\bar{\mathcal{A}}^l}^{\Delta t^l} \bar{R}^l$ is classical since it is equivalent to inverting a linear system.

### 3.3.2. Douglas-Rachford algorithm

The second method that we consider is the Douglas-Rachford algorithm [70, 72]. The original method was introduced as a first-order-in-time operator splitting to solve the evolution problem

$$ \left\{ \begin{array}{l} \partial_{\tau} U + (A+B)U \ni R, \\ U(\tau=0) = U^0 \in \mathcal{D}(A) \cap \mathcal{D}(B), \end{array} \right. \qquad (17) $$

where $A$ and $B$ are maximal monotone operators. The algorithm is described as follows: we choose $V^0 \in (\mathcal{I} + \Delta \tau B)U^0$ and for $m \ge 0$ define

$$ V^{m+1} = J_{\bar{A}}^{\Delta \tau} (2J_{\bar{B}}^{\Delta \tau} - \mathcal{I})V^m + (\mathcal{I} - J_{\bar{B}}^{\Delta \tau})V^m + \Delta \tau J_{\bar{A}}^{\Delta \tau} R. $$

The approximation of $U$ at time $m\Delta\tau$ is then $U^m = J_{\bar{B}}^{\Delta\tau} V^m$. A crucial property of this method is that as $m \to \infty$, $U^m$ converges to the steady-state solution $U^\infty$ of eq. (17) (i.e. $(A+B)U^\infty \ni R$) if the latter exists.

We adopt the Douglas-Rachford algorithm to solve the semi-discrete problem eq. (11) by considering it as a stationary problem. More precisely, we
search for the steady-state solution of the following pseudotime [73] problem,

$$
\begin{cases}
\partial_{\tau} U(\tau) + \left( \frac{\mathcal{I}}{\Delta t^l} + \mathcal{A}^l + \mathcal{B}_{\psi} \right) U(\tau) \ni \frac{1}{\Delta t^l} \mathbf{R}^l, \\
U(0) = \mathbf{U}^l,
\end{cases}
$$

where $\tau$ is the pseudotime. We apply the Douglas-Rachford algorithm for

$$
A = \frac{\mathcal{I}}{\Delta t^l} + \mathcal{A}^l, \quad B = \mathcal{B}_{\psi}, \quad R = \frac{1}{\Delta t^l} \mathbf{R}^l, \quad \Delta \tau = \Delta t^l,
$$

and when it converges, we put $\mathbf{U}^{l+1} = \mathbf{U}^{\infty}$.

The requirement that $A$ is monotone imposes a constraint on the time step $\Delta t^l$ which is related to the Townsend ionization rate $\alpha$. For simplicity, let us assume that $D_i$, $\mathbf{u}_i$, $\alpha$ are constant and other reaction coefficients are null. We define the scalar product on the space $\mathcal{H} \equiv H \times H \times H$ as

$$
(\mathbf{U}, \mathbf{V})_{\mathcal{H}} = (n_e, w_e) + (n_p, w_p) + (n_n, w_n),
$$

with $\mathbf{U} = (n_e \ n_p \ n_n)^t$, $\mathbf{V} = (w_e \ w_p \ w_n)^t$. Since $A$ is a linear operator it is sufficient to check the sign of $(A\mathbf{U}, \mathbf{U})_\mathcal{H}$. By some algebraic manipulation we have

$$
(A \mathbf{U}, \mathbf{U})_{\mathcal{H}} = \frac{1}{\Delta t^l} \sum_{i \in \mathcal{L}} |n_i|^2 + \sum_{i \in \mathcal{L}} D_i |\nabla n_i|^2 + \frac{1}{2} \sum_{i \in \mathcal{L}} \int_{\Gamma} (\mathbf{u}_i \cdot \mathbf{\nu}) n_i^2 d\mathbf{l} - \alpha N |n_e|(1 + |n_p|).
$$

For $(A\mathbf{U}, \mathbf{U})_{\mathcal{H}} \ge 0$, $\Delta t^l$ must satisfy

$$
\Delta t^l \le \left( \alpha N |n_e|(1 + |n_p|) - \sum_{i \in \mathcal{L}} D_i |\nabla n_i|^2 - \frac{1}{2} \sum_{i \in \mathcal{L}} \int_{\Gamma} (\mathbf{u}_i \cdot \boldsymbol{\nu}) n_i^2 \mathrm{d}\mathbf{l} \right)^{-1} \sum_{i \in \mathcal{L}} |n_i|^2. \quad (18)
$$

In case of very large $\alpha$ such as in streamers or microdischarges, the time step
is severely restricted in spite of implicit integration.

Finally, the discrete version of the Douglas-Rachford method is presented
in algorithm 1.

3.3.3. Gauss-Seidel algorithm

For a discrete linear system $AU = R$ where $A$ is a square matrix, the
solution $U$ can be obtained iteratively via

$$
U^{m+1} = (M + D)^{-1} (R - N U^m), \quad m \ge 0,
$$
**Algorithm 1** Douglas-Rachford

1: Let $\epsilon > 0$, $U^0 = \overline{U}^l$, $V^0 \in (\overline{\mathcal{I}} + \Delta t^l \mathcal{B}_{\bar{\psi}})U^0$, $\delta > \epsilon$

2: **while** $\delta > \epsilon$ **do**

3: Set $V^{m+1} = J_{\overline{\mathcal{I}}+\Delta t^l \overline{\mathcal{A}}^l}^1(2J_{\mathcal{B}_{\bar{\psi}}} - \overline{\mathcal{I}})V^m + (\overline{\mathcal{I}} - J_{\mathcal{B}_{\bar{\psi}}})V^m + J_{\overline{\mathcal{I}}+\Delta t^l \overline{\mathcal{A}}^l}^1 \overline{\mathcal{R}}^l$

4: Set $U^{m+1} = J_{\mathcal{B}_{\bar{\psi}}}V^{m+1}$ and let $(n_e^{(m+1)}, n_p^{(m+1)}, n_n^{(m+1)})^t = U^{m+1}$

5: Update $\delta = \max_{i \in \mathcal{L}} \left( \frac{|n_i^{(m+1)} - n_i^{(m)}|_{\infty}}{|n_i^{(m)}|_{\infty}} \right)$ where $| \cdot |_{\infty}$ is the $L^\infty$ norm defined on $H$

6: **end while**

7: Set $\overline{U}^{l+1} = U^m$

where $U^0$ is an initial guess, $D$, $M$, $N$ are resp. the diagonal, strictly lower and upper triangular components of $A$. By using the forward substitution, the elements of $U^m$ can be computed sequentially as follows,

$$U_p^{m+1} = \frac{1}{1 + A_{pp}} \left( R_p - \sum_{q>p} A_{pq} U_q^m - \sum_{q<p} A_{pq} U_q^{m+1} \right), \quad p = 1, 2, \dots,$$

where $U_p^m$ is the $p^{th}$ element of $U^m$ and $A_{pq}$ is the element of $A$ on the $p^{th}$ row and $q^{th}$ column.

For the nonlinear problem eq. (14), by analogy to the Gauss-Seidel method we propose the following iteration for forward substitution,

$$U^{m+1} = (\overline{\mathcal{I}} + \Delta t^l \mathcal{M}^l + \Delta t^l \mathcal{D}^l + \Delta t^l \mathcal{B}_{\bar{\psi}})^{-1} (\overline{\mathcal{R}}^l - \Delta t^l \mathcal{N}^l U^m),$$

where $\mathcal{D}^l$, $\mathcal{M}^l$, $\mathcal{N}^l$ are resp. the diagonal, strictly lower and upper triangular components of $\overline{\mathcal{A}}^l$. Its element-wise form reads as

$$(I + \frac{\Delta t^l}{1 + \Delta t^l \overline{\mathcal{A}}_{pp}^l} \mathcal{B}_{\bar{\psi}}) U_p^{m+1} = \frac{1}{1 + \Delta t^l \overline{\mathcal{A}}_{pp}^l} \left( \overline{\mathcal{R}}_p^l - \Delta t^l \sum_{q>p} \overline{\mathcal{A}}_{pq}^l U_q^m - \Delta t^l \sum_{q<p} \overline{\mathcal{A}}_{pq}^l U_q^{m+1} \right).$$

As long as $1 + \Delta t^l \overline{\mathcal{A}}_{pp}^l > 0^8$, $\forall p$ so that $\frac{\Delta t^l}{1 + \Delta t^l \overline{\mathcal{A}}_{pp}^l} \mathcal{B}_{\bar{\psi}}$ is monotone, the

<footnote n="8">⁸this imposes a restriction on $\Delta t^l$ which is a discrete and local (on each cell $K \in \mathcal{T}$) version of eq. (18)</footnote>
resolvent of $\mathcal{B}_{\bar{\psi}}$ is well defined as in eq. (15) and we have

$$U_p^{m+1} = \frac{J_{\mathcal{B}_{\bar{\psi}}}}{1 + \Delta t^l \mathcal{A}_{pp}^l} \left( \bar{\mathbf{R}}_p^l - \Delta t^l \sum_{q>p} \bar{\mathcal{A}}_{pq}^l U_q^m - \Delta t^l \sum_{q<p} \bar{\mathcal{A}}_{pq}^l U_q^{m+1} \right).$$

The complete description of the proposed algorithm is presented in algorithm 2.

**Algorithm 2 Gauss-Seidel**

1: Let $\epsilon > 0$, $U^0 = \bar{U}^l$, $\delta > \epsilon$

2: **while** $\delta > \epsilon$ **do**

3: Forward substitution

4: $U^{m+\frac{1}{2}} = (\bar{\mathcal{I}} + \Delta t^l \mathcal{M}^l + \Delta t^l \mathcal{D}^l + \Delta t^l \mathcal{B}_{\bar{\psi}})^{-1} (\bar{\mathbf{R}}^l - \Delta t^l \mathcal{N}^l U^m)$

5: Backward substitution

6: $U^{m+1} = (\bar{\mathcal{I}} + \Delta t^l \mathcal{N}^l + \Delta t^l \mathcal{D}^l + \Delta t^l \mathcal{B}_{\bar{\psi}})^{-1} (\bar{\mathbf{R}}^l - \Delta t^l \mathcal{M}^l U^{m+\frac{1}{2}})$

7: Let $(n_e^{(m+1)}, n_p^{(m+1)}, n_n^{(m+1)})^t = U^{m+1}$

8: Update $\delta = \max_{i \in \mathcal{L}} \left( \frac{|n_i^{(m+1)} - n_i^{(m)}|_{\infty}}{|n_i^{(m)}|_{\infty}} \right)$

9: **end while**

10: Set $\bar{U}^{l+1} = U^m$

**4. Comparison of algorithms in one-dimensional discharges**

Our first test case is the wire-to-wire corona discharge proposed in [23]. The sketch of the actuator is illustrated in fig. 3. It is composed of a small wire of radius $r_1 = 0.35mm$ and a bigger wire $r_2 = 1mm$, parallel to each other and $d = 40mm$ apart. The wires are $L = 16cm$ long. Between them we apply a voltage $V_G = 40kV$. The computational domain is the straight line $\Gamma_0$ in the figure, which is well meshed near the electrodes to capture correctly the ionization process.

If we suppose that the plasma is uniform along the wires then the three-dimensional problem can be easily transformed into a two-dimensional one. In order to transform it into a one-dimensional problem, we suppose that the plasma is contained in a layer whose thickness $S$ depends on the space coordinate $x$ (see fig. 3). $S(x)$ is compute from the relation $\frac{d(V_{ext}S)}{dx} = 0$,
<figure><img src="image_4.png" /><figcaption>Figure 3: Wire-to-wire actuator schematics</figcaption></figure>

where $V_{ext}$ is the analytic expression of the electrostatic potential generated by the electrodes, which can be found in [23]. We set $\max_x S = 5mm$.

We set $n_p^0 = n_e^0 = 10^9 m^{-3}$ accordingly to section 2.3. On the contrary, the seed negative ion density is considered negligible: $n_n^0 = 1m^{-3}$.

For this one-dimensional test, we suppose that $D_p = D_n = 0$. By affecting the variable change $v_i = Sn_i$, the quasi-2D model reads as follows,

$$ 
\left\{ 
\begin{array}{l}
\partial_t v_e + \nabla \cdot (\boldsymbol{u}_e^S v_e - D_e \nabla v_e) \in (\alpha - \eta) N v_e - k_{ep}^S v_e v_p + S_{S\psi}, \\
\partial_t v_p + \nabla \cdot (\boldsymbol{u}_p v_p) = \alpha N v_e - k_{ep}^S v_e v_p - k_{np}^S v_n v_p, \\
\partial_t v_n + \nabla \cdot (\boldsymbol{u}_n v_n) = \eta N v_e - k_{np}^S v_n v_p, \\
v_i(t=0) = Sn_i^0, \quad i \in \mathcal{L}, \\
-\nabla \cdot (\varepsilon S \nabla \phi) = q \sum_{i \in \mathcal{L}} z_i v_i,
\end{array}
\right. 
$$

with $\boldsymbol{u}_e^S = -\mu_e \boldsymbol{E} + D_e S'/S$, $k_{ep}^S = k_{ep}/S$ and $k_{np}^S = k_{np}/S$.

### 4.1. Conservation of steady state

We define the steady state of a corona discharge as when the value of the electric current does not change more than $10^{-4}\%$. The analytic expression of the total current $I$ is given in [23]. For the described test case, it happens that the steady state exists, but the capacity of the code to reproduce this steady state depends strongly on the algorithm used to solve the species conservation laws as we show in the next paragraphs. In this section, all the simulations are conducted on a 400-element grid with the smallest cell size $\Delta x_{min} = 5.5 \mu m$, and using the Scharfetter-Gummel scheme as the flux approximation method.
We compute the time step as $\Delta t = \min(c\Delta t_{ions}, \Delta t_{\phi})^9$. When the simulations occur between 0 and $0.04\text{ms}$, we fix $c = 10^{-2}$ since at the discharge onset, there is a rapid charge multiplication and as a result, the dielectric relaxation time $\Delta_{\phi}$ decreases quickly and so having a small time step $\Delta t$ ensures that the plasma dynamics are captured correctly. After $0.04\text{ms}$, we choose the value of $c$ between 1 and $10^2$. After around $0.4\text{ms}$, the discharge enters the ion collection phase and we have $c\Delta t_{ions} < \Delta t_{\phi}$.

*Lie splitting.* The computed electric currents with the Lie algorithm are shown in fig. 4. The discharges were simulated until $T = 7\text{ms}$ when we obtained steady states. We show that the Lie splitting does not conserve the steady state: indeed, the current changes when we vary the time step. We also show the electric current obtained with the forward Euler scheme as reference (black dashed line).

<figure><img src="image_5.png" /><figcaption>Figure 4: Lie splitting: electric current $I$ with different CFL numbers $c$</figcaption></figure>

*Douglas-Rachford algorithm.* We defined the stop criterion as $\epsilon = 0.0001\%$ (see algorithm 1). The electric current at steady state is conserved by the Douglas-Rachford algorithm as shown in fig. 5. This feature was remarked

<footnote>$^9c$ is the CFL number</footnote>
in [70]. At $T = 7\,\text{ms}$, the current of $\Delta t = 10^2 \times \Delta t_L$ differs only 0.007% from the current of $\Delta t = \Delta t_L$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Douglas-Rachford algorithm: electric current $I$ with different CFL numbers $c$</figcaption></figure>

*Gauss-Seidel algorithm.* We defined the stop criterion as $\epsilon = 10^{-4}\%$ (see algorithm 2). The electric current at steady state is also conserved by the Gauss-Seidel algorithm as shown in fig. 6. At $T = 7\,\text{ms}$, the current of $\Delta t = 10^2 \times \Delta t_L$ differs only $8 \times 10^{-7}\%$ from the current of $\Delta t = \Delta t_L$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Gauss-Seidel algorithm: electric current $I$ with different CFL numbers $c$</figcaption></figure>

## 4.2. Performance comparison of DR and GS algorithms

We have seen that the steady state current computed with the Douglas-Rachford (DR) and Gauss-Seidel (GS) algorithms is virtually independent of the time step. Therefore, these algorithms should be preferred over the Lie splitting to conduct implicit simulations of corona discharges. However, there is a huge difference between the two methods in terms of computation time as shown in table 1.

In order to understand why, we count the number of iterations $N(t^l)$ at each time level $t^l$ that is necessary for the algorithms to converge. Then we
evaluate the averaged number of iterations $\overline{N}$ by dividing it with the number of time levels. Furthermore, we only count the iterations from $T = 1ms$ to $7ms$, since before $1ms$ the dielectric time step $\Delta t_{\phi}$ can fall below $c\Delta t_{ions}$ so $N(t^l)$ would not depend on $c$. To sum up, we have,

$$ \overline{N} = \frac{\sum_{t^l=1ms}^{7ms} N(t^l)}{\text{Number of time levels}} $$

<table><caption>Table 1: Averaged number of iterations & CPU time (minutes) of Douglas-Rachford and Gauss-Seidel algorithms for corona discharge simulations on the 400-element grid</caption><thead><tr><th rowspan="2">c</th><th colspan="2">DR</th><th colspan="2">GS</th></tr><tr><th>N&#773;</th><th>CPU time</th><th>N&#773;</th><th>CPU time</th></tr></thead><tbody><tr><td>1</td><td>14</td><td>45</td><td>1</td><td>17</td></tr><tr><td>10</td><td>84</td><td>26</td><td>1.13</td><td>4.65</td></tr><tr><td>100</td><td>856</td><td>24</td><td>1.33</td><td>1.68</td></tr></tbody></table>

We remark from table 1 that the averaged number of iterations of the DR algorithm scales almost linearly with $c$. This explains why the DR simulation time is not shortened although we increase the time step. On the contrary, the averaged number of iterations of the GS algorithm is almost the same regardless of the time step and boosts significantly the simulation time. Therefore, we exclusively use the GS algorithm from now on.

### 4.3. Mesh convergence

In this section, we conduct a mesh convergence study to ensure that the Gauss-Seidel algorithm, combining with the SGCC flux schemes [68], is able to capture correctly the discharge dynamics. We perform the simulations on five grids with 400 cells ($\Delta x_{min} = 5.5\mu m$), 800 cells ($2.6\mu m$), 1600 cells ($1.3\mu m$), 3200 cells ($0.6\mu m$) and 6400 cells ($0.3\mu m$). We use the Scharfetter-Gummel scheme and the SGCC$_1$, SGCC$_2$ schemes [68] for flux approximation. Since the latter two are non-linear because of slope limiters, we use the fixed-point algorithm to iterate on the slopes$^{10}$.

<footnote>$^{10}$see Appendix B</footnote>
The electric currents of each flux scheme are shown in fig. 7, fig. 8 and fig. 9. They are compared to the current of the SGCC₂ scheme obtained on the 0.3µm-grid (the black dashed curve). We can conclude that the whole numerical scheme is able to perform corona discharge simulation with confidence by observing that the electric currents converged to the black reference curve as the mesh size decreased. We also remark the high-order flux schemes (SGCC₁, SGCC₂) increased significantly the numerical precision.

<figure><img src="image_8.png" /><figcaption>Figure 7: Scharfetter-Gummel flux scheme: electric current *I* with different mesh size $\Delta x_{min}$</figcaption></figure>

Finally, we note that a mesh convergence study for the SGCC schemes with explicit time integration was partially conducted in [68] but the simulation with the 3200-element grid was only launched until 300µs. The reason is that it would have taken *months* to obtain the steady state current. In this paper, the simulations are rapid. For example, a computation with the SGCC₂ scheme takes only *an hour* on the 3200-element grid and *4 hours* on the 6400-element grid<fnref n="11" />.

<footnote n="11">¹¹on a Dell Latitude 5410</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8: SGCC₁ flux scheme: electric current *I* with different mesh size Δ*x*<sub>min</sub></figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: SGCC₂ flux scheme: electric current *I* with different mesh size Δ*x*<sub>min</sub></figcaption></figure>
5. Two-dimensional simulations of DC corona discharges

The same wire-to-wire discharges in section 4 are now studied here. A
sketch of the actuator in two-dimension is presented in fig. 10. The electrodes
are connected to a resistor R. The presence of the circuit resistance modifies
of course the voltage V_a at the anode following the relation V_a = V_G - RI.
The computation of V_a is given in [23]. The zone inside the solid-line loop is
the computation domain. Since we assume that the wires are hanged in air,
the bottom boundary condition is symmetric. The domain height equals to
h. The distance between an electrode and its closest vertical border equals
to l.

<figure><img src="image_11.png" /><figcaption>Figure 10: Sketch of the computation domain</figcaption></figure>

The particular geometry of the wire-to-wire actuator intensifies the electrostatic field in the electrodes vicinity [74]. Therefore, all the chemical reactions are much stronger in these zones than they are in the region between the wires. As a result, the numerical grids that we use are strongly refined near the electrodes as illustrated in fig. 11. The minimal grid size is fixed to a value Δx<sub>min</sub>.

We set $n_p^0 = n_e^0 = 10^9 m^{-3}$ accordingly to section 2.3. On the contrary, the seed negative ion density is considered negligible: $n_n^0 = 1m^{-3}$.

The characteristics of the actuator and the numerical parameters are
grouped in table 2.
<figure><img src="image_12.png" /><figcaption>Figure 11: Example of grid refinement near the smaller wire generated by Gmsh [66]</figcaption></figure>

<table><caption>Table 2: Parameters of the wire-to-wire simulations</caption><thead><tr><th>r<sub>1</sub></th><th>0.1 to 0.5mm</th><th>V<sub>G</sub></th><th>9 to 37kV</th><th>Δx<sub>min</sub></th><th>5μm</th></tr></thead><tbody><tr><th>r<sub>2</sub></th><td>1mm</td><th>R</th><td>10kΩ</td><th>n<sub>e</sub><sup>0</sup></th><td>10<sup>9</sup>m<sup>-3</sup></td></tr><tr><th>L</th><td>20cm</td><th>h</th><td>8mm</td><th>n<sub>p</sub><sup>0</sup></th><td>10<sup>9</sup>m<sup>-3</sup></td></tr><tr><th>d</th><td>10 to 40mm</td><th>l</th><td>4mm</td><th>n<sub>n</sub><sup>0</sup></th><td>1m<sup>-3</sup></td></tr></tbody></table>

## 5.1. Steady-state solution

We present the numerical results for a case of $V_G = 13kV$, $r_1 = 0.1mm$, $d = 10mm$, $\psi = 10^{10}m^{-3}$ and $\Delta x_{min} = 5\mu m$. To avoid the formation of streamers that can reduce significantly the dielectric time step and prolong the simulation, we increase the voltage gradually on $2\mu s$ until it reaches the maximal value $V_G$.

We set the CFL number to $c = 10^3$ and run the simulation on four MPI processes$^{12}$. The waiting time of our implicit method is about three to four hours instead of a week for explicit methods in COPAIER.

Figure 12 shows the $x$-component of the EHD force density $F_x \equiv \rho E_x$ at $T = 4ms$ when the discharge reaches the steady state. The result shows that between the electrodes, there is a force density of $16kNm^{-3}$ which directs towards the cathode. On the left of the anode, the EHD reverses sign and reaches $-4.6kNm^{-3}$. Overall, the force density radiates from and concentrates near the anode.

<footnote>$^{12}$see Appendix C</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 12: EHD force density $F_x$ ($Nm^{-3}$) at $T = 4ms$ for $r_1 = 0.1mm$, $d = 10mm$, $V_G = 13kV$, $\psi = 10^{10}m^{-3}$ (figure reflected on symmetric boundaries)</figcaption></figure>

## 5.2. Influence of the floor density on numerical discharges

We conduct a parametric study in which we compute the steady state current for each (constant) floor density $\psi$. In fig. 13, we show the $\psi - I$ curves for fixed $r_1 = 0.1mm$ and we change the wire distance $d$. For each $d$, the circuit voltage $V_G$ is chosen closed to the maximal potential reached before the apparition of spark discharges, as documented in [7]. In fig. 14, we fix $d = 10mm$ and change the anode radius $r_1$.

In the first case, the minimal electron density, when there is no floor density, is more than $10^8 m^{-3}$ for all $d$, suggesting that the discharges are self-sustaining and do not need other sources of electrons other than impact ionization and secondary emission from the cathode surface. For this reason, we do not show in fig. 13 the values of $I$ for $\psi$ smaller than $5 \times 10^8 m^{-3}$. Although the discharges maintain themselves, the current is far from experimental value if there is no additional electron source (see section 5.3 for $d = 10mm$). The current $I$ clearly depends on $\psi$ as suggested in fig. 13, but the dependence seems to differ in each value interval of $\psi$. For $\psi = 10^9 - 10^{11}m^{-3}$, we have the estimated law $I \sim \psi^{0.07}$. For $\psi = 10^{12} - 10^{14}m^{-3}$, we have $I \sim \psi^{0.7}$. This suggests that $\psi$ should be around $10^{10}m^{-3}$ for this wire-to-wire discharge, since a too high $\psi$ may lead to a too high electric current and potentially a spark discharge.

In fig. 14, we highlight the fact that in certain cases, the discharge is not self-sustaining unless there is an additional source of electrons. The example of the case $r_1 = 0.5mm$ shows that the current $I$ is close to 0 when $\psi$ is small
<figure><img src="image_14.png" /><figcaption>Figure 13: $\psi - I$ curves with $r_1 = 0.1\text{mm}$ and $V_G = 13\text{kV}$ for $d = 10\text{mm}$, $V_G = 20\text{kV}$ for $d = 20\text{mm}$, $V_G = 37\text{kV}$ for $d = 40\text{mm}$</figcaption></figure>

$(1\text{m}^{-3})$, suggesting that there is no discharge at all though this is not true (see section 5.3). Therefore, the inclusion of $\psi$ is necessary.

### 5.3. Comparison with experimental data

Reproducing experimental data is the first step for a plasma solver to become a predictive simulation tool. In this section, we compare our numerical results to the data of Bérard et al. [7]. The simulations are conducted for $d = 10\text{mm}$, $r_1$ ranging from 0.1 to 0.5mm and $V_G$ ranging from 9 to 17kV. Figure 15 shows the numerical $V-I$ curves with filled markers and the experimental ones with hollow circle markers.

As an attempt to fit the given data, we use two values of $\psi$, $10^{10}\text{m}^{-3}$ (shown by filled square markers) and $10^{11}\text{m}^{-3}$ (filled triangles). For $r_1 = 0.1\text{mm}$, the numerical currents of $\psi = 10^{10}\text{m}^{-3}$ fit the data very well. For $r_1 = 0.175\text{mm}$, the currents of $\psi = 10^{11}\text{m}^{-3}$ only agree with the experiments for high voltages, and as $V_G$ decreases so does the slope of the numerical curve. The two smallest values seem to agree more with the $10^{10}\text{m}^{-3}$-curve but the curve itself deviates too much for large voltages. The situation is quite similar for $r_1 = 0.2\text{mm}$, although the numerical results differ even more from experimental data.
<figure><img src="image_15.png" /><figcaption>Figure 14: $\psi - I$ curves with $d = 10mm$ and $V_G = 13kV$ for $r_1 = 0.1mm$, $V_G = 16.5kV$ for $r_1 = 0.5mm$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: $V_G - I$ curves with $d = 10mm$. Numerical currents with $\psi = 10^{10}m^{-3}$ (filled squares) and $\psi = 10^{11}m^{-3}$ (triangles). Experimental data in hollow circles.</figcaption></figure>
For corona discharges, the $V_G - I$ characteristics is given by the law $I = kV_G(V_G - V_c)$ [38, Chapter 12] where $V_c$ is the breakdown voltage and $k$ is a coefficient. The fitting of experimental and numerical currents with this quadratic law gives $k (\mu A(kV)^{-2})$ and $V_c (kV)$ for each curve. These characteristics are laid out in Table 3.

<table><caption>Table 3: $V_G - I$ characteristics of experimental data ('data') and numerical results ($\psi = 10^{10}, 10^{11}m^{-3}$)</caption><thead><tr><th rowspan="2"></th><th colspan="3">r<sub>1</sub> = 0.1mm</th><th colspan="3">r<sub>1</sub> = 0.175mm</th><th colspan="3">r<sub>1</sub> = 0.2mm</th></tr><tr><th>data</th><th>ψ = 10<sup>10</sup></th><th>ψ = 10<sup>11</sup></th><th>data</th><th>ψ = 10<sup>10</sup></th><th>ψ = 10<sup>11</sup></th><th>data</th><th>ψ = 10<sup>10</sup></th><th>ψ = 10<sup>11</sup></th></tr></thead><tbody><tr><th>k</th><td>2.42</td><td>2.3</td><td>2.33</td><td>2.8</td><td>1.71</td><td>2.13</td><td>2.77</td><td>1.69</td><td>2.24</td></tr><tr><th>V<sub>c</sub></th><td>8.88</td><td>8.54</td><td>7.64</td><td>10.55</td><td>10.66</td><td>9.68</td><td>11.53</td><td>11.51</td><td>10.51</td></tr></tbody></table>

The results show that the breakdown voltages $V_c$ of $\psi = 10^{10}$ are closer to experimental data. $V_c$ decreases with increasing $\psi$ since the breakdown condition is favorable with large floor density. However, the proportional coefficient $k$ of $\psi = 10^{11}$ is closer to experimental data. Overall, we find that the proposed numerical model provides quite good estimates of experimental data, considering the simplicity of the model, the first-order discretization and other factors in the experiments that could change the current, such as impurities on electrode surface.

## 6. Conclusions

The main results of this study are summarized as follows.

1. The continuity equations were reformulated into the form of differential inclusions. The constraint on the electrons $n_e \ge \psi$ were integrated into the conservation laws via the subdifferential of a characteristic function of a convex set.

2. The conservation laws were discretized in time using the backward Euler scheme and ensued a system of inclusions that was nonlinear because of the subdifferential terms. Some numerical methods were studied to solve this system. The Douglas-Rachford and Gauss-Seidel algorithms demonstrated to conserve the steady state of the discharges. The Gauss-Seidel algorithm converged however much faster and thus was used by default for simulations.
3. The proposed time integration approach reduces drastically the computation time of DC corona discharges: an implicit simulation takes only *hours* or *days* to reach the steady-state, while an explicit simulation takes *weeks* or even *months*.

4. We have proposed a simplified discharge model that is capable to reduce significantly the simulation time but also reproduce with quite good agreement the experimental data of wire-to-wire DC corona discharges. Works are still needed to further improve the quality of numerical results.

## Appendix A. Implementation of the Scharfetter-Gummel scheme on two-dimension grids

<figure><img src="image_17.png" /><figcaption>Figure A.16</figcaption></figure>

Let us consider a conformal triangulation<fnref n="13" /> *T* of the domain Ω. Let *K*, *L* ∈ *T* be two neighbor triangles sharing the edge λ<sub>KL</sub> and *x*<sub>KL</sub> be the midpoint of λ<sub>KL</sub>. Let *ŷ*<sub>KL</sub> (resp. *ŷ*<sub>LK</sub>) be the intersection of the line passing through *x*<sub>KL</sub> orthogonal to λ<sub>KL</sub> and one of other edges of *K* (resp. *L*) such that the segment bounded by *x*<sub>KL</sub>, *ŷ*<sub>KL</sub> (resp. *ŷ*<sub>LK</sub>) is contained in *K* (resp. *L*) (see fig. A.16). Let $d_{KL} = \frac{1}{2} |\hat{x}_{KL} - \hat{x}_{LK}|$.

We use the following extension of the Scharfetter-Gummel scheme [67] on

<footnote n="13">¹³i.e. no hanging nodes</footnote>
two-dimensional grids to approximate the flux $f_{i,K,\lambda_{KL}} \cdot \nu_{K,\lambda_{KL}}$ in eq. (12),

$$ \bar{f}_{KL} = \frac{D_{i,KL}}{d_{KL}} (\mathcal{B}(\mathfrak{p}_{i,KL})\bar{n}_{i,K} - \mathcal{B}(-\mathfrak{p}_{i,KL})\bar{n}_{i,L}), \quad i \in \mathcal{L}, $$

where $\mathcal{B}(z) \equiv \frac{z}{\exp(z) - 1}$ is the Bernoulli function, $\mathfrak{p}_{i,KL} \equiv -\frac{u_{i,KL}d_{KL}}{D_{i,KL}}$ is the numerical Péclet number, $u_{i,KL}$ (resp. $D_{i,KL}$) is an approximation of $u_i(\mathbf{x}_{KL}) \cdot \nu_{K,\lambda_{KL}}$ (resp. $D_i(\mathbf{x}_{KL})$).

## Appendix B. Implicit integration of the SGCC schemes

The $(p+1)^{th}$-order SGCC flux scheme [68] reads as follows,

$$ \bar{f}_{KL}^p = \frac{D_{i,KL}}{d_{KL}} (\mathcal{B}(\mathfrak{p}_{i,KL})\bar{n}_{i,K}^p - \mathcal{B}(-\mathfrak{p}_{i,KL})\bar{n}_{i,L}^p), \quad i \in \mathcal{L}, $$

It can be decomposed as

$$ \bar{f}_{KL}^p = \bar{f}_{KL} + \tilde{f}_{KL}^p, $$

which is the sum of the Scharfetter-Gummel flux and a high-order correction which reads

$$ \tilde{f}_{KL}^p = \frac{D_{i,KL}}{d_{KL}} (\mathcal{B}(\mathfrak{p}_{i,KL})(\bar{n}_{i,K}^p - \bar{n}_{i,K}) - \mathcal{B}(-\mathfrak{p}_{i,KL})(\bar{n}_{i,L}^p - \bar{n}_{i,L})), $$

In the above equation, $(\bar{n}_{i,K}^p - \bar{n}_{i,K})$ is a weighted sum of reconstructed derivatives of $n_i$ on $K$. In the presence of slope limiters to ensure the stability of the scheme, it is a nonlinear term even if the problem is linear.

Following the steps in section 3.2, we can write the discretization of the conservation laws with the SGCC scheme as follows,

$$ (\bar{\mathcal{I}} + \Delta t^l \bar{\mathcal{A}}^l + \Delta t^l \tilde{\mathcal{A}}^{p,l} + \Delta t^l \mathcal{B}_{\bar{\psi}}) \bar{\mathcal{U}}^{l+1} \ni \bar{\mathcal{R}}^l, $$

where $\tilde{\mathcal{A}}^{p,l}$ is the nonlinear operator ensued from $\tilde{f}^p$. We adapt a fixed-point method to solve this problem (see algorithm 3). This algorithm has been only applied for one-dimensional simulations.
**Algorithm 3 Fixed-point**

1: Let $\epsilon > 0$, $U^{(0)} = \bar{U}^l$, $\delta > \epsilon$

2: **while** $\delta > \epsilon$ **do**

3: Set $V^{(q)} = \tilde{\mathcal{A}}^{p,l}U^{(q)}$

4: Use algorithm 1 or algorithm 2 to solve $(\bar{\mathcal{I}} + \Delta t^l \bar{\mathcal{A}}^l + \Delta t^l \mathcal{B}_{\bar{\psi}})U^{(q+1)} \ni \bar{\mathcal{R}}^l - \Delta t V^{(q)}$

5: Let $(n_e^{(q+1)}, n_p^{(q+1)}, n_n^{(q+1)})^t = U^{q+1}$

6: Update $\delta = \max_{i \in \mathcal{L}} \left( \frac{|n_i^{(q+1)} - n_i^{(q)}|_{\infty}}{|n_i^{(q)}|_{\infty}} \right)$

7: **end while**

8: Set $\bar{U}^{l+1} = U^q$

**Appendix C. MPI implementation of the Gauss-Seidel algorithm**

The sequential nature of the Gauss-Seidel algorithm makes the parallel-computing implementation not evident. In fact, if we imagine a simple partition of the simulation domain (see fig. C.17) and the variables on the second subdomain are enumerated after those on the first subdomain, then the density update on the second subdomain would have to wait for the update on the first subdomain to finish before getting started. For this reason, we rather send $U^{(m)}$ than $U^{(m+1)}$ of the border zone (yellow-filled triangles) from the first subdomain to the second and vice versa. Meanwhile, the interior variables (on red-filled triangles) are updated normally by the Gauss-Seidel algorithm. The overall method resembles a block-Jacobi method where each block is affiliated to a subdomain.

<figure><img src="image_18.png" /><figcaption>Figure C.17: In this figure the computation domain is partitioned into two subdomains, sharing the thick black line as interface. The triangles which share an edge with the interface are called border cells and are painted in yellow. The other triangles are called interior cells and are painted in red.</figcaption></figure>
References

[1] E. Moreau, Airflow control by non-thermal plasma actuators, Journal of physics D: applied physics 40 (3) (2007) 605.

[2] K.-S. Choi, J.-H. Kim, Plasma virtual roughness elements for cross-flow instability control, Experiments in fluids 59 (10) (2018) 1–15.

[3] Y. Zhu, Y. Wu, B. Wei, H. Xu, H. Liang, M. Jia, H. Song, Y. Li, Nanosecond-pulsed dielectric barrier discharge-based plasma-assisted anti-icing: modeling and mechanism analysis, Journal of Physics D: Applied Physics 53 (14) (2020) 145205.

[4] H. Xu, Y. He, K. L. Strobel, C. K. Gilmore, S. P. Kelley, C. C. Hennick, T. Sebastian, M. R. Woolston, D. J. Perreault, S. R. Barrett, Flight of an aeroplane with solid-state propulsion, Nature 563 (7732) (2018) 532–535.

[5] R. Van Dyken, T. McLaughlin, C. Enloe, Parametric investigations of a single dielectric barrier plasma actuator, in: 42nd AIAA aerospace sciences meeting and exhibit, 2004, p. 846.

[6] C. Enloe, T. E. McLaughlin, R. D. Van Dyken, K. Kachner, E. J. Jumper, T. C. Corke, Mechanisms and responses of a single dielectric barrier plasma actuator: plasma morphology, AIAA journal 42 (3) (2004) 589–594.

[7] P. Bérard, D. Lacoste, C. Laux, Measurements and simulations of the ionic wind produced by a dc corona discharge in air, helium and argon, in: 38th AIAA Plasmadynamics and Lasers Conference In conjunction with the 16th International Conference on MHD Energy Conversion, 2007, p. 4611.

[8] C. Porter, J. Baughn, T. McLaughlin, C. Enloe, G. Font, Plasma actuator force measurements, AIAA journal 45 (7) (2007) 1562–1570.

[9] C. Enloe, M. McHarg, T. E. McLaughlin, Time-correlated force production measurements of the dielectric barrier discharge plasma aerodynamic actuator, Journal of applied physics 103 (7) (2008) 073302.
[10] A. R. Hoskinson, N. Hershkowitz, D. E. Ashpis, Force measurements of single and double barrier dbd plasma actuators in quiescent air, Journal of Physics D: Applied Physics 41 (24) (2008) 245209.

[11] D. Opaits, A. Likhanskii, M. Edwards, S. Zaidi, M. Shneider, S. Macheret, R. Miles, Surface plasma induced wall jets, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2010, p. 469.

[12] A. Debien, N. Benard, L. David, E. Moreau, Unsteady aspect of the electrohydrodynamic force produced by surface dielectric barrier discharge actuators, Applied Physics Letters 100 (1) (2012) 013901.

[13] Y. Zhang, L. Liu, Y. Chen, J. Ouyang, Characteristics of ionic wind in needle-to-ring corona discharge, Journal of Electrostatics 74 (2015) 15-20.

[14] J. R. Roth, X. Dai, Optimization of the aerodynamic plasma actuator as an electrohydrodynamic (ehd) electrical device, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, 2006, p. 1203.

[15] T. Corke, B. Mertz, M. Patel, Plasma flow control optimized airfoil, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, 2006, p. 1208.

[16] T. C. Corke, M. L. Post, D. M. Orlov, Sdbd plasma enhanced aerodynamics: concepts, optimization and applications, Progress in Aerospace Sciences 43 (7-8) (2007) 193-217.

[17] M. Forte, J. Jolibois, J. Pons, E. Moreau, G. Touchard, M. Cazalens, Optimization of a dielectric barrier discharge actuator by stationary and non-stationary measurements of the induced flow velocity: application to airflow control, Experiments in fluids 43 (6) (2007) 917-928.

[18] F. O. Thomas, T. C. Corke, M. Iqbal, A. Kozlov, D. Schatzman, Optimization of dielectric barrier discharge plasma actuators for active aerodynamic flow control, AIAA journal 47 (9) (2009) 2169-2178.

[19] D. Opaits, M. Shneider, R. B. Miles, A. Likhanskii, S. Macheret, Surface charge in dielectric barrier discharge plasma actuators, Physics of Plasmas 15 (7) (2008) 073505.
[20] A. Cristofolini, C. A. Borghi, G. Neretti, Charge distribution on the surface of a dielectric barrier discharge actuator for the fluid-dynamic control, Journal of Applied Physics 113 (14) (2013) 143307.

[21] D. Hong, H. Rabat, Y.-K. Pu, A. Leroy, Measurement of the surface charging of a plasma actuator using surface dbd, Journal of Electrostatics 71 (3) (2013) 547-550.

[22] J. Boeuf, Y. Lagmich, T. Unfer, T. Callegari, L. Pitchford, Electrohydrodynamic force in dielectric barrier discharge plasma actuators, Journal of Physics D: Applied Physics 40 (3) (2007) 652.

[23] J. Matéo-Vélez, P. Degond, F. Rogier, A. Séraudie, F. Thivet, Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment, Journal of Physics D: Applied Physics 41 (3) (2008) 035205.

[24] P. Seimandi, G. Dufour, F. Rogier, An asymptotic model for steady wire-to-wire corona discharges, Mathematical and computer modelling 50 (3-4) (2009) 574-583.

[25] M. Mamunuru, T. Simon, D. Ernie, U. Kortshagen, Plasma actuator simulation: Force contours and dielectric charging characteristics, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2010, p. 1221.

[26] Y. Lagmich, T. Callegari, L. Pitchford, J. Boeuf, Model description of surface dielectric barrier discharges for flow control, Journal of Physics D: Applied Physics 41 (9) (2008) 095205.

[27] J. Boeuf, Y. Lagmich, L. Pitchford, Contribution of positive and negative ions to the electrohydrodynamic force in a dielectric barrier discharge plasma actuator operating in air, Journal of applied physics 106 (2) (2009) 023115.

[28] T. Unfer, J. Boeuf, Modelling of a nanosecond surface discharge actuator, Journal of physics D: applied physics 42 (19) (2009) 194017.

[29] J.-P. Boeuf, T. Unfer, L. Pitchford, Studies of the electrohydrodynamic force produced in a dielectric barrier discharge for flow control. report
no. 2, phase 3, Tech. rep., Université Paul Sabatier Toulouse (France)
(2010).

[30] H. Nishida, T. Nonomura, T. Abe, Characterization of electrohydrodynamic force on dielectric-barrier-discharge plasma actuator using fluid simulation, International Journal of Aerospace and Mechanical Engineering 6 (11) (2012) 2324-2328.

[31] R. Morrow, The theory of positive glow corona, Journal of Physics D: Applied Physics 30 (22) (1997) 3099.

[32] K. Adamiak, P. Atten, Simulation of corona discharge in point-plane configuration, Journal of electrostatics 61 (2) (2004) 85-98.

[33] L. Zhao, K. Adamiak, Ehd flow in air produced by electric corona discharge in pin-plate configuration, Journal of electrostatics 63 (3-4) (2005) 337-350.

[34] A. A. Martins, Simulation of a wire-cylinder-plate positive corona discharge in nitrogen gas at atmospheric pressure, Physics of Plasmas 19 (6) (2012) 063502.

[35] S. Wang, J. Qu, L. Kong, J. Zhang, Z. Qu, Numerical and experimental study of heat-transfer characteristics of needle-to-ring-type ionic wind generator for heated-plate cooling, International Journal of Thermal Sciences 139 (2019) 176-185.

[36] Y. Guan, R. S. Vaddi, A. Aliseda, I. Novosselov, Experimental and numerical investigation of electrohydrodynamic flow in a point-to-ring corona discharge, Physical Review Fluids 3 (4) (2018) 043701.

[37] S. Chen, J. Nobelen, S. Nijdam, A self-consistent model of ionic wind generation by negative corona discharges in air with experimental validation, Plasma Sources Science and Technology 26 (9) (2017) 095005.

[38] Y. P. Raizer, J. E. Allen, Gas discharge physics, Vol. 1, Springer, 1991.

[39] P. L. Ventzek, T. J. Sommerer, R. J. Hoekstra, M. J. Kushner, Two-dimensional hybrid model of inductively coupled plasma sources for etching, Applied physics letters 63 (5) (1993) 605-607.
[40] G. Hagelaar, G. Kroesen, Speeding up fluid models for gas discharges by implicit treatment of the electron energy source term, Journal of Computational Physics 159 (1) (2000) 1-12.

[41] G. Dufour, F. Rogier, Numerical modeling of dielectric barrier discharge based plasma actuators for flow control: the copaier/cedre example, Aerospace Lab (10) (2015).

[42] P. Belenguer, J. Boeuf, Transition between different regimes of rf glow discharges, Physical Review A 41 (8) (1990) 4447.

[43] B. Parent, M. N. Shneider, S. O. Macheret, Sheath governing equations in computational weakly-ionized plasmadynamics, Journal of Computational Physics 232 (1) (2013) 234-251.

[44] A. Villa, L. Barbieri, M. Gondola, R. Malgesini, An asymptotic preserving scheme for the streamer simulation, Journal of Computational Physics 242 (2013) 86-102.

[45] T. Unfer, J.-P. Boeuf, F. Rogier, F. Thivet, An asynchronous scheme with local time stepping for multi-scale transport problems: Application to gas discharges, Journal of Computational Physics 227 (2) (2007) 898-918.

[46] I. Kossyi, A. Y. Kostinsky, A. Matveyev, V. Silakov, Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures, Plasma Sources Science and Technology 1 (3) (1992) 207.

[47] J. Teunissen, U. Ebert, 3d pic-mcc simulations of discharge inception around a sharp anode in nitrogen/oxygen mixtures, Plasma Sources Science and Technology 25 (4) (2016) 044005.

[48] B. Bagheri, J. Teunissen, U. Ebert, M. M. Becker, S. Chen, O. Ducasse, O. Eichwald, D. Loffhagen, A. Luque, D. Mihailova, et al., Comparison of six simulation codes for positive streamers in air, Plasma Sources Science and Technology 27 (9) (2018) 095002.

[49] F. Rogier, G. Dufour, K. Kourtzanidis, Numerical simulation of sinusoidal dbd actuators and comparison with experiments, in: 45th AIAA Plasmadynamics and Lasers Conference, 2014, p. 2808.
[50] K. Kourtzanidis, G. Dufour, F. Rogier, Self-consistent modeling of a surface ac dielectric barrier discharge actuator: in-depth analysis of positive and negative phases, Journal of Physics D: Applied Physics 54 (4) (2020) 045203.

[51] K. Kourtzanidis, G. Dufour, F. Rogier, The electrohydrodynamic force distribution in surface ac dielectric barrier discharge actuators: do streamers dictate the ionic wind profiles?, Journal of Physics D: Applied Physics 54 (26) (2021) 26LT01.

[52] H. Borradaile, K. Kourtzanidis, F. Rogier, K.-S. Choi, X. Mao, Flow reversal in millimetric annular dbd plasma actuator, Journal of Physics D: Applied Physics 54 (34) (2021) 345202.

[53] A. Loehrmann, J. Rivenc, G. Peres, G. Dufour, F. Rogier, K. Kourtzandis, Numerical and experimental correlation on electro-hydrodynamic system for small aircraft propulsion, in: 2021 AIAA/IEEE Electric Aircraft Technologies Symposium (EATS), IEEE, 2021, pp. 1-12.

[54] G. J. M. Hagelaar, Modeling of microdischarges for display technology, Vol. 109, Technische Universiteit Eindhoven Eindhoven, 2000.

[55] S. Mahadevan, L. L. Raja, Simulations of direct-current air glow discharge at pressures 1 torr: Discharge model validation, Journal of Applied Physics 107 (9) (2010) 093304.

[56] G. Hagelaar, L. C. Pitchford, Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models, Plasma sources science and technology 14 (4) (2005) 722.

[57] J.-P. Aubin, A. Cellina, Differential inclusions: set-valued maps and viability theory, Vol. 264, Springer Science & Business Media, 2012.

[58] H. Brezis, H. Brézis, Functional analysis, Sobolev spaces and partial differential equations, Vol. 2, Springer, 2011.

[59] V. Barbu, Nonlinear differential equations of monotone types in Banach spaces, Springer Science & Business Media, 2010.

[60] H. Brezis, Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert, Elsevier, 1973.
[61] J.-L. Lions, Quelques méthodes de résolution des problemes aux limites non linéaires, Dunod, 1969.

[62] H. Brezis, Inéquations variationnelles paraboliques, Séminaire Jean Leray (1971) 1–10.

[63] K. Ito, K. Kunisch, Parabolic variational inequalities: The lagrange multiplier approach, Journal de mathématiques pures et appliquées 85 (3) (2006) 415–449.

[64] V. Gorin, A. Kudryavtsev, J. Yao, C. Yuan, Z. Zhou, Boundary conditions for drift-diffusion equations in gas-discharge plasmas, Physics of Plasmas 27 (1) (2020) 013505.

[65] S. Pancheshnyi, Role of electronegative gas admixtures in streamer start, propagation and branching phenomena, Plasma Sources Science and Technology 14 (4) (2005) 645.

[66] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities, International journal for numerical methods in engineering 79 (11) (2009) 1309–1331.

[67] D. L. Scharfetter, H. K. Gummel, Large-signal analysis of a silicon read diode oscillator, IEEE Transactions on electron devices 16 (1) (1969) 64–77.

[68] T. D. Nguyen, C. Besse, F. Rogier, High-order scharfetter-gummel-based schemes and applications to gas discharge modeling, Journal of Computational Physics 461 (2022) 111196.

[69] D. Breden, K. Miki, L. Raja, Self-consistent two-dimensional modeling of cold atmospheric-pressure plasma jets/bulletts, Plasma Sources Science and Technology 21 (3) (2012) 034011.

[70] P.-L. Lions, B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM Journal on Numerical Analysis 16 (6) (1979) 964–979.

[71] N. Parikh, S. Boyd, et al., Proximal algorithms, Foundations and trends® in Optimization 1 (3) (2014) 127–239.
[72] J. Eckstein, D. P. Bertsekas, On the douglas—rachford splitting method and the proximal point algorithm for maximal monotone operators, Mathematical Programming 55 (1) (1992) 293–318.

[73] A. Jameson, Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings, in: 10th Computational fluid dynamics conference, 1991, p. 1596.

[74] S.-Z. Li, H. S. Uhm, Investigation of electrical breakdown characteristics in the electrodes of cylindrical geometry, Physics of Plasmas 11 (6) (2004) 3088–3095.