User:Egm6322.s09.Three.nav/Project

 See my comments below. Egm6322.s09 14:29, 1 May 2009 (UTC)

About me
I am a doctoral candidate enrolled in the Mechanical & Aerospace Engineering Department at the University of Florida, Gainesville. I work in the Computational Plasma Dynamics Laboratory & Testing Facility under the supervision of Dr.Subrata Roy.

Plasma Physics


On earth, we know and have learnt to live with only three states of matter. However Sir William Crookes identified Plasma, the fourth state in 1879. It is by far the most common form of matter. It is present in the stars and in the tenuous space between them. It is present, though invisible, in almost 99.9% of the visible universe. Ideally plasma consists of neutrals, electrons and ions. However commercial plasma is ionized gas. In order to separate electrons and ions from atoms, energy is needed and this energy can be in any form, thermal, electric etc. Most commonly, a gas is passed through a pair of electrodes and an electric field. This electric field provides the energy needed to ionize the gas into atoms and ions. Plasma has helped create a better understanding of the world we live in. It provides us the prospect of abundant energy. It provides alternatives to surface cleaning, waste removal, in the medical industry and many more.

My research focuses on Dielectric Barrier Discharges(DBDs) & their applications. Commerically produced plasma is classified into two types- Volume Plasma (wherein the plasma is produced in the discharge space between two electrodes) & Surface Plasma (wherein both electrodes are embedded on the same surface). DBDs are a kind of surface plasma. They are generated between two electrodes with a dielectric barrier in between them. The electrodes are separated by a gap of a few millimetres. Typically DBDs are operated at voltages of 1-100 kV and frequencies of 50Hz-1MHz. When electric field is sufficiently high to cause breakdown of the discharge gas, a large number of micro-discharges can be observed emanating from the electrodes.

Current Research Focus- Plasma Sterilization Plasma Sterilization is a process in which low-temperature plasma is utilized to kill micro-organisms. Our research focus is on using Dielectric Barrier Discharges for this purpose. Plasma is produced usually by passing a gas through an electric field and ionizing it. Dielectric Barrier Discharges are a special case, in which two electrodes are embedded into the top and bottom of a dielectric surface (an insulator). An electric field is produced between these two electrodes. When a gas is passed through this electric field, it is ionized, producing ions, electrons, neutrals and UV photons. All of these species collide with the micro-organism, reacting with its essential elements and forming harmless compounds. The micro-organism, when devoid of its essential elements becomes harmless and is thus inactivated. More information about preliminary experiments on this subject is given here.

Plasma Physics, for the most part,deals with a lot of non-linear partial differential equations. The governing equations defining the behaviour of electrons, ions and neutrals inside a plasma are often expressed in differential forms and resemble closely the governing equations for fluid mechanics. Some of the differential equations are 1st order while some others are 2nd order but inevitably problems of non-linearity of PDEs are experienced. For a good explanation of the basic principles of Plasma Physics, refer to this site :Fundamentals of Plasma Physics

Abstract
One of the oldest problems that has sparked a lot of research in plasma physics is the plasma-wall problem. The sheath is a region of electric field that exists between the plasma and the boundary into which it runs. The basic problem of plasma running into a wall needs to be understood in order to better visualize some of the other aspects of plasma physics. In this paper, we attempt to understand the processes happening in the sheath region. There exist a set of governing equations, akin to the governing equations of fluid dynamics that outline the physical processes running in the plasma. These equations are discretized using a simple, first order finite difference scheme and then solved using the Newton-Raphson Method. The boundary conditions applied however are those pertaining to the sheath region existing between plasma and a solid wall. The objective in conducting this analysis is to investigate the accuracy of the numerical scheme used. If the numerical scheme is accurate, the electron density at the wall should go to zero. Further a probability density function is constructed using electron density data generated from the solver code.

Introduction
To completely describe the state of plasma, we would need to write down all the particle locations and velocities and describe the electromagnetic field in the plasma. But considering the zillions of particles, this would be an impossible task! In order to make our life easier, we construct models which simplify the processes taking place in plasma.

Models are of two types:
 * Fluid Model- This describes plasmas in terms of smoothed quantities like density and averaged velocity around each position. A general description is the two-fluid picture, where ions and electrons are described separately. Fluid models usually model the particles using a Maxwell-Boltzmann Distribution. However they are usually unable to capture accurately finer structures such as beams or double layers.
 * Kinetic Model- While being more accurate than the fluid model, it is usually more computationally intensive. It describes the particle velocity distribution function at each point in the plasma and therefore does not need to assume a distribution function like the fluid model. The more commonly used Particle-In-Cell (PIC) technique includes kinetic information by following the trajectories of a large number of individual particles

Understanding plasma sheaths is perhaps one of the oldest problems in plasma physics. A plasma sheath is the localized electric field that separates plasma from a material boundary. It confines the more mobile species (electrons) in the plasma and accelerates the less mobile species (ions) out of the plasma and toward the walls.This means that in the typical case, where the electrons are more mobile than the positively charged ions, the electric field in the sheath points towards the boundary.This behaviour can be further influenced by external processes such as variations in the applied electric field or the imposition of an external magnetic field.

Modeling of the behaviour of species in the Sheath Region
Many numerical models have been developed to describe Sheaths. The simplest visualization of this process is the behaviour of electrons and ions in the sheath region. Consider a planar slab of collisionless plasma, modeled in 1D space, using a fluid approach and assuming a Boltzmannian electron distribution. The governing equations for such a model are given below:

Continuity Equation for Ion flux:
 * $$ \frac{\partial (n_{i}v_{i})}{\partial x}= zn_{e}$$ ---(1)

Continuity Equation for Electron flux:
 * $$ \frac{\partial (n_{e}v_{e})}{\partial x}= zn_{e}$$ ---(2)

Momentum of ions: 
 * $$ \frac{\partial (Mn_{i}v_{i}^{2})}{\partial x}+en_{i}\frac{\partial(\psi)}{\partial{x}}= 0 $$ ---(3)

Momentum of electrons:
 * $$ \frac{\partial (mn_{e}v_{e}^{2}+kn_{e}T_{e})}{\partial x}-en_{e}\frac{\partial(\psi)}{\partial{x}}= 0$$(4)

Poisson Equation:
 * $$ \epsilon_{0}\frac{\partial^2 \psi}{\partial x^2}= -e(n_{i}-n_{e})$$ --(5)

where k is the Boltzmann Constant, m is the electron mass, M is the ion mass, Te is the electron temperature,ve is the electron velocity, viis the ion velocity, $$\psi$$ is the electric potential, z is the frequency of direct ionization, e is the electron charge, $$\epsilon_{0}$$ is the permittivity of free space, niis the ion density and ne is the electron density.

More on the Poisson Equation The Poisson's Equation given by $$\nabla^{2}\phi= E$$ gives way to the Laplace's Equation for vanishing 'E'. The Poisson's Equation is one of the cornerstones of Electrostatics. The derivation of Poisson's Equation from basic concepts of electrostatics is a fascinating exercise and is given below:


 * Gauss' Law of Electricity (in a differential control volume) $$\Rightarrow \nabla.D= \rho_{f}$$

where D= electric field displacement, $$\rho_{f}$$ describes the external charge density


 * For a linear, isotropic and homogenous medium, D= $$\epsilon E$$

where $$\epsilon$$ is the permittivity of medium and E is the electric field.


 * Substituting the above equation into Gauss' Law $$\Rightarrow \nabla.E= \frac {\rho_{f}}{\epsilon}$$


 * Faraday's Law of Induction(for a constant magnetic field of strength B)$$\Rightarrow \nabla \times E= -\frac{\partial B}{\partial t}$$

When the curl of a vector is zero, it can be expressed as the gradient of a scalar. Express
 * $$\Rightarrow E= -\nabla \phi$$, where $$\phi$$ is defined as electric potential(scalar)

$$\Rightarrow \nabla^{2}\phi= -\frac {\rho_{f}}{\epsilon}$$

Defining Charge $$\rho_{f}$$ as the product of the charge of a single electron multiplied by the number density of electrons inside the plasma, we get the form shown in (5).

More details on these equations and their non-dimensionalization can be obtained from the paper by Sternberg and Godyak. For our purpose, we take the normalized versions of Equations (1)-(5) used in the paper by


 * $$\frac{\partial{(y_{i}V_{i}})}{\partial{x}}= y_{e}= exp(-\phi)$$--(5)
 * $$ V_{i}\frac{\partial{V_{i}}}{\partial{x}}-E+V_{i}\frac {y_{e}}{y_{i}}=0$$---(6)
 * $$ \frac{\partial{E}}{\partial{x}}- \kappa^{2}(y_{i}-y_{e})= 0$$ --(7)
 * $$ E= \frac{\partial{\phi}}{\partial{x}}$$(8)

where $$y_{i}$$ is the normalized ion density, $$y_{e}$$ is the normalized electron density, $$V_{i}$$ is the normalized ion velocity, E is the electric field strength and $$\phi$$ is the electric potential.

Here the variables solved for will be $$y_{i}, V_{i}$$ and E or $$\phi$$ (E can be calculated from $$\phi$$ from (8)). Knowing $$y_{i}$$, $$y_{e}$$ can be calculated (from (5)).

 It is useful for the general readers to identify the unknowns to be solved for by numerical methods as described further below. Egm6322.s09 17:57, 1 May 2009 (UTC) Unknowns identified here Egm6322.s09.Three.nav 13:15, 4 May 2009 (UTC)

Boundary Conditions
It is assumed with sufficiently high accuracy that at the center of the plasma ni= ne. This is termed a quasi-neutral plasma.

Debye Length & Quasi-neutrality of Plasmas The characteristic length scale in a plasma is the electron Debye length, $$\lambda_{DE}$$. It is typically the distance scale over which significant charge-densities can exist.

However across a plasma of length L>>>$$\lambda_{DE}$$,

Eqn.(5) $$\Rightarrow \nabla^{2}\phi= -\frac{e(zn_{i}-n_{e})}{\epsilon}$$

We generally expect that $$\phi \leq T_{e}= \frac{e n_{e}\lambda_{DE}^{2}}{\epsilon}$$

$$\Rightarrow \frac {|zn_{i}-n_{e}|}{n_{e}} \leq \frac {\lambda_{DE}^{2}}{L^2}$$

For L>>>$$\lambda_{DE}$$, $$\frac {\lambda_{DE}^{2}}{L^2} \leq 1$$

$$\Rightarrow \frac {|zn_{i}-n_{e}|}{n_{e}}\leq 1$$

$$\Rightarrow |zn_{i}-n_{e}| \leq n_{e}$$

$$\Rightarrow zn_{i}= n_{e}$$

Which defines a quasi neutral plasma. This assumption helps simplify the solution of the macro-scopic equations listed from (1)-(5) greatly.

The boundary conditions at the plasma center are specified as follows:

$$ n_{i}(0)= n_{0}  = n_{e}(0) $$

$$ v_{i}(0)= 0= v_{e}(0)$$

$$\psi(0)= 0= \frac{\partial{\psi}}{\partial{x}}(0)$$

The normalized boundary conditions are then going to be as follows:

$$y_{e}(x=0)=0$$ - Defining the electron density at the start of the plasma

$$\phi(x=0)=0$$- Defining the electric potential at this location

$$y_{i}(x=0)=1$$- Defining the ion density at the start of the plasma

$$V_{i}(x=0)=0$$- Defining the ion-velocity at the start of the plasma.

Remember: 'Start of the plasma' implies the starting point, at which there is no plasma generated. Obviously, here the gas is still not fully ionized, wherein ideally, the ion density, electron density and ion velocity are zero.

Solution through Numerical Schemes
A finite difference scheme is applied to discretize the first derivatives in Equations(5)-(7). Equation(8) is just needed to calculate the electric potential $$\phi$$ Finite difference schemes are usually applied for complex, non-linear partial differential equations, wherein the solution space is divided into n number of grid points and the required value at each grid point is calculated. For the single derivative, there are typically three types of schemes: Forward difference, Backward difference and Central Difference. Each of these have a different order of accuracy and are applied depending on how accurate the numerical scheme applied has to be.

Forward Difference$$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n+1}-A_{n}}{\delta x}$$

Backward Difference$$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n}-A_{n-1}}{\delta x}$$

Central Difference $$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n+1}-A_{n-1}}{2\delta x}$$

where An+1 is the value of A(the quantity being calculated) at the n+1th grid point and An-1 is the value of A at the n-1th grid point.

The Central Difference Scheme usually has an order of accuracy higher than the other two schemes but this does not necessarily mean that using this scheme might give us the most accurate solution. The nature of the PDE being solved should also be taken into account. This is where we define terms like 'stability' and 'consistency' of a finite difference scheme. For more information, the book by Anderson is a very informative book to read. The basics of fundamental numerical schemes as well as more complex ones are outlined along with illustrative applications to common problems in Fluid Mechanics like the Burger's Equation.

Newton Raphson Method

Consider a typical problem that gives N functional relations to be zeroed involving variables xi, i=1,2....N.

f(xi)= 0, i=1,2....N

Each of these functions can be expanded into a Taylor Series i.e $$ f(x_{i}+\delta x_{i})= f(x_{i})+ \delta x_{j}\frac {\partial{f_{i}}}{\partial{x_{j}}}+O(\delta x^2)$$

Assuming $$f \to 0$$, we get
 * $$ \Sigma \alpha_{ij} \delta x_{ij}= \beta_{ij} \ where \ \beta_{ij}= -f_{i} $$

$$ \alpha_{ij} $$ is known as the Jacobian. Here Equations 5-7 are the equations that are discretized and contain three unknowns. So $$\alpha$$ is going to be a 3x3 matrix and $$\beta$$ will be a 3x1 matrix. $$\delta x_{ij}$$ is then solved for iteratively in MATLAB.


 * Hence once Equations (5)-(7) are discretized using forward difference, they reduce to a set of non-linear equations in which the unknowns are yi, Vi and E of the form


 * Defining Vi= u, yi= y and ye= ne (for simplicity), Eqns.(5)-(7) when discretized reduce to the following equations


 * $$\frac {y_{n}(u_{n+1}-u_{n})+u_{n}(y_{n+1}-y_{n})}{\delta x}= ne = exp(-\phi)$$


 * $$\frac {u_{n}(u_{n+1}-u_{n})-E+\frac {u_{n}ne}{y_{n}}}{\delta x}= 0$$


 * $$\frac {\phi_{n+1}- 2\phi_{n}+\phi_{n-1}}{\delta x^2}= \kappa^2(y_{n}-ne)$$

where the subscript 'n+1' implies the (n+1)th grid point, 'n' implies the nth grid point and 'n-1' implies the (n-1)th grid point.

This set of non-linear equations is solved using the Newton-Raphson Method,described above.

 You can also linearize the PDEs in Eqs.(5)-(7); the linearized PDEs would look close to an elliptic PDE, particularly after you replace the linearized Eq.(6) into Eq.(7). Egm6322.s09 17:57, 1 May 2009 (UTC)

Hence we get a matrix of type Ax= B where A= $$\alpha_{ij}$$ and B= $$\beta_{ij}$$ as described above. Once A and B are calculated, x= A-1B. This process is repeated for several iterations in order to calculate the necessary plasma parameters. Once these are known, ye can be calculated using Equation(5).
 * Now, once the electron density distributions and ion-density distributions are calculated, how do we estimate what happens down the line, say after like 100 iterations? Does the electron density in the Sheath Region remain zero always, as it should? Does the ion-density increase or decrease? This check achieves two purposes: The accuracy of the numerical scheme is verified and the sheath behaviour after a long time is understood. This check is performed with the help of a probability density analysis. A MATLAB code was written to simulate the entire process described above. This code furthe takes the electron-density data generated, divides it into 100 blocks and constructs a probability density function. This is repeated for every iteration so that we know the nature of the pdf after about 100 iterations.

Results
 Remember to put a license to each image that you uploaded to prevent their removal due to a lack of a license. Egm6322.s09 17:57, 1 May 2009 (UTC)

Licenses inserted for all images Egm6322.s09.Three.nav 13:46, 4 May 2009 (UTC)



Fig(1) shows the published results from the paper by Roy et.al(2004) In this paper, equations (1-7) to (1-10) are solved using a second-order time accurate, implicit finite element method using a computational domain extending from x=0 to the wall. Fig. 1(a) is a plot of ion density versus position x. Fig. 1(b) is the plot of ion velocity as a function of position, x. Fig. 1(c) is the plot of electric potential as a function of x, while Fig. 1(d) represents the plot of electric field versus x. As expected, since the more mobile species is confined within the plasma and the less mobile species is accelerated out, it is seen accordingly that at xw~0.71, the ion velocity increases. Since electrons are being confined, one would expect electric field at wall to rise (more charge being accumulated). This is seen in Fig. 1(c) and consequently in Fig. 1(d).

In Fig.(2), the plasma parameter distributions obtained from the MATLAB code used are given. In fig. 2(a), as can be seen, ion density prediction matches the trend as predicted in Fig.1(a) i,e ion density decreases in the direction of the wall. This is expected. Fig. 2(b), which shows the plot of ion velocity versus x also seems to match the trend predicted in 1(b). It makes senses too because one would expect ion velocity to increases as sheath approaches. Fig.2(c) makes perfect sense because electron density decreases to zero as wall is approached. The value of x at which ne goes to zero is a little earlier than expected (at x=0.2). Also to be noted is the important fact that ion density at wall does not go to zero while electron density does go to zero, which simply confirms the fact the electrons stay behind in the plasma while ions accelerate towards the wall. Fig. 2(d) is a plot of the potential versus x, which also matches the plot shown in fig.1(c). This tells us that the code is predicting the behaviour of the species in the sheath correctly to some extent. The linear trend present in fig. 2(a), 2(b), 2(d) is of slight concern as we are looking for more of an exponential trend. This could be because of an error in the code or due to slight errors introduced by the numerical scheme used. Note that the results in fig.(1) are the result of a finite element code while the results in fig. (2) are those computed from a finite difference scheme. Finite element schemes are much more robust and stable schemes and hence what needs to be done here is take a higher order finite difference scheme and apply it to the problem at hand.



Fig. 3 shows the pdf analysis for the ion density. Ion density data is taken from the code, divided into 100 blocks and a pdf is constructed. Fig.3(a) shows this pdf as a bar diagram. Further the convergence of this pdf is investigated by calculating moments as high as E[[x7]. Fig. 4 shows the results of a similar analyis for the electron density data. There are a couple of important points to note in these two figures. Fig.3(a) shows that the ion density does not go to zero at the wall. Fig. 4(a) shows that there is a 4% probability of electron density going to zero at the wall. Fig. 3(b) and Fig. 4(b) confirm convergence of the pdf for both electron density and ion density data. Note also that the electron density is steady after about 40 iterations while the ion density is steady after about 20 iterations

Conclusion & Future Work
The governing equations for the plasma wall problem were solved using a finite difference scheme to discretize them. The electron and ion density are calculated using a MATLAB code. The computed distributions of ion and electron density exhibit the trend shown in the published results as shown in Fig.(1). The linear trend may be due to accuracy errors in solver code. However the ion and electron density plots are very similar to ion and electron behaviour in sheath region. Also the PDF analysis shown in Fig. (3) and (4) for electron and ion density data respectively shows that there is convergence after 100 iterations, proving that the results predicted by the code are accurate. Also we see that there is a 4% probability of the electron density going to zero. This value might further reduce, provided the accuracy errors in the code are resolved.

Future work involves using more complex, higher order numerical schemes to discretize the governing equations. Also a simplistic case of the plasma-wall phenomenon was examined here. More complex phenomenon can be introduced into this simulation in order to further understand this phenomenon.

Since understanding Plasma Physics requires a thorough understanding of PDES, their types, solution methods used in solving them, this project was a good review of the PDE principles that were discussed in class. Principles, discussed in class, such as the Laplace's Equation and the significance of the 'gradient' operator, 'curl' and so on were good tools in providing a better understanding of the theory behind the vast subject of Plasma Physics.

Signature
Egm6322.s09.Three.nav 20:40, 30 April 2009 (UTC)