An implicit time integration approach for simulation of corona discharges

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 ﬂow 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/ﬂoor density imposed on the electrons. Therefore, simpliﬁed 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 identiﬁed the most suitable algorithm, we carry on to two-dimensional wire-to-wire discharges and compare the numerical results with experiment data for diﬀerent ﬂoor densities.


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 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 [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 proposed throughout the years.Morrow [31] 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 [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-toplane [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 negative 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 electron 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 layers near the electrodes.Dufour & Rogier developed a sub-cycling algorithm in their multiscale plasma solver COPAIER [41].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 [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 ψ 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 [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.

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 Ω ⊂ 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, 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 ∈ Ω, 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 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 The electric potential φ(t, x) (V ) is given by Gauss' law, 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 [46,55].In order to limit the complexity of the system, we adopt a simplified model from J.-P.Boeuf 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 → 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, 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 [56]; 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, 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.
The characteristic function of K is defined as The subdifferential [57, Proposition 3] of I K at v is a set of functions defined as 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 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 ⊂ H be a non empty closed convex.Then for all v ∈ H, there exists a unique u ∈ K (its projection) such that, 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 [59] which reads as Finally, the simplified LFA model reads as follows, 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 [60] and variational inequalities [61,62,63].

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.
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 [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 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 8 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

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 ψ can also be envisioned to include the naturally formed electrons.

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 ∆x and at the time level t l , the time step ∆t l ≡ t l+1 − t l must satisfy6 ∆t l ≤ min 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), with . 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 [69] except the nonlinearity only stemming from S ψ .
Let us denote as I the identity operator,  .Then at t l , the semi-discrete (in time) conservation laws of eq. ( 10) can be rewritten as 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.(10), where h is an element of S l+1 ψ .Now integrating this equation on K and dividing it by |K|, we have 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 x K,λ the center of λ ∈ E K and then n e n p by n e n p .Finally we have for every K ∈ T .Now for any smooth function v defined on Ω, we define a piecewise constant function v on Ω such that It is not evident how h with h ∈ S l+1 ψ relates to ψ and n e .Therefore, we make an assumption that 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, where 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 [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 ∂I K is maximal monotone on the Hilbert H.We recall that a (multivalued) operator A is monotone if 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 [70].In our case, if (I + ∆t l ∂I K )w y then from the definition eq. ( 4) we have, 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 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 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 [70,72].The original method was introduced as a first-order-in-time operator splitting to solve the evolution problem 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 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 [73] problem, where τ is the pseudotime.We apply the Douglas-Rachford algorithm for 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 For (AU , U ) H ≥ 0, ∆t l must satisfy 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 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, 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, 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 As long as 1 + ∆t l A l pp > 08 , ∀p so that B ψ is monotone, the resolvent of B ψ is well defined as in eq. ( 15) and we have The complete description of the proposed algorithm is presented in algorithm 2.

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 Γ 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 [23].We set max x S = 5mm.
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,

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 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.

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 (∆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 [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 .
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 [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 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 [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.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 min .
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.

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 [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 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 (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. [7].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 [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 (µ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.

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

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

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

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

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

Figure 10 :
Figure 10: Sketch of the computation domain

Figure 11 :
Figure11: Example of grid refinement near the smaller wire generated by Gmsh[66]

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

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

Table 2 :
Parameters of the wire-to-wire simulations