Tuan Nguyen

Christophe Dung

François Besse

Rogier

Christophe Besse

Nguyen Tuan Dung

François Rogier

An implicit time integration approach for simulation of corona discharges

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

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.

Introduction

In the past two decades plasma actuators have been extensively studied for their applications in airflow control [START_REF] Moreau | Airflow control by non-thermal plasma actuators[END_REF][START_REF] Choi | Plasma virtual roughness elements for cross-flow instability control[END_REF], anti-icing [START_REF] Zhu | Nanosecond-pulsed dielectric barrier discharge-based plasma-assisted anti-icing: modeling and mechanism analysis[END_REF] and flight sustaining [START_REF] Xu | Flight of an aeroplane with solid-state propulsion[END_REF], among many others. These devices are robust, light, low energy consuming and easy to install, thus they have a potential interest for industrial 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 colliding with neutral particles in air, creating an electro-hydrodynamic (EHD) force. A lot of experimental efforts have been made to characterize the EHD force [START_REF] Van Dyken | Parametric investigations of a single dielectric barrier plasma actuator[END_REF][START_REF] Enloe | Mechanisms and responses of a single dielectric barrier plasma actuator: plasma morphology[END_REF][START_REF] Moreau | Airflow control by non-thermal plasma actuators[END_REF][START_REF] Bérard | Measurements and simulations of the ionic wind produced by a dc corona discharge in air, helium and argon[END_REF][START_REF] Porter | Plasma actuator force measurements[END_REF][START_REF] Enloe | Time-correlated force production measurements of the dielectric barrier discharge plasma aerodynamic actuator[END_REF][START_REF] Hoskinson | Force measurements of single and double barrier dbd plasma actuators in quiescent air[END_REF][START_REF] Opaits | Surface plasma induced wall jets[END_REF][START_REF] Debien | Unsteady aspect of the electrohydrodynamic force produced by surface dielectric barrier discharge actuators[END_REF][START_REF] Zhang | Characteristics of ionic wind in needle-to-ring corona discharge[END_REF], to optimize the actuator geometry [START_REF] Roth | Optimization of the aerodynamic plasma actuator as an electrohydrodynamic (ehd) electrical device[END_REF][START_REF] Corke | Plasma flow control optimized airfoil[END_REF][START_REF] Corke | Sdbd plasma enhanced aerodynamics: concepts, optimization and applications[END_REF][START_REF] Forte | Optimization of a dielectric barrier discharge actuator by stationary and non-stationary measurements of the induced flow velocity: application to airflow control[END_REF][START_REF] Thomas | Optimization of dielectric barrier discharge plasma actuators for active aerodynamic flow control[END_REF] or to measure the dielectric surface charge [START_REF] Opaits | Surface charge in dielectric barrier discharge plasma actuators[END_REF][START_REF] Cristofolini | Charge distribution on the surface of a dielectric barrier discharge actuator for the fluid-dynamic control[END_REF][START_REF] Hong | Measurement of the surface charging of a plasma actuator using surface dbd[END_REF]. 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 [START_REF] Boeuf | Electrohydrodynamic force in dielectric barrier discharge plasma actuators[END_REF][START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF][START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF][START_REF] Mamunuru | Plasma actuator simulation: Force contours and dielectric charging characteristics[END_REF][START_REF] Lagmich | Model description of surface dielectric barrier discharges for flow control[END_REF][START_REF] Boeuf | Contribution of positive and negative ions to the electrohydrodynamic force in a dielectric barrier discharge plasma actuator operating in air[END_REF][START_REF] Unfer | Modelling of a nanosecond surface discharge actuator[END_REF][START_REF] Boeuf | Studies of the electrohydrodynamic force produced in a dielectric barrier discharge for flow control[END_REF][START_REF] Nishida | Characterization of electrohydrodynamic force on dielectric-barrier-discharge plasma actuator using fluid simulation[END_REF].

For simulation of corona discharges, many fluid models have been proposed throughout the years. Morrow [START_REF] Morrow | The theory of positive glow corona[END_REF] studied the discharge for concentricsphere electrodes with a set of evolution continuity equations for five species (electrons, positive ions, negative ions, neutrals) coupled with a Poisson equation for the electrostatic potential. This approach took into account also the electron source from photoionization but was not validated experimentally. Adamiak & Atten [START_REF] Adamiak | Simulation of corona discharge in point-plane configuration[END_REF] 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 [START_REF] Adamiak | Simulation of corona discharge in point-plane configuration[END_REF], needle-toplane [START_REF] Zhao | Ehd flow in air produced by electric corona discharge in pin-plate configuration[END_REF], wire-cylinder-plate [START_REF] Martins | Simulation of a wire-cylinder-plate positive corona discharge in nitrogen gas at atmospheric pressure[END_REF] or needle-to-ring [START_REF] Wang | Numerical and experimental study of heat-transfer characteristics of needle-to-ring-type ionic wind generator for heated-plate cooling[END_REF]. A similar model, but evolutionary, was proposed by Guan et al. [START_REF] Guan | Experimental and numerical investigation of electrohydrodynamic flow in a point-to-ring corona discharge[END_REF] but the boundary condition for the charge equation depends on the experimental electric current. Chen et al. [START_REF] Chen | A self-consistent model of ionic wind generation by negative corona discharges in air with experimental validation[END_REF] 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 negative corona discharges [START_REF] Raizer | Gas discharge physics[END_REF].

Simulation of corona discharges proves to be very challenging in terms of waiting time [START_REF] Morrow | The theory of positive glow corona[END_REF] 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 electron motion is only on picoseconds. Many numerical efforts have been made to reduce the simulation time. Ventzek et al. [START_REF] Ventzek | Twodimensional hybrid model of inductively coupled plasma sources for etching[END_REF] proposed a semi-implicit scheme which allowed to compute the field implicitly and separately from the specie densities. Hagelaar [START_REF] Hagelaar | Speeding up fluid models for gas discharges by implicit treatment of the electron energy source term[END_REF] proposed a scheme in the same category in the presence of the electron energy equation. Seimandi et al. [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF] studied some asymptotic models to reduce the problem dimension in the ionization layers near the electrodes. Dufour & Rogier developed a sub-cycling algorithm in their multiscale plasma solver COPAIER [START_REF] Dufour | Numerical modeling of dielectric barrier discharge based plasma actuators for flow control: the copaier/cedre example[END_REF]. The idea was to repeatedly "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 [START_REF] Belenguer | Transition between different regimes of rf glow discharges[END_REF][START_REF] Parent | Sheath governing equations in computational weakly-ionized plasmadynamics[END_REF][START_REF] Villa | An asymptotic preserving scheme for the streamer simulation[END_REF] or asynchronous time stepping [START_REF] Unfer | An asynchronous scheme with local time stepping for multi-scale transport problems: Application to gas discharges[END_REF], 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 [START_REF] Kossyi | Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures[END_REF]. 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 [START_REF] Teunissen | 3d pic-mcc simulations of discharge inception around a sharp anode in nitrogen/oxygen mixtures[END_REF][START_REF] Bagheri | Comparison of six simulation codes for positive streamers in air[END_REF], but not for the glow regime in which we are interested [START_REF] Raizer | Gas discharge physics[END_REF] 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 [START_REF] Dufour | Numerical modeling of dielectric barrier discharge based plasma actuators for flow control: the copaier/cedre example[END_REF] with explicit time integration. Simulations with floor density have been contributing to the study of plasma actuators [START_REF] Rogier | Numerical simulation of sinusoidal dbd actuators and comparison with experiments[END_REF][START_REF] Kourtzanidis | Self-consistent modeling of a surface ac dielectric barrier discharge actuator: in-depth analysis of positive and negative phases[END_REF][START_REF] Kourtzanidis | The electrohydrodynamic force distribution in surface ac dielectric barrier discharge actuators: do streamers dictate the ionic wind profiles?[END_REF][START_REF] Borradaile | Flow reversal in millimetric annular dbd plasma actuator[END_REF]. Recently, simulation results from [START_REF] Loehrmann | Numerical and experimental correlation on electro-hydro dynamic system for small aircraft propulsion[END_REF] 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 ψ on the charged particle density n can be seen as a constraint: n ≥ ψ, 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 investigates 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 [START_REF] Bérard | Measurements and simulations of the ionic wind produced by a dc corona discharge in air, helium and argon[END_REF]. The studied case is a corona discharge generated by a wire-to-wire actuator and the simulations are conducted in two dimensions with COPAIER.

Modeling of gas discharges with floor density

Simplified fluid model

We use the Local Field Approximation (LFA) model [START_REF] Hagelaar | Modeling of microdischarges for display technology[END_REF] 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 Ω ⊂ R d (d = 1, 2, 3) (see fig. 1) and evolves on a time interval (0, T ). Let us denote as L the set of particle species that are considered in the discharge and as n i (t, x) (measured in m -3 ) the particle density of the specie i ∈ L. The densities are grouped in a vector denoted as U ≡ (n i ) i∈L . The electrostatic field is denoted as E(t, x) (V m -1 ).

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

∂ t U + ∇ • F (E, U ) = S(E/N, U ), (t, x) ∈ (0, T ) × Ω, (1)
where

F ≡ (f i ) i∈L , S ≡ (S i ) i∈L , E ≡ |E| and N (m -3
) is the air density. More precisely, f i (E, U ) is the particle flux density of specie i and S i (E/N, U ) is the creation/decay rate per unit volume of i via plasmachemical processes. For x ∈ Ω,

f i = u i n i -D i ∇n i where u i (E) (ms -1 ) is the drift velocity of i and D i (E/N ) (m 2 s -1
) is its diffusion rate. The charged species are subject to a mobility law which reads as

u i = sign(z i )µ i E,
where z i is the charge number of i and µ i (E/N ) (m 2 V -1 s -1 ) is its mobility coefficient. The field is derived from the electric potential φ as E = -∇φ.

The electric potential φ(t, x) (V ) is given by Gauss' law,

-∇ • (ε∇φ) = ρ, (2)
where ε(x) (CV -1 m -1 ) is the electric permittivity of the domain Ω, ρ ≡ q i∈L z i n i (Cm -3 ) is the space charge density and q ≈ 1.6 × 10 -19 C is the elementary charge. We assume for simplicity that n i , φ, E, S i are smooth functions defined on (0, T ) × Ω.

x Ω Γ

x 1

x 2 The number of species and physicochemical reactions is generally very numerous [START_REF] Kossyi | Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures[END_REF][START_REF] Mahadevan | Simulations of direct-current air glow discharge at pressures 1 torr: Discharge model validation[END_REF]. In order to limit the complexity of the system, we adopt a simplified model from J.-P. Boeuf et al.'s work [START_REF] Boeuf | Contribution of positive and negative ions to the electrohydrodynamic force in a dielectric barrier discharge plasma actuator operating in air[END_REF] 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 → 2e + p, electron attachment e + N → n, electron-ion recombination e + p → N and ion-ion recombination p + n → N . In the following, we suppose that the neutral particle density, still denoted as N , is constant and N = 2.5 × 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, L = { e, p, n }. The specie source terms are written as,

S e = (α -η)N n e -k ep n e n p , S p = αN n e -k ep n e n p -k np n n n p , S n = ηN n e -k np n n n p ,
where α, η, k ep and k np (m 3 s -1 ) are resp. the ionization, electron attachment, electron-ion recombination and ion-ion recombination coefficients, which are positive. In the LFA model, these coefficients depend only on the reduced field strength E/N . In the simulations, α, η are extracted from look-up tables provided by the BOLSIG+ application [START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF]; k ep = 2 × 10 -13 and k np = 1.7 × 10 -13 . As discussed in section 1, we introduce a positive and smooth function on Ω, ψ(x) ≥ 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, x) ≥ ψ(x), ∀(t, x) ∈ (0, T ) × Ω. ( 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 K as the closed convex set of smooth functions defined on Ω that satisfy eq. ( 3), i.e.

K = { v | v ≥ ψ } 2 .
The characteristic function of K is defined as

I K (v) = 0 if v ∈ K, +∞ if v / ∈ K.
The subdifferential [57, Proposition 3] of I K at v is a set of functions defined as

∂I K (v) =      ∅ if v / ∈ K, { 0 } if v > ψ, { u | ∀w ∈ K, (u, w -v) ≤ 0 } otherwise, (4)
where (•, •) is a scalar product defined on a Hilbert space H such that for every t ∈ (0, T ), n i (t, .) ∈ H, i ∈ L. For example, H = L2 (Ω) and (u, v) = Ω u(x)v(x)dx.

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, •) ∈ K ∀t ∈ (0, T ), then we would have from eq. ( 1) that S e -∂ t n e -∇ • f e + n e = n e ≥ ψ, ∀(t, x) ∈ (0, T ) × Ω, 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 [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF].

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

|v -u| = min w∈K |v -w|, where | • | is the norm induced from (•, •). Furthermore, u is characterized by the property (v -u, w -u) ≤ 0, ∀w ∈ K.
Thus, we have (S e -∂ t n e -∇ • f e + n e -n e , w -n e ) ≤ 0 for all w ∈ K and t ∈ (0, T ). This and eq. ( 4) imply that S e (t, •) - •) is an element of the set ∂I K n e (t, •) for all t ∈ (0, T ). Let us denote S ψ (t) ≡ -∂I K n e (t, •) . The electron conservation law which integrates the floor density ψ now transforms into a differential inclusion [START_REF] Barbu | Nonlinear differential equations of monotone types in Banach spaces[END_REF] which reads as

∂ t n e (t, •) -∇ • f e (t,
∂ t n e (t, •) + ∇ • f e (t, •) ∈ S e (t, •) + S ψ (t) 3 , ∀t ∈ (0, T ). (5)
Finally, the simplified LFA model reads as follows,

               ∂ t U + ∇ • F ∈ S +    S ψ 0 0    , U (0, •) = U 0 , ∇ • (εE) = ρ, (6)
where U 0 = (n 0 e n 0 p n 0 n ) t and n 0 i are smooth functions defined on Ω. We remark that the above system is not charge conservative, i.e. ∂ t ρ + ∇• j = 0 4 . However, since the floor density ψ that we choose in the simulations 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 ∂ t n e + ∇ • f e in eq. ( 5) is not arbitrarily chosen (we could have otherwise S e ∈ ∂ t n e + ∇ • f e + S ψ ). The existence and uniqueness of solution of the nonlinear problem eq. ( 5) are subject of classical studies based on the theory of monotone operators [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] and variational inequalities [START_REF] Lions | Quelques méthodes de résolution des problemes aux limites non linéaires[END_REF][START_REF] Brezis | Inéquations variationnelles paraboliques[END_REF][START_REF] Ito | Parabolic variational inequalities: The lagrange multiplier approach[END_REF].

Boundary conditions

The boundary conditions for the conservation laws are imposed on the fluxes f i . Different boundary conditions exist and each is defined on a nonoverlapping boundary portion Γ k . The union of these portions is the domain boundary Γ: Γ = ∪ k Γ k (see fig. 10 for example).

1. For free-flow boundaries Γ 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.

(E • ν) |Γ f = 0, f i|Γ f ≡ (f i • ν) |Γ f = (n i u i • ν) |Γ f = 0, i ∈ L, (7)
where ν(x) is the unit outward normal of Ω at x ∈ Γ. 2. For wall boundaries Γ w , which are physical interfaces between the discharge domain and the electrodes 5 , we determine the specie fluxes using the two-stream approach [START_REF] Gorin | Boundary conditions for drift-diffusion equations in gas-discharge plasmas[END_REF]. 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

φ |Γw = φ 0 , f i|Γw = n i max (u i • ν, 0) |Γw , i ∈ L \ { e } , f e|Γw = 1 2 8k B T e πm e 1 2
n e|Γw -2γf p|Γw .

(8) 5 we do not consider interfaces with dielectric material: they do not appear in the discharges studied in this article where φ 0 is the electrode voltage, k B is the Boltzmann constant, T e is the electron temperature, m e is the electron mass and γ is the secondary emission coefficient. We set T e = 2eV and γ = 10 -4 . 3. For symmetric boundaries Γ s , which are fictive interfaces between symmetric discharge domains, the boundary conditions are

(E • ν) |Γs = 0, f i|Γs = 0, i ∈ L. (9)

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 [START_REF] Pancheshnyi | Role of electronegative gas admixtures in streamer start, propagation and branching phenomena[END_REF]. 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 ψ can also be envisioned to include the naturally formed electrons.

Numerical methods

Key features

A complete description of the COPAIER solver is found in [START_REF] Dufour | Numerical modeling of dielectric barrier discharge based plasma actuators for flow control: the copaier/cedre example[END_REF], 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 [START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF] 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 [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF].

The cell-centered finite volume approach is employed to discretized the continuity equations/inclusions. The Scharfetter-Gummel scheme [START_REF] Scharfetter | Large-signal analysis of a silicon read diode oscillator[END_REF] as well as the SGCC scheme [START_REF] Nguyen | High-order scharfetter-gummel-based schemes and applications to gas discharge modeling[END_REF] 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 ∆x and at the time level t l , the time step ∆t l ≡ t l+1 -t l must satisfy6 ∆t l ≤ min

x ∆x |u e E(t l , x) | tanh |u e E(t l , x) |∆x 2D e E(t l , x) ≡ ∆t l e , ∆t l ≤ min i,x ∆x |u i | tanh |u i |∆x 2D i ≡ ∆t l ions , i ∈ L \ { e } , ∆t l ≤ min x ε(x) σ(t l , x) ≡ ∆t l φ .
In corona discharges, the dielectric relaxation time ∆t l φ is much larger than the CFL condition ∆t l e (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.

Discrete plasma model

For l ≥ 0, let U l (x) be an approximation of U (t l , x) at discrete time t l (and similarly for φ, u i , D i , α, η, k ep , k np ) and ∂ t U (t l+1 ) be approximated by U l+1 -U l ∆t l , we consider the following time discretization of eq. ( 6),

         U l+1 -U l ∆t l + ∇ • F (E l , U l+1 ) ∈ S(E l /N, Ũ ) +    S l+1 ψ 0 0    , ∇ • (ε∇E l+1 ) = ρ l+1 , (10)
with

S(E l /N, Ũ ) =   (α l -η l )N n l+1 e -k l ep n l+1 e n l p + S l+1 ψ α l N n l+1 e -k l ep n l+1 e n l p -k l np n l n n l p η l N n l+1 e -k l np n l n n l p   and S l ψ ≡ S ψ (t l ) = -∂ K n e (t l , •)
. We remark that the coefficients u e , D i , α, η, k ep , k np are computed from the potential φ 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 [START_REF] Breden | Self-consistent two-dimensional modeling of cold atmospheric-pressure plasma jets/bullets[END_REF] except the nonlinearity only stemming from S ψ .

Let us denote as I the identity operator,

A l : U →   ∇ • f e (E l , U ) -(α l -η l )N n e + k l ep n e n l p ∇ • f p (E l , U ) -α l N n e + k l ep n e n l p ∇ • f n (E l , U ) -η l N n e + k l np n l n n l p   ≡   A l U e A l U p A l U n   , B ψ : U →   -S ψ 0 0   and R l =   n l e n l p -∆t l k l np n l n n l p n l n -∆t l k l np n l n n l p   ≡   R l e R l p R l n 
 . Then at t l , the semi-discrete (in time) conservation laws of eq. ( 10) can be rewritten as

(I + ∆t l A l + ∆t l B ψ )U l+1 R l . ( 11
) K x K x K,λ ν K λ Ω Figure 2: A triangulation T of a polygonal domain Ω in R 2
Furthermore, let us consider a polygonal domain Ω in R 2 and a mesh T consisting of closed polygonal convex subsets of Ω (see fig. 2). For an element K ∈ T , we denote as E K the set of its edges. Let us consider the electron conservation law of eq. [START_REF] Hoskinson | Force measurements of single and double barrier dbd plasma actuators in quiescent air[END_REF],

n l+1 e -n l e ∆t l + ∇ • f e (E l , U l+1 ) -(α l -η l )N n l+1 e + k l ep n l+1 e n l p = h,
where h is an element of S l+1 ψ . Now integrating this equation on K and dividing it by |K|, we have

n l+1 e,K -n l e,K ∆t l + 1 |K| λ∈E K λ f e • ν K dl -(α l -η l )N n l+1 e K + k l ep n l+1 e n l p K = h K , where v K ≡ 1 |K| K v(x)
dx for any function v, dl is the unit length element and ν K (x) is the unit outward normal of K at x ∈ λ ∈ E K . We approximate α l (x) for x ∈ K (and similarly for η, k ep ) by

α l (x K ) ≡ α l K with x K the center of K, λ f e • ν K dl by |λ|f e,K,λ • ν K,λ ≡ |λ|f e (x K,λ ) • ν K (x K,λ ) with
x K,λ the center of λ ∈ E K and then n e n p by n e n p . Finally we have

n l+1 e,K -n l e,K ∆t l + 1 |K| λ∈E K |λ|f e,K,λ • ν K,λ -(α l K -η l K )N n l+1 e,K + k l ep,K n l+1 e,K n l p,K ≈ h K , (12)
for every K ∈ T . Now for any smooth function v defined on Ω, we define a piecewise constant function v on Ω such that

v(x) = v K , ∀x ∈ K \ ∪ λ∈E K λ.
It is not evident how h with h ∈ S l+1 ψ relates to ψ and n e . Therefore, we make an assumption that

h ∈ S ψ (n l e ) ≡ -∂I K (n l e ), K = v v ≥ ψ . (13)
Finally, using the Scharfetter-Gummel scheme 7 [67] to approximate the fluxes f i,K,λ • ν K,λ , i ∈ L yields a linear relation of (n i,K ) K∈T . The left-hand 7 see Appendix A side of eq. ( 12), combining with similar equations for other species, can be written in the same way of eq. ( 11) as follows,

(I + ∆t l A l + ∆t l B ψ )U l+1 R l , (14)
where

U l ≡ (n l e,K n l p,K n l n,K ) t K∈T ∈ R |L||T | , I is the identity matrix, A l is a |L||T | × |L||T | matrix which is the discrete counterpart of A l and R l is
the discrete counterpart of R l . We remark that writing B ψ U l+1 is a notation abuse since 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).

Solving the conservation laws

The complex structure of the operator I + ∆t l A l + ∆t l B ψ makes the evaluation of (I + ∆tA + ∆tB) -1 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 [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] 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 [START_REF] Parikh | Proximal algorithms[END_REF].

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 ∂I K is maximal monotone on the Hilbert H. We recall that a (multivalued) operator A is monotone if

(y -u, w -v) ≥ 0, ∀v, w ∈ D(A) ⊂ H, u ∈ A(v), y ∈ A(w),
where D(A) is the domain of A; and is maximal monotone if its resolvent J ∆t A ≡ (I + ∆tA) -1 is a contraction defined on H [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]. In our case, if (I + ∆t l ∂I K )w y then from the definition eq. ( 4) we have,

y -w ∆t l ∈ ∂I K (n) ⇐⇒ ∀v ∈ K, (y -w, v -w) ≤ 0.
It is shown that w is the projection of y on K, i.e. w(x) = max(y(x), ψ(x)), which is unique according to theorem 1 and easily shown to be a contraction. Therefore, if Y = (y e y p y n ) t then

J ∆t l B ψ Y ≡ (I + ∆t l B ψ ) -1 Y = max(y e , ψ) y p y n t . (15)
Hence, from eq. ( 11) the electron density at t l+1 is n l+1 e = max(R l e -∆t l (A l U l+1 ) e , ψ) which implies that n l+1 e = max(R l e -∆t l (A l U l+1 ) e , ψ). By this point, making the assumption eq. ( 13) is in fact to approximate the righthand side with max R l e -∆t l (A l U l+1 ) e , ψ . Further approximations of flux and chemical coefficients lead finally to

n l+1 e = max R l e -∆t l (A l U l+1 ) e , ψ ,
which is the electron density solution of eq. ( 14).

Another remark is that J ∆t l B ψ is independent of ∆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 B ψ instead of J ∆t l B ψ .

Lie splitting

The simplest operator splitting method is the Lie algorithm which is firstorder in time. At t l , we compute the numerical solution, still denoted as U l+1 , by

U l+1 = J B ψ J ∆t l A l R l . (16)
The evaluation of J ∆t l A l R l is classical since it is equivalent to inverting a linear system.

Douglas-Rachford algorithm

The second method that we consider is the Douglas-Rachford algorithm [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]. The original method was introduced as a first-order-in-time operator splitting to solve the evolution problem

∂ τ U + (A + B)U R, U (τ = 0) = U 0 ∈ D(A) ∩ D(B), (17)
where A and B are maximal monotone operators. The algorithm is described as follows: we choose V 0 ∈ (I + ∆τ B)U 0 and for m ≥ 0 define

V m+1 = J ∆τ A (2J ∆τ B -I)V m + (I -J ∆τ B )V m + ∆τ J ∆τ A R.
The approximation of U at time m∆τ is then U m = J ∆τ B V m . A crucial property of this method is that as m → ∞, U m converges to the steady-state solution U ∞ of eq. ( 17) (i.e. (A + B)U ∞ 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 [START_REF] Jameson | Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings[END_REF] problem,

   ∂ τ U (τ ) + I ∆t l + A l + B ψ U (τ ) 1 ∆t l R l , U (0) = U l ,
where τ is the pseudotime. We apply the Douglas-Rachford algorithm for

A = I ∆t l + A l , B = B ψ , R = 1 ∆t l R l , ∆τ = ∆t l ,
and when it converges, we put U l+1 = U ∞ . The requirement that A is monotone imposes a constraint on the time step ∆t l which is related to the Townsend ionization rate α. For simplicity, let us assume that D i , u i , α are constant and other reaction coefficients are null. We define the scalar product on the space H ≡ H × H × H as (U , V ) H = (n e , w e ) + (n p , w p ) + (n n , w n ), with U = (n e n p n n ) t , V = (w e w p w n ) t . Since A is a linear operator it is sufficient to check the sign of (AU , U ) H . By some algebraic manipulation we have

(AU , U ) H = 1 ∆t l i∈L |n i | 2 + i∈L D i |∇n i | 2 + 1 2 i∈L Γ (u i • ν)n 2 i dl -αN |n e |(1 + |n p |).
For (AU , U ) H ≥ 0, ∆t l must satisfy

∆t l ≤ αN |n e |(1 + |n p |) - i∈L D i |∇n i | 2 - 1 2 i∈L Γ (u i • ν)n 2 i dl -1 i∈L |n i | 2 . (18)
In case of very large α 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.

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 ) , m ≥ 0, Algorithm 1 Douglas-Rachford 1: Let > 0, U 0 = U l , V 0 ∈ (I + ∆t l B ψ )U 0 , δ > 2: while δ > do 3: Set V m+1 = J 1 I+∆t l A l (2J B ψ -I)V m + (I -J B ψ )V m + J 1 I+∆t l A l R l 4: Set U m+1 = J B ψ V m+1 and let n (m+1) e n (m+1) p n (m+1) n t = U m+1 5: Update δ = max i∈L |n (m+1) i -n (m) i | ∞ |n (m) i | ∞ where | • | ∞ is the L ∞ norm defined on H 6: end while 7: Set 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 m+1 p = 1 1 + A pp R p - q>p A pq U m q - q<p A pq U m+1 q , p = 1, 2, . . . ,
where U m p 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 = (I + ∆t l M l + ∆t l D l + ∆t l B ψ ) -1 R l -∆t l N l U m ,
where D l , M l , N l are resp. the diagonal, strictly lower and upper triangular components of A l . Its element-wise form reads as

I + ∆t l 1 + ∆t l A l pp B ψ U m+1 p = 1 1 + ∆t l A l pp R l p -∆t l q>p A l pq U m q -∆t l q<p A l pq U m+1 q .
As long as 1 + ∆t l A l pp > 08 , ∀p so that

∆t l 1 + ∆t l A l pp
B ψ is monotone, the resolvent of B ψ is well defined as in eq. ( 15) and we have

U m+1 p = J B ψ 1 + ∆t l A l pp R l p -∆t l q>p A l pq U m q -∆t l q<p A l pq U m+1 q .
The complete description of the proposed algorithm is presented in algorithm 2.

Algorithm 2 Gauss-Seidel

1: Let > 0, U 0 = U l , δ > 2: while δ > do 3:
Forward substitution 4:

U m+ 1 2 = (I + ∆t l M l + ∆t l D l + ∆t l B ψ ) -1 R l -∆t l N l U m 5:
Backward substitution 6:

U m+1 = (I + ∆t l N l + ∆t l D l + ∆t l B ψ ) -1 R l -∆t l M l U m+ 1 2 7: Let n (m+1) e n (m+1) p n (m+1) n t = U m+1 8: Update δ = max i∈L |n (m+1) i -n (m) i | ∞ |n (m) i | ∞ 9: end while 10: Set U l+1 = U m

Comparison of algorithms in one-dimensional discharges

Our first test case is the wire-to-wire corona discharge proposed in [START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF]. 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 Γ 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 threedimensional 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 where V ext is the analytic expression of the electrostatic potential generated by the electrodes, which can be found in [START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF]. We set max x S = 5mm.

d(V ext S) dx = 0, S(x) r 2 r 1 V G x Γ 0
We set n 0 p = n 0 e = 10 9 m -3 accordingly to section 2.3. On the contrary, the seed negative ion density is considered negligible: n 0 n = 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,

                 ∂ t v e + ∇ • (u S e v e -D e ∇v e ) ∈ (α -η)N v e -k S ep v e v p + S Sψ , ∂ t v p + ∇ • (u p v p ) = αN v e -k S ep v e v p -k S np v n v p , ∂ t v n + ∇ • (u n v n ) = ηN v e -k S np v n v p , v i (t = 0) = Sn 0 i , i ∈ L, -∇ • (εS∇φ) = q i∈L z i v i ,

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 [START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF]. 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 specie 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 ∆x min = 5.5µm, and using the Scharfetter-Gummel scheme as the flux approximation method.

We compute the time step as ∆t = min(c∆t ions , ∆t φ ) 9 . When the simulations occur between 0 and 0.04ms, we fix c = 10 -2 since at the discharge onset, there is a rapid charge multiplication and as a result, the dielectric relaxation time ∆ φ decreases quickly and so having a small time step ∆t ensures that the plasma dynamics are captured correctly. After 0.04ms, we choose the value of c between 1 and 10 2 . After around 0.4ms, the discharge enters the ion collection phase and we have c∆t ions < ∆t φ .

Lie splitting. The computed electric currents with the Lie algorithm are shown in fig. 4. The discharges were simulated until T = 7ms 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). Gauss-Seidel algorithm. We defined the stop criterion as = 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 = 7ms, the current of ∆t = 10 2 × ∆t L differs only 8 × 10 -7 % from the current of ∆t = ∆t L .

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 term 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 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 ∆t φ can fall below c∆t ions so N (t l ) would not depend on c. To sum up, we have, 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.

N = 7ms t l =1ms N (t

Mesh convergence

In this section, we conduct a mesh convergence study to ensure that the Gauss-Seidel algorithm, combining with the SGCC flux schemes [START_REF] Nguyen | High-order scharfetter-gummel-based schemes and applications to gas discharge modeling[END_REF], is able to capture correctly the discharge dynamics. We perform the simulations on five grids with 400 cells (∆x min = 5.5µm), 800 cells (2.6µm), 1600 cells (1.3µm), 3200 cells (0.6µm) and 6400 cells (0.3µm). We use the Scharfetter-Gummel scheme and the SGCC 1 , SGCC 2 schemes [START_REF] Nguyen | High-order scharfetter-gummel-based schemes and applications to gas discharge modeling[END_REF] 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 .

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 2 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 1 , SGCC 2 ) increased significantly the numerical precision. Finally, we note that a mesh convergence study for the SGCC schemes with explicit time integration was partially conducted in [START_REF] Nguyen | High-order scharfetter-gummel-based schemes and applications to gas discharge modeling[END_REF] 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 2 scheme takes only an hour on the 3200-element grid and 4 hours on the 6400-element grid 11 .

11 on a Dell Latitude 5410

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 [START_REF] Matéo-Vélez | Modelling wire-to-wire corona discharge action on aerodynamics and comparison with experiment[END_REF]. 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. The particular geometry of the wire-to-wire actuator intensifies the electrostatic field in the electrodes vicinity [START_REF] Li | Investigation of electrical breakdown characteristics in the electrodes of cylindrical geometry[END_REF]. 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 min .

Γ s Γ s Γ s Γ w Γ w
We set n 0 p = n 0 e = 10 9 m -3 accordingly to section 2.3. On the contrary, the seed negative ion density is considered negligible: n 0 n = 1m -3 . The characteristics of the actuator and the numerical parameters are grouped in table 2. Table 2: Parameters of the wire-to-wire simulations

Steady-state solution

We present the numerical results for a case of V G = 13kV , r 1 = 0.1mm, d = 10mm, ψ = 10 10 m -3 and ∆x min = 5µ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µ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 ≡ ρ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 16kN m -3 which directs towards the cathode. On the left of the anode, the EHD reverses sign and reaches -4.6kN m -3 . Overall, the force density radiates from and concentrates near the anode.

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 ψ. In fig. 13, we show the ψ -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 [START_REF] Bérard | Measurements and simulations of the ionic wind produced by a dc corona discharge in air, helium and argon[END_REF]. 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 selfsustaining 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 ψ smaller than 5×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 ψ as suggested in fig. 13, but the dependence seems to differ in each value interval of ψ. For ψ = 10 9 -10 11 m -3 , we have the estimated law I ∼ ψ 0.07 . For ψ = 10 12 -10 14 m -3 , we have I ∼ ψ 0.7 . This suggests that ψ should be around 10 10 m -3 for this wire-to-wire discharge, since a too high ψ 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 ψ is small 10 2 (1m -3 ), suggesting that there is no discharge at all though this is not true (see section 5.3). Therefore, the inclusion of ψ is necessary.

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. [START_REF] Bérard | Measurements and simulations of the ionic wind produced by a dc corona discharge in air, helium and argon[END_REF]. The simulations are conducted for d = 10mm, 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 ψ, 10 10 m -3 (shown by filled square markers) and 10 11 m -3 (filled triangles). For r 1 = 0.1mm, the numerical currents of ψ = 10 10 m -3 fit the data very well. For r 1 = 0.175mm, the currents of ψ = 10 11 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 m -3 -curve but the curve itself deviates too much for large voltages. The situation is quite similar for r 1 = 0.2mm, although the numerical results differ even more from experimental data. For corona discharges, the V G -I characteristics is given by the law I = kV G (V G -V c ) [START_REF] Raizer | Gas discharge physics[END_REF]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 (µA(kV ) -2 ) and V c (kV ) for each curve. These characteristics are laid out in Table 3. The results show that the breakdown voltages V c of ψ = 10 10 are closer to experimental data. V c decreases with increasing ψ since the breakdown condition is favorable with large floor density. However, the proportional coefficient k of ψ = 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.

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 ≥ ψ 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. We use the following extension of the Scharfetter-Gummel scheme [67] on 13 i.e. no hanging nodes

Figure 1 :

 1 Figure 1: A two-dimensional open domain Ω with boundary Γ ≡ Ω \ Ω. The coordinates of an arbitrary point is denoted as x.

Figure 3 :

 3 Figure 3: Wire-to-wire actuator schematics

  with u S e = -µ e E + D e S /S, k S ep = k ep /S and k S np = k np /S.

Figure 4 :

 4 Figure 4: Lie splitting: electric current I with different CFL numbers c

Figure 5 :

 5 Figure 5: Douglas-Rachford algorithm: electric current I with different CFL numbers c

Figure 6 :

 6 Figure 6: Gauss-Seidel algorithm: electric current I with different CFL numbers c

Figure 7 :

 7 Figure 7: Scharfetter-Gummel flux scheme: electric current I with different mesh size ∆x min

Figure 9 :

 9 Figure 9: SGCC 2 flux scheme: electric currentI with different mesh size ∆x min

Figure 10 :

 10 Figure 10: Sketch of the computation domain

Figure 11 :

 11 Figure11: Example of grid refinement near the smaller wire generated by Gmsh[START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]

r 1 0

 1 .1 to 0.5mm V G 9 to 37kV ∆x min 5µm r

Figure 12 :

 12 Figure 12: EHD force density F x (N m -3 ) at T = 4ms for r 1 = 0.1mm, d = 10mm, V G = 13kV , ψ = 10 10 m -3 (figure reflected on symmetric boundaries)

Figure 13 :

 13 Figure 13: ψ -I curves with r 1 = 0.1mm and V G = 13kV for d = 10mm, V G = 20kV for d = 20mm, V G = 37kV for d = 40mm

10 1 10 3

 3 10 5 10 7 10 9 10 11 10 13

Figure 14 : 1 =

 141 Figure 14: ψ -I curves with d = 10mm and V G = 13kV for r 1 = 0.1mm, V G = 16.5kV for r 1 = 0.5mm

Figure 15 :

 15 Figure 15: V G -I curves with d = 10mm. Numerical currents with ψ = 10 10 m -3 (filled squares) and ψ = 10 11 m -3 (triangles). Experimental data in hollow circles.

r 1

 1 = 0.1mm r 1 = 0.175mm r 1 = 0.2mm data ψ = 10 10 ψ = 10 11 data ψ = 10 10 ψ = 10 11 data ψ = 10 10 ψ = 10

Figure A. 16 2 |

 162 Figure A.16

Table 1 :

 1 Averaged number of iterations & CPU time (minutes) of Douglas-Rachford and Gauss-Seidel algorithms for corona discharge simulations on the 400-element grid

				l )
			Number of time levels	.
	c		DR		GS
		N	CPU time	N	CPU time
	1	14 45	1	17
	10 84 26	1.13 4.65
	100 856 24	1.33 1.68

Table 3 :

 3 V G -I characteristics of experimental data ('data') and numerical results (ψ = 10 10 , 10 11 m -3 )

	11

z e = -1, z p = +1, z n = -1

without ambiguity, all inequalities are point-wise, i.e. v(x) ≥ ψ(x) ∀x ∈ Ω

S e (t, •) + S ψ (t) ≡ { S e (t, •) + v | v ∈ S ψ (t) }

j ≡ q(f p -f n -f e ) is the conduction current density

σ(t, x) ≡ q i∈L |z i |µ i E(t, x) n i (t, x) is the electrical conductivity

this imposes a restriction on ∆t l which is a discrete and local (on each cell K ∈ T ) version of eq.[START_REF] Thomas | Optimization of dielectric barrier discharge plasma actuators for active aerodynamic flow control[END_REF]

see Appendix B

r 2 r 1 V G R Ω Γ f Γ f Γ f

see Appendix C

two-dimensional grids to approximate the flux f i,K,λ KL • ν K,λ KL in eq. ( 12),

where B(z) ≡ z exp(z) -1 is the Bernoulli function, p i,KL ≡ -

Appendix B. Implicit integration of the SGCC schemes

The (p + 1) th -order SGCC flux scheme [START_REF] Nguyen | High-order scharfetter-gummel-based schemes and applications to gas discharge modeling[END_REF] reads as follows,

It can be decomposed as

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

In the above equation, n p i,K -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,

where Ãp,l is the nonlinear operator ensued from 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

Set V (q) = Ãp,l U (q)

4:

Use algorithm 1 or algorithm 2 to solve (I + ∆t l A l + ∆t l B ψ )U (q+1) R l -∆tV (q)

5:

Let

MPI implementation of the Gauss-Seidel algorithm

The sequential nature of the Gauss-Seidel makes the parallelcomputing 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.