Keywords: Numerická Simulace Propagace Streameru na Dynamicky Adaptovaných Sítích numerická simulace, pohyb streameru, rozvětvení streameru, minimální model, dynamická adaptace Numerical Simulation of Streamer Propagation numerical simulation, streamer motion, streamer branching, minimal model, dynamical grid adaptation

Anotace

Tato dizertační práce se zabývá numerickou simulací propagace streameru (elektrický výboj ve vysokonapět'ovém elektrickém poli). Pro popis streameru je použit minimální model, který se skládá ze soustavy transportních rovnic pro elektricky nabité částice spárovaných s Poissonovou rovnicí pro elektrický potenciál. V práci simulujeme obecný pohyb streameru ve 3D. Tento obecný pohyb je prezentován rozvětvením streameru, kterého se dosáhne pomocí lokálních poruch

List of Symbols

A [cm -2 Torr -2 ] fitting parameter D e [cm 2 •s -1 ] diffusion coefficient D σ ij diamond cell E [Vcm -1 ]
intensity of electric field I ion emission N

[cm -3 ] neutral gas density n e

[cm -3 ] electron density n i (n p ) [cm -3 ] positive ion density n n

[cm -3 ] negative ion density p 0 p ratio of the atmospheric pressure to a gas pressure

q [C] elementary charge v e [cm• s -1 ] electron velocity v p [cm• s -1 ] positive ion velocity v n [cm• s -1 ] negative ion velocity V [V ] electric potential t [s] time △t [s] time step T i mesh cell x, y, y [cm] coordinates α [cm -1 ] ionization coefficient ǫ 0 [Fm -1 ] permittivity β ep [cm

Introduction

When a neutral gas is exposed to high intensity electric field, non-equilibrium ionization processes (discharges) occur. The ionization causes that neutral particles are decomposed to ions and free electrons. Meanwhile the gas temperature remains close to ambient temperature, the electron temperature increases to a high level (T e > 10000 K). The discharges can appear in various forms depending on spatiotemporal conditions such as the intensity of the electric field, volume and pressure of the medium etc. One can distinguish stationary and unsteady discharges. The dark, glow or ark discharges belong to the stationary ones. Leaders and streamers are transient unsteady discharges. We are interested in streamers that are filamentary discharges which dynamics are controlled by highly localized and nonlinear space charge regions. The streamers can be used in several technical applications, e.g.: streamers are used for treatment of contaminated media like exhaust gases, polluted water or bio gas because the streamers emit the reactive radicals . The streamers are also used in the processing of electronics, especially for etching and deposition of thin films. They can be applied in surface modification (e.g. hardening, corrosion resistance) too.

The atmospheric pressure plasma jets belong among other utilizations of streamers. They can be used for biomedical applications because they can produce various biocidal agents such as UV radiation, reactive species and charged particles. The atmospheric plasma jets can be easily handled and can treat a target from a distance without thermal damage. Therefore they have a potential to be used in diverse fields of biomedicine. Another utilization of the weakly ionized gas is to control a boundary layer in airflow at high speeds because achievements in aerospace allow to create hypersonic air vehicles at a Mach number greater than 6. The traditional flow control mechanism is often ineffective at such high velocities due to significant physical-chemical modification of the airflow. The streamers can be undesirable too. They can cause e.g. transmission line outgoes at power system when a streamer appears in a disturbed electric field around an insulator string.

The interest is how to protect the insulator string from the streamer.

The generation of streamer discharges is described by coupling of electrostatic with motion of charged particles (electrons, positive and negative ions). The electrostatic is represented by Poisson's equation for the electric potential and the particles motion is described by a set of convection-diffusion-reaction equations. The charged particles are considered as continuum.

It means that we deal with the particles in a form of densities and we can use continuum mechanics for them. The motion of charged particles causes a local nonlinearity in the electric field which is an essential element for creation of the streamer discharge. Another important phenomenon in the streamer motion is photoionization which decompose neutral particles to charged ones and it also locally changes the electric field. The photoionization effect is very often described by Helmholtz or Eddington models or it is replaced by background densities of the electrons and the positive ions in many studies. The charged particles react with surrounding environment in chemical reaction during the streamer propagation. It produce chemically active species which are interesting for many applications.

Governing equations

The streamer motion can be describe by different models which differ in a complexity. The streamer propagation is essentially formulated in a form of fluid models that couple unsteady transport equations for charged species and Poisson's equation for the electric potential.

One of these approaches is two-fluid Euler-Poisson system which is considered e.g. in [START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF]. This approach consist of a continuity equation and a convection-diffusion-reaction equation for the electrons and the positive ions (one pair for each specie)

(n i ) t + ∇ • (n i u i ) = 0, m i [(n i u i ) t + ∇ • (n i u i ⊗ u i )] + ∇p i (n i ) = qn i E,
(n e ) t + ∇ • (n e u e ) = 0, m e [(n e u e ) t + ∇ • (n e u e ⊗ u e )] + ∇p e (n i ) = -qn e E, ∇ 2 V = -q ǫ 0 (n i -n e ) , E = -grad(V ),

where n i,e denotes positive ion (electron) density, u i,e is its velocity. V is the electric potential, E is the intensity of the electric field q is the elementary charge and ǫ 0 is the dielectric constant. p i,e denotes particle pressure of ions (electrons). Another approach, which is considered in many papers, deal with three convection-diffusionreaction equations for charged species coupled with Poisson's equation for the electric potential

+ ∇ • (n n v n ) = n e β att | v e | -n n n p β np , ∇ 2 V = - q ǫ 0 (n p -n n -n e ) , E = -grad(V ), ∇ 2 S ph -(λp O 2 ) 2 S ph = -Ap 2 O 2 I,
where n e , n p , n n are electron, positive and negative ion densities, v e , v p , v n are drift velocities of the particles and D e is a diffusion coefficient. V is the electric potential and E denote the intensity of the electric field. The source terms on the right hand side of convectiondiffusion-reaction equations represent ionization (terms multiplied by a constant α), attach-ment (terms multiplied by a constant β att ), electron-positive ion recombination (β ep ) and negative-positive ion recombination (β np ). The source term S ph denotes the photoionization and is expressed by a Helmholtz equation (see above). A and λ in the Helmholtz equation are fitting parameters and we can find them e.g. in [START_REF] Papageorgiou | Three-dimensional numerical modelling of gas discharges at atmospheric pressure incorporating photoionization phenomena[END_REF]. p O 2 is a partial pressure of the oxygen and I represents photon emission which can be described by relation I = A • n e • v e • λ. We can see the precise expression for drift velocities, diffusion coefficient and coefficients in source terms in relations (B.1) -(B.8). There is the simplest approach, so-called minimal model, among the approaches. It's the simplest model which is able to capture the streamer propagation and consists of a convectiondiffusion-reaction equation for the electron density and an ordinary differential equation for the positive ion density. These two equations are coupled with the Poisson's equation for the electric potential [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF], [START_REF] Nikandrov | Streamer simulations with dynamically adaptive cartesian mesh[END_REF] ∂n e ∂t + ∇ • (n e v e ) -∇ • (D e ∇n e ) = n e α| v e |, ∂n p ∂t

= n e α| v e |, ∇ 2 V = -q ǫ 0 (n p -n e ) , E = -grad(V ).

The motion of the positive ions is neglected in the minimal model because the velocity of the electrons is much higher than the positive ion velocity. Only the ionization is considered in source term and the photoionization is replaced by background densities for electrons and positive ions. All approaches have the same relation between transport equations of charged particles and the electric field. The electric potential depends on charged particles because their densities have an influence on the right hand side of Poisson's equation. On the other hand the material properties (drift velocities, diffusion coefficient etc.) depend on the intensity of the electric field (see appendix B) which is a negative gradient of the electric potential. It leads to the coupled system of equations which can't be solved separately.

Numerical simulation

Because of the streamer nature which is a filament of a plasma, the first attempts how to simulate streamer motion were 1D simulations where the 3D nature of the electric potential (Poisson's equation) was solved by a disk method. This approach gives sufficient information about distribution of electric field, electron drift velocity, charge etc. along the streamer axis but it isn't able to tells us the distributions in a space. Therefore there is a need to simulate streamer in more dimensions. A plane to plane electrode system or a pin to plane electrode system is very often considered in streamer simulations. The computational domain is very simple when the plane to plane system is taken into account. It allows to simulate the streamer as a 2D axisymmetric case. The system of governing equations is then written in a cylindrical coordinates in such a case. The computational domain is in this case a rectangle and the initial nonlinearity in the electric field is introduced by distribution of the electron and the positive ion densities. The cartesian grid is suitable for such simple computation domain (see e.g. [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF], [START_REF] Nikandrov | Streamer simulations with dynamically adaptive cartesian mesh[END_REF], [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF]) and the finite difference method is mainly used. The finite difference method has advantages that the cell vertex schemes are simple and mostly come from Taylor expansion. To avoid the high numerical viscosity in these schemes, they can be easily extended to obtain higher orders of accuracy. The influence of higher orders schemes on streamer motion is studied e.g. in [START_REF] Bourdon | Influence of numerical schemes on positive streamer propagation. 27 ICPIG[END_REF] where the first order upwind scheme, the second order MUSCL scheme family and QUICKEST scheme family (3 rd order and 5 th order) are compared.

When we would like to solve streamer motion in more general geometries, it is better to use unstructured grids. The simplest case of a such geometry is the pin to plane electrode system where the shape of the pin electrode is necessary to capture. The finite volume method (FVM) is better choice for computation on unstructured grids but on the other hand it's difficult to extend these schemes to higher orders of accuracy (higher than the second order). Another possibility is to use the finite element method (FEM) on the unstructured grids. The comparison of FVM and FEM methods is discussed e.g. in [START_REF] Ducasse | Critical analysis on two-dimenional point-to-plane streamer simulation using finite element and finite volume methods[END_REF]. It seems that the computational cost is higher for the FEM than for the FVM. It's necessary to skip to 3D if we want to simulate the streamer motion in more complex geometries (e.g. different shapes of electrodes) or more general motion which is caused by local disturbances in the electric field. The streamer propagation in the plane to plane electrode system computed on 3D cartesian grid is e.g. in [START_REF] Nikandrov | Streamer simulations with dynamically adaptive cartesian mesh[END_REF]. The attempts of general streamer motion (streamer branching), which is introduced by local seeds of plasma spots, are discussed e.g. in [START_REF] Papageorgiou | Three-dimensional numerical modelling of gas discharges at atmospheric pressure incorporating photoionization phenomena[END_REF] and an unstructured mesh is considered at this approach.

Mesh adaptation

The difficulty in a numerical simulation of the streamer propagation comes from a very narrow region where the variables governing the physics of the motion vary in several orders of magnitude. Therefore we need a fine grid in this region. A mesh size would be too large if the mesh will be fine in a whole computational domain and the computational time and memory cost would be too high for a such mesh. A way how to avoid this disadvantage in the simulation is a dynamical mesh adaptation. Almost all of streamer simulations use adapted meshes now. We can see the mesh adaptation on 2D cartesian grids e.g. in [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF], [START_REF] Nikandrov | Streamer simulations with dynamically adaptive cartesian mesh[END_REF], [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF]. In [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF], they treat independently with refinement procedures for continuity equations and for Poisson's equation to capture the streamer head front in a cartesian grid. The 3D mesh adaptation is used e.g. in [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF] where the streamer motions is computed in a plane to plane electrode system. The streamer motion is also simulated on 2D unstructured grids where the finite-element method is applied together with using of an adaptive mesh generator [START_REF] Min | A study on the streamer simulation using adaptive mesh generation and FEM-FCT[END_REF], [START_REF] Papadakis | New high quality adaptive mesh generator utilized in modelling plasma streamer propagation at atmospheric pressures[END_REF]. The adaptation criterion deals several variables, ratio of the gradient of the electron density and the magnitude of the electron density|gradn e |/n e , the magnitude of the net charge density ρ = e(n i -n e ) and the magnitude of the ionization source term. Advantage of the mesh adaptation on a structured grid is that we can use relatively simple schemes easily extended to higher orders of accuracy but on the other hand we obtain non-conform grids. The mesh adaptation on unstructured grids gives conform grids because a conformity step is applied after cells refinement. The mesh adaptation is especially necessary in the 3D computation where the mesh size could oversize several millions of cells. That leads to high memory cost and the computation takes several weeks to obtain results.

There are still lots of challenges in the streamer simulation regardless of the great progress in the streamer simulation. Most of the streamer simulations published in the literature take into account the axisymmetric approach and geometries which aren't able to capture the structure observed in practice simulations or in the nature. The full 3D simulation is needed to solve the streamer branching or the streamer propagation in complex geometries.

Goals of the Thesis

The main aim of the thesis is to simulate the streamer branching. Therefore we have to develop fully 3D method with dynamical grid adaptation. The physical model can be simplified to the minimal model of the streamer. The unstructured grid is a proper choice even with the disadvantage of lower order schemes to simulate a general motion of the streamer. 3D computation is very time consuming. Numerical method have to be developed and tested on simplified cases. The proper intermediate step is a 2D case (planar discharge) even if it's a different type of problem than the filamentary discharge. We want to test properties of the scheme (numerical tests of convergence), the adaptation algorithm and proper formulation of conditions. Computing code for an unstructured grid needs series of support utilities for description and manipulation with a grid (e.g. description of connectivity). Mainly the algorithm of the grid adaptation needs a special care. Modern computer languages (e.g. C++) offer suitable tools for this purposes. Applicability of such programs exceeds range of the streamer motion simulation. They can be used in many applications of finite volume method for different mathematical models (Euler equations, Navier-Stokes equation, flow in porous media, shallow water equations etc.) and different type of problems (combustion, flow with shock waves, multi-face flow etc.) where the dynamical mesh adaptation can increase the quality of results.

We therefore formulate following goals of the thesis:

• Simulation of streamer branching by introducing temporal disturbances;

• Development of a numerical method for a 3D streamer motion calculation described by the minimal model and based on dynamic adaptation of unstructured grids.

Intermediate steps:

• Development of 2D method (in a form which allows simple transition to 3D) for numerical experiments to verify the properties of the method;

• To develop a general platform based on object programming in C++ language for simulation of various physical problems on unstructured grids with a dynamical mesh adaptation in 2D and 3D.

Chapter 2

Governing equations

The modeling of streamer propagation mechanism in atmospheric pressure discharges received a large attention and was the subject of a large research effort during the last two decades. Until very recently, streamer propagation was essentially described through fluid models that couple non-stationary transport equations for charged species (electron and ions) with a Poissons (elliptic) equation for the electrical potential. Usually the transport equations of charged species consist of time-dependent convection-diffusion-reaction equations. In principle, when using advanced physical model, the source terms in these equations include a non-local component induced by the photoionization effect that was often invoked as the main propagation driver at least for streamer in field-free space. In this case, a radiation transport model has to be used and a non local radiation transfer equation has to be coupled to the transport and Poisson's equations. This results in a very high computational costs and recently an advanced physical model based on Eddinton's approximation allowed a very efficient and physically sound description of the photoionization phenomenon. It remains that very often a background seed electron field was used to simulate the electrons generated by photoionization phenomena and responsible for the avalanches that takes place at the streamer tip and that are responsible for streamer propagation. The photoionization effect is deeply studied e.g. in [START_REF] Célestin | Etude de la dynamique des streamers dans l'air a la pression atmospherique[END_REF] where several approximation is taken into account. There are also more complex boundary conditions for the Poisson's equation in [START_REF] Célestin | Etude de la dynamique des streamers dans l'air a la pression atmospherique[END_REF].

The simplest minimal model of discharge motion [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF], [START_REF] Nikandrov | Streamer simulations with dynamically adaptive cartesian mesh[END_REF] (

The diffusion coefficient D e is a function of the electron drift velocity and the electric field [START_REF] Morrow | Streamer propagation in air[END_REF]]

D e =   0.3341 • 10 9 • E N 0.54069   • v e E . (2.5)
The source terms depend on the electron drift velocity and the electron density

S e = α N • v e • n e • N, (2.6)
S + i = S e (2.7)
where the ratio

Initial conditions

The electric field is established between two planar electrodes. High voltage 25000 V is applied at the anode and the cathode is grounded. The disturbance in the electric field between these two planar electrodes is given by an initial Gaussian pulse. The photo-ionization effect, which is necessary to run streamer, is substitute by a background density:

2 dimensions n e (x, y, 0) = 10 16 • e -(x-0.2) 2 +(y-0.25) 2 σ 2 + 10 9 [cm -3 ], σ = 0.01, n i (x, y, 0) = n e (x, y, 0) , (2.9)
3 dimensions n e (x, y, z, 0) = 10 + 10 8 [cm -3 ], σ = 0.01, n i (x, y, z, 0) = n e (x, y, 0) .

(2.10)

Boundary conditions

The Neumann's homogeneous condition is considered at all boundaries for the electron density

∂n e ∂ n = 0. (2.11)
Following Dirichlet's conditions are prescribed for the electric potential at the anode and the cathode

V = 25000 [V], for anode, V = 0 [V],
for cathode.

(2.12)

The Neumann's homogeneous condition is prescribed at the rest of boundaries

∂V ∂ n = 0. (2.13)
Chapter 3

Numerical Approximation

The streamer propagation phenomenon is described through the time-space evolution of physical characteristics such as electron density, ion density, electric field and space charge. All these characteristics experience a several orders of magnitude variation in a very narrow region located at the microdischarge tip (the so called streamer region). Further, the transport equations and the electric field equation are strongly coupled through source-terms, transport fluxes and space charge. The propagation dynamics intimately depends on the space-time distribution of the microdischarge characteristics in this tinny region with a characteristics length of few micrometers. Researchers therefore concentrate on the development of high order schemes, e.g., FCT, QUICK, (W)ENO, MUSCLE, improved Sharfeter-Gummel, etc., for 1D geometries and 2D-axisymmetrical geometries. The use of highly non-symmetrical electrode system, along with the non-planar nature of the streamer propagation motivated some groups to adapt some of these numerical schemes to the case of non-structured grids.

Nevertheless the streamer simulation remains costly and time consuming. That's the reason why a simple upwind scheme extended by a Van Leer's type MUSCL algorithm together with Barth-Jespersen limiter is used in this work. The dynamical mesh adaptation improve the accuracy of the simple scheme. There is possible to use another limiters such as Van Albada limiter (see section 4.1.2). The Barth-Jespersen limiter is chosen because it's a quite simple limiter and doesn't depend on parameters which should be tuned.

As it was mentioned before, streamer propagation is represented by unsteady convection diffusion reaction equations for electric particles coupled with Poisson's equation for electrical potential. The unsteady equations are solved by an explicit scheme on unstructured grid. The discretized Poisson's equation leads to a system of algebraic linear equations which is solved with a direct solver at each time step.

Numerical Method

The system of convection-diffusion-reaction equations for the charge particle densities (electrons, positive ions) and the Poisson's equation for the electric potential are treated by the finite volume method on unstructured meshes. The algorithm is

• new values of the electron density n e and the ion density n i are computed by an explicit scheme;

• the electrical potential is established from the new values n e , n i by numerical solution of the Poisson's equation;

• the electric field E is computed from the electric potential;

• the electron drift velocity v e and the diffusive coefficient D e is computed as functions of the electric field E;

• a new time layer is computed => new values of n e , n i ;

• a new mesh is made by a dynamical mesh adaptation after several iterations;

• values of unknowns are interpolated form the old mesh to the new one.

Discretization of Electron Density Equation

The finite volume method is used for discretization of the electron density equation. The domain Ω is divided in a system of subdomains µ i , Ω ≡ i µ i , which do not overlap each other and have common boundaries. Integration over volume µ i gives (the treatment is the same in 3D and leads to the same general form)

dn e dt = - 1 µ i m j=1 n n e ij v e ij n ij |σ ij | -D e ij ∇n n eij n ij |σ ij | + S N e i , (3.3)
where m is number of faces of volume µ i , n ij is the unit normal vector of the face σ ij (face between volumes µ i and µ j ) and |σ ij | is its length. Other variables denoted by subscript ij represent variables on the face σ ij . At the end, the time derivative is approximated by the three steps Runge -Kutta method

n (0) e i = n n e i , n (k) e i = n n e + α (k) Res n (k-1) e i , (3.4
)

n n+1 e i = n (3) e i ,
with k = 1, 2, 3 and the coefficients α (1) = 0.5, α (2) = 0.5, α (3) = 1. The residual in equation (3.5) is given by Res(n (k-1)

e i ) = - △t µ i m j=1 n (k-1) e ij v e ij n ij |σ ij | -D e ij ∇n (k-1) eij n ij |σ ij | + △t • S (k-1) e i , (3.5)
where △t is a time step.

Discretization of Convective Terms

The convective flux in electron density equation is computed by simple upwind scheme

n e ij = n e i if ( v e ij • n ij ) ≥ 0, n e j in other case, (3.6)
with assumption that the normal vector n ij is oriented from the cell T i to the cell T j . Computation of the interface drift velocity v e ij will be explained later.

Higher order approximation in space

The upwind scheme has a first order of accuracy. It leads to a high numerical dissipation which is necessary to compensate by higher order algorithm. The scheme is therefore extended by a Van Leer's type MUSCL algorithm together with Barth-Jespersen limiter.

One introduce n e i , n e j in the equation (3.6) instead of n e i , n e j

n e i = n e i + ψ i • ∇n e i • r i , (3.7
)

n e j = n e j + ψ j • ∇n e j • r j ,
where ∇n e i , ∇n e j are gradients of electron density in cells T i , T j . We can compute the gradients assuming that the electron density is piecewise linear function and its value n e i is in the center of gravity of the cell T i . This linear function is computed by the least square method which include all neighbors of the cell T i (neighbors with at least one common vertex). r i is a vector coming from the center of gravity of the cell T i to the center of gravity of a face σ ij . In the same way, r j is a vector coming from the center of gravity of the cell T j to the center of gravity of the face σ ij . ψ i and ψ j are Barth-Jespersen limiter function.

The gradient ∇n e i = ((n e ) x , (n e ) y ) from equation (3.7) is computed by following formulas in 2D

(n e ) x =

I yy • J x -I xy • J y D , (n e ) y = I xx • J y -I xy • J x D , (3.8)
where I xx , I yy , I xy , J x , J y , D are given by equation (A.1).

We obtain a little bit more complicated relations in 3D for components of the gradient ∇n e i = ((n e ) x , (n e ) y , (n e ) z )

(n e ) x = I yy I zz -I 2 yz J x + (I xz I yz -I xy I zz ) J y + (I xy I yz -I xz I yy ) J z D , (3.9)
(n e ) y = (

I xz I yz -I xy I zz ) J x + (I xx I zz -I 2 xz ) J y + (I xy I xz -I yz I xx ) J z D , (3.10)
(n e ) z = (I xy I yz -

I xz I yy ) J x + (I xy I xz -I yz I xx ) J y + I xx I yy -I 2 xy J z D , (3.11)
where I xx , I yy , I zz , I xy , I xz , I yz , J x , J y , J z , D are defined by equations (A.2)-(A.4).

The Barth-Jespersen limiter is given by following formula

ψ i =        min 1, ne max -ne i △ 2 , △ 2 > 0, 1, △ 2 = 0, min 1, ne min -ne i △ 2 , △ 2 < 0, (3.12)
where n emax is maximum of electron density from values in the cell T i and all neighboring cells with common face, n e min is minimum of electron density from values in the cell T i and all neighboring cells with common face. △ 2 is given by relation △ 2 = ∇n e i • r with r which comes from the cell's center of gravity to the center of gravity of one of faces of the cell T i . Final ψ i is the minimum over the faces.

Discretization of Dissipative Terms

The curvilinear integral of the diffusive term in the equation (3.2) is approximated by numerical integration around a cell boundary, so one can write

∂µ i D e ∇n n e nds = m j=0 D eij ∇n n eij n ij |σ ij |, (3.13)
where ∇n n eij represents gradient of the electron density on the face σ ij , n ij is a unit normal vector of the face σ ij and |σ ij | is its measure. The constant m indicate the total number of a cell's faces. Computation of the diffusion coefficient D eij is described in the equation (3.32).

Gradient Computation in 2D

Approximation of the electron density gradient ∇n eij in 2D is done on the diamond cell (see Fig. One can assume that the electron density gradient is constant on the co-volume D σ ij . According to Green-Gauss theorem the approximation leads to

∇n eij = 1 µ D σ ij 1 2 (n e (A) + n e (R)) n AR |σ AR | + 1 2 (n e (R) + n e (B)) n RB |σ RB | + + 1 2 (n e (B) + n e (L)) n BL |σ BL | + 1 2 (n e (L) + n e (A)) n LA |σ LA | , (3.14)
where n e (A) (n e (B), When we factor one half out of the square bracket and modify the expression inside the square brackets, we obtain

• • • ) represents value of the electron density in the point A (B, • • • ), n AR is
∇n eij = 1 2µ D σ ij ( n AR |σ AR | + n LA |σ LA |) n e (A) + ( n RB |σ RB | + n BL |σ BL |) n e (B) + + ( n AR |σ AR | + n RB |σ RB |) n e (R) + ( n BL |σ BL | + n LA |σ LA |) n e (L) . (3.15)
The following equalities are obvious

n LA |σ LA | + n AR |σ AR | = n LR |σ LR |, n AR |σ AR | + n RB |σ RB | = n ij |σ ij |, (3.16
)

n RB |σ RB | + n BL |σ BL | = -n LR |σ LR |, n BL |σ BL | + n LA |σ LA | = -n ij |σ ij |.
If we substitute relations (3.16) in the equation (3.15), it gives

∇n eij = 1 2µ D σ ij (n e (A) -n e (B)) n LR |σ LR | + (n e (R) -n e (L)) n ij |σ ij | . (3.17)
The electron density values n e (A), n e (B) in the face σ ij end points are computed by the least square method

n e (A) = N (A) p=1 α p (A)n ep , n e (B) = N (B) p=1 α p (B)n ep , (3.18)
where n ep is value of the electron density in the cell T p , N (A) (N (B)) is number of cells including vertex A (B), α p (A) (α p (B)) are weights coming from the least square method. One can write (for vertex A)

α p (A) = 1 + λ x (x p -x A ) + λ y (y p -y A ) N (A) + λ x R x + λ y R y , (3.19)
where the weights parameters are given by formulas (A.5).

Gradient approximation in 3D

The diamond cell D σ ij in 3D is constructed by the centers of gravity L, R of cells T i , T j which contain the common face σ ij . The rest of vertices of the diamond cell are vertices of the face σ ij (see Fig. 3.2). The electron density gradient ∇n eij is approximated by following formula which coming from the equations (A.6) -(A.9)

∇n eij = 1 3µ D σ ij (n e (A) -n e (C)) n BRDL |σ BRDL | + (3.20) + (n e (B) -n e (D)) n ALCR |σ ALCR | + (n e (R) -n e (L)) n ij |σ ij | .
The electron density values n e (A), n e (B), n e (C), n e (D) are computed by the least square method

n e (A) = N (A) p=1 α p (A)n ep , n e (B) = N (B) p=1 α p (B)n ep , n e (C) = N (C) p=1 α p (C)n ep , n e (D) = N (D) p=1 α p (D)n ep . (3.21)
The weights α p (A) is given by

α p (A) = 1 + λ x (x p -x A ) + λ y (y p -y A ) + λ z (z p -z A ) D , (3.22)
where λ x , λ y , λ z , D are expressed in relations (A.10) -(A.12).

Discretization of the Positive ion Density

The equation for ion density is in the minimal model only an ordinary differential equation and its numerical approximation is represented by three steps Runge -Kutta method (equivalent of equation (3.5))

n (0) i i = n n i i , n (k) i i = n n i + α (k) Res n (k-1) i i , (3.23
)

n n+1 i i = n (3)
i i , where the coefficients α (k) are the same as in the equation (3.5) and the residual is expressed by relation

Res(n

(k-1) i i ) = △t • S (k-1) e i .
(3.24)

Discretization of Poisson's Equation

Discretization of Poisson's equation leads to a system of linear equations with a sparse matrix. Such system of equation can be solved by direct methods or by iterative methods. Iterative methods gives possibilities as Jacobi iterative method, Gauss-Seidel iterative method or super relaxation method which is improved Gauss-Seidel iterative method. There are also gradients method such as Biconjugate gradients method or gmrs method (generalized minimum residual method). The gmrs method together with proper preconditioning is very often used in numerical simulations (solving of implicit scheme, etc.). Gauss elimination or LU factorization belong among direct methods. The LU factorization is more suitable when right hand side of the system of equation is changed. Because of dynamic mesh adaptation, the matrix of the system of equations is changed after a given number of iterations and the right hand side is changed each iteration. This leads to "quasi-stationary" matrix of equations and therefore LU factorization is the proper method in our case. We use MKL library from Intel to solve the system of liner equations by LU factorization. The Poisson's equation is discretized by a central type approximation which leads to a system of linear equations

A • V n+1 = b n+1 . (3.25)
A is a matrix of coefficients, V is a vector of unknowns (its dimension is equaled to the total number of cells) and b is a vector of right hand side. The updated values n n+1 e , n n+1 i are used for the evaluation of the source term in equation (2.2). A row i in the matrix A corresponds to the cell T i . We use a similar finite volume method approximation as for diffusive terms in the equation for electron density (3.13).

1 µ i Q = 1 µ i µ i △V dV = 1 µ i ∂µi ∇V nds ≈ 1 µ i m j=1 ∇V ij n ij |σ ij |. (3.26)
An approximation of the gradient ∇V ij is also made by the diamond cell (see Fig. 3.1, Fig. 3.

2).

According to the equation (3.17), one can write for 2D case

∇V ij = 1 2µ(σ ij ) [(V A -V B ) • n LR |σ LR | + (V R -V L ) • n ij |σ ij |] . (3.27)
When we substitute the equation (3.27) to a relation

△V = 1 µ i m j=1 ∇V ij n ij |σ ij |, (3.28)
we obtain the final form of approximation of the Poisson's equation in 2D

△V = m j=1 1 2µ i µ D σ ij • n LR |σ LR | • n ij |σ ij | • V A - - 1 2µ i µ D σ ij • n LR |σ LR | • n ij |σ ij | • V B + + 1 2µ i µ D σ ij • n ij |σ ij | • n ij |σ ij | • V R - - 1 2µ i µ D σ ij • n ij |σ ij | • n ij |σ ij | , (3.29)
where electrical potential values V A , V B are computed by the least square method (see equation (3.18)) with weights computed according to equation (3.19).

For 3D case, the electrical potential gradient is approximated according to equation (3.20), so one can write

∇V eij = 1 3µ D σ ij (V A -V C ) n BRDL |σ BRDL | + (3.30) + (V B -V D ) n ALCR |σ ALCR | + (V R -V L ) n ij |σ ij | .
If we substitute the equation (3.30) to the equation (3.28), the numerical approximation of the electrical potential in 3D leads to

△V = m j=1 1 3µ i µ D σ ij • n BRDL |σ BRDL | • n ij |σ ij | • V A - - 1 3µ i µ D σ ij • n BRDL |σ BRDL | • n ij |σ ij | • V C + + 1 3µ i µ D σ ij • n ALCR |σ ALCR | • n ij |σ ij | • V B - - 1 3µ i µ D σ ij • n ALCR |σ ALCR | • n ij |σ ij | • V D + + 1 3µ i µ D σ ij • n ij |σ ij | • n ij |σ ij | • V R - - 1 3µ i µ D σ ij • n ij |σ ij | • n ij |σ ij | • V L . (3.31)
The electrical potential values V A , V B , V C , V D come from the least square method (equation (3.21)). The least square weights are computed by equation (3.22).

Discretization of Intensity of Electric Field, Electron Drift Velocity and Diffusion Coefficient

The intensity of electric field is given as a negative gradient of the electrical potential E n = -grad (V n ). The intensity of electric field is computed on cell faces by diamond scheme (see equation (3.17) and (3.20)).

The electron drift velocity is a function of the intensity of electric field v n e ij = f E n ij and depends on the ratio

E n ij N [V cm 2 ],
where N is the neutral gas density (N = 2.5 • 10 19 cm -3 ). The precise formula is expressed in (B.1).

Finally we can compute the diffusion coefficient D e as a function of the intensity of electric field and the electron drift velocity

D n e ij =   0.3341 • 10 9 • E n ij N 0.54069   • v n e ij E n ij . (3.32)
Values at the cell center of gravity are estimated as average value of the values on the cell faces.

Discretization of Source Terms

The source terms depends on the electron drift velocity and the electron density, so one can write for the cell T i

S n e i = α N • v n e i • n n e i • N, S n i i = S n e i , (3.33)
where the ratio α N [cm 2 ] is computed by the following formula

if E N > 1.5 • 10 -15 , α N = 2 • 10 -16 • exp -7.248 • 10 -15 E n i /N , else, α N = 6.619 • 10 -17 • exp -5.593 • 10 -15 E n i /N . (3.34)
3.2 Some Results of Numerical Analysis in 2D

Convergence

The analysis is describe in more details and is taken from [START_REF] Martin | Élaboration de solveurs volumes finis 2D/3D pour résoudre le problème de l'élasticité linéaire[END_REF]. We suppose a case when elements of a mesh satisfied classic hypotheses of regularity:

There are constants α 1 , α 2 , β 1 , β 2 > 0 such that for all T j ∈ T h and for all e ∈ ∂T j

α 1 h 2 j ≤ µ(T j ) ≤ α 2 h 2 j and β 1 h j ≤ mes(e) ≤ β 2 h j (3.35)
where h j is a diameter of the triangle T j and µ(T j ) is its volume.

We consider a elliptic problem

-div A ∇u = f in Ω ⊂ R 2 u = g on Γ = ∂Ω (3.36)
where A(x, y) is a symmetric positive definite matrix with coefficients

a ij in C 1 (Ω), f ∈ C 0 (Ω) and g ∈ C 2 (Γ).
We can construct co-volume D σ ij around the face σ ij (see Fig. 3.1). We can discretize

∂u ∂x | σ ij by following formula ∂u ∂x | σ ij ≈ 1 µ(D σ ij ) ε∈Dσ ij 1 2 (u N 1 + u N 2 ) ε n xε dσ (3.37)
where ε is a face of the co-volume D σ ij and n xε is an axial component of its unit normal vector. N 1 and N 2 are starting and ending points of the face ε and u N 1 , u N 2 represent values of the state u in the points N 1 , N 2 .

For a node P of a mesh, we can use use linear approximation v of u in all cells which contain the node P . We can write

u P = K∈V (P ) α K (P )u K , (3.38)
where V (P ) are all triangles containing the node P , a K is a state in the center of a triangle K and α K (P ) are weights of interpolation.

Let's label u h an approximation of the exact solution u. u j is an approximation of an average value in a triangle T j .

We define exact and approached fluxes for a face e of a triangle K

Φ Exact (u) = 1 mes(e) e
A(s)p(s)

• n e ds with p = ∇u,

Φ App (u h ) = (A e p e ) n e , (3.39)
where A e is an average value of A on the face e and it can be express in local base directions n e , t e by following way

A = λ e µ e µ e ν e . (3.40)
Approximation of the gradient p e on the face e can be express in the base n e , t e like (see Fig. 3.1)

p e = u R -u L h e -α e u B -
Φ App (u h ) = λ e u R -u L h e + β e u B -u A mes(e) , (3.42)
where β e = µ e -α e λ e . Remark: if A = Id. the flux Φ App is coherent with the formula (3.17).

Definition 3.2.1 For u ∈ H 2 (Ω), we define the consistency error as a difference between the exact and approached fluxes

R e (u) = Φ Exact (u) -Φ App (π h u) (3.43)
where π h u is a L 2 projection of u on the piecewise constant functions on cells T h . The scheme is called weakly consistant if

R h (u) L 2 (Ω) -→ h→0 0 (3.44) where R h|µ(Dσ ij ) (u) = R e (u).
The finite volume discretization of the equation (3.36) can be written in following discrete form, which depends on the parameter h and tends to the exact problem

Let's find u h such that A h u h = f h in Ω and u h = g on Γ = ∂Ω where f h,i = 1 µ(T i ) T i f (x, y)dxdy.
A h is a discrete operator defined as

T i A h dxdy = e∈∂T i ǫ T i e mes(e)Φ App e (u h ) (3.45)
where ǫ T i e = ±1 if the normal vector of the face e is interior or exterior in the triangle T i .

Let's denote

|u h | 1 = e∈A h mes(e) u R -u L d e 1 2 (3.46)
where u L and u R are upstream and downstream values of u h in the direction of the normal of the face e and d e = dist(X R , X L ), and let's denote

A h u h , u h = T i ∈T h (A h u h ) i u h,i . (3.47) Definition 3.2.2 Scheme is said coercive if ∃α > 0 such as A h u h , u h ≥ α|u h | 2 1 . (3.48) Definition 3.2.

(Discrete Poincaré inequality)

We say that discrete Poincaré inequality is verified on a regular triangulation

T h if ∃β > 0 independent of h such that u h 2 L 2 (Ω) ≤ β|u h | 2 1 .
(3.49) Theorem 3.2.4 If the scheme is weakly consistant, coercive and verifies discrete Poincaré inequality, then it is convergent

u -u h 2 L 2 (Ω) -→ h→0 0. (3.50)
The following lemma gives a sufficient condition, based on the interpolation weights, for the weak consistency.

Lemma 3.2.5 Whether the interpolation weights verify for all nodes P of a mesh

(i) α K (P ) ≥ 0, (ii)
K∈V (P ) α K (P ) = 1, (iii)
K∈V (P ) α K (P ) --→ P K = 0,
then the scheme (3.37) is weakly consistant.

The least square approximation respects this conditions. More details about the analysis in [START_REF] Martin | Élaboration de solveurs volumes finis 2D/3D pour résoudre le problème de l'élasticité linéaire[END_REF].

Discrete Maximum Principle for Diamond Type scheme

The algorithm deals with discrete maximum principle based on repair technique for diamond type scheme for diffusion problem. We consider equation (3.36) and T is a set of all cells in a domain.

The maximum principle says if f (x) ≤ 0 for all x ∈ Ω than the maximum of u(x) is on the boundary u(x) ≤ min The discrete maximum principle for all cells K says if f K ≤ 0 for all K ∈ T then

U K ≤ min J∈∂Ω g J for all K ∈ T , (3.53) if f K ≥ 0 for all K ∈ T then U K ≥ max J∈∂Ω g J for all K ∈ T , (3.54)
where J is boundary midpoint and the Dirichlet boundary condition is prescribed g J = g(J).

There are discussed two repair techniques in [START_REF] Wang | Discrete Maximum Principle based on Repair Technique for Diamond Type Scheme of Diffusion Problems[END_REF]. The first one is a local repair technique and the second one is a global repair technique.

The total energy of of the diffusion equation in a discrete form is defined as

E(U ) = K∈T U K V K , (3.55)
where V K is a volume of the cell K. It is important to keep the conservation of the total energy when the discrete solution is modified because of the satisfaction of the maximum principle. We define U min = min J∈∂Ω g j as the minimum on the boundary. We require the repaired solution to have the same total energy E(U r ) = E(U ).

Local Repair Technique

The solution U is repaired cell by cell in this technique. If U K ≤ U min in the cell K, we have to add a needed energy

E K = -(U K -U min )V K to the cell K. It leads to U r K = U min . S(K)
is a set of all cell which have a common edge with the cell K. The available energy which can be transfered from a cell L (L ∈ S(k)) to the cell K is

E a L = max(0, (U L -U min )V L ). (3.56)
The total available energy coming from all neighboring cells is

E a = L∈S(K) E a L . (3.57)
There is enough available energy to correct the solution in the cell K if E a ≥ E K . In such a case we set U r K = U min and the needed energy E K from neighboring cells by following formula

U r L = U L - E a L E K E a V L . ( 3

.58)

In case that E a < E K , we enlarged the set S(K) by adding of all neighboring cells with cells in S(K).

Global Repair Technique

If we divide the set T into subsets T p = {K; U K ≥ U min }, T n = {K; U K < U min }, we can define a total available energy in the whole computational domain end a total need energy in the whole computational domain. The total available energy is

E p = K∈Tp (U K -U min )V K (3.59)
and the total needed energy is

E n = K∈Tn -(U K -U min )V K . (3.60)
We can write for each K ∈ T

E K = (U K -U min )V K K ∈ T p -(U K -U min )V K K ∈ T n . (3.61)
It is necessary to be satisfied E p ≥ E n in order to we can use the repair technique. In such a case the total energy is positive and we have enough energy which can be transport from T p to T n . When the needed energy E n is compensated, the total energy in T n is

K∈Tn U r K V K = K∈Tn U K V K + E n = = K∈Tn U K V K + - K∈Tn (U K -U min )V K = = K∈Tn U min V K . (3.62)
Because the repaired solution has to satisfy U r K ≥ U min we can see that

U r K = U min for all K ∈ T n . (3.63)
We prescribe

U r K = U K - E k En Ep V K for all K ∈ T p . (3.64)
Because E p ≥ E n , we obtain

U r K ≥ U K - E K V K = U min for all K ∈ T p . (3.65)
More details in [START_REF] Wang | Discrete Maximum Principle based on Repair Technique for Diamond Type Scheme of Diffusion Problems[END_REF].

Dynamic Mesh Adaptation

Because of the increase of the efficiency of the presented scheme, the dynamic mesh adaptation is performed to construct a mesh which is able to capture the stiff gradient of variables in areas of streamer head propagation.

There are several possibilities how to adapt a mesh. One can distinguish mesh adaptation for steady cases and on the other hands for the unsteady cases.

A progressive refinement is often use in steady computations to catch significant phenomenons such as shock waves in fluid flows. It means that we refine (split cells) a mesh in the direction of a shock wave for a once time and let the computation run. If its necessary we can refine the mesh in the same area again. Another possibility is a deformation of cells -moving of vertices and if it necessary omission of some nodes.

When an unsteady simulation run, the dynamical mesh adaptation can increase the efficiency of the computation because we need less cells. It leads to decrease of CPU time and used memory. A mesh is therefore adapted in each time step if it necessary in such way that cells are reduced in the significant area and enlarge in areas where it isn't necessary. It means when we want to reduce a cell we move its vertices closer to each other, to enlarge a cells we move its vertices apart or we omit a vertex. We can also use different way how to reduce cells. We can divide cells into sub-cells which are similar to the original ones and unrefine them when the significant phenomenon moves from this area. A quasi-stationary adaptation is also a possibility how to increase the efficiency of the code. We use a reference grid where we compute a criterion to know which cells should be refined.

Then we refine the cells, satisfy the conformity and the computation run on the refined grid for given time steps. After that we create new mesh and interpolate values of unknowns from the last refined grid to the new one.

We choose the last type of mesh refinement because the refined area is sufficient large and the streamer moves in this area for some time. The mesh is again refined before the streamer head leaves the refined area.

Criterion

The algorithm starts by computing a criterion on a fixed coarse base mesh. The criterion consists of several parameters (e.g. gradient of source terms). The criterion is scaled to an interval 0, 1 . Then the criterion is smoothed by ten iterations of diffusion equation.

After the smoothing the criterion once more scaled to the interval 0, 1 . Next the interval split on subintervals which say how many times each cell of the reference grid will be divided.

The algorithm satisfies that the level of refinement among neighboring cells cannot be higher than one. The grid conformity is applied when each triangle is subdivided. Moreover in 2D if a cell has two neighbors with a common face which should be divided on four sub-triangles, the cell will be also divided on four triangles. Mores about the adaptation criterion is in the section 5.3. After a given number of iteration, the new criterion is computed on reference grid and we create new refined mesh. The values from old refined grid a interpolated to the new one.

Mesh refinement in 2D

The refinement in 2D is performed for meshes which consist only of triangular cells. First the reference triangle is subdivided in four new triangles. Then the conformity satisfies that two neighboring triangles have only one common face.

Cell dividing

The faces midpoints are used for the cell dividing. They are used together with the cell vertices to create four new triangles from the original one. The new triangles have a level of refinement decreased of one. The algorithm is repeated until all cells have level of refinement equaled to zero (see Fig. 3.3).

Conformity

If the levels of refinement among neighboring cells are different, we have to satisfy the mesh conformity. For each cell, we check the level of refinement of neighboring cells with common face. In case that at least two neighboring cells has higher level of refinement, we increase the level of refinement for given cell and it's divided by standard algorithm. In case that only one neighboring cell has higher level of refinement, we split given cell on two triangles (see Fig. 3.3, part Level 2). This algorithm slightly enlarges the refined area.

Mesh refinement in 3D

The mesh refinement algorithm in 3D is very similar to the algorithm in 2D. The conformity step is a little bit different and will be explained later. The presented mesh refinement algorithm works only with tetrahedrons.

Cell dividing

The level of refinement for each cell is set according to a criterion as in 2D case. In one step of refinement a given cell is subdivided into eight new tetrahedrons using edges midpoints and vertices from the original cell. New tetrahedron is created at each vertex. When we imaginary cut these tetrahedrons out of the original cell, a double-pyramid remains. This double-pyramid is subsequently split into four tetrahedrons (see Fig. 3.4). • mp1, v2, mp2, mp4,

• mp6, mp4, mp5, v4,

• mp3, mp2, v3, mp5,

• mp1, mp2, mp3, mp6.

• mp1, mp4, mp6, mp2,

• mp3, mp5, mp6, mp2,

• mp2, mp4, mp5, mp6.

Conformity

If a given cell shouldn't be divided, we have to satisfy conformity of mesh. We check all neighboring cells common at least on edge if there are cells which has been divided. There are several cases how to split the cell according number of edges which contains divided neighboring cells. We split the given cell into two new tetrahedrons if only one edge contains divided neighbors (see Fig. 3.5). • mp1, v2, v3, v4.

If each of six edges of a given cell should be divided, the cell is split as it is regular cell described in section 3.3.3.1.

In case that more than one but less then six edges should be split, we firstly divide the given cell into four tetrahedrons using the cell center of gravity C (see Fir. 3.6). It leads to set of tetrahedrons

• v1, v2, C, v4,
• v1, v2, v3, C,

• v2, v3, C, v4,

• v3, v1, C, v4. • v1, mp1, C, mp2,

• mp2, mp1, C, v4,

• mp1, v2, C, v4. Last there is possibility that three edges of the tetrahedron v1, v2, C, v4 should be split. This gives four new tetrahedrons (see Fig. 3.9)

• v1, m1, C, mp2,

• mp2, mp1, C, mp3,

• mp1, v2, C, mp3,

• mp2, mp3, C, v4. This way of satisfying of conformity doesn't expand the area of refinement. On the other hand, we obtain worse shape of cells. The conformity could be improved by the proposal in [START_REF] Bey | Tetrahedral Grid Refinement[END_REF] and [START_REF] Liu | Quality Local Refinement of Tetrahedral Meshes Based on 8-Subtetrahedron Subdivision[END_REF].

Chapter 4

Simple Tests of Grid Adaptation Algorithm

Linear Advection Equation, Grid Adaptation

These section deals with testing of the advection together with the mesh adaptation on a simple scalar equation

u t + a • ∇u = 0, (4.1)
in one, two and three dimensions. The velocity vector a is a = 1 in an one dimensional case, a = (1, 0) in a two dimensional case and a = (1, 0, 0) in a three dimensional case. The numerical test are as for 1 st order of accuracy as for 2 nd order of accuracy scheme. The simple upwind scheme (see eq. (3.6)) is used in the 1 st order scheme. The extension of the upwind scheme on the second order in space is made by Van Leer's type MUSCL algorithm together with Barth-Jespersen limiter (see eq. (3.7)), more in [START_REF] Blazek | Computational Fluid Dynamics: Principles and Applications[END_REF]. The limiter is also used e.g. in [START_REF] Dobes | Numerical solution of single and multi-phase internal transonic flow problems Int[END_REF]. The second order of accuracy in time is due to three steps Runge-Kutta method.

The testing domains are

• an interval x ∈ 0, 1 for 1D case,

• a rectangle 0, 1 × 0, 0.5 for 2D case,

• a brick 0, 1 × 0, 0.5 × 0, 0.5 for 3D case.

Initial condition

The initial condition is represented by Gaussian pulse with maximum magnitude of 10 u(x, y, z, 0) = 10 • e -(x-0. The initial condition 4.2 is written for 3D case where the initial Gaussian pulse is a sphere. The initial condition for 2D and 1D cases are derived from 3D initial condition by omitting z variable in 2D (y, z variables in 1D). The 3D case is also tested for cylindrical initial Gaussian pulse and in such case the initial condition is same as for 2D case.

Boundary conditions

For all test cases, the homogeneous Neumann's conditions are used at all boundaries.

Test Results

All presented results are at physical time t = 0.6 s. We can see results obtained by 1D computation in Figs. 4.1-4.2. Fig. 4.1 shows results computed by 1 st order scheme and 2 nd order scheme with CF L number equaled to 0.4. There is more numerical diffusion in the 1 st order scheme which leads to decrease of maximum value of the initial Gaussian pulse and its spreading during transportation. When the 2 nd order scheme is used, the maximum value almost doesn't decrease and the initial pulse doesn't spread. shows comparison of 1 st order scheme and 2 nd order scheme along the axis of computational domain y = 0.25. We can again see that the first order scheme is more dissipative than the second order scheme. Fig. 4.4 presents isolines of the same results. almost the same results. The Barth-Jespersen limiter is unlike the Van Albada limiter independent on chosen parameters. The Bart-Jespersen limiter is therefore more suitable to use because we don't have to "tune" any parameter if a computation condition is changed.

Chapter 5

Numerical Results for 2D Planar Discharge (5.1)

The initial Gaussian pulse for the electron and ion densities creates disturbance in the electric field which is necessary for initialization of ionization process in the wave propagation.

The background electron and ion densities 10 9 cm -3 substitute the photoionization effect which is neglected in this simple streamer model. The computation is stopped at prescribed time when the head of electron avalanche is still inside the computational domain. In further computation, the discharge goes out through the right boundaries and the source terms increase to infinity but this effect doesn't represents physical behavior of the electric discharge.

The second domain (domain with parabolic left boundary) has different initial conditions n e (x, y, 0) = 10 9 [cm -3 ], n i (x, y, 0) = n e (x, y, 0) .

(5.

2)

The disturbance in the electric field is created by the shape of the left boundary so we don't need any initial Gaussian pulse.

Boundary Conditions

The left boundaries of the domains are anode, the right boundaries are plane cathodes. There are prescribed homogeneous Neumann conditions for all unknowns for upper and lower boundaries. Following boundary conditions are prescribed for the anode and the cathode ∂n e ∂ n = 0, for anode and cathode,

V = 25000 [V], for anode, V = 0 [V],
for cathode.

(5.3)

Discharge in Rectangular Domain

Criterion Variables

The first row in Fig. 5.5 shows results, when the magnitude of the source term is chosen as criterion function crit = S e (case 1). The second row shows results of case 2, when the criterion is based on gradient of the electron density n e , crit = |gradn e |. Finally we present results of case 3 , when both previous criterion are used simultaneously crit = max (|grad n e |, S e ). We use this expression in two different algorithms V1 and V2.

In the version V1 (row 3 in Fig. 5.5), both criteria are scaled to one (dividing by the maximum value respectively). Then for each cell is chosen maximum of these two scaled criteria. The new criterion is smoothed using diffusion equation and scaled once more to one and it gives the final new criterion.

In the version V2 (see Fig. 5.6), each criterion is smoothed and scaled to one separately and the new criterion is maximum of these two values.

Discussion of Results

We can observe a significant influence of choice of criterion variable on position of head of the planar wave in fixed time T and consequently on discharge propagation velocity. The The table 5.2 shows, that the highest number of grid cells is also in case 3, version V2.

Comparing to case 1 it uses twice more cells, the number of time steps is almost the same and CPU time is 50% increasing. The maximum of grid points (about 22 000) confirms high efficiency of dynamic grid adaptation algorithm for solved problem. When all triangles of reference grid will be adapted up to fourth level, the number of volumes overcomes 1 200 000.

Distribution of Subintervals of Adaptation Levels

The results in first row in Fig. 5.6 as well as all results presented in Fig. 5.5 were computed by the same basic definition of subintervals of 0, 1 of different levels of adaptation -level 4 for crit ∈ 0.64, 1 , level 3 for crit ∈ 0.36, 0.64), level 2 for crit ∈ 0.16, 0.36), level 1 for crit ∈ 0.04, 0.16). The results presented in the second row in Fig. 5.6 were achieved by the same algorithm, but with enlarged interval of highest level of adaptation -level 4 for crit ∈ 0.36, 1 , level 3 for crit ∈ 0.16, 0.36), level 2 for crit ∈ 0.08, 0.16), level 1 for crit ∈ 0.04, 0.08).

Discussion of Results

We can observe higher density of grid in plotted regions, but also very similar distribution of n e as well as position and shape of source term. The number of grid volumes is two times higher as for basic distribution, number of time steps is almost the same and CPU is 1.7 time higher.We can conclude, that the quality of achieved result is not too sensitive on choice of subinterval of highest level of adaptation.

Grid Density and Numerical Convergence Test

The influence of the grid density is studied for both 1 st and 2 nd order of approximation in space variables. The test cases are constructed with the same way as in previous subsection.

Test Description

The solution at time t in = 3.14 • 10 -8 s (see Fig. 5.4) is once more chosen as initial approximation. The computational grid is obtained as uniformly refined reference grid, i.e. all triangles of reference grid are split into the same number of 4 lr triangles in dependence on level of refinement lr = 2, • • • , 5. The computation starts with generation of grid and conservative interpolation of unknowns at starting time into volumes of grid. Then the time evolution is calculated up to the prescribed time. We use the same CFL number for all cases.

Uniformly Refined Grid Fig. 5.7 shows the isolines of net charge density ρ = e (n i -n e ) obtained by first order scheme. The isolines of net charge density computed by second order scheme are plotted in Fig. 5.8. 1D plots of the electron density n e and the net charge density ρ along the axis of discharge (y = 0.25) document the influence of the grid density (see Fig. 5.9 and Fig. 5.10). Presented results approve that the first order scheme over predicts the electric discharge propagation speed. This property is in accordance with the conclusion in [START_REF] Célestin | Etude de la dynamique des streamers dans l'air a la pression atmospherique[END_REF]. We can observe in both Fig. 5.7, and Fig. 5.9 that the position of the discharge head as well as the magnitude of the electron density decreases with higher grid density. The magnitude of net charge density is increasing. The distribution of net charge also shows, that the difference between results on neighboring coarse and fine grid slowly decreases with increasing of gird density (position and magnitude of maximum). The method with second order scheme in space variables has different behavior. The position of the discharge head and magnitude of electron density increases and magnitude of net charge decreases when the grid density is increasing. We can also observe much smaller difference between results (net charge distribution in Fig. 5.10) on neighboring coarse and fine grid as in case of first order scheme. This property also in some sense document higher accuracy and convergence of second order scheme. We can see that the first order scheme is still too dissipative and inaccurate even if very fine mesh is used. The Fig. 5.10 shows that 4 levels of mesh refinement is sufficient to obtain satisfactory accurate results with second order scheme.

Comparison of Uniformly and Dynamically Adapted Grid

Finally the accuracy and efficiency of proposed dynamic mesh adaptation algorithm is tested.

The results are compared with the numerical results obtained on uniformly refined grid with 4 levels of refinement. These comparisons are plotted for both the fist and the second order scheme in form of the electron density n e and the net charge density ρ along the axis (Fig. 5.12 and Fig. 5.13). The comparison is also shown in form of isolines for the net charge density in Fig. 5.7 and Fig. 5.8 where the third column is uniformly refined grid and the fourth column is dynamically refined grid. We can see that the results are almost identical. It confirms that the accuracy of uniformly refined grid is conserved. The table 5.3 shows the high efficiency of the dynamic mesh adaptation algorithm. Figs. 5.14-5.18 present time evolution planar wave in 2D domain where the anode has a parabolic shape. The disturbance in the electric field is given by the shape of the anode and there is no initial Gaussian pulse. We can see in the figures evolution of the electron density, the net charge density and the source terms.

Chapter 6

Numerical Results for 3D Streamer

Description of Test Cases

The best method from 2D simulations is adopted to the 3D experiments. It means treatment with the choice of parameters for the adaptation criterion, smoothing of criterion, using of the second order scheme. The algorithm V2 (see section 5.3) is used in adaptation criterion with 3 levels of refinement -level 3 fro crit ∈ 0.5625, 1 , level 2 for crit ∈ 0.25, 0.5625), level 1 for 0.0625, 0.25). The lower level of refinement is because of the CPU time and also the computer memory. The 3D simulations take almost 3 weeks to obtain results even with the dynamcal mesh adaptation. There are two different domains for experiments in3D. The first of them is a simple brick 1 × 0.5 × 0.5 cm (see Fig. 6.1). + 10 8 [cm -3 ], σ = 0.01, n i (x, y, z, 0) = n e (x, y, 0) . (

The meaning of initial Gaussian pulse and the background electron and ion densities is same as in 2D. The increase of background electron and ion densities comes from axisymmetrical simulations [START_REF] Trdlička | Application of Finite Volume Method for Numerical Simulation of Streamer[END_REF] where a low magnitude of the densities have made blow up of the solution.

Boundary Conditions

The left boundary is an anode with high voltage, right boundary is a cathode. The homogeneous boundary conditions are prescribed for all unknowns on upper, lower front and back boundaries. Boundary conditions for the anode and the cathode are given by following formulas ∂n e ∂ n = 0, for anode and cathode,

V = 25000 [V], for anode, V = 0 [V], for cathode. (6.2)
The second domain ("pyramid") has half voltage on the anode

V = 12500 [V]
, for the anode. (6.3)

Streamer Propagation

Figures in this section shows results coming from computation in the domain where the initial disturbance of the electric field is due to the initial Gaussian pulse. The results are similar to the 2D computations. It demonstrate that we can trust to the method also in 3D simulations.

The mesh has been generated by GMSH mesh generator. The generator produced a mesh where the ratio between the smallest and the biggest cells is 1 : 10. The total number of cells in used mesh is around 150000. During the tuning of the code, it has been found out that the streamer simulation is sensitive on a bad quality of an used mesh. The old version of the GMSH generator produces grids where the ratio between the smallest and the biggest cells is 1 : 10 4 . When such grid is used, the method diverges.

The computation takes three weeks of CPU time. The number of iterations needed to obtain the results at the time t = 2.52 • 10 -8 s is around 20500. The dynamical mesh adaptation leads to a mesh with around 450000 cells at the final time t = 2.52 • 10 -8 s.

We can see a time evolution of the electron density in Figs. The results seem to be a little bit more oscillatory in the electron density n e . 2D test presents that the numerical method has such behavior on coarse grids. The oscillatory can be also caused by a slightly different refinement algorithm in 3D and there is also different conformity step in the 3D. The cell dividing is different in 2D and 3D from the nature of the geometrical shapes of cells. A triangle is divided on the four same sub-triangles similar to the original one. The tetrahedron is divided on eight sub-tetrahedrons where four of them (tetrahedrons at vertices) are similar to the original one but the rest of the sub-tetrahedrons have different shape. This property together with the conformity step can deform the mesh. These properties will be tested in the future. The domain with the "pyramid" cathode is next used for simulations of the streamer branching where a plasma spot is added to disturb the electric field (see section 6.4). The 2D results aren't included to the thesis.

Streamer Branching

The tested streamer branching is caused by an instability in the electric field. The effect of plasma spot [START_REF] Papageorgiou | Three-dimensional numerical modelling of gas discharges at atmospheric pressure incorporating photoionization phenomena[END_REF] is applied near the head of the streamer in this work. We add neutral plasma spot which is introduced into the model as another source term into the continuity equations for the electron and the positive ions densities. The plasma spot is canceled after short time (less than nanosecond).

Figs. 6.16-6.107 show results of four test versions when one or two plasma spots were added. Tests are computed in the "pyramid" domain. The figures are zoomed around the streamer paths.

The plasma spot is added at time t = 1.26 • 10 -8 s and canceled after the "duration time". In Version 3 and Version 4 , the second plasma spot is added when the first one is canceled. An evolution of the streamer branching, where the plasma spot is introduced in a double distance from the streamer filament axis, is presented in Figs. [START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].37 and 6.56. The evolution of the electron avalanche is in Figs. [START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].37-6.46. Figs. 6.47-6.56 show the evolution of the net charge density. As in the previous case, the first row in the figures show the streamer just in the moment when the plasma spot is applied to the electric field. The moment when the plasma spot is abolished is presented in the second row. The electrons from the streamer filament are attached to the plasma spot area and we can see that the electrons also continue in the straight direction. We can observe a branch coming from the main filament in the third row. The fourth and fifth rows show that the main streamer discharge still propagates toward the cathode while the propagation in the branch almost stops. An adapted mesh in the planar cut at the time t = 1.89 • 10 -8 s is plotted in Fig. 6.57. The mesh size is 660000 cells and the computation has made 31000 time steps. The results, where two plasma spots are sequentially added, are plotted in Fig. 6.58-6.81. Both plasma spots lie above the streamer filament. The electron avalanche evolution is presented in Figs. present a moment when the first plasma spot is voided and the second plasma spot is added to the electric field. Figs. 6.60, 6.66, 6.72 and 6.78 display a moment when the second plasma spot is abolished. We can see the electron motion toward to the first plasma spot place. We observe the progression toward the second spot and the main filament slightly deflects from the straight direction. The last figures in the sets show following the streamer evolution. Two plasma spots on the same side from the main streamer filament produce a longer branch than the case where only one plasma spot is added. On the other hand the propagation velocity of the main filament is slower than in the case with one plasma spot. The adapted mesh around the branched streamer at the time t = 1.89 • 10 -8 s is plotted in Fig. 6.82. There is 570000 cells in the mesh. There has been needed 28900 time steps to obtain the results at the time t = 1.89 • 10 -8 s. The last test simulate conditions when two independent plasma spots are added to the electric field in two different times. The first spot is placed above the main streamer filament and the second plasma spot is placed under the streamer filament. Results of the test are plotted in Figs.

Conclusions

The goals of the thesis have been fulfilled. The work deals with the streamer propagation which is an electric discharge in a high voltage electric field. The streamer branching has been successfully simulated by introducing of temporal disturbances in the electric field. We have simulated a change of streamer direction, a singular branching and a multiple branching, see the section 6.4. The shape and position of streamer branches depend on position of plasma seeds (local disturbances) and also on the time in which the disturbances influence the electric field. We have referred about the streamer branching in [START_REF] Benkhaldoun | Trdlička: A Full 3D dynamically adaptive unstructured grid mesh finite volume approach to simulate multiple branching in streamer propagation[END_REF] and we haven't seen successful simulation of the streamer branching in full 3D coordinates yet.

Results in chapter 6 demonstrate the applicability of the developed method. The results present a streamer propagation in a simple 3D domain. Computation of a one 3D test takes around three weeks and the computer memory costs are also very high even if the dynamical mesh adaptation is used. That's the reason why only three levels of refinement are taken into account in the calculations although 2D tests show that at least four levels of the mesh refinement should be considered. The properties of the method have been tested on 2D problems (see chapter 5). The tests of the adaptation criterion in the section 5.3 reveal that the grid density around the streamer (planar discharge) head and the streamer path has influence on the velocity of the propagation.

Therefor we need a fine mesh in the computational domain at these places. We also found an optimal adaptation level (as compromise of the grid size and quality of results) and confirm, that results achieved on the adapted grid are the same as results of a static fine grid in the whole domain (see Figs. 5.12, 5.13). The dynamical mesh adaptation decreases the CPU time almost 60 times. We have released findings about the influence of the dynamical mesh adaptation and the proper choice of the refinement criterion in [START_REF] Benkhaldoun | Dynamical grid adaptation method for model plasma problem[END_REF], [START_REF] Benkhaldoun | Simulation of Planar ionization wave Front Propagation on Unstructured Adaptive Grid[END_REF] and [START_REF] Benkhaldoun | Trdlička: Numerical Study of Planar Discharge Motion[END_REF]. We don't deem that the development of the method is finished. The parallelization should be improved and maybe more effective algorithm for solving of a liner systems should be included.

A general platform based on C++ programing have been used and successfully tested. It can solve different physical problems. The streamer motion has been solved by the platform in this thesis. It's been used for a shallow water system at University Paris 13 in thesis [START_REF] Sari | Modélisation mathématique et numérique de systèmes multicouches avec transport de sédimets pour les problèmes à surface libre[END_REF].

The platform is also utilized at Czech Technical University for fluid mechanics computation where an implicit method has been included.

A.2 Relations for Dissipative Terms

Weights parameters for the 2D gradient computation in diffusion terms: (A.12)

Appendix B B.1 Material Functions for the Air

The relations for the material functions are taken from [START_REF] Morrow | Streamer propagation in air[END_REF]. η 2 is two-body attachment coefficient, η 3 is three-body attachment coefficient.

Electron
Recombination coefficients β ep , β np (there is only β in [START_REF] Morrow | Streamer propagation in air[END_REF]):

β = 2 • 10 -7 [cm 3 • s -1 ] (B.

  ∂n e ∂t + ∇ • (n e v e ) -∇ • (D e ∇n e ) = S ph + n e α| v e | -n e β att | v e | -n e n p β ep , ∂n p ∂t + ∇ • (n p v p ) = S ph + n e α| v e | -n e n p β ep -n n n p β np , ∂n n ∂t

  3.1). The diamond cell is constructed by connection of centers of gravity (L, R) of cells T i , T j which share the face σ ij and its endpoints A, B. One obtains the co-volume D σ ij by this construction.

Figure 3 . 1 :

 31 Figure 3.1: Diamond cell in 2D

  a unit normal vector of the co-volume face σ AR and |σ AR | is its length. The others co-volume faces and their normal vectors are labeled analogically (see Fig. 3.1).

Figure 3 . 2 :

 32 Figure 3.2: Diamond cell in 3D

  x∈∂Ω g(x), for all x ∈ Ω, (3.51) if f (x) ≥ 0 for all x ∈ Ω than the minimum of u(x) is on the boundary u(x) ≥ max x∈∂Ω g(x), for all x ∈ Ω. (3.52)

Figure 3 . 3 :

 33 Figure 3.3: Mechanism of multi-level refinement in 2D

Figure 3 . 4 :

 34 Figure 3.4: Mechanism of refinement in 3D

Figure 3 . 5 :

 35 Figure 3.5: Conformity when one edge should be split

Figure 3 . 6 :

 36 Figure 3.6: Splitting of tetrahedron using the center of gravity

Figure 3 . 7 :

 37 Figure 3.7: One edge to be split

Figure 3 . 8 :

 38 Figure 3.8: Two edge to be split

Figure 3 . 9 :

 39 Figure 3.9: Three edge to be split

2 0.05 2 .

 22 2) 2 +(y-0.25) 2 +(z-0.25)(4.2)

Figure 4 . 1 : 4 Figure 4 . 2 :Fig. 4 .

 414424 Figure 4.1: 1D computation: 1 st order scheme (left), 2 nd order scheme (right), CF L = 0.4

Figure 4 . 3 : 4 Figure 4 . 4 : 4 Fig. 4 .

 4344444 Figure 4.3: 2D computation: 1 st order scheme (left), 2 nd order scheme (right), CF L = 0.4

Figure 4 . 5 : 1 Fig. 4 .

 4514 Figure 4.5: 2D computation: 1 st order in time, 2 nd order in space, CF L = 0.1

Figure 4 . 6 :

 46 Figure 4.6: 3D computation: 2 nd order scheme CF L = 0.4, structured grid 1200000 cells (left), unstructured grid 60000 cells (right)

Figure 4 . 7 :Figure 4 . 8 :

 4748 Figure 4.7: 3D computation: 2 nd order scheme CF L = 0.4, dynamical mesh adaptation, transportation of a sphere, isolines (left), mesh (right)

Figure 4 . 9 :

 49 Figure 4.9: Comparison of results from different dimensions

4. 1 . 2

 12 Comparison of limitersThis section compare results obtained by 2D computations where Barth-Jespersen and Van Albada limiters are used. The initial condition (4.2) is slightly modified u(x, y, z, 0) = 10 16 • e -(x-0.2) 2 +(y-0.25) 2 +(z-0.25) 2 0.05 2

( 4 . 3 )

 43 to verify that we can use such orders of magnitude.

Figure 4 . 10 :

 410 Figure 4.10: Comparison of limiters

5. 1

 1 Description of Test CasesWe consider the propagation of 2D problem of discharge in homogeneous electric field described by equations (2.1)-(2.2). The first domain is a simple rectangle 1 × 0.5cm, the second domain has a left boundary partly replaced by a parabola. The parabola goes through points [0, 0.15] , [0.2, 0.25] and [0, 0.35] (see Fig. 5.1).

Figure 5 . 1 : 2 + 10 9

 51210 Figure 5.1: 2D computational domains

Figs. 5 .

 5 Figs. 5.3 and 5.2 present results obtained with the second upwind scheme on dynamically adaptive grid in rectangular domain. The shape and magnitude of velocity vectors agrees with physical mechanism of the electric discharge propagation. The results in front of the scale in Fig. 5.3 are rotated about 90 • in positive sense.

Figure 5 . 4 :

 54 Figure 5.4: Reference mesh, mesh, electron density, electric field at t in = 3.14 • 10 -8

Figure 5 . 5 :

 55 Figure 5.5: Criterion based on source terms (first row), criterion based on gradient of electron density (second row), criterion based on source terms and gradient of electron densityalgorithm V1 (third row)

Figure 5 . 6 :

 56 Figure 5.6: Criterion based on source terms and gradient of electron density -algorithm V2, basic distribution of level of adaptation(first row), enlarged interval of higher level of adaptation (second row)

Figure 5 . 7 : 1 Figure 5 . 8 : 2 Figure 5 . 9 : 1

 571582591 Figure 5.7: 1 st order, comparison of net charge density computed on different uniformly refined meshes at time t = 5.25 • 10 -8 s. From left to right: 2 levels, 3 levels, 4 levels, adaptation (4 levels)

Figure 5 . 2 Figure 5 . 11 :

 52511 Figure 5.10: 2 nd order, comparison of electron density (left) and net charge density (right) computed on different coarse meshes at time t = 5.25 • 10 -8 s.t: 2 levels, 3 levels, 4 levels

Figure 5 .

 5 Figure 5.12: 1 st order, comparison of electron density (left) and net charge density (right) computed on fine mesh and mesh with dynamic adaptation

Figure 5 . 14 :

 514 Figure 5.14: Electron density (left), net charge density (middle) and source term (right) at time t = 1.155 • 10 -8 s.

Figure 5 . 15 :

 515 Figure 5.15: Electron density (left), net charge density (middle) and source term (right) at time t = 1.47 • 10 -8 s.

Figure 5 . 16 :

 516 Figure 5.16: Electron density (left), net charge density (middle) and source term (right) at time t = 1.89 • 10 -8 s.

Figure 5 . 17 :

 517 Figure 5.17: Electron density (left), net charge density (middle) and source term (right) at time t = 2.1 • 10 -8 s.

Figure 5 . 18 :

 518 Figure 5.18: Electron density (left), net charge density (middle) and source term (right) at time t = 2.31 • 10 -8 s.

Figure 5 . 19 :

 519 Figure 5.19: Evolution of electron density (left) and net charge density (left).

Figure 5 . 20 :

 520 Figure 5.20: Evolution of source term.

Figure 6 . 1 :

 61 Figure 6.1: 3D computational domain

Figure 6 . 2 :

 62 Figure 6.2: 3D computational domain -pyramid

6 . 3 - 6 . 5 .

 6365 The figures present isolines of the electron density in a cut plane y = 0.25 (left part of the figures) and full 3D view (right part of the figures).A time evolution of the net charge density is shown in Figs. 6.6-6.8. As in previous set of figures, the left part of the figures presents isolines in the cut plane y = 0.5 and the right part of the figures shows full 3D view.

Figure 6 . 3 :Figure 6 . 4 :Figure 6 . 5 :Figure 6 . 6 :Figure 6 . 7 :Figure 6 . 8 :

 636465666768 Figure 6.3: Electron density at time t = 1.155 • 10 -8 s

Figure 6 . 9 :Figure 6 . 10 :Figure 6 . 11 :Figure 6 . 12 :Figure 6 . 13 :Figure 6 . 14 :

 69610611612613614 Figure 6.9: Electron density at time t = 0.945 • 10 -8 s

Figure 6 . 15 :

 615 Figure 6.15: Electron density at time t = 1.995 • 10 -8 s (left), electron density at time t = 1.89 • 10 -8 s in "pyramid" domain(right)

Version 1

 1 Figs.[START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].16-6.35 present an evolution of the electron density and the net charge density where the neutral plasma spot is introduced too close the streamer head. Figs. 6.16-6.20 and 6.26-6.30 show isolines in a cut plane y = 0.25, Figs. 6.21-6.25 and 6.31-6.35 show full 3D view. Figs. 6.16, 6.21, 6.26 and 6.31 show the state when the plasma spot is added to the electric field. Figs. 6.17, 6.22, 6.27 and 6.32 captures a moment when the plasma spot is canceled. We can observe that the electron from the streamer head are attached to the area of the plasma spot. Figs. 6.18, 6.23, 6.28 and 6.33 show that the streamer motion is deflected to the area of the plasma spot. Figs. 6.19, 6.20, 6.24, 6.25, 6.29, 6.30, 6.34 and 6.35 display that the streamer avalanche next moves through the the area of the plasma spot toward the cathode. We can see that introducing of the plasma spot too close the streamer head causes only a deflection of the streamer motion. Fig.6.36 shows an adapted mesh around the streamer path in the cut plane y = 0.25 at the time t = 1.89•10 -8 s. The mesh has around 440000 cells at the final time and it's been needed 28500 time steps to obtain the results.

Figure 6 . 16 :Figure 6 . 17 :Figure 6 . 18 :Figure 6 . 19 :Figure 6 . 20 :Figure 6 . 21 :Figure 6 . 22 :Figure 6 . 23 :Figure 6 .Figure 6 .Figure 6 .Figure 6 . 2

 61661761861962062162262366662 Figure 6.16: Electron density at times t 1 = 1.26 • 10 -8 s, isolines in a cut plane y = 0.25

Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .

 66666666 Figure 6.37: Electron density at times t 1 = 1.26 • 10 -8 s, isolines in a cut plane y = 0.25

Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 . 4 1

 6666666664 Figure 6.58: Electron density at times t 1 = 1.26 • 10 -8 s, isolines in a cut plane y = 0.25

Fig. 6 .

 6 107 shows the adapted mesh at the time t = 1.89 • 10 -8 s. The mesh contains around 540000 cells. We have needed 28300 time steps to obtain the results at the final time.

Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .Figure 6 .

 6666666666666 Figure 6.83: Electron density at times t 1 = 1.26 • 10 -8 s, isolines in a cut plane y = 0.25

  -x A ) , R y = N (A) p=1 (y p -y A ) , I xx = N (A) p=1 (x p -x A ) 2 , I yy = N (A) p=1 (y p -y A ) 2 , I xy = N (A) p=1 (x p -x A ) (y p -y A ) , λ x = I xy R y -I yy R x D , λ y = I xy R x -I xx R y D , D = I xx I yy -I 2 xy .(A.5)3D computation of the gradient in the diffusion terms is given by following relation∇n eij = 1 µ D σ ij1 3 (n e (A) + n e (B) + n e (R)) n ABR |σ ABR | + + 1 3 (n e (B) + n e (C) + n e (R)) n BCR |σ BCR | + + 1 3 (n e (C) + n e (D) + n e (R)) n CDR |σ CDR | + + 1 3 (n e (D) + n e (A) + n e (R)) n DAR |σ DAR | + + 1 3 (n e (A) + n e (B) + n e (L)) n BAL |σ BAL | + ((B) + n e (C) + n e (L)) n CBL |σ CBL | + + 1 3 (n e (C) + n e (D) + n e (L)) n DCL |σ DCL | + + 1 3 (n e (D) + n e (A) + n e (L)) n ADL |σ ADL | ,where n e (A) (n e (B),• • • ) is the electron density value in the point A (B, • • • ), n ABR ( n BCR , • • • ) is an unit outward vector of the face σ ABR (σ BCR , • • • ) and |σ ABR | (|σ BCR |, • • • ) is its measure.When we factor one third out of the square bracket and modify expressions inside the bracket, we obtain∇n eij = 1 3µ D σ ij ( n ABR |σ ABR | + n DAR |σ DAR | + n BAL |σ BAL | + n ADL |σ ADL |) n e (A) + + ( n ABR |σ ABR | + n BCR |σ BCR | + n BAL |σ BAL | + n CBL |σ CBL |) n e (B) + + ( n BCR |σ BCR | + n CDR |σ CDR | + n CBL |σ CBL | + n DCL |σ DCL |) n e (C) + (A.7) + ( n CDR |σ CDR | + n DAR |σ DAR | + n DCL |σ DCL | + n ADL |σ ADL |) n e (D) + + ( n ABR |σ ABR | + n BCR |σ BCR | + n CDR |σ CDR | + n DAR |σ DAR |) n e (R) + + ( n BAL |σ BAL | + n CBL |σ CBL | + n DCL |σ DCL | + n ADL |σ ADL |) n e (L) .The vector equilibrium impliesn ABR |σ ABR | + n DAR |σ DAR | + n BAL |σ BAL | + n ADL |σ ADL | = n BRDL |σ BRDL |, n ABR |σ ABR | + n BCR |σ BCR | + n BAL |σ BAL | + n CBL |σ CBL | = n ALCR |σ ALCR |, n BCR |σ BCR | + n CDR |σ CDR | + n CBL |σ CBL | + n DCL |σ DCL | =n RBDL |σ RBDL |, n CDR |σ CDR | + n DAR |σ DAR | + n DCL |σ DCL | + n ADL |σ ADL | =n ALCR |σ ALCR |, n ABR |σ ABR | + n BCR |σ BCR | + n CDR |σ CDR | + n DAR |σ DAR | = n ij |σ ij |, (A.8) n BAL |σ BAL | + n CBL |σ CBL | + n DCL |σ DCL | + n ADL |σ ADL | =n ij |σ ij |.When we substitute the equations (A.8) in the equation (A.7), we can write for the electron density gradient ∇n eij∇n eij = 1 3µ D σ ij (n e (A) -n e (C)) n BRDL |σ BRDL | + (A.9) + (n e (B) -n e (D)) n ALCR |σ ALCR | + (n e (R) -n e (L)) n ij |σ ij | .Parameters in the weights for gradient computation are given by following relationsI xx = N (A) p=1 (x p -x A ) 2 , I xy = N (A) p=1 (x p -x A ) (y p -y A ) , R x = N (A) p=1 (x p -x A ) , I yy = N (A) p=1 (y p -y A ) 2 , I xz = N (A) p=1 (x p -x A ) (z p -z A ) , R y = N (A) p=1 (y p -y A ) , I zz = n(A) p=1 (z p -z A ) 2 , I yz = N (A) p=1 (y p -y A ) (z j -z v ) , R z = N (A) p=1 (z p -z A ) ,(A.10) D = I xx I yy I zz + 2I xy I xz I yz -I xx I 2 yz -I yy I 2 xz -I zz I 2 xy , (A.11) 107 λ x = I 2 yz -I yy I zz R x + (I xy I zz -I xz I yz ) R y + (I xz I yy -I xy I yz ) R z D , λ y = (I xy I zz -I xz I yz ) R x + (I 2 xz -I xx I zz ) R y + (I yz I xx -I xz I xy ) R z D , λ z = (I xz I yy -I xy I yz ) R x + (I yz I xx -I xz I xy ) R y + I 2 xy -I xx I yy R z D .

8 )

 8 The neutral gas density N = 2.5 • 10 19 [cm -3 ] is use in relations (B.1), (B.3) -(B.7).

  is the time, n e denotes the electron density, n i the positive ion density, v e the electron drift velocity, D e the diffusion coefficient. S e = S +

	is coupled with the Poisson's equation for the electric potential V
								∆V = Q,	(2.2)
	where Q = -constant, e the electron charge. The intensity of the electric field is computed as a negative e ǫ (n i -n e ) (2D tests were computed with Q = -e ǫ (n i -n e ) 10 4 ), ǫ is the dielectric
	gradient of the electric potential
								E = -∇V.	(2.3)
	The electron drift velocity v e is a function of the electric field E and depends on the ra-
	tio || E||/N where N is the neutral gas density (N = 2.5 • 10 19 cm -3 ) [5]
	for	E N	> 2 • 10 -15 ,		v e = -7.4 • 10 21 •	E N	+ 7.1 • 10 6 •	E E	,
	for 10 -16 <	E N	≤ 2 • 10 -15 ,	v e = -1.03 • 10 22 •	E N	+ 1.3 • 10 6 •	E E	,
	for 2.6 • 10 -17 <	E N	≤ 10 -16 , v e = -7.2973 • 10 21 •	E N	+ 1.63 • 10 6 •	E E	,
	for	E N	≤ 2.6 • 10 -17 ,	v e = -6.87 • 10 22 •	E N	+ 3.38 • 10 4 •	E E	.
									is taken into account in our work .
	It consists of a convection-diffusion-reaction for the electron density, an ordinary differential
	equation for the positive ion density coupled by the Poisson's equation for the elecric potential
							∂n e ∂t	+ div n e v e -D e ∇n e = S e ,
									∂n i ∂t	= S + i ,	(2.1)
	where t i are source terms. The system (2.1)

Table 4 .

 4

	1: 1D computation, 1000 cells
		max	min
	1 st order	8.176	0
	2 nd order	9.9433 -3e -8
	1/2 order CF L = 0.1 10.021 -4e -8

Table 4 .

 4

	2: 2D computation, 775x387 cells
		max	min
	2 nd order, grid 775x387x4	9.970 -1e -5
	2 nd order, unstructured grid 60000 cells,	5.1348 -0.19
	2 nd order, unstructured grid 60000 cells, sphere 3 levels 9.7774 -0.002

Table 4 .

 4 3: 3D computation Tables 4.1-4.3 present maximum and minimum values of the results at time t = 0.6 s. 1/2 order in tables means 1 st order in time and 2 nd order in space.

  6.58-6.69 and the net charge evolution is in Figs. 6.70-6.81. The introducing of the first plasma spot is shown in Figs. 6.58, 6.64, 6.70 and 6.76. Figs. 6.59, 6.65, 6.71 and 6.77

  [START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].83-6.106. Figs.[START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].83, 6.89, 6.95 and 6.101 have the same meaning as in previous cases. Figs.[START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].84, 6.90, 6.96 and 6.102 show the streamer when the first spot is canceled and the second spot is introduced. There is a state when the second plasma spot is abolished in Figs.[START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].85, 6.91, 6.97 and 6.103. We see that the first branch begins here. The beginning of the second branch is plotted in Figs.[START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF].86, 6.92, 6.98 and 6.104. The last figures in the sets present the progression of the main streamer filament and the branches.

  drift velocity v e :

	Ionization coefficient α:
		if else, N E	> 1.5 • 10 -15 ,	α N α N	= 2 • 10 -16 • exp = 6.619 • 10 -17 • exp -7.248 • 10 -15 E /N E /N -5.593 • 10 -15	[cm 2 ] [cm 2 ]	(B.5)
	Attachment coefficient β att , (labeled η in [5]):
			if else, N E	> 1.05 • 10 -15 ,	η 2 N η 2 N	= 8.889 • 10 -5 • = 6.089 • 10 -4 •	E N N E	+ 2.567 • 10 -19 [cm 2 ] -2.893 • 10 -19 [cm 2 ]	(B.6)
						η 3 N 2 = 4.7778 • 10 -59 •	E N	-1.2749	[cm 5 ]	(B.7)
	for	E N	> 2 • 10 -15 ,			v e = -7.4 • 10 21 •	E N	+ 7.1 • 10 6 •	E E	[cm• s -1 ]
	for 10 -16 <	E N	≤ 2 • 10 -15 ,	v e = -1.03 • 10 22 •	E N	+ 1.3 • 10 6 •	E E	[cm• s -1 ]
	for 2.6 • 10 -17 <	E N	≤ 10 -16 , v e = -7.2973 • 10 21 •	E N	+ 1.63 • 10 6 •	E E	[cm• s -1 ]
	for	E N	≤ 2.6 • 10 -17 ,		v e = -6.87 • 10 22 •	E N	+ 3.38 • 10 4 •	E E	[cm• s -1 ]
									(B.1)
	Positive ion drift velocity v p :
									v p = 2.34 • E •	p 0 p	[cm• s -1 ]	(B.2)
	Negative ion drift velocity v n :
						for for		E N E N	> 5 • 10 -16 , v n = -2.7 • E • < 5 • 10 -15 , v n = -1.86 • E • p 0 p p p 0	[cm• s -1 ] [cm• s -1 ]	(B.3)
	p 0 p	is a ratio of the atmospheric pressure to a gas pressure.
	Diffusion coefficient D e : D e =   0.3341 • 10 9 •	E N	0.54069	  •	v e E	[cm 2 • s -1 ]	(B.4)

Acknowledgments I would like to thank to the supervisors Professor Ing. Jaroslav Fořt, CSc. (Czech Technical University in Prague, Faculty of Mechanical Engineering) and Professor Fayssal Benkhaldoun (Université Paris 13, LAGA). I also would like to thank Professor Laurence Halpern for welcoming me at LAGA and Professor Khaled Hassouni (Université Paris 13, LSPM) and Doc. Ing. Jiří Fürst, PhD. (CTU in Prague, Faculty of Mechanical Engineering) for valuable advices and for participation in the jury. I also would like to thank Professor Philippe Helluy, Professor Richard Liska and Director of Research at CNRS Anne Bourdon for they acceptation to be reviewers of the work and also for their participation in the jury. And finaly I also would like to thank to Doc. Ing. Jan Halama, PhD for his participation in the jury. The work has been partly supported by the Marie Curie Actions of the European Commission in the frame of the DEASE project (MEST-CT-2005-021122) and project of the Grant of Czech Technical University in Prague No SGS10/244/OHK2/3T/12.

Grid Adaptation Test

The properties of the grid adaptation algorithm are the subject of the tests. It is important to have sufficiently high density grid in critical areas of the computational domain. These areas are the streamer (planar wave in 2D) path and the head of the discharge. The grid adaptation up to level 4 is therefore used. We need on the other hand to avoid too large extension of region with high level of adaptation because of almost explosive increase of number of cells.

Test Description

All computations are done on the same reference grid with 4700 cells (see Fig. 5.4). The first order scheme is used in the tests. As the initial condition we use the temporal result achieved on adapted grid at physical time t in = 3.14 • 10 -8 s (see Table 5.1) calculated from initial conditions (5.3).

We can see ( Fig. 5.4), that in this time the locally modified electric field is developed. It further initiated the growth and motion of electron avalanche. The computation starts with grid adaptation step and next the time evaluation is calculated. It means, that the new grid is created according to prescribed specific version of grid adaptation algorithm and initial solution is conservatively interpolated into finite volumes of new grid. After that we start to compute the time evaluation of solution. We use the same criterion for calculation of time step, the same CFL number, fixed number of computed time steps between grid adaptation. The computed cases differs only in grid adaptation algorithm or its parameters. We compare results achieved at the same physical time T = 5.25 • 10 -8 s. The number of triangles and nodes in Table 5.2 means the numbers on adapted grid at time t = T . The CPU time is the relative CPU time necessary for computation of solution in time interval t in , T . cells nodes reference mesh 4700 2424 mesh at t in = 3.14 • 10 -8 9154 4651 We can also observe, that only simultaneous use of both variables (case 3) in grid adaptation criterion produce acceptable results. The position of planar wave head and shape of source term as in case 1 is similar only in version V2.

Future Work

Despite the achieved results, there is still lots of work to do. It's necessary to improve the dynamical mesh adaptation algorithm in three dimension. Especially the mesh conformity step needs changes to produce better quality of a mesh. One of the possibilities is described e.g. in papers [START_REF] Bey | Tetrahedral Grid Refinement[END_REF], [START_REF] Liu | Quality Local Refinement of Tetrahedral Meshes Based on 8-Subtetrahedron Subdivision[END_REF].

It is also important to implement massive parallelization to the developed code because the 3D computation is very time consuming. The simple minimal model of streamer motion needs to be extended to more complex model (adding equation for other charged species, more realistic time depending boundary conditions, chemical reactions etc.).

Appendix A A.1 Relations Used in Reconstruction

Relations coming from the least square method which are used in 2D reconstruction

Relations coming for the least square method which are used in 3D reconstruction