Numerical Analysis Project

Moshe Damkani 309935815 24/01/2012

Abstract This project involves the analysis of plasma, which is essentially an ionized gas composed of positive ions and free electrons, as it comes into contact with a metal wall. Due to its conductivity, the lighter electrons rush towards the neutral wall and give it an electric charge and creating an electrostatic field which repels any additional electrons and attracts the nearby positive ions. Giving the plasma kinetic energy causes ions that were far away from the wall to be forced past a line called the sheath edge into the sheath area, where they now begin to feel the electric field and are pulled towards the negatively charged wall. For a critical value of kinetic energy given to the plasma, the density of ions becomes greater than that of the electrons for the whole length of the sheath area. A major focus of this project is the examination of this kinetic energy and how it relates to density distributions and the behavior of the electric field in the sheath area but, more importantly, is the exposition and explanation of the numerical methods used to reach such conclusions.

Table of Contents Physical Problem Mathematical Model Numerical Methods Numerical Results Physical Results Summary and Conclusion Flowchart Code 4 5 6 10 14 18 19 20

Physical Problem The goal of this project is to analyze the behavior of plasma as it come into contact with a metal wall. Plasma is a state of matter, similar to gas, whose particles are ionized. It is thus composed of both positive ions and negative electrons in equal proportions as to maintain a uniformly neutral charge. As plasma has a significant amount of charged particles, it is electrically conductive and responds strongly to differences in electric potential. When the plasma comes into contact with the metal wall, its free electrons, being lighter than the ions, rush towards the wall, giving the wall a negative charge. This now negatively charged wall creates a steady electrostatic field which both repels any additional electrons and attracts nearby positive ions. The effective force created from this field is significant in a range of plasma near the wall called the sheath layer, and ions that are able to cross the sheath edge into this range will be affected and pulled towards the wall. By giving the ions in the pre-sheath plasma a certain kinetic energy (expressed in velocity), they are able to cross the sheath edge into the electric field. As this kinetic energy increases, the density of ions pushed towards the wall also increases and will surpass the density of electrons there for a certain critical value. The focus of this project will be directed towards examining the relationships between this kinetic energy, the resulting density distribution of the plasma, and the behavior of the electrostatic forces in the sheath area.

Mathematical Model In order to create a working mathematically model, three unit-less parameters are suggested and defined as follows:

x 0 kTe ne e2 e0 kTe

normalized position; increases with the distance from the sheath edge ( x = 0 ) towards the wall. normalized kinetic energy given to the ions. electric potential energy, normalized and negative; increases as potential decreases (receiving a larger value near the negatively charged wall).

0 =

e kTe

0 = 8.85 10 12
e = 1.61 10 19 k = 1.38 10 Te
23

dielectric constant electron charge Boltzmann constant electron temperature

F m

[C ]
m 2 kg s2 K

[K ]

These parameters simply the expression of the governing differential equations:

F6

d2 = 1 + 2 d 0

exp ( )
1/2

F7

1/2 d = 4 0 1 + + 2 exp ( ) 2 ( 2 0 1) , ( = 0 ) = 0 d 0

where F7 represents the electric field and F6 represents its gradient.

Numerical Methods The classical fourth order Runge-Kutta method was used in order to solve the given differential equation for the normalized electric potential (6). The full method for solving an initial value problem of two variables is given below:

k1 = hf ( xn , yn ) h k k2 = hf xn + , yn + 1 2 2 h k k3 = hf xn + , yn + 2 2 2 k4 = hf ( xn + h, yn + k3 ) yn +1 = yn + 1 ( k1 + 2k2 + 2k3 + k4 ) 6

y' = f ( t, y ) y ( t 0 ) = y0

Each value, yn +1 , is determined by the present value, yn , plus the weighted average of 4 deltas, where each delta is the product of the size of the interval, h , and the slope, f = y' . As our differential equation (6) is expressed as

d = f ( ) , and is not a function of independently, the method can be d

formulated as follows, as a function of one variable:

k1 = hf ( n ) k k2 = hf n + 1 2 k k3 = hf n + 2 2 k4 = hf ( n + k3 ) 1 ( k1 + 2k2 + 2k3 + k4 ) 6

n +1 = n +

This simplification makes the RK4 method especially suitable for this analysis, and as it is an explicit iterative method, it is also easy to code. That is to say that each step is a function solely of the value preceding it, with the first value being built off of the initial condition. The accuracy of this method is extremely high, having a local error of order O h 5 and an accumulated error of order O h 4 .

( )

( )

The False Position Method is an iterative method used to find the zero of a function f in a certain domain [ m1 , m2 ] , and it is utilized frequently in this project due to its efficiency. The function is evaluated at both the upper and lower bounds, and a line is made between two points ( m1 , f ( m1 )) and ( m2 , f ( m2 )) whose intersection with the independent axis becomes an approximation for the functions root. This method is formulated by equating the slopes y / x between two pairs of points that lie on the line and then imposing f ( M ) = 0 such that M may be the functions root.

f ( m1 ) f ( M ) f ( m2 ) f ( m1 ) = m1 M ( m2 m1 )

f ( m1 ) 0 f ( m1 ) ( m2 m1 ) m M = 1 M = m1 f ( m2 ) f ( m1 ) ( m2 m1 ) f ( m2 ) f ( m1 )

The bounds are then re-evaluated as follows: if f ( m1 ) and f ( M ) have different signs, then we can be assured that the function has a zero between m1 and M . If their signs are the same, then they are on the same side of the axis and there is no zero between them, so the zero must be between M and m2 . This algorithm is coded as follows:
if f(m1)*f(M)<0 m2=M; else m1=M; end

The new domain [ m1 , m2 ] will now give a better approximation for M , and this process is repeated until the method converges. Convergence is satisfied as each successive approximation becomes more similar to the previous approximation, until finally M M prev error. becomes smaller than a demanded

The method of Least Squares is used for fitting a mathematical model to a given set of data. The models constants are chosen as to minimize the sum of the squares of the errors between the actual values i and the models approximations ( i ) at each point i . This is done by creating such an expression, S , and setting its partial derivatives with respect to each constant equal to zero. The model given in this project and its subsequent derivations are exposited below:

( CH )3/2 = A 2
( ) = A 1/2 ( CH )
S = ( ( i ) i )
i =1 n +1 2 3/ 4

S= A
i =1

n +1

1/2

( i CH )

3/ 4

0=

S A

0 = 2 A 1/2 ( i CH )
i =1

n +1

3/ 4

3/ 4 1 i ( i CH ) A 3/2 2

0 = A 1/2 ( i CH )
i =1

n +1

3/ 4

i ( i CH )

3/ 4

A1/2

0 = ( i CH ) i =1 A1 =
i =1 n +1 i =1

n +1

3/2

A1/2i ( i CH )
1.5

3/ 4

(
i

n +1

CH )

CH )

0.75

0=

S CH

0 = 2 A 1/2 ( i CH )
i =1

n +1

3/ 4

1/ 4 3 i A 1/2 ( i CH ) ( 1) 4

0 = A 1/2 ( i CH )
i =1

n +1

3/ 4

i ( i CH )

1/ 4

A1/2

0 = ( i CH ) i =1 A2 =

n +1

1/2

A1/2i ( i CH )
0.5

1/ 4

( (
i =1 n +1 i =1 i

n +1

CH )

CH )

0.25

i
8

The two equations for

A are then equated and expressed as f ( CH ) = 0 . This single variable

function can now be solved using the False Position Method, and the CH

plugged into either expression for A to find A . The two functions used in this project for these calculations are arbitrarily named and are defined as follows:

gadget(E,X,Xch)

A1 A 2 - gadget(E,X,Xch)

monster(E,X,Xch) =

Another noteworthy numerical method for data fitting involves formulating the model as a linear equation (with a quick change of variables) as follows,

= A 2 / 3 4 / 3 + CH
y = ax +b

and solving a handful of determinants to arrive at each constant. The writer, however, feels that this method is far less intuitive and, not caring to derive it within the space provided, has used the previous method in its stead.

Numerical Results Runge-Kutta method The step value h for this project was chosen based on the value final function, which begins to stabilize after log ( h ) = 2.5 given by the RK4

Figure 1

In this project n = 800 :

h=

(upperbound lowerbound ) = ( 2 0 ) = 0.0025

log ( 0.0025 ) = 2.6021 > 2.5

and so, all sub-sequential utilization of the RK4 function is valid.

False Position Method The False Position method is used in this project for approximating a critical normalized kinetic energy 0C and is used again within the Least Squares Method to approximate a constant CH . The accuracy of this method is dependent upon the error required to satisfy convergence.
10

Approximating 0C

The ion and electron densities ni and ne are first normalized using the given parameters:

ni = n0 1 0

1/2

F1

ni ( ) = 1 + n0 0 ni ( ) = exp ( ) n0

1/2

e ne = n0 exp kTe

F2

For part C it is asked to find a certain critical value 0C above which F1 will always be greater than F2 on the interval 0 < < 2 . To do this a function is created:
critdif(X,X0) = F1(X,X0) - F2(X);

which should return positive values for all 0 > 0C . Now, noting that F2 is not a function of 0 , and rearranging F1 as follows...
1/2

+ F1 = 0
0

0 0 +

...it can be seen that, for a given , critdif(X,X0) increases with 0 . Therefore, any 0 greater than the 0 for which critdif(X,X0) equals zero will return a critdif(X,X0) greater than zero, and thus satisfies the requirements for 0C . This value can be found using the false position method, solving
critdif(X,X0)=0

as a function of 0 with a fixed value of . It can also be seen above that F1

increases as decreases. Therefore, F1 is maximized as approaches 0

critdif(X,X0) is maximized as approaches 0 0C is maximized as approaches 0, as 0 increases with critdif(X,X0)

11

As seen in the graph, the method converges and converges to 0C = 0.5 .

Figure 2

12

Least Squares

Approximating CH The accuracy of the least squares method used in this project is determined again by the error required to satisfy convergence in the False Position Method, as solving CH is done by finding when the following function equals zero, as derived earlier in the Numerical Methods section:

n +1 n +1 1.5 0.5 i CH ) ( ( i CH ) n +1i =1 n +1i =1 =0 ( )0.75 ( )0.25 CH i CH i i i =1 i =1 i

As seen in the graph, the method converges and converges to CH = 0.9375 .

Figure 8

13

Physical Results The main focus of this project is to examine the behavior of the plasma in the sheath layer next to the negatively charged wall. As the kinetic energy given to the plasmas ions increases, more and more ions are able to cross the sheath edge into the electric field. The resulting distributions of the ion and electron densities in this area are shown below as functions of the normalized potential for the kinetic energies 0 = 0.1 and 0 = 2 . Remember that the normalized potential increases as potential decreases and receives a larger value near the negatively charged wall. So for 0 = 0.1 the normalized ion density F1 only surpasses the normalized electron density F2 close to the wall, while for 0 = 2 the ion density is always greater. The density distribution of the electrons in the sheath layer is not effected by changes in the kinetic energy given to the plasma, as any further electrons crossing the sheath edge are being repelled by the negative potential difference, while the ions are attracted and able to accumulate.

Figure 3

Figure 4 14

The critical value 0C above which the ion density F1 will always be greater than the electron

density F2 was shown to converge to 0.5 and is graphed below.

Figure 5

15

The Electric Field is graphed below for each 0 . For all three values there is a tendency for the field to become more and more negative (as the normalized potential increases), due to the pull exerted by the negative potential difference at the wall. However, for 0 values below 0C , there also exists a force in the opposite direction because for these values there still exists portions in the beginning of the sheath layer with a greater density of electrons, which are pulled away from the wall. For the critical value 0C the densities of electrons and ions are just equal at the sheath edge, and so the forces cancel momentarily, giving the initial constant slope seen below.

Figure 6

The second derivative of the normalized potential is the Electric Field Gradient. It shows again that for 0 values below 0C , the electric field begins with a negative slope.

Figure 7

16

In the last section of the project, a model for the normalized position as a function of the normalized electric potential was fitted to the actual values given by the RK4 function. The optimal value for one of the models constants, CH , was shown to converge to 0.9375 . The second constant,

A , can now be solved as a function of CH using

A=

(
i =1

i =1 n +1

(
i

n +1

CH )

1.5

CH )

0.75

and is found to equal 1.6942 . The model

( ) = A 1/2 ( CH )

3/ 4

= 1.6942 1/2 ( 0.9375 )

3/ 4

is graphed below along with the actual data and is indeed shown to be an appropriate approximation.

Figure 9

17

Summary and Conclusion This project involved the analysis of the electric fields and potential differences created in plasma due to varying density distributions of the its ions and electrons. The numerical methods used include the Runge-Kutta method for solving ordinary differential equations, the False Position method for finding a functions zeros, and the Least Squares method for optimizing a given mathematical model to fit real world data. These methods demonstrate the power and usefulness of numerical approaches in solving complex physical problems.

18

19

Code Structure main.m RK4.m falseposition.m Runge-Kutta method for solving the ordinary differential equation, F7 False Position iterative method for finding the root of any given function between two given bounds for a given error of convergence formula for calculating the normalized ion density formula for calculating the normalized electron density formula for calculating the electric field gradient formula for calculating the electric field equals the difference between F1 and F2 and used along with falseposition.m to find the critical normalized kinetic energy 0C a function of CH which equals A 0.5 a function of CH which should equal zero and is used along with falseposition.m to find CH twenty.m used along with falseposition.m for finding the value of which corresponds with = 20

F1.m F2.m F6.m F7.m critdif.m

gadget.m monster.m

20

clc clear all; close all; n=800; upperbound=2; lowerbound=0; %[ Stability of RK4 ] figure(1); set(gcf,'color','w') k=1:0.05:3.5; h = 10.^-k; for i=1:length(h) P(i) = max(RK4(0,0.1,(upperbound-lowerbound)/h(i),h(i))); end plot(k,P,'o-') title('Stability of the Runge-Kutta method') xlabel('-log(h)'); ylabel('\chi_f_i_n_a_l'); grid on h=(upperbound-lowerbound)/n; x=lowerbound:h:upperbound; %[ Finding X0c ] k=1:0.1:6; for i=1:length(k) error = 10^-k(i); fhandle=@critdif; X0c(i) = falseposition(fhandle,0.1,2,error,0,0.00001); end figure(2); set(gcf,'color','w') plot(k,X0c,'o-') title('Convergence of \chi_0_C Appoximations with Decreasing Error') xlabel('-log(error)'); ylabel('\chi_0_C'); grid on X0c=X0c(i-1); %[ Graphing normalized ni and ne ] for X0=[0.1 X0c 2;3 4 5]; figure(X0(2,:)) set(gcf,'color','w') plot(x,F2(x,X0(1,:)),'k',x,F1(x,X0(1,:)),'r') title(['\chi_0 = ',num2str(X0(1,:))]) legend('n_e(\chi)/n_0','n_i(\chi)/n_0') xlabel('\chi'); ylabel('normalized ion and electron densities'); grid on end

21

%[ figure 6: graphing the derivative (gradient) for each X0 value ] %[ figure 7: graphing the 2nd derivative (e-field) for each X0 value for i=6:7 figure(i);set(gcf,'color','w') if i==6; fhandle=@F7; else; fhandle=@F6; end

D1=fhandle(RK4(0,0.1,n,h),0.1); D2=fhandle(RK4(0,X0c,n,h),X0c); D3=fhandle(RK4(0,2.0,n,h),2.0); plot(x,D1,x,D2,'m',x,D3,'r'); xlabel('\xi'); grid on legend('\chi_0 = 0.1',['\chi_0_,_c = ',num2str(X0c)],'\chi_0 = 2.0') if i==6; ylabel('d\chi/d\xi'); title('Electric Field') else; ylabel('d^2\chi/d\xi^2'); title('Electric Field Gradient') end end

%[

Best Fit Function

%finding Emax, such that X(Emax)=20, as we have no formula solved for E(X) fhandle=@twenty; Emax=falseposition(fhandle,2,20,0.001,X0c,n); h=(Emax-0)/n; X=RK4(0,X0c,n,h); E=0:h:Emax; %Creating vectors X6 and E6 satisfying 6<X<20, to use for fitting X6=X(find(X>6,1):length(X)); E6=E(find(X>6,1):length(X)); %Approximating Xch, demanding successively smaller error... fhandle=@monster; k=1:(5-1)/20:5; for i=1:length(k) error = 10^-k(i); Xch(i)=falseposition(fhandle,0,5.9,error,E6,X6); end figure(8); set(gcf,'color','w') plot(k,Xch,'o-')% ...and plotting the results. title('Convergence of \chi_C_H Appoximations with Decreasing Error') xlabel('-log(error)'); ylabel('\chi_C_H'); grid on %Solving A, creating Fitted Equation A=gadget(E6,X6,Xch(i))^2; E1=(((X-Xch(i)).^1.5)/A).^0.5; j=find(imag(E1)==0,1);%for plotting real values only figure(9); set(gcf,'color','w') plot(X,E,X(j:length(X)),E1(j:length(X)),'--r') legend('actual \xi','fitted \xi(\chi)') title('Comparing actual and fitted distances as functions of potential') xlabel('\chi'); ylabel('\xi','rot',0); grid on

22

%fourth order Runge-Kutta method function [x]= RK4(x,x0,n,h) %input f(1) = initial condition for i=1:n k1 = h*F7(x(i),x0); k2 = h*F7(x(i)+0.5*k1,x0); k3 = h*F7(x(i)+0.5*k2,x0); k4 = h*F7(x(i)+k3,x0); x(i+1) = x(i) + (k1+2*k2+2*k3+k4)/6; end end

function [M] = falseposition(fhandle,m1,m2,error,E,X) M=m1; Mprev=M+error+1; while abs(M-Mprev) > error; Mprev = M; M = m1-fhandle(E,X,m1)*(m2-m1)/(fhandle(E,X,m2)-fhandle(E,X,m1)); if fhandle(E,X,m1)*fhandle(E,X,M)<0; m2=M; else m1=M; end end end

23

%normalized ion density function [ans] = F1(X,X0) ans=(1+X./X0).^-0.5; end

%normalized electron density function [ans] = F2(X,X0) ans=exp(-X); end

%electric field gradient function [ans] = F6(X,X0) ans=(1+X/X0).^-0.5-exp(-X); end

%electric field function [ans] = F7(X,X0) ans=(4*X0*(1+X/X0).^0.5+2*exp(-X)-2*(2*X0-1)).^0.5; end

function [ans] = critdif(dummyvar,X,X0) ans = F1(X,X0)-F2(X,X0); end

function [ans] = gadget(E,X,Xch) ans = sum((X-Xch).^1.5)/sum(E.*(X-Xch).^0.75); end

function [ans] = monster(E,X,Xch) ans= sum((X-Xch).^0.5)/sum(E.*(X-Xch).^-0.25) - gadget(E,X,Xch); end

function[ans]=twenty(X0c,n,Emax) h=(Emax-0)/n; X=RK4(0,X0c,n,h); ans=max(X)-20; end

24