## 1. Introduction

For the analysis of collisional plasmas one needs to analyze the different constituent particles which characterize these discharges, mainly the charged particles such as electrons, positive and negative ions as well as the neutral gas particles. Due to the presence of the charged particles, one needs to calculate also the electromagnetic fields within the plasma which are characterized by the Maxwell’s equations. The general Maxwell equations can be characterized by substituting a scalar and vector potential equations using Lorentz Gauge transformation to calculate the electric and magnetic fields. In the case of electrostatic fields, Maxwell equations reduce to Poisson equation that characterizes the electric field in the absence of a magnetic field, whereas in the magnetostatic case, a steady current exists invariant in time with the magnetic field related quantities considered constant. At high frequencies, the electromagnetic wavelength is small and if one is outside this single wavelength, the electric and magnetic field are directly coupled to each other such that if one of the two parameters is calculated, then the other is known. In the case of low frequencies, where the wavelength is high, if you are within the near field region, the electric and magnetic fields are completely independent, therefore the solution of both electric and magnetic field equations is necessary to calculate the electric and magnetic field distribution within the collisional plasma.

The charged particles behaviour is characterized using the continuity conservation equations of mass, momentum and energy for electrons, positive and negative ions including convective, diffusive and source term phenomena. Regarding the behaviour of the neutral gas particles within collisional plasmas, the Navier-Stokes equations, which are the conservation equations for mass, momentum and energy for the neutral gas particles need to be solved including convective, diffusive as well as source terms phenomena such as shear stresses, momentum transfer by elastic collisions, Lorentz forces, Joule heating and thermal conduction. Having established the necessary equations, the implementation of the solution procedure for solving Poisson, charged particle continuity and Navier-Stokes equations is presented. Thereafter, the Finite Element-Flux Corrected Transport method (FE-FCT) formulation follows in two-dimensional Cartesian, two-dimensional cylindrical axisymmetric and three-dimensional Cartesian coordinates, comprising of the predictor-corrector step based on the Taylor-Galerkin finite element method to calculate the high and low order schemes. The validation of the fluid flow equations using the FE-FCT is performed using the shock tube type problem, the shock wave incident on a wedge test case and the energy source term that result in sound and shock wave generation. Having validated the fluid model thoroughly, an adaptive mesh generation technique is discussed which reduces computational needs significantly and at the same time guaranteeing the stability and accuracy of the results. Finally, different collisional plasma configurations are analyzed including avalanche, primary and secondary streamer propagations and finally heating effects in constant voltage Dielectric Barrier Discharges and normal and abnormal glow discharges in ambient atmospheric air.

## 2. Model description

The complete plasma model in its multidimensional form consists of the Maxwell equations to account for the electromagnetic field, the continuity equations for charges to account for the charged particles (electrons, positive and negative ions) and the Navier-Stokes equations to account for the neutral gas charges.

### 2.1. Maxwell’s equations

The Maxwell equations consist of the following four differential equations in macroscopic form:

where **D** is the electric field strength, **B** is the magnetic flux density, **E** is the electric field strength, **H** is the magnetic field strength, ρ_{c} is the net charge density, **J** is the current density and t is the time.

Equation (1) is known as the Gauss’s law equation for electricity, equation (2) the Gauss’s law for magnetism, equation (3) as the Michael-Faraday’s or Faraday’s law of induction and finally equation (4) shows the Ampere’s circuital law with Maxwell’s correction.

where:

where ε_{c} and μ_{c} are the dielectric and magnetic permeability constants. One can define a vector potential **A** that **B** is the curl of:

which satisfies equation (2) above of Maxwell equations which basically states that no magnetic monopoles exist such that:

since the divergence of a curl is zero satisfying the magnetic monopole constraint.

Now by substituting equation (8) into Faraday’s equation (3), one gets:

By rewriting equation (10) above, one gets:

The term in the parenthesis of equation (11) above has no curl present implying that a potential V exists such that:

which gives the electric field to be:

Substituting equation (13) into Amperes law of electromagnetism equation (4) gives:

where :

with c being the speed of propagation within the medium. From vector calculus, one has that:

Substituting equation (16) above and equation (13) into equation (14) gives:

which gives:

Rearranging equation (18) above gives:

In order to make the choice of vector potential **A** not arbitrary, one imposes the Lorentz gauge which sets the parenthesis of equation (19) equal to zero:

with equation (20) above becoming:

Equation (21) above is Maxwell’s equation for the vector potential. On another note, by combining equations (1) and (13), one gets:

which gives:

Now using the Lorentz Gauge tranformation from equation (20) into equation (23) above gives:

Finally, equations (21) and (25) are respectively the vector and scalar potential equations of the Maxwell’s equations which when solved by using the equations (26) and (27):

and:one can calculate the electric and magnetic field of time-dependent Maxwell’s equations.

#### 2.1.1. Electrostatic case

In the electrostatic case, one assumes that ∂/∂t, **J**, **H** and **B** are all zero, and equation (25) becomes:

which is the Poisson equation. In the presence of no free charges, ρ = 0 and the above equation (29) becomes the Laplace equation:

#### 2.1.2. Magnetostatic case

In the magnetostatic case, there are steady currents in the system under consideration, which generate magnetic fields (ferromagnetic media are ignored), therefore ones sets ∂/∂t=0 and assumes that **J**, **H** and **B** are constant vectors. Substituting these relations into the scalar and vector potential equations (21) and (25) of Maxwell’s gives:

which is again the Poisson equation:

and:which is the magnetic equivalent of the electrostatic Poisson equation.

### 2.2. Continuity equations for charges

The continuity equations which comprise of the conservation of mass for electrons, positive and negative ions are described in equations (36) to (38) below:

where N_{e}, N_{p}, N_{n} are respectively the electron, positive and negative ion densities, **W _{e}**,

**W**and

_{p}**W**the corresponding velocity vectors and S

_{n}_{e}, S

_{p}, S

_{n}are the source terms for the electrons, positive and negative ions respectively, which are calculated according to equations (39) to (41) below:

where α the ionisation coefficient, η the attachment coefficient, β_{ep} the recombination coefficient between electrons and positive ions, β_{pn} the recombination coefficient between positive and negative ions and D_{e} the electron diffusion coefficient. The positive and negative ion diffusion coefficients are omitted since their effect is negligible compared to electron diffusion for the time scales considered here.

The conservation of momentum and energy density for electrons, positive and negative ions equations are described below in equations (42) and (43) and have the general form:

Subscripts s, p and n represent the electron, positive and negative ions, and the various terms are similar to the ones used for the neutral gas particle conservation equations that are explained below.

### 2.3. Navier-Stokes equations

#### 2.3.1. Navier-Stokes general form

The Navier-Stokes equations which account for the neutral gas particles are described by equations (44) to (46) below:

where ρ the neutral gas density, **v** the neutral gas velocity vector, S the neutral gas density source term, P the pressure, **τ** the shear stress tensor which comprises of the τ_{ii}_{} and τ_{jj} components, **MT _{s}** the momentum transfer of charged to the neutral particles due to elastic collisions, ε the neutral gas thermal energy density,

**Q**is the thermal conductivity term, f

_{ts}the percentage of energy density of charged particle with subscript s that is transferred as thermal energy to the neutral particles due to inelastic collisions and

**J**the current density,

**JxB**is the Lorentz force and

**(JxB).v**is the magnetic Lorentz force acting on the energy of the flow.

The neutral gas density source term is calculated by:

where m is the neutral gas particle mass which is constant.

The thermal conductivity, current density and pressure are calculated respectively by equations (48) to (50):

where k is the thermal conductivity, T is the neutral gas temperature, q_{s}, N_{s,}
**W _{s}** the charge, density per unit volume and velocity vector of charged particle s, N the number of neutral gas particles per unit volume and K the Boltzmann constant.

The shear stress tensors are calculated by:

where the viscosity coefficient, and subscript i and j represent the different degrees of freedom in the different space directions.

The momentum transfer of charged to the neutral particles due to elastic collisions **MT _{s}** is calculated by:

where μ_{s} is the mobility of charged particle s.

The numerical algorithm and its implementation are presented below in two-dimensional Cartesian, two-dimensional axisymmetric cylindrical and three-dimensional Cartesian coordinates in the following sections.

#### 2.3.2. Two-dimensional Cartesian coordinates of Navier-Stokes

The Navier-Stokes equations in two-dimensional Cartesian (x, y) coordinates are formulated as follows:

#### 2.3.3. Two-dimensional cylindrical axisymmetric coordinates of Navier-Stokes

The Navier-Stokes equations in two-dimensional cylindrical (r, z) axisymmetric coordinates are formulated as follows:

#### 2.3.4. Three-dimensional Cartesian coordinates of Navier-Stokes

The Navier-Stokes equations in three-dimensional Cartesian (x, y, z) coordinates are formulated as follows:

(66) |

The two-step Lax-Wendroff technique comprising of the predictor-corrector steps is used for time stepping. The charged particle continuity and Navier-Stokes are discretised using Taylor-Galerkin Finite Elements [1], whereas Poisson’s equation using Galerkin Finite Elements. Accuracy and efficiency for these long calculations are a crucial factor. The Flux Corrected Transport method (FCT) [2-4] ensures that accurate and efficient results are obtained free from inaccuracies building up and from non-physical oscillations.

(67) |

## 3. Implementation of the solution procedure

The solution procedure follows the pattern shown in Figure 1 for a single time step. It assumes that an initial distribution of electrons and positive ions at a given neutral gas temperature, velocity, density and applied voltage exists. The charged particle densities at time n (ρ^{n}) are fed to Poisson’s equation (PO) to obtain the electric field distribution at time n (E^{n}). The electric field and neutral gas parameters (N^{n}) are used to calculate the transport parameters (TP) at time n (TP^{n}). Then, the transport parameters are passed to the charge particle continuity equation (CON) to calculate the charge densities at time n+1/2 (ρ^{n+1/2}). Next, the charge densities (ρ^{n}), electric field (E^{n}) and transport properties (TP^{n}) at time n are passed to the Navier-Stokes solver (NS) to calculate the neutral gas properties at time n+1/2 (N^{n+1/2}). This completes the half time step solution. Subsequently, the charge particle densities at time n+1/2 are passed to the Poisson solver to calculate the electric field at time n+1/2 (E^{n+1/2}). The electric field and neutral gas parameters at time n+1/2 are used to calculate the transport properties at time n+1/2 (TP^{n+1/2}). Consequently the transport properties are fed to the charged particle continuity solver to calculate the charge densities at time n+1 (ρ^{n+1}). Finally, the Navier-Stokes solver uses the transport parameters, the electric field and the charge densities at time n+1/2 to calculate the neutral gas parameters at time n+1 (N^{n+1}) and this process is repeated continuously to proceed forward in time.

In fluid analysis, the collision coefficients and drift velocities of the charged particles can be calculated either using experimental or numerical techniques. Experimentally they are calculated as a function of the ratio of the electric field to the number of neutral gas particles per unit volume (E/N), depending also on the density and pressure of the neutral gas that has a direct effect on the frequency of collisions between the charged and neutral gas particles.

Numerically, commercial software packages calculate the transport properties by solving the electron distribution Boltzmann equation by utilizing the collisional cross sectional data and such software examples are ELENDIF [5], BOLSIG [6], and BOLSIG+ [7] which are used to calculate the above transport parameters [8].

The charged particle continuity equations are coupled to Poisson’s equation via the electric field strength and to the Navier-Stokes equations via the number of neutral gas particles per unit volume and their pressure. Poisson’s equation is coupled to the continuity equations via the net charge density. As far as the Navier-Stokes equations are concerned, they are coupled to Poisson’s equation via the electric field strength and to the charged particle continuity equations via production and loss processes, and via momentum and kinetic energy exchange between charged and neutral particles.

## 4. Finite element formulation of continuity and Navier-Stokes fluid equations

### 4.1. Two-dimensional Cartesian coordinates

The fluid transport equations take the general form in the two-dimensional Cartesian coordinates as:

where

The Taylor Galerkin scheme is used to develop the FE-FCT scheme, with the time discretization preceding the space discretization. The time stepping is performed in two steps, the advective predictor and the corrector step.

#### 4.1.1. Advective predictor step

The advective predictor step is formulated using the Taylor series expansion for the unknown variable

where A_{e} is the area of a triangular element, Δt is the time step, b_{i} and c_{i} are the linearly interpolated piecewise shape functions. The velocities at time (n+1/2) are assumed to be the average value of the nodal velocities at time n.

#### 4.1.2. Corrector step

The corrector step utilizes the Taylor series expansion for the full time step to the second order approximation to get:

(71) |

where D_{c} is the consistent mass matrix, _{e} is the boundary surrounding element e and

where _{l}) is used for computational purposes, resulting in:

where α is an integer representing the iteration number.

To calculate the low order scheme, mass diffusion is added of the form:

to get the low order scheme:

where c_{d} is the variable diffusion coefficient and is dependent on the mesh size, time step and speed of an element. A detailed analysis of the above formulation procedure is found in [9].

### 4.2. Two-dimensional cylindrical axisymmetric coordinates

The fluid transport equations take the general form in the two-dimensional cylindrical axisymmetric coordinates as:

where

The Taylor Galerkin scheme is used to develop the FE-FCT scheme, with the time discretization preceding the space discretization. The time stepping is performed in two steps, the advective predictor and the corrector step.

### 4.2.1 Advective predictor step

The advective predictor step is formulated using the Taylor series expansion for the unknown variable

where A_{e} is the area of a triangular element, Δt is the time step, b_{i} and c_{i} are the linearly interpolated piecewise shape functions.

where

where

The velocities at time (n+1/2) are assumed to be the average value of the nodal velocities at time n.

#### 4.2.2. Corrector step

The corrector step utilizes the Taylor series expansion for the full time step to the second order approximation to get:

(80) |

The above discretization results in the following equations:

where:

andwhere _{c} and lumped matrix D_{l} need to be multiplied by

### 4.3. Three-dimensional Cartesian coordinates

The fluid transport equations take the general form in the three-dimensional Cartesian coordinates as:

The Taylor Galerkin scheme is used to develop the FE-FCT scheme, with the time discretization preceding the space discretization. The time stepping is performed in two steps, the advective predictor and the corrector step.

#### 4.3.1. Advective predictor step

The advective predictor step is formulated using the Taylor series expansion for the unknown variable

where V_{e} is the volume of a tetrahedral element, b_{i}, c_{i} and d_{i} are the linearly interpolated piecewise shape functions. The velocities of the charged species at time (n+1/2) are assumed to be the average value of the nodal velocities at time n.

#### 4.3.2. Corrector step

The corrector step utilizes the Taylor series expansion for the full time step to the second order approximation to get:

(86) |

where A_{e} is the boundary surrounding element e. Writing the above equation in the general form results in:

Instead of using a consistent mass matrix, a lumped mass matrix is used for computational purposes, resulting in:

To calculate the low order scheme, mass diffusion is added of the form:

to get the low order scheme:

## 5. Fluid flow validation of the above FE-FCT formulation

In order to validate the fluid flow equations and the methodology developed, the Euler equations which are the Navier-Stokes equations without diffusion and source terms were tested using the shock tube type problem, where different parameters, such as density, diffusion coefficient and mesh size were considered. Furthermore, a shock wave incident on a wedge was also tested in the Cartesian case. Finally, energy source terms are introduced into the two-dimensional Cartesian and axisymmetric cylindrical geometries, which result in sound and shock wave generation, resembling the heating expected during the development of a spark and an arc plasma gas discharge.

### 5.1. Riemann test case

The Euler equations are validated in one-dimensional Cartesian coordinates using the shock tube or Riemann test case in air with the ratio of the specific heats, γ, equal to 5/3, where analytical results are available. Figure 2 shows a one-dimensional comparison between the analytical and numerical solution of density for three different meshes at times t = 0.1 and 0.5 s. It is shown that in the case of the finer mesh, the overshoot observed is reduced. Figure 3 shows the analytical and numerical solution for density for the finest mesh (Mesh 3) at three different instants in time of t = 0.1, 0.3 and 0.5. A direct comparison shows clearly that the Euler solver developed is capable of simulating the propagation of shock waves and that the results are of adequate accuracy. This is also verified by calculating the maximum percentage error for the three different meshes as shown in Table 1 below.

### 5.2. Heating source term test case

The two-dimensional Cartesian Euler equations in air with γ equal to 1.4 are tested by introducing a heating source term of magnitude 1x10^{12} J/m^{2} for a duration of 3x10^{-8} s within a square block of size 0.002 m x 0.002 m, with the simulation being run for 5000 steps at a constant time step increments of 3x10^{-8} s. The contour plot of the momentum in the x direction is shown in Figure 4, where at the exterior of the wave front, the wave is moving outwards in the x direction with maximum momentum. At the interior of the wave front, there is a sudden decrease and a change in the direction of the momentum. This is due to the fact that the inertial effects of the sudden explosion cause an overexpansion of the wave, therefore a rarefraction wave is created moving inwardly. Figure 5 shows a contour plot of the energy at a time of 9x10^{-5} s. High energies occur at the shock front of the wave with low energies forming behind the shock front, again due to the overexpansion of the wave, which are in agreement with blast wave theory Baker (1973) [10].

### 5.3. Shock wave incident on a wedge test case

A shock wave of Mach number 2 in air with γ equal to 1.4 is incident on a wedge angle of 46^{o} degrees. The neutral gas is at a temperature of 300 ^{o}K and a pressure of 30 kPa. The simulation is run for 13000 time steps of 1x10^{-7} s duration. The simulation is terminated 0.1 m before the shock wave hits the boundary on the right. Figure 6 shows the Cartesian contour plots of the density at a time of 1.3x10^{-3} s with the results being compared with the benchmark test case of a shock wave impinging on a wedge Takayama (1997) [11]. The results are found to be in good agreement capturing the shock wave created at (0.9 m, 0.7 m) coordinates.

### 5.4. Micro-blast waves test case

The two-dimensional Euler solver is validated in the cylindrical axisymmetric case using numerical blast wave tests that involve the release of sudden burst of energy by focusing a laser beam within small volumes in air with γ equal to 1.4, leading to the development of shock waves. Jiang et al (1998) [12] has analyzed both experimentally and numerically the propagation of micro-blast waves in ambient air. The duration of the pulse was of the order of 18 ns and the total amount of energy released was measured to be 1.38 J. This sudden release of energy raises instantaneously the pressure, density and temperature to values many times higher than those of the ambient conditions, thereby creating highly compressive acoustic waves that will in turn generate strong shock waves.

Figures 7 and 8 compare respectively one-dimensional plots of the density against time along the axis of symmetry obtained by Jiang and the authors for the same mesh and different diffusion coefficients and for different mesh and the same diffusion coefficient. Given the diversities of the two numerical schemes, the results show to be in good agreement.

## 6. Adaptive mesh algorithm

In order to reduce computational needs considerably, an adaptive mesh generator has been developed, rendering possible the analysis of heating effects both in short and long gaps. For an adaptive mesh generator, an error indicator is necessary to decide on the amount of refinement. The error indicator used by the author is the one used by Lohner [13], which in multidimensional form is calculated as 999999follows:

where N^{i}_{k} is the shape function of node i in element k, N^{j}_{l} is the shape function of node j in element l, U_{j} is the value of the variable chosen to be used as error indicator at node j, E^{i} is the error indicator value at node i, and is a factor varying from 0 to 1. The value of is added as noise filter, so that any loss of monotonicity such as wiggles or ripples are not refined. The above error indicator is dimensionless, fast to calculate, varies from 0 to 1 such that prefixed tolerances and many variables as error indicators can be used at the same time.

In this case, an adaptive mesh algorithm in two-dimensional Cartesian, and two-dimensional cylindrical axisymmetric coordinates has been developed by the author in order to analyze plasmas much faster and more accurately. This is achieved by constructing a computer resource efficient algorithm, which is automated in providing the necessary results. Generally, the numerical solution of time dependent differential equations is classified into the categories of static and dynamic. For static methods, any addition of nodes, and edge swap and topological movement of nodes is performed at a fixed time. In the dynamic case, or moving mesh methods, a mesh equation is introduced that involves node velocities such that a fixed number of nodes are moved in such a way that the nodes are always concentrated near regions of rapid variation of the solution. Thereby the simultaneous solution of the differential and mesh equations is necessary, having the advantage that no interpolation between existing and future meshes is necessary [14]. The author of this paper has adopted the static method, which is more widely used and tested in the literature.

The algorithm incorporates an innovative element quality improvement procedure that has the ability to guarantee the generation of new meshes, as well as the treatment of existing bad element quality meshes, to nearly ideal standards, guaranteeing a minimum element quality of 0.85, and above 0.90 of average element quality in uniform and non-uniform geometric domains. Furthermore, another novelty of the paper is also the utilisation of an interpolation method between meshes which is very fast, making the adaptive meshing more attractive. The re-meshing times are decided optimally according to the maximum speeds of the fastest particle within the simulation to guarantee that the desired mesh resolution is in place at all times, whereas the numerical diffusion error due to the interpolation between meshes is minimized by interpolating between meshes of similar size.

Figure 9 shows a flow chart of the implementation procedure of the adaptive mesh developed by the author. The first step is to create a reference coarse mesh that will be used as a reference mesh for future refinement and coarsening of the meshes. The second step is to apply the initial conditions on the reference coarse mesh, and calculate the amount of refinement by multiplying the error by a constant factor in each element. The refinement factor is decided on a basis of experience, and trial and error. Then a newly adapted initial mesh is developed using a freely available two-dimensional software package which uses the divide and conquer method to refine elements [15]. Since the elements that are created from this package are of bad element equality, they are first treated using adaptive mesh techniques. Specifically, the h-refinement/coarsening technique is firstly used, which includes the edge swap and node additions/removals to improve the interconnectivity between adjacent nodes. Secondly, the r-refinement technique, which involves jiggling of the mesh, follows, i.e. the movement of nodes around the geometry in a controlled way, such that the overall element quality is improved.

Once the adapted initial mesh is created that can be used for the solution of the differential equations, an interpolation of the results from the reference coarse mesh to the adapted initial mesh is performed. Having defined the mesh to be used and calculated the corresponding values that all the variables have at the nodes of this mesh, the simulation proceeds forward in time, until an indicator signifies that the solution has reached the outer boundary of the refined region, and a re-meshing operation needs to be performed. The indicator calculates the maximum distance that the fastest particle travels, and ensures that it does not exceed the geometric tolerance of the initial coarse mesh, thereby it is guaranteed that the correct resolution is provided at all times. Then in order to decide on the appropriate mesh to run the simulations, the results are interpolated from the adapted initial mesh back to the reference coarse mesh, and then the error is calculated to create a final mesh. Then this mesh is treated as above, using edge swap and node additions/removals operations to improve the interconnectivity between adjacent nodes, and finally mesh jiggling operations to improve the overall element quality, creating the adapted final mesh. Then one needs not to interpolate from the reference coarse mesh to the created adapted final mesh, but instead from the created adapted initial mesh to the created adapted final mesh. This operation is performed due to the fact that interpolation is a source of numerical diffusion on the results, and by interpolating just once, instead of twice, makes the results during the interpolation processes more accurate. Then having decided on the adapted final mesh and having interpolated the results, then one runs the simulation forward in time and the procedure is repeated many times. This completes the implementation of the adaptive mesh algorithm. For the implementation of the above adaptive mesh algorithm, three tools are necessary, which are (a) the error calculation, (b) the element quality improvement algorithm, and (c) the algorithm for interpolation between meshes tools, and are thoroughly explained below.

## 7. Application of electromagnetic and fluid code on microplasmas

### 7.1. RF applicator at 40 MHz

In this section, results are presented showing the avalanche and streamer propagation in a gap of 1 cm in an RF applicator at 40 MHz. Figures 10 and 11 refer to times t_{1} = 50, t_{2} = 55, t_{3} = 62.5, t_{4} = 63.5 ns. Figure 10 shows the radial field along the symmetry axis from times t_{1} to t_{4.} It is shown that during the time interval t_{1} to t_{2}, its magnitude increases to a value of approximately 0.5 x 10^{5} V/m due to the net charge that exists and that it stays fairly constant at the bottom electrode end, but gradually increases at the upper electrode end, with both radial fields always extending closer towards the electrodes as time progresses. It is interesting to note that the radial field is positive during this time at the upper end and negative towards the bottom end. This initiates the propagation of streamers towards the electrodes. At time t_{3}, the radial field reverses at both ends, when compared with the time interval of t_{1} to t_{2}, with the radial field magnitude being larger than before and this is due to the streamer impinging on the two electrodes. At time t_{4,} the radial field increases even further at the bottom electrode side, but changes sign at the upper electrode side. This is because the streamer that hits the upper electrode forms earlier than the one that hits the lower electrode, as the discharge is overall exposed to a larger negative voltage cycle. The change of sign of the radial field at the upper electrode is due to the absorption of the electrons into the electrode and the start of the formation of a cathode fall region of net positive charge near it. Figure 11 shows the axial electric field along the symmetry axis from times t_{1} to t_{4.} At time t_{1,} at the middle of the gap, the axial electric field gets distorted from the Laplacian field which initiates the streamer propagation. During the time interval t_{1} to t_{2}, the distortion increases in magnitude and spreads out further towards the two electrodes. At time t_{3}, there is a sudden decrease of the axial electric field at the two electrodes, leaving a nearly constant plateau in the middle of the gap due to the streamer impinging on the two electrodes. At a later stage of t_{4},_{,} the ripple profile reverses shape and starts to increase again at an ever higher rate, which corresponds to the times when the streamers charge is accumulated on the two electrodes.

### 7.2. Heating effects of normal and abnormal glow discharge

The development of a glow discharge in atmospheric pressure air and its associated Joule heating are analyzed by solving the Poisson, charged particle continuity and Navier-Stokes equations using the FE-FCT method in two-dimensional cylindrical axisymmetric coordinates. An applied direct current voltage of 20% above the breakdown voltage is applied at the anode in a 1 mm air gap between two parallel plate electrodes and a normal glow discharge is shown to consist of the cathode fall, negative glow, positive column and anode regions. Neutral gas heating occurs with the initiation of the glow discharge, with the temperature at the anode shown to increase by only a few degrees ^{o}K as shown in Figure 12 due to the low electron densities and axial fields at the adjacency of the anode, whereas at the proximity of the cathode the temperature increases by approximately 180 ^{o}K as shown in Figure 13, with the temperature maximum at the cathode fall region, reducing abruptly at almost ambient temperature just outside it, where it remains almost constant.

The numerical results for the development of a glow discharge and its associated Joule heating by modeling its transition from a single electron, in a uniform applied electric field in atmospheric air, to a fully fledged state which consists of the cathode fall, negative glow, positive column and anode regions are presented. It is shown that the positive column travels towards the anode in the form of a return wave and maximum heating by approximately 180 ^{o}K occurs at the cathode fall region during the glow discharge development.

### 7.3. Secondary streamers

A constant voltage of 130 kV is applied on a dielectric barrier discharge configuration with the two metallic electrodes placed 4 mm apart, and with each one covered by alumina dielectric of 1 mm thickness, leaving an air gap of 2 mm for the discharge to develop. The two-dimensional neutral gas temperature distribution at time t = 0.3263 μs is shown in Figure 14, with most of the heating within the microplasma occurring along the symmetry axis where most of the activity takes place [8]. The column of heated air extends to a radial distance of 0.2 mm, and the temperature is raised approximately 120 ^{o}K. It is also observed that striations occur between the anode and the cathode. Along the dielectric electrodes, the temperature is larger than in the inter-electrode gap, but smaller than that at the symmetry axis. On average, higher neutral gas temperatures are observed at the cathode due to the arrival of the primary and secondary streamer charges on the dielectric cathode [16]. The two-dimensional axial electric field distribution at time t = 11.8 ns is shown below in Figure 15. Along the symmetry axis, there exists a cathode fall, negative glow, positive column and anode regions, similar to those of a normal glow discharge. A secondary streamer is shown to form at a radial distance of 1.5 mm that travels from the anode towards the cathode, extending radially up to the outer boundaries of the discharge and at later stages, when it hits the cathode, the sheath region spreads throughout the cathode dielectric surface.

### 7.4. Adaptive mesh results on DC avalanche and streamers in ambient air

The adaptive mesh generator has been tested in a geometric configuration that comprises of two metallic parallel-plates that are placed 1 mm apart in ambient atmospheric air. With a single electron released at the cathode as initial condition, the development of the avalanche and the primary streamer are analyzed as shown in Figures 16a-f, which show the mesh plots and the corresponding electron densities along the symmetry axis at three different instances in time. Figures 16a and 16b show respectively the initial mesh refinement and electron density distribution along the symmetry axis, whereas Figures 16c and 16d show the mesh and electron density during cathode directed primary streamer propagation, and Figures 16e and 16f show the mesh and electron density distributions just before the streamer hits the cathode. In the simulations, the photoemission effect is included in the calculation using the model developed by the authors [17], whereas photoionization and impact ionization phenomena are excluded. It has been shown that the avalanche and streamer propagation are captured by the adaptive mesh generator in an optimum way.

## 8. Conclusions

The electromagnetic and fluid analysis of collisional plasmas has been thoroughly discussed including the necessary differential conservation equations to characterize such plasmas. The finite-element formulations were presented for the solution of these equations and the implementation procedure to couple the above set of equations was discussed. Thereafter, the FE-FCT algorithm was validated against theoretical and experimental fluid flow results and it was then used in a variety of collisional plasma configurations to study different plasma phenomena.