Open access peer-reviewed chapter

Electromagnetic and Fluid Analysis of Collisional Plasmas

By Antonis P. Papadakis

Submitted: November 22nd 2011Reviewed: May 3rd 2012Published: October 10th 2012

DOI: 10.5772/48328

Downloaded: 1176

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:

D=ρcE1
B=0E2
×E=BtE3
×H=J+DtE4

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:

J=σ(E+uxB)E5
D=εcEE6
B=μcHE7

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

B=xAE8

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

.B=.(xA)=0E9

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:

(xA)t=xEE10

By rewriting equation (10) above, one gets:

x(E+At)=0E11

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

E+At=VE12

which gives the electric field to be:

E=AtVE13

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

Etc2(xB)Jεc=0E14

where :

c2=1μcεcE15

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

×B=××A=(A)2AE16

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

(AtV)tc2((A)2A)Jεc=0E17

which gives:

2At2(Vt)+c22Ac2((A)Jεc=0E18

Rearranging equation (18) above gives:

2At2c22A+c2(A+1c2Vt)Jεc=0E19

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:

(A+1c2Vt)=0E20

with equation (20) above becoming:

2At2c22AJεc=0E21

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

(VAt)=ρcεcE22

which gives:

2V(.A)t=ρcεcE23

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

2V+(1c2Vt)t=ρcεcE24
2Vt2c22Vc2ερc=0E25

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):

B=xAE26
and:
E=VAtE27

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:

c22Vc2ερc=0E28
2V=ρcεcE29

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

2V=0E30

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:

c22Vc2ερc=0E31

which is again the Poisson equation:

2V=ρcεcE32
and:
c22AJεc=0E33
2A=Jc2εcE34
2A=μJE35

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:

Net+.(NeWe)=SeE36
Npt+.(NpWp)=SpE37
Nnt+.(NnWn)=SnE38

where Ne, Np, Nn are respectively the electron, positive and negative ion densities, We, Wp and Wn the corresponding velocity vectors and Se, Sp, Sn are the source terms for the electrons, positive and negative ions respectively, which are calculated according to equations (39) to (41) below:

Se=αNe|We|ηNe|We|βepNeNp+.(DeNe)E39
Sp=αNe|We|βepNeNpβpnNpNnE40
Sn=ηNe|We|βpnNpNnE41

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 De 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:

(ρsWs)t+.(ρsWsWs)=Ps.τssMTsE42
εst+.(εsWs)=.(Qs).(Ws(PsI+τs))sMTs.vsftsJs.EE43

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:

ρt+.(ρv)=SE44
(ρv)t+.(ρvv)=P.τ+sMTs+JxBE45
εt+.(εv)=.(Q).(v(PI+τ))+sMTs.v+sftsJ.E+(JxB).vE46

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, MTs 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, fts 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:

S=m(αNe|We|ηNe|We|+2βpnNpNn+βepNeNp)E47

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):

Q=kTE48
J=qsNsWsE49
P=NKTE50

where k is the thermal conductivity, T is the neutral gas temperature, qs, Ns, Ws 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:

τii=μ(2vixi23.v)E51
τij=μ(vixj+vjxi)E52

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 MTs is calculated by:

MTs=Nsmm+msqsμs(Wsv)E53

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:

Qt+Kxx+Kyy+Lxx+Lyy=ME54
Q=[ρρvxρvyε]E55
Kx=[ρvxρvxvx+Pρvyvx(ε+P)vx]Ky=[ρvyρvxvyρvyvy+P(ε+P)vy]E56
Lx=[0τxxτyxKTx+τxxvx+τxyvy]Ly=[0τxyτyyKTy+τyyvy+τyxvx]E57
M=[m(αNe|We|ηNe|We|+βepNeNp+2βpnNpNn)MTs{x}+(JxB){x}MTs{y}+(JxB){y}sftsJ.E+sMTs{x}vx+sMTs{y}vy+(JxB){x}vx+(JxB){y.}vy]E58

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:

Qt+1r(Krr)r+Kzz+1r(Lrr)r+Lzz=ME59
Q=[ρρvrρvzε]E60
Kr=[ρvrρvrvr+Pρvzvr(ε+P)vr]Kz=[ρvzρvrvzρvzvz+P(ε+P)vz]E61
Lr=[0τrrτzrKTr+τrrvr+τrzvz]Lz=[0τrzτzzKTz+τzzvz+τzrvr]E62
M=[m(αNe|We|ηNe|We|+βepNeNp+2βpnNpNn)MTs{r}+Pr+(JxB){r}vrMTs{z}+(JxB){z.}vzsftsJ.E+sMTs{r}vr+sMTs{z}vz+(JxB){r}vr+(JxB){z.}vz]E63

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:

Qt+Kxx+Kyy+Kzz+Lxx+Lyy+Lzz=ME64
Q=[ρρvxρvyρvzε]E65
Kx=[ρvxρvxvx+Pρvxvyρvxvz(ε+P)vx]Ky=[ρvyρvyvxρvyvy+Pρvyvz(ε+P)vy]Kz=[ρvzρvzvxρvzvyρvzvz+P(ε+P)vz]E66

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.

Lx=[0τxxτyxτzxKTx+τxxvx+τxyvy+τxzvz]Ly=[0τxyτyyτzyKTy+τyxvx+τyyvy+τyzvz]Lz=[0τxzτyzτzzKTz+τzxvx+τzyvy+τzzvz]E67
M=[m(αNe|We|ηNe|We|+βepNeNp+2βpnNpNn)MTs{x}+(JxB){x}vxMTs{y}+(JxB){y.}vyMTs{z}+(JxB){z}vzsftsJ.E+sMTs{x}vx+sMTs{y}vy+sMTs{z}vz+(JxB){x}vx+(JxB){y.}vy+(JxB){z}vz]E68

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 (En). The electric field and neutral gas parameters (Nn) are used to calculate the transport parameters (TP) at time n (TPn). 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 (En) and transport properties (TPn) at time n are passed to the Navier-Stokes solver (NS) to calculate the neutral gas properties at time n+1/2 (Nn+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 (En+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 (TPn+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 (Nn+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].

Figure 1.

Solution procedure of the complete model

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:

Qt+j=12Kjxj+j=12Ljxj=ME69

where Qis the unknown or independent variable, Kand Lare respectively the convective and diffusive fluxes, Mis the source term.

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 Qat a first order approximation, thereby resulting in the half time values of the independent variableQn+1/2:

Qn+1/2=Ae3i=13Qin+Δt4i=13Kinbi+Δt4i=13KinciΔt6Aei=13MinE70

where Ae is the area of a triangular element, Δt is the time step, bi and ci 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:

DcΔQn=Δtbie2Kxn+1/2Δtcie2Kyn+1/2Δtbie2Lxn+1/2Δtcie2Lyn+1/2+ΔtMn+1/2Ae3+ΔtΓe(Knn+1/2+Lnn)NiedΓE71

where Dc is the consistent mass matrix, Knn+1/2and Lnnare the convective term at time n+1/2 and diffusive term at time n fluxes outward in the normal direction of the boundary, Γe is the boundary surrounding element e and Nieis the interpolation function. Writing the above equation in the general form results in:

DcΔQn=BnE72

where ΔQnis the vector of nodal increments, Bnis the vector of added element contributions to the nodes. Instead of using a consistent mass matrix, a lumped mass matrix (Dl) is used for computational purposes, resulting in:

DlΔQa,n(DcDl)ΔQa1,n=BnE73

where α is an integer representing the iteration number.

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

cd(DcDl)QnE74

to get the low order scheme:

DlΔQa,n(DcDl)ΔQa1,ncd(DcDl)Qn=BnE75

where cd 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:

Qt+1r(Krr)r+(Kz)z+1r(Lrr)r+(Lz)z=ME76

where KrKzLrLrare respectively the radial and axial convective and diffusive fluxes, and r is the radial coordinate in cylindrical axisymmetric coordinates.

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 Qat a first order approximation, thereby resulting in the half time values of the independent variableQn+1/2:

Qn+1/2=112i=13(3+ri')(Qin+Δt2Mn)Δt4Aei=13((Kzn)ici+(Krn)ibi)Δt6rei=13(Krn)iE77

where Ae is the area of a triangular element, Δt is the time step, bi and ci are the linearly interpolated piecewise shape functions.

where ri'is given by:

ri'=rireE78

where reis the r-coordinate of the centroid of the element e given by:

re=i=13ri3E79

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:

reDcΔQn=Δtrebie2Krn+1/2Δtrecie2Kzn+1/2Δtrebie2LrnΔtrecie2Lzn+ΔtreAe3Mn+1/2+ΔtreAe(Knn+1/2+Lnn)NiedAE80

The above discretization results in the following equations:

D'cΔQn=B'nE81

where:

D'c=rDcE82
and
B'n=rBnE83

where rdenoting a matrix of reentries everywhere, Dc'is the consistent mass matrix and B'nis the vector of added element contributions to the nodes in the cylindrical axisymmetric case. In order to implement the Cartesian and axisymmetric cases together, the half time step values Qn+1/2have to be calculated slightly differently, and the consistent mass matrix Dc and lumped matrix Dl need to be multiplied by rein the cylindrical axisymmetric case, Finally, during the limiting procedure, the antidiffusive element contribution, the consistent and lumped matrices need to be multiplied by refor each element.

4.3. Three-dimensional Cartesian coordinates

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

Qt+j=13Kjxj+j=13Ljxj=ME84

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 Qat a first order approximation, thereby resulting in the half time values of the independent variableQn+1/2:

Qn+1/2=Ve4i=14Qin+Δt12i=14Kinbi+Δt12i=14Kinci+Δt12i=14KindiΔt8Vei=14MinE85

where Ve is the volume of a tetrahedral element, bi, ci and di 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:

DcΔQn=Δtbie3Kxn+1/2Δtcie3Kyn+1/2Δtdie3Kzn+1/2Δtbie3Lxn+1/2Δtcie3Lyn+1/2Δtdie3Lzn+1/2+ΔtMn+1/2Ve4+ΔtAe(Knn+1/2+Lnn)NiedAE86

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

DcΔQn=BnE87

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

DlΔQa,n(DcDl)ΔQa1,n=BnE88

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

cd(DcDl)QnE89

to get the low order scheme:

DlΔQa,n(DcDl)ΔQa1,ncd(DcDl)Qn=BnE90

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.

Mesh numberMaximum density error
11.6177 %
21.0895 %
30.2345 %

Table 1.

Maximum density error for three different meshes

Figure 2.

Comparison of analytical and numerical results using three different meshes at two different instants in time t = 0.1 and 0.5 s. □ line: Analytical solution, ◊ line: Mesh 1, 500 nodes, x line: Mesh 2, 1000 nodes and + line: Mesh 3, 5000 nodes

Figure 3.

Comparison of analytical and experimental results using the finer mesh at three different instants in time t = 0.1, 0.3 and 0.5 s. □ line: Analytical solution, + line: Mesh 3, 5000 nodes

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 1x1012 J/m2 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 46o degrees. The neutral gas is at a temperature of 300 oK 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.

Figure 4.

Contour plot of the x-momentum at time t = 9x10-5 s

Figure 5.

Contour plot of the energy at time t = 9x10-5s

Figure 6.

Contour plot of the density at a time of 1.3x10-3 s

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.

Figure 7.

Non-dimensionalized density (ρ*) in radial direction at t1 = 5.36 μs, t2 = 6.65 μs, t3 = 10.05 μs, t4 = 12.58 μs. - line: Jiang, ◊ line : Mesh 2 : 27570 nodes, Diff. Coeff., 0.0012, x line: Mesh 2:, 27570 nodes, Diff. Coeff., 0.002

Figure 8.

Non-dimensionalized density (ρ*) in the radial direction at t1 = 5.36 μs, t2 = 6.65 μs, t3 = 10.05 μs, t4 = 12.58 μs. - line: Jiang, ◊ line: Mesh 2, Diff. Coeff., 0.002, x line: Mesh 1, Diff. Coeff., 0.002

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:

Ei=k,l(ΩNkiNljdΩ.Uj)2k,l(Ω|Nki|[|NljUj|+ε(|Nlj||Uj|)]dΩ)2E91

where Nik is the shape function of node i in element k, Njl is the shape function of node j in element l, Uj is the value of the variable chosen to be used as error indicator at node j, Ei 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.

Figure 9.

Flow chart of the implementation procedure of the adaptive mesh generator

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 t1 = 50, t2 = 55, t3 = 62.5, t4 = 63.5 ns. Figure 10 shows the radial field along the symmetry axis from times t1 to t4. It is shown that during the time interval t1 to t2, its magnitude increases to a value of approximately 0.5 x 105 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 t3, the radial field reverses at both ends, when compared with the time interval of t1 to t2, with the radial field magnitude being larger than before and this is due to the streamer impinging on the two electrodes. At time t4, 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 t1 to t4. At time t1, 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 t1 to t2, the distortion increases in magnitude and spreads out further towards the two electrodes. At time t3, 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 t4,, 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.

Figure 10.

One-dimensional plots of the radial field at times t1 to t4 along the symmetry axis

Figure 11.

One-dimensional plots of the axial field at times t1 to t4 along the symmetry axis

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 oK 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 oK 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 oK occurs at the cathode fall region during the glow discharge development.

Figure 12.

One-dimensional plots of the temperature along the symmetry axis close to the cathode at t = 6.75 ns at V= 5600 V

Figure 13.

One-dimensional plots of the temperature along the anode at four instances in time at V = 5600 V

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

Figure 14.

Plot of the two-dimensional cylindrical axisymmetric neutral gas temperature at time t = 0.3263 μs

Figure 15.

Plot of the two-dimensional axial electric field distribution at a time t = 11.8 ns

Figure 16.

a) Two-dimensional cylindrical mesh at time t = 0 s at a voltage of 5600 V, (b) One-dimensional electron distribution along the symmetry axis at time t = 0 s at a voltage of 5600 V, (c) Two-dimensional cylindrical mesh at time t = 3.59 ns at a voltage of 5600 V, (d) One-dimensional electron distribution along the symmetry axis at time t = 3.59 ns at a voltage of 5600 V, (e) Two-dimensional cylindrical mesh at time t = 4.54 ns at a voltage of 5600 V, (f) One-dimensional electron distribution along the symmetry axis at time t = 4.54 ns at a voltage of 5600 V

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.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Antonis P. Papadakis (October 10th 2012). Electromagnetic and Fluid Analysis of Collisional Plasmas, Finite Element Analysis - Applications in Mechanical Engineering, Farzad Ebrahimi, IntechOpen, DOI: 10.5772/48328. Available from:

chapter statistics

1176total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Finite Element Analysis of Functionally Graded Piezoelectric Spheres

By A. Ghorbanpour Arani, R. Kolahchi, A. A. Mosalaei Barzoki, A. Loghman and F. Ebrahimi

Related Book

First chapter

Finite Element Analysis in Dental Medicine

By Josipa Borcic and Alen Braut

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us