Finite Volume Method for Streamer and Gas Dynamics Modelling in Air Discharges at Atmospheric Pressure

Electrical discharges in air at atmospheric pressure like corona or dielectric barrier discharges are generally crossed by thin ionized filament called streamers (about 100μm diameter). The streamer develops and propagates inside the background gas with a high velocity (around 106 m/s) higher than the electron drift velocity (around 105 m/s). During the transport of charged particles within the filaments under the action of the electric field, the energetic charged particles undergo many collisions with the background gas (neutral particles). The interactions between charged and neutral particles generate in turn a gas dynamics characterized by gas temperature and density gradients. The variation of density, momentum transfer and energy of the different particles, present within the ionized filaments, are governed by the fluid conservation laws (or continuity equations) coupled, in the charged particles case, to the electric field or Poisson equation.


Introduction
Electrical discharges in air at atmospheric pressure like corona or dielectric barrier discharges are generally crossed by thin ionized filament called streamers (about 100µm diameter).The streamer develops and propagates inside the background gas with a high velocity (around 10 6 m/s) higher than the electron drift velocity (around 10 5 m/s).During the transport of charged particles within the filaments under the action of the electric field, the energetic charged particles undergo many collisions with the background gas (neutral particles).The interactions between charged and neutral particles generate in turn a gas dynamics characterized by gas temperature and density gradients.The variation of density, momentum transfer and energy of the different particles, present within the ionized filaments, are governed by the fluid conservation laws (or continuity equations) coupled, in the charged particles case, to the electric field or Poisson equation.
It is very important to well known the electro-dynamics characteristics of these atmospheric pressure non thermal plasma generated by streamer or micro-discharge dynamics for an efficient use in the associated applications such as the pollution control of flue gases (Kim, 2004;Marotta et al., 2007), the combustion and ignition improvement (Starikovskaia, 2006), the airflow control (Eichwald et al., 1998;Moreau, 2007) and the biomedical fields (Laroussi, 2002;Fridman, 2008).
In fact, up to now, the optimal use of the atmospheric non thermal plasma sources needs further experimental research works and also modelling investigations in order to better understand the electro-dynamics processes and phenomena induced by the microdischarges (Ebert&Sentman, 2008;Eichwald et al., 2008).In the frame work of the microdischarge modelling, the obtained results (the streamer morphology and velocity, the production of charged and radical particles, the dissipated power) depend on the hydrodynamics physical model (Eichwald et al., 2006;Li et al., 2007) the discretisation method (Finite Difference Method: FDM, Finite Element Method: FEM, or Finite Volume Method: FVM) and the numerical solver used (Ducasse et al., 2007(Ducasse et al., , 2010;;Soria-Hoyo et al., 2008).Indeed, the solution of the micro-discharge fluid models requires high resolution numerical schemes in order to be able to consider the strong coupling between both the Finite Volume Method -Powerful Means of Engineering Design 284 transport and field equations and the steep gradients of the charged particles evolution in a sharp and very fast ionizing wave (Soria-Hoyo et al., 2008).Therefore, the streamer dynamics modelling does not only depend on the selected physical model but also on the accuracy and the stability of the numerical algorithm.Furthermore, the parametric analysis of non thermal plasma discharge requires less time consuming and optimized numerical algorithms.
The present work is dedicated to the use for Finite Volume Method through the streamer discharge simulation and the gas dynamics simulation.We start with an overview on streamer and gas dynamics modelling followed by the model and numerical algorithms for streamers and gas dynamics; in this main part, we explicitly discuss how the model equations are discretized with the help of FVM.Finally, some results about both the streamer discharge and the gas dynamics simulations are shown; in the case of the streamer discharge we also discuss the validation of the present models from comparison between the experiment and the simulation.

Bibliographic overview on streamer discharge and gas dynamics modelling
Initials attempts at the numerical treatment of the electro-hydrodynamic model, in the case of gas discharges at atmospheric pressure, began in the 1960's with Davies et al. (Davies et al., 1964) and Ward (Ward, 1971).They used a first order method of characteristics in the context of the Finite Difference (FD) method.However, poor spatial resolution restricted their study to qualitative results.Towards the late 1970's, Davies improved the method of characteristics by introducing an iterative counterpart that increased the overall accuracy of the algorithm to second order (Davies et al., 1971(Davies et al., , 1975(Davies et al., , 1977)).This method was adopted by several research teams (Kline, 1974), (Yoshida & Tagashira, 1976) and (Abbas & Bayle, 1980) and, as a result, it became the dominant method until the early 1980's.In 1981, Morrow and Lowke (Morrow & Lowke, 1981) presented a work that numerically integrated the system of continuity equations with the two-step Lax-Wendroff method of Roach (Roach, 1972).However, due to numerical dispersion and numerical instability, calculations were restricted to low density plasmas.Such restrictions were overcame by the introduction of the Finite Difference (FD) -Flux Corrected Transport (FCT) technique, originally developed by Boris and Book (Boris & Book, 1973) and extended to two dimensions by Zalezak (Zalesak, 1979).The FCT technique adds an optimal amount of diffusion and is remarkably stable in presence of sharp density gradients.In this context, FCT has become the most frequently used numerical method in streamer discharge modelling since Morrow introduced it for the first time to the gas discharge community in the 1980's (Morrow, 1981(Morrow, , 1982)).Thus, Morrow was the first to offer an analysis (Morrow, 1982) for high density plasmas (up to electron density of 10 13 cm -3 ) and particular attention was paid to the selection of the Courant-Friedrichs-Lewy condition (CFL) (Courant et al., 1928).According to Morrow (Morrow, 1982) the CFL has to take values lower than 0.1 in order to quickly damp out any numerical oscillations resulting from steep density gradients.However, until the early 1980's, the FD based models did not exceed 1.5D description (1D for the continuity equations, and 2D for the electric field calculation in order to take into account the filamentary structure of the streamer: radial extension).By the middle of the 1980's, Dhali and Williams (Dhali & Williams, 1985) had launched the first fully two-dimensional Finite Volume Method for Streamer and Gas Dynamics Modelling in Air Discharges at Atmospheric Pressure 285 simulation, using the FCT technique to solve the continuity equations.Thus, they elucidated several aspects of both positive and negative streamer phenomena.Subsequently, Kunhardt and Wu (Kunhardt & Wu, 1987) improved the FCT method and described a self-consistent numerical simulation of the formation and propagation of streamers in electropositive (N 2 ) and electronegative (N 2 -SF 6 ) gases.Finally, an implicit version of FCT for gas discharge problems was presented by Steinle and Morrow (Steinle & Morrow, 1989).This new algorithm gave a threefold increase in the overall simulation speed because it was able to use a CFL ~ 1 while maintaining the scheme accuracy.However, this new method never became popular among the scientific community.In the mid-1990's, using the same model as Dhali and Williams (Dhali & Williams, 1985), Vittelo et al (Vittelo et al., 1993(Vittelo et al., , 1994) ) reported a more accurate analysis of the negative streamer in N 2 with a quite small spatial resolution (2.5µm -10µm).They also made the first simulations for streamer propagation in non-uniform gaps (point-to-plane electrode configuration) using a fully two-dimensional model.More systematic work on non-uniform gaps was also performed by Babaeva and Naidis (Babaeva & Naidis, 1996) (using the FCT technique), Kulikovsky (Kulikovsky, 1995a(Kulikovsky, , 1995b) ) (using an optimized second-order Shurfetter-Gummel scheme) or Pancheshnyi and Starikovskii (Pancheshnyi & Starikovskii, 2003) (using a first-order upwind scheme).Another efficient second-order numerical scheme was introduced by Van Leer (Van Leer, 1979) and named the second order Monotonic Upwind-Centered Scheme For Conservation Laws (MUSCL) scheme.This algorithm was used through the Finite Volume Method (FVM) in the 3D modelling of high pressure micro-discharges in micro-cavities (Eichwald et al., 1998) in the 1.5D (Eichwald et al., 2006) and the 2D (Ducasse et al., 2007) modelling of the positive streamer propagation.At the beginning of the century (2000) a new approach to gas discharge modelling was presented in the works of Georghiou et al. (Georghiou et al., 1999(Georghiou et al., , 2000) ) and Min et al. (Min et al., 2000, 2001) in which they used the Finite Element Method (FEM) to solve the electro-hydrodynamic model for parallel plate and wire-plate gaps.Based on a Finite Element Flux Corrected Transport algorithm (FEM-FCT) the simulations maintain the ability to handle steep gradients through the use of FCT, but also allow for the use of unstructured triangular cells.This method significantly reduces the number of unknowns, consequently reducing the computing time.Fine resolution is used only where it is necessary, enabling the model to be extended to a fully two-dimensional form and thereby making it possible to model complex geometries.Moreover, the first works that adapted the FEM to the charged carrier conservation equations were the papers of Yousfi et al. (Yousfi et al., 1994) and Novak and Bartnikas (Novak & Bartnikas, 1987) but the former restricted its application to 1.5D problems.Finally, a few works on 3D streamer discharge modeling using either FVM or FEM numerical schemes have been performed by Kulikovsky (Kulikovsky, 1998), Park et al. (Park et al., 2002), Akyuz et al. (Akyuz et al., 2003) Georghiou et al. (Georghiou et al., 2005) Pancheshnyi (Pancheshnyi, 2005) and Papageorghiou et al. (Papageorghiou et al., 2011).

Model and equation discretisations with finite volume method (FVM) for streamer discharge and gas dynamics
In this work, streamer formation and propagation (streamer dynamics) is modelled using the first order electro-hydrodynamic model in the framework of the drift-diffusion approximation (Eichwald, 2006).Moreover, neutral dynamics is not taken into account; only the charged particle dynamics is considered.Thus, the equations involved in this model are the following: In (1) (Maxwell-Gauss) E  is the total electric field (due to geometry and space charge) e the absolute value of the electronic charge, n e , n p and n n the electron, positive and negative ion densities,  c the space charge density and  0 and  r the free space and relative permittivities; here, the relative permittivity is equal to 1.With regard to the continuity equations (3) of the density, subscript "s" stands for the electrons (e) or the positive (p) or negative (n) ions.Moreover, s v  is the velocity and s  the source term; both are functions of the reduced electric field E/n g (n g is the gas density) according to the local electric field approximation.The electron velocity is calculated using the classical drift-diffusion approximation (4), where µ s and D s are the mobility and the diffusion coefficients respectively.
The Maxwell-Gauss equation ( 1) is discretised with finite volume method (FVM).Thus, the equation is integrated on an elementary volume ij V (cell) of the three-dimensional (3D) space.In addition, the Gauss-Ostrogradsky theorem (or divergence theorem) is used to transform the volume integration in surface integration.The Gauss-Ostrogradsky theorem is a mathematical statement of the physical fact that, in the absence of the creation or destruction of matter, the density within a region of space can change only by having it flow into or away from the region through its boundary.After resolution, the following equation is obtained: In the case of a z-axis of symmetry (Fig. 1), the elementary volume and surfaces in cylindrical coordinates are: In equation ( 5), the electric field components are expressed in function of the potential ( 7) to obtain equation ( 8) (FVM discretised Poisson equation); ( 7) is determined by first order finite difference of equation ( 2).Finally, ( 8) is ordered to generate the linear equation ( 9) for the (i,j) cell. where and For the whole elementary volumes of the 3D space, the equations rearrange in a matrix way as Ax=b, with A, x and b, respectively of dimension n r 2 ×n z 2 , n r ×n z and n r ×n z ; x is the solution of the linear equation system.As regards the boundary conditions, we applied a Part of the elementary volume (cell) in the (Orz)plane; localised by the (i,j) discretised coordinates.
system at a later time from the state of the system at the current time.Thus, the discretised equation is the following: The above continuity equations ( 10) and Poisson equation ( 9) are coupled through the source terms and the reduced electric field, which depends on the space charge density.Thus, the algorithms used have to be robust in order to prevent the development of nonphysical oscillations or diffusion phenomena within the electro-hydrodynamic model.
The streamer development is described in a two-dimensional cylindrical (Orz) geometry, where (Oz) is associated with the streamer propagation axis, and (Or) with its radial extension.The next section is devoted to validating and comparing the efficiency of the algorithms to solve the continuity equations ( 10) and Poisson equation ( 9).
As regards the gas dynamic model, the system of equation bellow is used.
The mean energy source term translation processes (depending on the reduced electric field) and   T j (r, t) the current density vector within the streamer discharge simulation.The mean energy is evaluated (for each cell of the calculation domain) on 150ns; the duration is negligible compared with the neutral gas dynamic duration process >1µs (shock wave).The thermal flux density is expressed in function of the thermal conductivity coefficient and temperature gradient of the gas (air).The gas viscosity is not taking into account, there is no reactivity of the gas, the electrodes are at ambient temperature, and gliding conditions are www.intechopen.comFinite Volume Method for Streamer and Gas Dynamics Modelling in Air Discharges at Atmospheric Pressure 289 applied to the electrodes.Finally, if the thermal transfer is not taking into account, the simulation diverges and there is no propagation wave.
One can notice the three equations ( 11) to ( 13) are based on the same structure as (3), the transport term in one side (left term) and source term on the other side (right term).Thus, the algorithm used to solve the equations is the same as (10).

Algorithm tests to solve the energy, momentum, density continuity and poisson equations
Several kinds of algorithms are presented to solve the two equation types of the model: continuity and Poisson equations.The algorithms have been tested in accuracy and computing time on special tests.For the continuity equations (or transport equations) we used the Davies' test (Davies & Niessen, 1990;Davies, 1992;Yousfi et al., 1994;Ducasse et al., 2010) in two directions of a cylindrical coordinate system (r, z).Six algorithms are tested: Upwind, Superbee Monotonic Upstream-centred Scheme for Conservation Law (MUSCL Superbee), Piecewise Parabolic Method PPM, ETBFCT, and Zalesak Peak Preserver (Ducasse et al., 2010).For the Poisson equation we compare the analytic solution to the numeric ones given by MUMPS (Direct method; not iterative) and SOR (Iterative method) (Amestoy, 2001(Amestoy, , 2006;;Fournié, 2010; Press, 2 nd edition).
The Davies' test was first introduced by Davies and Niessen in order to compare the algorithm efficiency to solve one-dimensional continuity equations ( 14) (( 15) is the FVM discretised form) without source term; what we write here is valid for energy, momentum and density continuity equations.The test is interesting for streamer modelling since it reproduces mathematically the behaviour of the streamer head propagation which is a fast ionizing wave that propagates steep density gradients in a sharp velocity field.Thus, the test is performed along a normalized z-axis [0, 1] divided into N=100 regular cells, and consists in propagating a square density profile (wave) n(z,t) in a stationary oscillating velocity field, as shown in Fig. 2. Equations ( 16) and ( 17) give respectively the mathematical expressions for the density and the velocity profiles: This gives a velocity peak value at z0 .5  au (arbitrary unit) which is ten times greater than the values at the beginning and the end of the domain.The initial density of the square wave distribution n(z,t 0)  is enclosed between z 0.05  and z 0.25  with a constant value of 10au.In addition, the time step is chosen equal to 10 -5 au which corresponds to a Courant-Friedrich-Levy number (CFL) equal Finally, the boundaries at z=0 and z=1 are periodic in the sense that any particle leaving the right side boundary enters at the left side; so that, after the period , the transport density profile solution should be identical to the initial distribution n(z,t 0)  .
The comparison of the transported density profile with the exact solution at time t=T will determine the accuracy of the algorithms in handling discontinuities and steep gradients (for t grater than T see (Ducasse, 2010)).In order to quantify the algorithm accuracy, the mean absolute error ( 18) is calculated.Fig. 3 shows the numerical results obtained after one period, whereas Table 1 quantifies the performances of the algorithms in term of accuracy and time consumption.Moreover, with a CFL number equal to 10 -2 , one period T of the square wave evolution corresponds to 59 070 iterations or time steps, which means that the flux correction is applied 59 070 times at the edges of each cells.
MUSCL Superbee, PPM, ETBFCT, Zalesak without Peak Preserver (ZNOPP) and Zalesak Peak Preserver (ZPP) algorithms generate similar results and nearly preserve the solutions from numerical diffusion and dispersion (Fig. 3).Nevertheless, after one period, the results clearly indicate that the PPM algorithm generates the most accurate solution since both the www.intechopen.comsteep gradients and the floor of the square wave are better reproduced (Fig. 3a).Furthermore, the PPM-AE is roughly two times lower in comparison with the other tested algorithms (Table 1); this is in accordance with the previous observations.For Upwind, we observe on Fig. 3a it introduces a large amount of numerical diffusion (comparable to physical diffusion) and its AE is ten times higher than PPM-AE.

Finite Volume Method -Powerful Means of Engineering Design 292
The conservation criterion of the algorithms has been tested too.We observed that the particle conservation is verified for all the algorithms except for ZPP.Indeed, the particle number associated with ZPP increases of 1.1% as it was already emphasized by Morrow (Morrow, 1981).
The last two columns of Table 1, specify the absolute and relative computation time of the six algorithms.The processors used for this comparison are a 3GHz Intel® Pentium® IV with 512Ko of cache memory (768Mo of RAM) and a 2.8GHz Intel® quad-core Nehalem® EX with 8 Mo of cache memory for each processor (18Go of RAM).The results indicate ETBFCT is the fastest but also the less accurate (if Upwind is omitted).Therefore, by taken the ETBFCT values as the reference, it becomes possible to compare the gain of precision relatively to the computing time rise.For example, PPM is 2.32 times more accurate than ETBFCT but the computing time is multiplied by a factor 2.6 (2.5 with Nehalem®).In addition, ZPP is the less efficient since the precision increases by a factor 1.05 only, while the computation time increases by a factor 3.6 (4.3 with Nehalem®).MUSCL and ETBFCT show similar behaviours in term of computing time and accuracy.Moreover, we notice an important computing time fall from Pentium® IV to Nehalem® with a factor 5 for ETBFCT and about 3.5 to 4 for the others.In the particular case of Upwind we see a time consumption divided by 4 compared to ETBFCT, but more than 4 times less accurate than ETBFCT; some author still use the Upwind algorithm with a high space resolution to compensate the numerical diffusion (Pancheshnyi & Starikovskii, 2003;Urquijo et al., 2007).Afterwards, the MUSCL Superbee algorithm is selected to be tested on the Kreyszig radial test (Kreyszig, 1999).Indeed, MUSCL with ETBFCT is the most interesting in terms of computing time and accuracy, with boundaries conditions simpler to implement than ETBFCT.

www.intechopen.com
Finite Volume Method for Streamer and Gas Dynamics Modelling in Air Discharges at Atmospheric Pressure 293 The Kreyszig radial test (Kreyszig, 1999) was used to observe the behaviour of the MUSCL Superbee algorithm for a physical quantity movement along the r-axis, both towards and away from the z-axis (symmetry axis).The test consists of the advection of a normalized square profile along the normalized radial direction, with a constant speed of 10 8 cm.s -1 , a mesh of 100 uniform cells, and a time step fixed at 10 -11 s, which corresponds to a CFL equal to 0.1.In these conditions, Fig. 4  solutions, in a positive constant velocity field and a negative constant velocity field respectively.We note the correct behaviour of the algorithm, which introduce relatively small amounts of diffusion.The sharp corners are determined quite well despite the low spatial resolution.Thus, the solutions can be considered as satisfactory, all the more satisfactory the numerical solution tend to the analytic one with 1000 points and more accurate if a CFL=10 -2 is added.
We conclude the MUSCL algorithm gives interesting results for both the absolute error of the solution and computing time.Moreover, we noted that the numerical solutions of the continuity equation tend to the analytic solutions in both axial and radial directions when the mesh step and (or) the time step are decreased (CFL).
At this stage we examine the numerical behaviour of SOR and MUMPS algorithms (Amestoy et al., 2001(Amestoy et al., , 2006;;Fournié et al., 2010; Press, 2 nd edition) used to solve the Poisson equation (elliptic partial differential equation) without charge (i.e.Laplace's equation).We adopted a hyperbolic point-to-plane configuration for which the analytical potential field is known; the analytical solution was initially proposed by Eyring and first used by Morrow for streamers simulation (Eyring et al., 1928;Morrow & Lowke, 1997).Thus, the analytical solution is compared to the numerical one and the algorithm efficiency is quantified thanks to the relative error.
The curvature radius of the tip is 20µm and the inter-electrode space is 10mm; the applied voltage on the tip is equal to 9kV.The computational domain consists of a structured grid with none constant space cells in each direction.The limits of the domain are 19×19mm in z and r directions (cylinder of 19mm height and 19mm radius) and the number of nodal points in this domain is rz n n 307 1186 364102    (r min =1µm and z min =1µm).The spatial resolution along the z-axis is z=10µm from the plan until 200µm of the point, decrease down to z=1µm on the point, and increase again until the upper boundary; the radial step is progressively increased from r=1µm at the centre until the lateral boundary.
Because of the symmetrical axis (Oz), the radial derivatives along the z-axis are set equal to zero.To perform the potential field comparisons, we set at each nodal point of the open boundaries (r=19mm and z=19mm) the analytical solution (Dirichlet conditions); we also performed the comparisons with a zero Neumann condition, since the simulation use this boundary condition.
The isopotential maps in Fig. 5 compare the analytic field with the MUMPS and SOR solutions.Very good agreements are observed between all results, as well as around the tip than in the whole domain; it can be quantitatively discussed using Fig. 6 for the solution along the z-axis.Thus, the relative error shows the direct method MUMPS gives the closest solution to the analytical one with a value less than 0.1%.For SOR, the error depends on the tolerance chosen for the convergence and the spectral radius (sr influence the speed convergence and solution accuracy); the convergence tolerance is defined as the maximum of the relative difference between the solution at iteration k and k-1.Indeed, with a 10 -5 tolerance and an optimised convergence for the same tolerance (specific sr), we observe an error distribution lower than 0.5%, totally different from the MUMPS one, whereas with a 10 -5 tolerance and an optimised convergence at 10 -10 tolerance (other sr) the error tends to error distribution is constant on a large part of the inter-electrode space and decreases the MUMPS one; the result with a 10 -10 tolerance and an optimised convergence for the same www.intechopen.comFig. 6.Relative error between the analytic, and the SOR and MUMPS solutions along the z-axis.
www.intechopen.comtolerance is superimposed to the MUMPS one.In addition, for the three accurate results the reaching the point; it is due to the mesh, constant at the beginning and that starts to decrease close to the point.The contrary is observed if a constant step mesh is used: the error increases from the plan to the point (Kacem et al., 2011).
Table 2 gives the mean relative error calculated in the whole domain, the computing time and the number of iterations performed to satisfy the convergence tolerance fixed at 10 -5 ; these quantities are given for each method, several domain dimensions, several tolerance conditions, Dirichlet and Neumann Thus, SOR shows the highest mean relative error (0.30%) even if the value is acceptable.Moreover, one can notice MUMPS has the smallest mean relative error (4.2 10 -2 %) since this direct method gives the nearest solution compared to the analytic one; SOR generate the same solution accuracy with both a 10 -10 tolerance, or a 10 -5 tolerance with an optimised spectral radius (convergence speed optimised for a 10 -10 tolerance).Concerning the computing time, the SOR method needs between 17 and 34s to reach the specific tolerance criterion.MUMPS needs 830s to perform the direct calculation; thus, 830s are needed to construct the main matrix, analyse and performed the LU decomposition (Kacem et al., 2011) Afterwards, the point used is not hyperbolic anymore, so the potential at the boundaries of the domain is unknown (a Dirichlet condition is not possible) that is why we impose Neumann conditions.Thus, it is important to check the impact of a Neumann condition on the previous hyperbolic point and compare the numerical solution to the analytic one; the test is performed with the SOR method.Table 2 shows the boundary position in the previous tests (Dirichlet) is 2x2mm².If we choose a 1x2mm² domain, than we observe the mean relative error is too high (40%).But with a 2x2mm² domain, the error decreases down to 11%; the value is acceptable.With a 3x3mm² domain the error is still improved (3.9%) but at 4x4mm² the improvement starts to slow down (2.6%).
For the boundary positions of the next simulation we use the 2x2mm² domain with a tolerance of 10 -5 and a spectral radius obtain for an optimised convergence at a tolerance of 10 -10 ; it is a compromise between the solution accuracy and the computing time.
To finish, the code was compared to another one developed with finite element method.We found a very good agreement with less than 10% of difference (Ducasse et al., 2007).

Simulation results
The first part presents the results obtain with the streamer discharge simulation and the second the results obtained with the gas dynamic simulation.The streamer ionizing wave and the shock wave involved by the streamer discharge are simulated via a PRHE MPI parallelised streamer code; the simulator is able to reproduce both phenomena thanks to efficient algorithms we previously studied (no commercial software is able to do it).Both the streamer and gas dynamics simulations are 3D simulations with axial symmetry; cylindrical coordinates are used.Thus, only half space of a 2D domain (plane) is solved.
The electric discharge is obtained with a point-to-plane electrode system (see the algorithm test part above).The tip curvature radius is 20µm, the gap is 10mm, and the discharge occurs in air at atmospheric pressure; a time varying positive potential (reaches a 9kV maximum on 60ns about) is applied to the point (Fig. 9).The transport of charged particles, their reactivity, their influence on the electric field and the photoionisation phenomenon are taken into account; the air neutral particles are fixed.Moreover, the reaction scheme is composed of 28 reactions the reader can find in (Bekstein et al., 2010) plus 10 ionic recombination reactions and 15 reactions with metastables (Table 3); so the reaction scheme is composed of 53 reactions.Finally, the reaction scheme is composed of electrons, two negative ions (O -and O 2 -) seven positive ions (N 2 + , O 2 + , N + , O + , N 4 + , O 4 + and N 2 O 2 + ) two radical atoms (O, N) and three metastables (N 2 (A), N 2 (a') and O 2 (a)).
Electronic density and electric field are shown; the simulation results are compared with experimental ones for several inter-electrode space values.The simulation code has been parallelised and the computing time is given varying several parameters.Fig. 7 and Fig. 8 show the electric field and electronic density evolutions from 20ns up to 200ns.We observe the primary streamer starts its propagation at 20ns about, i.e. when the applied voltage reaches 3.8kV (Fig. 9).At this time, the space charge formation near the positive point is responsible for the little current peak which appears in the current curve in Fig. 9.The streamer arrives at the plane at 87.5ns about, and corresponds to the maximum calculated current peak.Moreover, a secondary front propagates from the point (the mechanism is different from an ionising wave) in the same time than the streamer propagation and after, during the relaxation phase (beyond 87.5ns); but the speed of this second front is definitely slower.More details about the streamer mechanism, the radical production are available in (Eichwald et al., 2008(Eichwald et al., , 2011)).
A first comparison between simulation and experiment (through the current) was done at one specific potential and inter-electrode space (Eichwald at al., 2008).Here, the simulated currents are compared with the experimental ones at 9, 10, 11mm inter-electrode space and a 9kV DC applied voltage (the power is increased by hand to reach 9kV, afterward we observe a streamer discharge quite periodically).Thus the experimental applied potential shape is different, but the streamer channel is filiform in both cases (simulation and experiment; no branching phenomenon is observed).Fig. 10 shows there is a difference on the main peak of 10mA about at 9mm, and increases with the inter-electrode space; nevertheless the orders of magnitude are the same.Moreover, the bump observed on the experimental curve at 9mm (due to the secondary streamer) is not visible at 10 and 11mm; it seems the phenomenon is of the same amplitude as the primary streamer, but not at all visible on the simulated curves (even at 9mm); it could explain the

Conclusion
In this chapter we have shown how the Finite Volume Method can be used for the discretisation of the transport and Poisson equations, allowing the simulations of streamer discharge development and the associated gas dynamics.It is clear that FVM is attractive since we work directly on elementary volumes that make sense from a physical point of view.Moreover, through a very carefully study, we showed important results as regards the algorithm accuracy, the algorithm convergence (SOR iterative method) the boundary conditions, and the computing time.In the case of the algorithms tested in our research group we have shown first that MUSCL Superbee is the most efficient to treat the conservation laws (Energy, momentum and density); we have also shown that SOR and MUMPS are both interesting in term of computing time.In fact SOR is efficient from a time step to another when the space charge varies slowly (it is the case at relatively low applied potential); MUMPS is efficient if the space charge varies rapidly (the computing time remaining practically the same).In addition, the PRHE-MPI-parallelised code is efficient with 16 processors on the Hyperion HPC system (2.8GHzIntel-Nehalem®; 8Mo cache memory).At this stage, it is possible to do parametric studies since the calculation is fast enough (around three days at 10mm inter-electrode distance and 200ns for the duration).Nevertheless, some improvements on the discharge model still have to be performed from a physical point of view.Indeed, we have shown the experimental current behaves differently when varying the inter-electrode space, whereas the simulation always showed the same shape.May be one of the improvement would be to take into account the local modifications of the air gas properties due to the streamer discharge by the direct coupling of gas dynamics and streamer discharge dynamics.
Fig. 1.Schematic representation of an elementary cell in the three-dimensional space.Due to the z-axial symmetry, the calculation on the half space is sufficient; thus computing time and memory size are saved.
Dirichlet or Neumann condition, following an electrode or an open space is considered.The continuity equation for density of each charge species (3) is also discretised with FVM (like Poisson equation) and an explicit scheme.The explicit scheme calculates the state of a

Fig. 2 .
Fig.2.Initial conditions in the Davies-test case for the density and the velocity field with the densities at times t=T.

FiniteFig. 3 .
Fig. 3. Solutions obtained after one period (59 070 iterations) in the case of (a) Godunov-type schemes and (b) FCT technique.PPM is the most efficient.

Fig. 4 .
Fig. 4. Radial density solution given by the FVM-MUSCL algorithm, (a) for a positive radial velocity, and (b) for a negative radial velocity.

FiniteFig. 5 .
Fig.5.Isopotentials given by the analytic solution (left side) and the numerical solutions (right side) which are identical for the two methods (SOR and MUMPS).NB: Zoom with a size of 5mm×18mm centred over the symmetry axis.The tip isopotential is 9000V and the interval between each isopotential is 900V

Fig. 7 .
Fig. 7. Reduced electric field (Td) distribution as a function of time.

Table 1 .
Absolute Error after one period (59 070 iterations) and mean computing time per iteration (no compilation option) for six numerical schemes.

Table 2 .
, whereas only 0.23s are necessary to calculate the final product of matrices (if the LU decomposition is known, than only 0.23s is necessary to know the potential field).As regards the iteration number, SOR needs between 2400 and 5900 iterations in order to converge; the number depending on the tolerance.Quantitative criteria to compare the method efficiency for Dirichlet and Neumann conditions on the open boundaries (r=rmax and z=zmax).