N O T I C E 

THIS DOCUMENT HAS BEEN REPRODUCED FROM 

MICROFICHE. ALTHOUGH IT IS RECOGNIZED THAT 

CERTAIN PORTIONS ARE ILLEGIBLE, IT IS BEING RELEASED 

IN THE INTEREST OF MAKING AVAILABLE AS MUCH 

INFORMATION AS POSSIBLE 
FUNDAMENTAL STUDIES ON A

HEAT DRIVEN LAMP

FINAL REPORT FOR NASA GRANT NAG 3-437

(NASA-Ch- 176361 ) F UbLAMEhIAL STUDIES CN A h86-15122 HEAT DRIVEN LAME Final ficicrt, 1 Sep. 1983

- 1 Sep. 1984 1 CaLL cg j e — M el1cc UDiv.) 40 p hL AU3/Mf A01 CSLL 201 Unclds

G3/75 04699

September 1. 19 F 3 - September 1. 1984

Submitted by the

Department of Mechanical Engineering

Carnegie-Mellon University

Pittsburgh, Pennsylvania 15213

John L. Lawless

Principal Investigator

Y, 4 - y^

DEC 1985

RECEIVED

MASH Sn FACIt^ -DAL -

i
'14

i

Table of Contents

Abstract

1. Introduction

II. Papers and Presentations

III. Student Participation

APPENDIX A: IEEE Conference Abstract

APPENDIX B: F. Stefani's Master's Thesis

1

2

2

2

3

4

Y
I

Abstract

A detailed theoretical study of a heat-driven lamp has been performed. This lamp uses a

plasma produced in a thermionic diode. The light is produced by the resonance transition of

cesium. An important result of this study is that up to 30% of the input heat is predicted to be

converted to light in this device. This is a major improvement over ordinary thermionic energy

converters in which only —1% is converted to resonance radiation. Efficiencies and optimum

inter-electrode spacings have been found as a function of cathode temperature and the radiative

escape factor. The theory developed explains the operating limits of the device.
2

I. Introduction

The topic of this work is the direct conversion of heat to light. The method investigated is

to use the cesium plasma of a thermionic discharge to produce resonance radiation. The

motivation for this study is to determine the feasibility of the tbermionic-photovoltaic energy

converter proposed by D. L. Chubb. The reason cesium is used is that the frequency of its

resonance radiation is closely matched to the band-gap of several semiconductors. This match

permits efficient operation of photovoltaic cells in strong contrast to the inefficient operation of

photovoltaics when powered by broadband solar radiation.

The major result of this work is that up to 30% of the input beat to the discharge may be

converted to resonance radiation even with very low escape factors for the radiation. A brief

description of the technical work is included in Appendix A. The details of the model, the

calculations. and the results are given in Appendix B.

I1. Apers and Presentations

• A paper on this work was presented at the 1985 IEEE International Conference on

Plasma Science. The abstract of this paper is included as Appendix A.

• F. Stefani completed his Master's Thesis on this work. It is included as Appendix B.

• A paper by F. Stefani and J. L. Lawless will be submitted to IEEE Journe/ of

P/ssms Science (in progress).

III. Student Participation

Student: F. Stefani

Period: 1983-1985

Degree: Master of Science

Thesis: See Appendix B
OFWINAL PAGE IS

QE FOOR QUALITY APPENDIX A: IEEE Conference Abstract

J

Presented at the 1983 IEEE International

Conference on Plums Science

Pittsburgh. Pennsylvania

June, 1983

Plasma Dynamies

In a Thermionio-Photovoltaic

Energy Converters

F. Stefant

J.L. Lawlew

Gruesk-Mellon University

pitlrbnfdb. P*MWIVOaia

D.L. Cbabb

NASA Lewis Research

Cleveland, Ohio

The naoaanee radiation from a tbermlosic discharge can be

sued to generate electricity by employing pbotoroltaie calls

to collect this eneW2. Since the A resosaace tine is well

matched to at band-gap of many uemi000duebn, light-to.

electricity sonversios effleiencis of 110% or more are

possible. This paper will Inrertlpte bow aMently a

tbermbaie discharge aw be made to aoaves beat to light.

Normally, an optimised dwrmtowe soup eoaverter

produces very little roo naoae radiation. Our theoretical

An" show that at higher currents and Iooger Inter

doetrode gaps, a eignMeont portion of the Input beat can be

enverted to line radiation erea If the radlatiee escape factor

Is a loo a 1%. Neglecting lad lows and optical collection

Amcultia, efficiencies of up to 30% for the cos atioa of

bat to radiation In the discharge are predicted.

The tbeaetieal studies model the electron transport, ion

trassport, and isel"c eollidow lading to iosisation,

naoeabluation, sad ankh" at rmonance radiation. The

transport of "nuance radiation Is treated miss an oeape

factor. The emakdoo rate le coupled to the loulsetion and

recombination rates through the a welted state kinetics.

DIffudon and energy balances determine the aleetron

temperature and density distribution. The model makes

dear an mdntuitiee result that the maxtmam aMeieaey with

which Use radiatios is amlued is cey weakly dependent on

the optical escape factor. 1W efficiencies we typically

found at law-electrode dKanew an the order of A em

An apartment has has sonavocled to measure the

efficleney for as.erting dermal sheep to neasanee

radiation and aiso the maximum discharge length.

Wrursments are In progrew. Available esperimeatal

romans will be preasated.

1D.L. saw. n F1 1 0 ig /wow ewr.we.., raw M.

J11W iwasm" asferaw w pl.raa alma, p.s DJres, t]diferai.

(IOU)
APPENDIX R- F. Stdani's Master's Thesis
LINE RADIATION EMISSION FROM A CESIUM

THERMIONIC CONVERTER

Francesco Stefani

Department of Mechanical Engineering

Carnegie-Mellon University
TABLE OF CONTENTS

Introduction .....................................................................................................................1

Anoverview of the Model ............................................................................................1

Ambipolar Diffusion Equation .....................................................................................2

TransportEquation .......................................................................................................4

Energy equation ............................................................................................................6

Solution Procedure .......................................................................................................7

Results............................................................................................................................ B

Discussion......................................................................................................................9

References.....................................................................................................................11

Figures............................................................................................................»............12

Appendix A: Reaction Rate Constants for Excited State IGnetics

Appendix B: Boundary Conditions

Appendix C: Program Listing
INTRODUCTION

it been suggested (t) that if line radiation from a thermionic device could be

emitted with sufficient intensity, then a device employing photo cells to collect

this energy could be made to generate electricity efficiently. This follows from

the observation that the band gap energies of several common photo-voltaics

are well matched to the One radiation energies of cesium diode converters.

Resonance Radiation

6P112 - 6S12

1.39 eV

6P3,2 - 6S12

1.46 ev

Excitation Band Gaps

Silicon 1.12 eV

Indium Phosphide 1.27 eV

Gallium Aresenide 1.35 eV

Table 1 Radiation and band gap energies.

This study investigates the efficiency with which resonance radiation can be

generated from a thermionic devico and identifies the operating conditions for

which radiation output and the radiation efficiency are maximized. h has been

tentatively established that radiation efficiencies on the order of 30% are

possible even for optically dense plasmas. This result, if confirmed by

experiment, means that the prospects for the proposed device are indeed

promising.

The first part of this report is a discussion of the model that was employed. This

Is followed by a brief outline of the solution procedure. The last part of this

report is a presentation of the results and an attempt to relate them to the

physical events occurring in the plasma. A listing of the computer progran

developed for this study is Included in the appendix.

AN OVERVIEW OF THE MODEL

The work presented herein is based on an analytical model of a cesium plasma

diode. This model, described in detail in reference 2, is summarized here

primarily to review the assumptions made in the context of the problem and to

provide some insights that help make the 'results more understandable. .

A control volume is taken- as the repion of the inter-electrode space extending

from the top of the emitter sheath to the top of the collector sheath. Three

governing equations are written for the plasma in the control volume. These

are conservation equations for mass, momentum, and energy. A soluton

procedure is employed that permits the optical thickness of the gas to be varied
2

as a parameter. Since the problem is 1-D in the direction of current flow it Is

possible to reduce the governing equations to algebraic expressions. When

these are solved simultaneously for the electron temperature, the electron

density, and the operating current, a solution is obtained.

From the onset, it was assumed that a thermionic device capable a generating

high levels of radiation would be operating in the ignited mode. This means

that the ions and electrons In the plasma are created predominantly by electronatom collisions in the control volume rather than by contact ionization at the

walls. In this case emission effects are of relatively little importance. Of greater

Importance are the assumptions made about the plasma in the inter•electrods

space. This work has assumed the following:

1)The plasma is singly ionized and is composed of electrons, cesium

atoms, and cesium ions.

2) The plasma exists at two temperatures, one for the heavier particles and

one for the electrons.

3) Both temperatures are isothermally distributed in space.

The two temperature plasma assumption is based en the fact that most

collisions are elastic, and very little energy Is transferred from the fast moving

electrons to the slower moving heavier particles. The Isothermal assumption is

In effect a statement that the thermal conductivity of the plasma is high. This

assumption has been confirmed experimentally as well as through numerical

tests [3].

AMBIPOLAR DIFFUSION EQUATION

The Ambipolar Diffusion Eqn. is a mass balance for the electrons and the ions

ML in the control volume. It Is derived by writing the conditions for mass

conservation for the steady state plasma in terms of Ionization and

recombination rate constants S, and a.

VT. =vri = SN9.n - an3

(1)

To reflect conditions in the device transport equations for the electrons and ions

must include the effects of electric fields as well as diffusion effects. In ore

dimension they are written as follows:

r, •-μ,( n d'P/dx + (kT,/e) do/dx) (2)

ri =*,(-n cPF/dx + (kT^ /e) dn/dx) (3)
3

where: r. , T+ = electron and ion fluxes

μ, , μi = electron and ion mobilities

e = elementary charge

n = electron density = ion density

IF = electron motive

k = Boltzmann's constant

T9, Ti = electron and ion temperatures

Substituting the transport equation for the ions (eqn. 3) into the conservation

equation (eqn. 1) yields a diffusion equation for the ions.

- d/dx((kT,/e) do/dx - n dP'/dx) = S Nos n - an3

(4)

Substuting the electron transport (eqn. 2) into this expression eliminates the

electric field term and results in a conservation equation in terms of the coupled

motion of electrons and ions:

-D, d2n/dx2 +T.d(tL /μ.)Jdx-= (1 + μ i/μ.)(S No, n - an3) (5)

where D,=μi(kTi+kT,)/e (6)

This equation is called the Ambipolar diffusion equation and D. is the ambipolar

diffusivity. The ambipolar diffusion equation can be simplified by neglecting the

ratio of the ion and electron mobilities IL /μ. (since it is much less than one) and

by making the reasonable assumption that this ratio does not vary much across

the inter-electrode distance. Assuming also that D, Is constant reduces the

ambipolar diffusion equation to the following form:

- D. d2n/dx2 = S Nos n - an3

(7)

Solving the ambipolar diffusion equation subject to the appropriate boundary

conditi ons (Appendix B) yields the following solution and associated

eigenvalue condition:'

SNos = (4K2)Ds,/d2 + (1/2) a n2

(9)

where n(,,) Is the maximum electron density.

Both these equations are of interest and reveal different things about the

plasma. The solution to the ambipolar diffusion equation gives the distribution

of the electrons in the inter-electrode space. The eigenvalue condition gives

the correspondence between electron temperature and the maximum density
4

for which ionization and recombination are such to balance the loss of particles

tv '-e walls by diffusion. This correspondence is termed the ignited condition.

The function an that appears in the solution of the ambipolar diffusion equation

is the Jacobi elliptic function (4). It Is characterized by the quarterperiod K,

whose value can range from x/2 to infinity. Both limits represent Important

cases in the problem we are solving. In the limit where Ks x/2, the an function is

exactly the sin function and the electron distribution is similar to that shown in

fig. 1. This case corresponds to conditions where recombination Is negligible,

usually at lower values of the electron density. The other limiting case is that

where the electron density is very close to the Saha value. This condition is

referred to as the Saha limit. In this case the electron distribution departs

significantly from the sin profile, being much steeper at the boundaries and flat

at the center, as in fig. 2.

In non-Saha conditions the quarter period Is calculated from the ratio of the

electron density to the Saha density; first by calculating the modulus m, and

then by taking the arithmetic geometric mean to find K 15J.

m•1/(2SNs,/(an2 ) — 1) (10)

In the Saha limit, the modulus approaches an asymptotic value of unity and

round-of error in the a.g.m. calculation of K becomes appreciable. In this case

K is calculated from the following expression:

K•nd(a/emD,)1R , m:*1 (11)

This is derived by solving eqn. 10 for SN,,, substituting this expression into the

eigenvalue condition (eqn. 9) and solving for K.

TgANSPORT EQUATION

The second governing equation is a momentum balance for the electrons. It can

be written as a generalized Ohm's law and is referred to in this work as the

transport equation. The transport equation is obtained by solving the electron

flux (eqn. 2) for the electric field and then integrating this over the length of the

Inter-electrode distance to obtain the voltage drop across the plasma.

r, - -;A.,( n d`Y/dx + (kT,/e) do/dx) (12)

dT/dx . - r,/μ, - (kT,/e) (1/n) drVdx (131

WO - W, • - j r, (1/μ,n)dx - (kT,/e) In (nd/no)

(14)
5

vo - IFS . J R - (kTe/e) In (nd/no) (15)

where no and na are the electron densities at x • 0 and x . d.

The first term in the preceding expression is the voltage drop across the plasma

due to ohmic heating. The second term represents the voltage drop due to

concentration gradients (diffusion drop). To calculate the resistance in the

plasma R, the first term in eqn. 14 must be integrated over the length of interelectrode space. If the resistance is to be valid over a broad range of conditions

then the electron mobility that appears in the integral must include the effects of

both electron-atom collisions and election-ion collisions.

In evaluating this integral it is convenient to make use of an approximate

superposition rule for mobilities:

1111e • 1/Pea + 1/I*ef (16)

where μace and psi are the mobilities resulting from electron-atom

collisions and from electron-ion collisions.

Thus: R • Rea + Rei

R• j 1/(e1Len)dx . j 1/(eμean)dx +- j 1/(eμein)dx (17)

Since conductivity due to electron-Ion effects (Qei • eμein) is a function of T. and

Independent of the electron density it may be approximated as a constant and

removed from the integration. For the conditions of interest to us the resistance

from electron ion effects may be calculated from the following expression [6].

Rei • d/Qei

where Qei - .15085 Te13 M(A) (18)

and A • 12389 T.1 a(n/2)-1R

The electron-atom mobility is also a function of temperature and independent of

the electron density so it too may be taken out of the integration. It is:

pea - v 1Le/3KT (19)

In the above expression the mean velocity of the atoms is taken from kinetic

theory as v.4(8KT/x;n) and the moon free path is based on cross-sections

found in the literature. This gives the following expression for the resistance:

7
R w i /μ., f (1 /n) dx + d/a.i

(20)

If the appropriate substitutions are made the last remaining integral can be

evaluated and the following expression is obtained.

R • d/(2A..n) In 102 / nond(1-m)] + d/Q,i

(21)

The preceding discussion shows that the electron-atom resistance is inversely

proportional to the electron density while the electron-ion resistance is

practically independent of it. This being the case, at high electron densities

resistance depends only on temperature. The general behavior of resistance

versus electron density is plotted in figure 3. This behavior accounts for some

Important trends that appear in the results.

Often times the voltage across the electrodes is a more meaningful parameter

than the voltage drop across the plasma. For this reason it ic more useful to

express the transport equation in terms of the former quantity. This requires that

the voltage drops across the electrode sheaths be added to the voltage drop

across the plasma. h is assumed that positive (electron retaining) sheaths are

formed at the electrodes and that the flux at the sheaths is equal to the random

thermal flux times a Boltzmann factor.

coll: J = endv 

0/4 a (♦va/kT•) Va s kT./e In(endve/4J) (22)

emit: J - JE - endve/4 8 (♦Va'kT.) Vk ATO/e In [endv0/4(JE-J)] (23)

When the sheath voltages are added to the expression for the plasma drop

derive' above several of the terms combine and the result is the transport

equation in its final form:

Vd - JR + (kT,/e) In (J/ PEA] (24)

The first term still represents the ohmic drop while the second term now

incorporates the combir,id effects of the sheath and diffusion drops.

ENERGY EQUATION

The third necessary equation is an energy balance for the system. In its final

form this equation states that the electrical power loss in the plasma is caused

by three separate processes. Energy is lost as thermionically emitted electrons

deposit onto the collector, as electrons and ions that are created in the plasma

diffuse to the walls and las.ly as radiation out:

Vd J a 2KJE(T; TEYe + Pde + Prd (25)

_ _ 7-s,..
7

The fate of thermionically emitted electrons is illustrated in fig 4. These enter

the control volume at the emitter temperature and at a rate corresponding to the

saturation current given by the Richardson-Dushmann equation. Neglecting

back-emission from the collector, the operating current is equal to the emitted

current minus the backscattered current. Each of these currents has an energy

flux associated with it consisting of potential and kinetic energy terms. The net

sum of these fluxes is the first term in the energy equation. This term is only a

function of the electron temperature and is referred to as the kinetic energy term.

The sum of diffusion and radiation losses is the power required to maintain the

plasma ignited. This is calculated from the net rate of ionization established by

the the ambipolar diffusion equation and can tie expressed as the sum of the

terms Pdrtt and Prad.

PdM - 4 Eo K D. Ne/d where Eo is the ionization energy (26)

Pmd = E6p N6p gA (27)

As one might expect the diffusion losses are inversely proportional to the interelectrode distance. The power lost by radiation is equal to the product of the

number of atoms in the first excited level, the frequency at which they emit and

the energy of this radiation:

The escape factor g, is used to modify the Einstein A-coefficient so that (for

example) if the gas is optically dense, g is close to zero and very little energy is

lost in the form of radiation. To calculate the total number of atoms in the 6p

state the expression for the local 6p population (eqn. 12) is integrated over the

control volume. This integral cannot be eva!uated analytically and so it is

necessary to perform the integration numerically. In contrast to diffusion losses

radiation losses are proportional to the inter-electrode distance. Also it should

be understood that Prad represents radiation emitted in all directions not just the

radiation which is available for collection by photo-panels.

SOLUTION PROCEDURE

In summary of the proceeding sectiona, the model reduces to three nonlinear,

algebraic equations in three unknowns. These are:

The ambipolar diffusion equation.

Da den/dx2 = S N., n - acn3 (7)

The transport equation.

Vd s JR + (kT./e) in (J/ (JE-41 (24)

And the energy equation.
8

Vd = (J1JE)2K(T; TE)/e + Pd;n/JE + PmWJE (25)

A solution of model requires that these be solved simultaneously. Although any

method of solution will do, the method outlined below is the one that yields the

most insight into the problem. This method takes advantage of the fact that the

ambipolar diffusion equation is a function of only two of the three variables (T,

and n) and as such can be solved first for the ignited condition which the other

two equations must satisfy (fig. 5). This reduces the problem to that of solving

two equations in two unknowns. These can be solved simultaneously, however

it is useful to solve them independently for the variable J and then to search for

solution points (fig 6). This allows one to better understand why solutions exist

where they do.

RESULTS

In the study the escape factor, the emitter temperature, and the inter-electrode

distance were varied as parameters while other conditions were fixed at the

following values.

Ion Temperature =1500 K

Emitter work function = 2.4 V

Plasma drop - 1 V

Electrode area (planar) - 1 cm2

Number density of plasma - 1016/Cm3

Briefly, the rationale that went into this selection was the following. 1) The

escape factor is necessarily a parameter since the optical density is

independent of the macroscopic quantities calculated in the model. For this

reason it is also true that some of the curves, especially those for higher values

of the escape factor, may not correspond to physically achievable solutions. 2)

Varying the emitter temperature is equivalent to varying to power of the device.

At the onset it was not clear what the relationship between power and efficiency

would be. 3) Finally, the inter-electrode distance has a direct effect on the area

available for radiation. Although distance does not appear explicitly in the

plotted insults it is perhaps the most important variable in determining the

radiation output. One important result is that some optimal solutions involve

Inter-electrode spacings which are significantly greater than normally

encountered.

Regarding the other conditions, the only one that warrants any mention is the

plasma drop of one volt. It is this characteristic that sets the proposed device

apart from conventional thermiornc generators and is largely responsible for the

higher radiaton eficiencies . This is discussed further in the next section.

The results of the study are summarized in figures 7,8, and 9. Three degrees of

optical density were investigated corresponding to escape factor values of g=1,

9-A, g=.01. For all three cases the emitter temperature was varied from 1500K
9

to 2000K. And for each combination of g and TE the inter-electrode distance

was varied from the minimum value for which a solution existed to greatest

possible value for which a solution existed. Since there is no simple criteria for

determining beforehand where solutions exist the method used was in effect

trial and error, letting the computer determine where solutions were

impermissible. An important part of understanding the results Involves

understanding why under certain conditions solutions fail to exist.

in each case studied the three governing equations were solved for the electron

temperature, electron density and operating current by the method outlined

above. From this information the radiation output and the efficiencies were

calculated. Radiation output is defined as Prad•and the overall efficiency is

defined as the ratio of the radiation output to the total power consumed by the

device.

tI s Pr APrad + Pdin + J (#A + 2kTa /e )] (27)

where *A is the work functon of the anode.

The volumetric efficiency was also calculated. This is a figure of merit defined

as the ratio of the radiation output to the power dissipated inside the plasma.

t1vW - Prae/JVd

(28)

The data points that are plotted correspond to the solutions for which the

efficiency was the greatest. In each case these coincide with the solutions for

which the interelectrode spacing was also the greatest. That this is hardly a

coincidence is taken up in the following section.

DISCUSSION

Two results of this study are of particular interest. These are that radiation

efficiencies on the order of thirty percent are attainable and that this maximum is

independent of the optical density of the gas. High radiation efficiencies result

from the greater plasma drop that characterizes the proposed device. A greater

plasma drop means that more energy is dissipated inside the plasma.

additionally it happens that a larger fraction of this energy is dissipated as

radiation. That this should be the case 'is not obvious and is explained in the

following section.

One important consequence of a higher plasma drop is that solutions exist at

greater inter-electrode spacings. This is because larger plasma -1 ro►;;; can

overcome a greater resistance. Ultimately it is the larger inter-e,'Ijctrode

spacings that account for the higher radiation efficiencies. As mentioned

earlier energy is lost from the plasma in three seperate processes. These are

Illustrated in figure 10 which shows the contributions of the various terms to the

energy equation. The kinetic energy term describes the energy that is lost as
10

emitted electrons are deposited onto the collector. Pd"f represents diffusion

losses to the walls and Prad is the energy lost as radiation. Since the electron

temperature is determined by the ignited condition, the kinetic-energy term is

independent of the plasma drop. Radiation and the diffusion losses on the

other hand depend strongly on the inter-electrode spacing. At short distances

diffusion losses are very high. At large distances radiation losses dominate.

From this discussion one can understand why at higher voltage drops a greater

percentage of the plasma losses will occur as radiation and consequently

radiation efficiencies will be greater.

This does not yet explain why the maximum radiation efficiency is independent

of the optical thickness of the gas, for Intuitively, one would expect optically

dense gasses to be less efficient at producing radiation. The explanation lies

two parts. First, optically dense solutions exist at greater distances then

optically thin solutions. This can be understood by referring once again to figure

10 and noting that since Prad is proportional to the escape factor g, a smaller

slope on the Prad line will cause the energy equation to intersect the tranport

equation at a greater distance. The consequence of this is that net radiation

output is not proportional to the escape factor. In some cases studied, reducing

the escape factor by a factor of 100 had the effect oT reducing the radiation

output by only a factor of 5. In addition to this there is the fact that optically

dense gasses produce their maximum radiation at lower value of the emitted

current (fig 9) and hence at lower input powers.

It is interesting to compare the proposed device with a Camot engine operating

between the two electrode temperatures. ff we take the collector temperature to

be such that the emitted current from the collector is one percent of the emitted

current from the emitter, then the device is operating between 1500K and 745K.

The Camot efficiency between these two temperatures is 50.3% and so our

device is operating at 56% of Camot efficiency.

The results then indicate that high radiation efficiencies can be obtained from

cesium thermionic devices provided that there is an adequate votage drop

across the plasma. The value of one volt that was used in the study is a

reasonable value in that it can be generated without any electrical power input

to the device. High efficiencies come by virtue of the higher inter-electrode

distances which the solutions will accomodate. One particularly positive result

is that efficiencies on the order of 30% can be obtained even with optically

dense gasses. Since this is most likely to be the case under conditions of

interest, this study lends some theoretical support to the proposed device.
11

REFERENCES

1. D.L. Chubb.'Therrnionic•Photovoltaic Energy Converter", 1983 IEEE

International Conference on Plasma Science. San Diego. .

2. J.L. Lawless. "An Analytical Model of Thermionic Discharges". To be

published.

3.- Ibid.

4. M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions,

Dover: New York, 1972. Ch 16.

5. Ibid. section 16.4

6. L.Sptizer, Jr. and R. HSrm, Physical. Review, Vol. 89, 977 (1953)

7. J.F. Nolan and A.V. Phelps. Physical. Review, Vol. 145,279 (1966)

8. D.W. Norcross and P. M. Stone, "Recombination Radiative Energy Loss and

Level Populations in Nonequilibrium Cesium Discharges" J. Quant.

Spectrosc. and Radiat. Transfer, Vol. 8, 655 (1968)
0.6

x 0.6

<

0.4

0.2

.r

^^ F

t

FIGURES

ELECTRON DENSITY VS. DISTANCE

0.8

x 0.6

0.4

0.2

0

0 0.2 0.4 _ 0.6 0.8 1

x/d .

Figure 1. Electron distribution under non-Saha conditions.

ELECTRON DENSITY VS. DISTANCE

1,

12

0

0 0.2 0.4 0.6 0.6 1 x/d'

Figure 2. Electron distribution under Saha conditions.
0.05

N 0.04

E

i

C 0.03

WU

Q 0.02

W- 0.01

I

I

1

1S

,old'

Legend

R

REA

REI

13

RESISTANCE VS. ELECTRON DENSITY

0.00

0 S 10

ELECTRON DENSITY cm's

Figure I Plasma resistance vs. electron density.

T

• ^^ + 2kTE/*) Vd

^t1A + 2kT,/e)

•t1^ + 2kT./e)

#A

Figure 4. Enemy fluxes for emitted electrons.
3200

Y

W

3000

2500

F

L

-r

14

TEMPERATURE VS. ELECTRON DENSITY

3100

2`00

0 s 10 ts 20

ELECTRON DENSITY cm -3

std,

Figure 5. Solution to the ambipolar dNusion equation.

CURRENT VS. ELECTRON DENSITY

100

soGo40.

U

20

/ lop*nd

s^^srotr ten

0 _"1

0 5 to is 20

ELECTRON DENSITY cm -3

•113

Figure 6. Solutions to the transport and energy equations.

•.r
15

EFFICIENCY VS. EMITTER TEMPERATURE

30

U 20

10

0

1500

so

40

\

lk. \

1600 floc 1800 1900 2000

EMITTER TEMPERATURE K '

Legend

0=1

rpz.ol

Figure 7. Efficiency vs. emitter temperature.

VOLUMETRIC EFFICIENCY VS. EMITTER TEMPERATURE

Legend

a_'^ _.

0=.0^

1500 Roo /700 1800 »00 2000

EMITTER TEMPERATURE K

Figure S. Volumetric efficiency vs. emitter temperature.

100

so

F

U so

U

40

Q 20

0
16

RADIATION DENSITY VS. EMITTER TEMPERATURE

Legend

(I--1

C=.1

Gs.01_._.

0 L= ---

---

t;00

--- --- __ _

ts00 1600 fl00 1800 1900 2000

EMITTER TEMPERATURE j

Figure 9. Radiation output vs. emitter temperature.

POWER LOSS VS. DISTANCE FOR nxED ELECTRON DENSITY

Energy Equollal

KE

Pdill

0

0 0.2 0.4 0.6 0.6 1

DISTANCE cm

Figure 10. Energy losses as a function of inter-electrode spadng-
ELECTRODE DISTANCE VS. EMITTER TEMPERATURE

0.6

C \

Z 0.4-

Q \.

W0 `\

0.2

W \-

Itzz"`

0.0

1500 1600 n00 1800 1200 2000

EMITTER TEMPERATURE K

Legend

G-- 1

c_.1`..

GZ.01

Figure 11. Maximum Inter-electrode spacing vs. emitter temperature.

J
APPENDIX A: REACTION RATE CONSTANTS FOR EXCITED STATE

KINETICS

To derive expressions for the ionizaton and recombination rate coefficients S

and a, it is necessary to know the atomic states that are present in the plasma

and the reactions that occur between these states. In this work three states and

three reactions (shown below) were considered important in determining the

ionization/recombination kinetics of the plazi ia.

CS6a + 6- b CS6p + 0-

CS6p CS6s + by

CS6p + e- a Cs* + 2eBased on these reactions the rate equations for the atomic states are written as

follows:

dN6Vdt = -N6snK6a-6p + N6pnV*,6t + N6pgA

dN6p/dt = -nN6pK" + nN6.K6.4p - N6pgA - nN6pK6p.c + n3Kc4p

where K6s .6p Is the foreward rate constant for the reaction between

the ground state and the 6p state, and K6p-s. the backwards rate

constant. The subscript c refers to the continuum.

Under quasi-steady conditions the populations of the excited states change

slowly in time. With this assumption the rate equation for the 6P state reduces

• to an algebraic expression that can be solved for the 6P population:

N6p = R, +R, N66

where: Ro = n3K04jP / (nKep4k + nK6p c + gA)

R, =nN64K4.4. / [nK4p-a. + nK6pt + gA]

Substituting the expression for N; p into the rate equation above tho,i yields the

desired form from which the coefficients S and scan be factored:

dN6./dt = —nSN6a + an3

I - — --- _
where: S- Ki..Gp-R1(K6p.0+gA/n)

a-RO(nKs,,.0+gA) /n3

To evaluate these expressions it Is necessary to know the toreward and

backward rate constants for the reactions between adjacent states.

K6,.6p Is the rate constant for collisional excitation from the ground state to first

excited level. h is found by integrating the product of the cross-section for the

reaction [ref. 7] and the velocity distribution of the electrons (assumed

Maxwellian) over the range 0 energies. Assuming that the cross-section

Increases linearly from the threshold the result is:

K6i-ep - exp (-CV02) [B (VO4 + V02/C + 1JC2 ) + A (V02 + 1/C )]/2C

where A - (m.120T.)2,24na0

B - (m,/2nkT.)"2nm,Aa/AE

C - m,/2kT9

and VO - 42m./E, E -1.33 eV

(YO-1.02X10-14 Cm2

AG/AE - 71 X10-16 cm2/eV

From this expression it is possible to obtain an expression for the backward

reaction, KWU. The quasi-equilibrium assumption states that:

P7lrwK6..ap-N6pV,Gp_"_

If equilibrium with the continuum is also assumed then the ratio of the

populations can be gotten from the Saha equation.

N6s/N6p -(1/3) exp(E6..G,*T,,)

This gives the following expression for K4,.U:

Kgp.0 - KN.6p (1/3) uxp(E ft.%ATs).

Kc.4 is an effective rate constant for collisional de-exdtation from tho continuum

to the 6P state. This Is given In reference 8 as:
Kc•6p .46 (e"It47LtOK r `) j" (okT^ ^^^T^^)tr1

Invoking once again quasi-equilibrium and the Saha equation yields:

Ksp.c • Kc-sp (202/mkTO) exp (Er AT•)/(2Lsp + 1)

rn
APPENDIX B: BOUNDARY CONDITIONS

In this work two different sets of boundary conditions are used for the electron

de g ity at the electrodes. In the lowest order approximation, the walls are

considored absorbing so that:

no=0 nd =0

where the subsctipts 0 and d are used to denote the emitter are collector

sheaths.

This assumption works well in finding the solution to the ambipolar diffusion

equation. h is not adequate however for expressions in which terms on the

order of lVd are important. This is the case when evaluating the resistance of

the plasma. For here the assumption of zero electron density at the walls will

lead to an infinite resistance.

Appropriate nonzero boundary conditions are derived in reference 2 by

developing expressions for electron and ion fluxes at the walls. These are the

following:

no = 0/vi 4TJO (D„ dn/dx + (μ,/μ,)I',)

nd = 0/vi 4T/To ( -D, dn/dx - (μ;/μ,)I',)

Where: A is the unit vector normal to the electrodes. On the cainode surface 6 is

directed into the plasma. On the anode it is directed out of the plasma. I'. is the

electron flux emitted thermionically from the cathode, given by the RichardsonDushmann equation
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCI'CCCCCCCCCCCCCCCCCCC

C C

C PROGRAM PL1 C

C

Wk1TTEN BY F. STEFANI, 1964.

C

C

C C

C THIS PROGRAM CALCULATES THE OPERATING CONDITIONS FOR A THERMIONIC C

C DEVICE BASED ON AN ANALYTICAL MODEL DEVELOPED BY J.L. LAWLESS. C

C SOLUTIONS ARE FOUND BY FIRST CALCULATING THE TEMPERATURE THAT C

C SATISFIES THE IGNITED CONDITION AT A GIVEN ELECTRON DENSITY, AND C

C THEN SOLVING THE TRANSPORT AND ENERGY EQUATIONS FOR THE OPERATING C

C CURRENT BASED ON THAT TEMPERATURE. THIS IS DONE FOR A RANGE OF C

C ELECTRON DENSITIES AND THE POINT WHERE THE THREE EQUATIONS ARE C

C SATISFIED SIMULTANEOUSLY CAN BE FOUND BY INSPECTION OF THE OUTPUT. C

C ''THE PROGRAM ALSO CALCULATES THE RADIATION EFICIENCY OF THE DEVICE. C

C C

C [REF: AN ANALYTICAL MODEL OF THERMIONIC DISCHARGES, J.L. LAWLESS, C

C CARNEGIE-MELLON UNIVERSITY] C

C C

C SOME ADDITIONAL COMMENTS: C

C C

C * THE FOLLOWING QUANTITIES ARE ASSIGNED FIXED VALUES IN THE PROGRAM: C

C PLASMA VOLTAGE DROP VD = 1.0 [V] C

C ION TEMPERATURE TI = 1500. [K] C

C # OF NUCLEI NNUC = 1.OE16 [CM**-3] C

C ELECTRODE AREA AELEC = 1.0 [CM-*2] C

C C

C * ANODE WORK FUNCTION: WFA = WFC - VD (WFC; CATHODE WK. FTN.) C

C C

C * RECOMBINATION IS ASSUMED THROUGHOUT. C

C C

C * N6P IS COMPUTED BY NUMERICALLY INTEGRATING THE LOCAL 6P POPULATION C

C OVER THE INTER-ELECTRODE SPACE. C

C C

C * RADIATION EFFICIENCY IS CALCULATED IN TWO WAYS: C

C EFFI=PRAD/(PRAD+PDIFRJ(VA+2KTE/CHG)) C

C EFF2=PRAD/(RJ-VD) C

C C

C * SIMPSONS RULE HiLS BEEN MODIFIED TO TAKE ADVANTAGE OF ,SYMMETRY. C

C C

C * OUTPUT IS WRITTEN TO FILE PLI.DAT C

C C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C LIST OF VARIABLES

C

C INPUTS:

C

C DIST INTER-ELECTRODE GAP [CM]

C TEMIT EMITTER TEMPERATURE [K]

C WFC EMITTER WORK FUNCTION [V]

C G RADIATION ESCAPE FACTOR

C

C OUTPUTS:

C

C NE ELECTRON DENSITY [CM**-3]

C XN1 ELECTRON DENSITY AT THE ANODE [CM**-3]

C TE ELECTRON TEMPERATURE [K]

C SMI 1 - MODULUS

C XJEN OPERATING CURRENT BASED ON ENERGY EQN. [A/CM-*2]

r
C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

RJ OPERATING CURENT BASED ON TRANSPORT EQN. [A/CM**2]

R PLASMA RESISTANCE [OHMS]

REA RESISTANCE DUE TO ELECTRON--ATOM COLLISIONS [OHMS]

REI RESISTANCE DUE TO ELECTRON-ION COLLISIONS [OHMS]

VK CATHODE SHEATH VOLTAGE [V]

VA ANODE SHEATH VOLTAGE IV]

EFF1 RADIATION EFFICIENCY (SEE HEADER)

EFF2 VOLUMETRIC EFFICIENCY "

SAHA FLAG FOR SAHA CONDITION:

SARA=O => PLASMA IS RECOMBINING

SAHA=1 => RECOMBINATION NEGLIGIBLE

SARA=2 => SM1 = 0

PHYSICAL CONSTANTS:

ME ELECTRON MASS [KG/1.602E-15]

MI ION MASS [KG/1.602E-15]

CHG FUNDAMENTAL CHA.IGE [C/1.602E-19]

K BOTZMANNS CONST [eV/K]

E6P ENERGY Cs6p-Cs6s [eV]

ESP ENERGY Cs6s-Cs6p [eV]

ACOEFF EINSTEIN A COEFFICIENT [1/S]

LS DEGENERACY OF THE GROUND STATE

LP DEGENERACY OF THE FIRST EXCITED LEVEL

SIGIA CROSS-SECTION FOR ION-ATOM COLLISIONS [CM**2]

SIGEA CROSS-SECTION FOR ELECTRON-ATOM COLLISIONS [CM**2]

SIGO Y-INTERCEPT OF CROSS-SECTION VS. ENERGY CRV. [CM**2]

DLTXCT SLOPE OF CROSS-SECTION VS. ENERGY CURVE. [CM**2/EV]

EXTERNAL F

EXTERNAL SIMPER

EXTERNAL TMIN

EXTERNAL XJT

EXTERNAL TDIFF

REAL KSP,KPS,KSPA,KSPB,isCP,KPC,NGS,NE,N6P,MUI,MUE,MFPI,MFPE

REAL MI,KSPCON,ME,H2AF,PI,LS,LP,K,NNUC,JE,JNET,QP,MOD,JION

DIMENSION DDIST(15),DTEMIT(15),DATWFC(15),DATG(15)

COMMON/C1/ME,MI,CHG,EPSO,K,HBAR,E6P,ESP,ACOEFF,PI,LS,

*LP,SIGIA,SIGEA,DLTXCT,VD,JE,TI,TEMIT,AELEC,DIST,G,NE,

*SIGSP,SIGO,NNUC

COMMON/C3/KPS,KPC,KSP,KCP,NGS

COMMON/C4/MOD,QP,SM1

COMMON/C2/TE,DA,MUE

COMMON/C5/REA,REI,R

COMMON/C6/SAHA

OPEN(UNIT=7,DEVICE='DSK',FILE='PL1')

C

C PHYSICAL CONSTANTS

C

ME=5.686E-16

MI=1.377E-10

CHG=1.0

EPSO=5.527E5

K=8.6176E-5

HBAR=6.5764E-16

E6P=2.4625

ESP=1.4319

EO=E6P+ESP

ACOEFF=3.57E7

PI=3.14159265
LS=0.0

LP=1.0

SIGIA=1.2E-13

SIGEA=4.E-14

SIGSP=50.E-16

SIGO=-1.01E-14

C •

C SLOPE OF CROSS-SECTION FOR 6S-6P TRANSITION IN CESIUM (DLTXCT)

C OBTAINED FROM NORCROSS & STONE: PP. 667, FIGURE 1, CURVE #4; REF. 29;

C NOLAN AND PHELPS, 1965.

C

DLTXCT=71E-16

C

C OTHER INPUTS

C

VD=1.0

T?=1500.

AELEC=1.

NNUC=1.E16

C

C READ AND STORE RUN PARAMATERS

C

ICOUNT=O

100 ICOUNT=ICOUNT+1

WRITE(5,475)ICOUNT

475 FORMAT(2X,'INPUT DATA:DIST,TEMIT,WFC,G',12X,'CASE',2I)

READ(5,476)DDIST(ICOUNT),DTEMIT(ICOUNT),L)ATWFC(ICOUNT),

$DATG(ICOUNT)

476 FORMAT(4G)

WRITE(50700)DDIST(ICOUNT),DTEMIT(ICOUNT),DATWFC(ICOUNT),

$DATG(ICCtP T)

700 FORMAT(2X,4(2X,E12.5))

IF(DTEMIT(ICOUNT).NE.0.0)GO TO 100

C

C MAIN PROGRAM DO LOOPP.

C

ICOUNT=ICOUNT-1

DO 10 NX=I,ICOUNT

DIST=DDIST(NX)

TEMIT=DTEMIT(NX)

WFC=DATWFC(NX)

G=DATG(NX)

JE=120.*TEMIT**2*EXP(-WFC/(K*TEMIT))

C

C WRIT; RELEVANT INPUT DATA AND HEADER TO OUTPUT FILE.

C

WRITE(7,4BO)NX

WRITE(7,481)DIST

WRITE(7,482)TEMIT

WRITE(7,483)G

WRITE(7,484)WFC

WRITE(7,485)JE

WRITE(7,486)VD

WRITE(7,487)

480 FORMAT(1H1,' INPUTS PL1

4B1 FORMAT(2X,'DIST ='F5.3)

482 FORMAT(2X,'TEMIT='FB.2)

483 FORMAT(2X,'G ='F5.3)

484 FORMAT(2X,'WFC ='F5.2)

485 FORMAT(2X,'JE -'F6.2)

CASE',I2,/)
486 FORMAT(2X,'VD ='F5.2)

487 FORMAT(' 1 ,2X,2X,'NE',BX,' XN1'

*,9X,'TE',12X,'SM1 ',

*BX,IXJENI,SX,IRJI,6X,IRI,BX,'RET',BX,IVKI.BX,'VAI,6X,

*'EFF',2X,'PRAD/PIN',1X,'SAHA')

C

C MAIN DO LOOP INDEXING NE

C OPTION 1: SCALING ORDERS OF MAGNITUDE

C

C DO 11 J=1,3

C DO 12 I=1,9

C NE=.1E12*FLOAT(I)*10.**FLOAT(J)

C+

C OPTION 2: FIXED INCREMENT

DO 11 I=1,30

NE=.05E15*FLOAT(I)

C

C ***SOLVE THE AMBIPOLAR DIFFUSION EQUATION***

C

C

C FIRST FIND 'TMIN'. WHICH IS THE SOLUTION TO THE AMBIPOLAR DIFFUSION

C EQUATION NEGLECTING RECOMBINATION.

C THIS TEMPERATURE IS THE LOWER BOUND IN THE ROOTFINDING

C FOR THE AMBIPOLAR DIFFUSION EQUATION WITH RECOMBINATION

C (WHICH IS EVALUATED BY THE FUNCTION TDIFF).

C

A=1500.

B=5500.

SAHA=O.

CALL BISECT(TMIN,A,B,.S,IFLAG)

IF(IFLAG.LT .2)GO TO 46

WRITE(7,333)IFLAG

333 FORMAT(2X,'ROOTFINDING FOR TMIN BOMBED',2X,2I)

GO TO 9

46 rONTINUE

:SAHA=(A+B)/2.

C

C SOLVE THE AMBIPOLAR DIFFUSION EQUATION AND CALCULATE THE FOLLOWING

C QUANTITIES IN THE PROCESS.

C

C

C If ROOTFINDING BOMBS THEN THE SOLUTION IS TOO CLOSE TO THE

C LOWER LIMIT CALCULATED BY TMIN. CONDITIONS ARE CALCULATED

C AT THE SAHA LIMIT. SAHA FLAG IS SET TO 1.

C

A=B

B=5500.

CALL BISECT(TDIFF,A,B,.01,IFLAG)

IF(IFLAG.LT .2)GO TO 47

C WRITF(7,334)IFLAG

C334 FORKI:T(2X,'ROOTFINDING FOR DIFF. EQN. BOMBED:

C enEFAULT TO SAHA VALUES',2X,2I)

SARA=1.

TE=TSAHA

GO TO 212

47 CONTINUE

TE=(A+B)/2.

ERROR-ABS(A-B)/2.

C

C CALCULATION OF THE RATE CONSTANTS.
C

212 CC=ME/(2.*K*TE)

DD=ESP/(K*TE)

VMSQ=2.*ESP/ME

C

C FORMULATION FOR KSP BASED ON VARIABLE CROSS-SECTION.

AA=(CC/PI)**(3./2.)*4.*PI*SIGO

KSPA=(AA/(2.*CC))*EXP(-CC*VMSQ)*(VMSQ+1./CC)

BB=(CC/PI)**1.5*2.*PI*(ME*1.E20)*DLTXCT

KSPB=BB/(2.*CC)*EXP(-CC*VMSQ)*((VMSQ/l.E20)*VMSQ+2.*VMSQ

*/(CC*1.E20)+2./(CC*l.E20*CC))

KSP=KSPA+KSPB

KPS=KSP*EXP(DD)*(2.*LS+1.)/(2.*LP+1.)

KCP=.46*EXP(.5*(ALOG(B.*TE/PI)+ALOG(K)-ALOG(ME))-5.*(ALOG(4.*PI*TE

*)+ALOG(EPSO*K)))

KPC=KCP/(( 2.*LP+1.)*(2.*PI*(HBAR/ME)*(HBAR/K)/TE)**(3./2.)

**EXP(E6P/(K*TE)))

R1=KSP/(KPC+KPS+G*ACOEFF/NE)

S=KSP-R1*(KPS+G*ACOEFF/NE)

RO=NE**2.*KCP/(KPS+KPC+G*ACOEFF/NF)

ALPHA=RO/NE*(NE*KPS+G*ACOEFF)/NE**2.

NGS=NNUC-NE

MFPI=1./(NNUC*SIGIA)

MFPE=1./(NNUC*SIGEA)

MUI=2.*CHG*MFPI/(3.*(K*TI*PI*MI/2.)**.5)

MUE=2.*CHG*MFPE/(3.*(K*TE*PI*ME/2.)**.5)

DA=MUI*K*(TE+TI)/CHG

IF(SAHA.EQ.0.0)GO TO 214

QP=SQRT(ALPHA/(B.*DA))*NE*DIST

IF(QP.LE.20.)GO TO 213

C

C IN THE CASE WHERE SM1 (SM1=1-MODULUS) IS TOO LOW TO BE CALCULATED BY

C THE A.G.M. LADDER (QP>20) 'THEN SM1 IS SET TO ZERO AND SAHA FLAG IS 2.

C

SARA=2.

SM1=0.

MOD=1.

N6P=(R0+R1*NGS)*AELEC*DIST

GO TO 215

C

C ***CALCUALTE OUTPUT CURRENT BY THE ENERGY EQN.***

C

C

C OUTPUT CURRENT CALCULATION (VARIABLE XJEN):

C N6P CALCULATED BY INTEGRATION OVER THE LENGTH.

C THE INTEGRAND, N6P(X), IS CALCULATED BY THE EXTERNAL FUNCTION F.

C

213 SM1=16.*EXP(-2.*QP)

MOD=1.-SM1

214 CALL SIMPER(F,AREA,DIST)

N6P=ARA

215 PDIFF=4.*EO*AELEC*DA*QP*NE/DIST*1.602E-19

PRAD=G*ACOEFF*ESP*N6P*1.602E-19

Q=PDIFF+PRAD

C

XJEN1=2.*JE *K *(TE-TEMIT)/(VD*CHG)

XJEN2=Q/VD

XJEN=XJENI+XJEN2

C

^. ***SOLVE THE TRANSPORT EQUATION-
C

C

C CALCULATE RJCRIT, THE UPPER LIMIT FOR ROOTFINDING, BASED ON THE CONDITION

C THAT THE ELECTRON DENSITY AT THE ANODE XNI, MUST BE POSITIVE.

C

RJCRIT=(SIGIA/SIGEA)*SQRT((TI/TE)*(MI/ME))

**(CHG*AELEC)*DA*NE*(PI/DIST)*1.602E-19

RJ=RJCRIT—RJCRIT*(.00001)

r IF(RJ.GE .JE)RJ=JE—JE*(.001)

C

r: C ROOTFINDING FOR OUTPUT CURRENT BASED ON THE TRANSPORT EQN.

C

C TRANSPORT EQUATION IS SOLVED BY EXTERNAL FUNCTION XJT.

C THE VARABLE USED FOR THE OUTPUT CURRENT IS RJ.

C

A=RJ/100.

B=RJ

CALL BISECT(XJT,A,B..00I,IFLAG)

IF(IFLAG.LT .2)GO TO 99

WRITE(7,717)NE,IFLAG

717 FORMAT(' ',X,E7.1,' ROOTFINDING XJT BOMBED',2X,'IFLAG=',2I)

GO TO 9

99 CONTINUE

s RJ=(A+B)/2.

ERROR=ABS(A—B)/2. —

`t C

C CALCULATION OF SHEATH VOLTAGES

C

VEBAR=SQRT(8.*K*TE/(PI*ME))

VIBAR=SQRT(8.*K*TI/(PI*MI))

A%FA=SQRTITE*2.*PI/TI)

XNA=(4. /(ALFA*VIBAR))*DA*2.*QP*NE/DIST

XNB=(SIGEA/SIGIA)*SQRT((TE/TI)*(ME/MI))*RJ/(1.602E-19*AELEC)

k*(4./(ALFA*VIBAR))

XNI=XNA—XNB

XNO=XNA+XNB

ARGVK=XNO*VEBAR*CHG*1.602E-19/(4.*(JE—RJ))

! VK=K*TE/CHG*ALOG(ARGVK)

ARGVA=XN1*VEBAR*C'iiG*1.602E-19/,4.*RJ)

VA=K*TE/CHG*ALOG(ARGVA)

C

C WFC IS THE CATHODE WORK FUNCTION

C WFA IS THE ANODE WORK FUNCTION.

C

I; WFA=WFC—VD

C

C CALCULATE SYSTEM EFFICIENCIES

C

EFFI=PRAD/(PRAD+RJ*(WFA+2.*K*TE/CHG))

EFF2=PRAD/(RJ*VD)

C

C ESTABLISH ION CURRENT

C

JION=PDIFF/EO

RATIO=JION/RJ

C

C WRITE CALCULATED QUANTITIES

C

WRITE(7,761)NE,XNI,TE,SMI,XJEN,RJ ,R,F:EI,VK,VA,EF'F1.,EFF2,SAHA

761 FORMAT(' ,X,E8.3,2X,E9 .3,3X,FIO.4,4X.E9.3,4X,F7.2,4X
*,F7.2,3X,F5.3,3X,F7.4,3X,F7.3,3X,F7.3,3X,F5.3,3X,F6.3,3X,F2.0)

9 CONTINUE

12 CONTINUE

11 CONTINUE

10 CONTINUE

WRITE(5,752)

752 FORMAT(2X,' OUTPUT WRITTEN TO FILE: PLI.DAT ')

CLOSE(UNIT=S)

STOP

END

FUNCTION TMIN(TE)

C

C CALCULATES IGNITED CONDITION BASED ON NO RECOMBINATION

C

REAL KSP,KPS,KSPA,KSPB,KCP,KPC,NGS,NE,MUI,MUE,MFPI,MFPE

REAL MI,ME,HBAR,PI,LS,LP,K,NNUC,JE,JNET,MOD

COMMON/Cl/ME,MI,CHG,EPSO,K,HBAR,E6P,ESP,ACOEFF,PI,LS,

*LP,SIGIA,SIGEA,DLTXCT,VD,JE,TI,TEMIT,AELEC,DIST,G,NE,

*SIGSP,SIGO,NNUC

CC=ME/(2.*K*TE)

DD=ESP/(K*TE)

VMSQ=2.*ESP/ME

AA=(CC/PI)**(3./2.)*4.*PI*SIGO

KSPA=(AA/(2.*CC))*EXP(-CC*VMSQ)*(VMSQ+1./CC)

BB=(CC/PI)**1.5*2.*PI*(ME*l.E'20)*DLTXCT

KSPB=BB/(2.*CC)*EXP(-CC*VMSQ)*((VMSQ/1.E20)*VMSQ+2.*VMSQ

*/(CC*1.E20)+2./(CC*l.E20*CC))

KSP=KSPA+KSPB

KPS=KSP*EXP(DD)*(2.*LS+l.)/(2.*LP+1.)

KCP=.46*EXP(.5*(ALOG(B.*TE/PI)+ALOG(K)-ALOG(ME))-5. n (ALOG(4.*PI*TE

*)+ALOG(EPSO*K)))

KPC=KCP/((2.*LP+1.)*(2.*PI*(HBAR/ME)*(HBAR/K)/TE)**(3./2.)

**EY..P(E6P/(K*TE)))

R1=2',SP/(KPC+KPS+G*ACOEFF/NE)

S=KSP-Iil*(KPS+G*ACOEFF/NE)

RO=NE*12.*KCP/(KPS+KPC+G*ACOEFF/NE)

ALPHA=RO/NE*(NE*KPS+G*ACOEFF)/NE**2.

NGS=NNUC-NE

TMIN=S *NGS-NE*ALPHA*NE

RETURN

END

4

FUNCTION TDIFF(TE)

C

C CALCULATES IGNITED CONDITION W/ RECOMBINATION ASSUMED.

C

REAL KSP,KPS,KSPA,KSPA,KCP,KPC,NGS,NE,MUI,MUE,MFPI,MFPE

REAL MI,ME,HBAR,PI,LS,LP,K,NNUC,JE,JNET,MOD

COMMON/C1/ME,MI,CHG,EPSO,K,HBAR,E6P,ESP,ACCEFF,PI,LS,

*LP,SIGIA,SIGEA,DL-TXCT,VD,JE,TI,TEMIT,AELEC,DIST,G,NE,

*SIGSP,SIGO,NNUC

COMMON/C4/MOD,QP,SM1

CC=ME/(2.*K*TE)

DD=ESP/(K*TE)

VMSQ=2.*ESP/ME

AA=(CC/PI)**(3./2.)*4.*PI*SIGO

KSPA=(AA/(2.*CC))*EXP(-CC*VMSQ)*(VMSQ+1./CC)

BB=(CC/PI)**1.5*2.*PI*(ME n 1.E20)*DLTXCT

KSPB=BB/(2.*CC)*EXP(-CC* VMSQ)*((VMSQ/l.E20)*VMSQ+2.*VMSQ
*/(CC*1.E20)+2./(CC*1.E20*CC))

KSP=KSPA+KSPB

KPS=KSP*EXP(DD)*(2.*LS+1.)/(2.*LP+1.)

KCP=.46*EXP(.5*(ALOG(B.*TE/PI)+ALOG(K)-ALOG(ME))-5.*(ALOG(4.*PI*TE

*)+ALOG(EPSO*K)))

KPC=KCP/((2.*LP+1.)*(2.*PI*(HBAR/ME)*(HBAR/K)/TE)**(3./2.)

**EXP(E6P/(K*TE)))

R1=KSP/(KPC+KPS+G*ACOEFF/NE)

S=KSP-R1*(KPS+G*ACOEFF/NE)

RO=NE**2.*KCP,/(KPS+KPC+G*ACOEF'F/NE)

ALPHA=RO/NE*(NE*KPS+G*ACOEFF)/NE**2.

NGS=NNUC-NE

MFPI=1./(NNUC*SIGIA)

MUI=2.*CHG*MFPI/(3.*(K*TI*PI*MI/2.)**.5)

DA=MUI*K*(TE+TI)/CHG

C

C A.G.M. CALCULATION OF QUARTER PERIOD; QP

C FOR MORE INFO SEE SECTION 16.4 IN THE HANDBOOK OF MATH. FUNCTIONS.

C MODULUS (YOD) BASED ON RATIO OF NE(MAX) TO N(SAHA).

C

MOD=1./k2.*S*NGS/(NE*ALPHA*NE)-l.)

SM1=1.-MOD

A0=1.

BO=SQRT(SM1)

CO=SQRT(MOD)

15 ANEW=(AO+BO)/2.

BNEW=SQRT(AO*BO)

CHEW=(AO-BO)/2.

ADIF=AO-ANEW

IF(ADIF.LT..00001)GO TO 16

AO=--ANEW

BO=BNEW

CO=CNEW

GO TO 15

16 CONTINUE

• QP=PI/(2.*ANEW)

X1=S*NGS

X2=4.*QP**2.*DA/DIST**2.

X3=.5*NE*ALPHA*NE

TDIFF=XI-X2-X3

28 RETURN

END

C

FUNCTION XJT(RJ)

C

C CALCULATES THE TRANSPORT EQUATION

C

REAL KSP,KPS,KCP,KPC,NGS,NE,MUI,MUE,MFPI,MFPE

REAL MI,ME,HBAR,PI,LS,LP,K,NNUC,JE,JNET,MOD

COMMON/C2/TE,DA,MUE

COMMON/Cl/ME,MI,CHG,EPSO,K,HBAR,E6P,ESP,ACOEFF,PI,LS,

*LP,SIGIA,SIGEA,DLTXCT,VD,JE,TI,TEMIT,AELEC,DIST,G,NE,

*SIGSP,SIGO,NNUC

COMMON/C4/MOD,QP,SM1

COMMON/C5/REA,REI,R

COMMON/C6/8AHA

VIBAR=SQRT(B. *K*TI/(PI*MI))

ALFA=SQRT(TE*2.*PI/TI)

XNA=(4. /(ALFA ► VIBAR))*DA*2.*QP*NE/DIST

XNB=(SIGEA/SIGIA)*SQRT((TE:/TI)*(ME/MI))*RJ/(1.602E-19*AELEC)
0*(4./(ALFA*VIBAR))

XN1=XNA-XNB

XNO=XNA+XNB

ARGI = NE**2./(XNO* XN1*4.)

TERM1 - 1./QP*ALOG(ARGI)

IF(SAHA . EQ.0.0 ) GO TO 100

TERM2=2

GO TO 101

100 TERM2=ALOG(16./SM1)/QP

101 REA= DIST/( 2.*NE*MUE * 1.602E-19 )*( TERMI +TERM2)

XLAM= 12389. * TE**1.5/SQRT ( NE/2.)

REI=DIST * 100/AELEC * ALOG ( XLAM)/(.015085*TE**1.5)

R=REA+REI

Cl-RJ*R

C2=RJ/( JE-RJ)

C3=ALOG(C2)

C4=C3*K*TE/CHG

XJT=(C1+C4)-VD

RETURN

END

C

SUBROUTINE BISECT(F,A,B,XTOL , IFLAG)

C

C NO-FRILLS BISECTION ROOT-FINDING SUBROUTINE

C

IFLAG=O

N--1

FA=F(A)

FB=F(B)

C WRITE ( 7,790)FA,FB

C790 FORMAT ( 2X,'ENDPOINTS:' , 3X,'FA=',E15.7, 7X,'FB=',E15.7)

IF (FA*F ( B).LE.O.) GO TO 5

I FLAG=2

C

C WRITE ( 7,601)A,B

C601 FORMAT ( 43H F ( X) IS OF SAME SIGN AT THE TWO ENDPOINTS

C 2F15.7)

RETURN

5 ERROR=ABS(B-A)

6 ERROR=ERROR/2.

C

C CHECK FOR SUFFICIENTLY SMALL INTERVAL

C

IF(ERROR.LE.XTOL) RETURN

XM=(A+B)/2.

C

C CHECK FOR UNREASONABLE ERROR REQUIREMENTS

C

IF(XM+ERROR.EQ . XM) GO TO 20

FM=F(XM)

C

C TEMPORARY PRINTOUT

C

N=N+1

C

C CHANGE TO NEW INTERVAL

C

IF(FA * FM.LE . O.) GO TO 9

A=XM

FA=FM

^..^- 
GO TO 6

9 B=XM

GO TO 6

20 IFLAG=1

MGTUmN

END

C

FUNCTION F(X)

C

C CALCULATES N6p(X) - RO + R1*NGS SO THAT IT CAN BE INTEGRATED BY

C THE SUBROUTINE SIMPER.

C

COMMON/C1/ME,MI,CHG,EPSO,K,HBAR,E6P,ESP,ACOEFF,PI,LS,

*LP,SIGIA; SIGEA, DLTXCT,VD,JE,TI,TEMIT,AELRC,DIST,G,NE,

*SIGSP,SIGO,NNUC

COMMON/C3/KPS,KPC,KSP,KCP,NGS

COMMON/C4/MOD,QP,SM1

DIMENSION A(20),B(20),C(20),PSI(20)

REAL KPS,KPC,KSP,KCP,NGS,NE,MOD

C

C CALCULATE AND STORE A.G.M. SCALE

C

A(1) =1.

B(1)=SQRT(SM1)

C(1)-SQRT(MOD)

N=2

20 A(N)=(A(N-1)+B(N-1))/2.

B(N)-SQRT(A(N-I)*B(N -1))

C(N)-(A(N-1)-B(N-1))/2.

ADIF-A(N-1)-A(N)

IF(ADIF.LT..0001)00 TO 26

N=N+1

GO TO 20

26 PSI(N)-2**(N-1)*A(N)*2.*QP*X/DIST

C

C RECURSION CALCULATION

C

DO 13 I=1,N

NN=N-I+l

TARG=C(NN)/A(NN)*SIN(PSI(NN))

IF(TARG.LT.1.0)GO TO 29

C WRITE(5,68)

C WRITE(5,69)C(N),A(N),PSI(NN),TARG

68 FORMAT(2X,'SUB:F(SAHA) BOMB: C(N),A(N),PSI(NN),TARG')

69 FORMAT(2X,E15 .9,4X,E15.9,4X,E15.9,4X,E15.9)

TARG=1.

29 PSI(NIA- 1)- (PSI(NN)+ASIN(TARG))/2.

13 CONTINUE

SN=SIN(PSI(1))

P-(SN**3.*NE*KCP*NE**2.+NERS,N-KSP*NGS)

*/(SN*NE*(KPS+KPC)+G*ACOEFF)

RETURN

END

C

SUBROUTINE SIMPFR(FX,AREA,DIST)

C

C SIMPSONS RULE: NSEG MUST BE EVEN, FX MUST BE SYMMETRICAL ABOUT MIDPOINT.

C

REAL LSEG

NSEG=20
LSEG=DIST/FLOAT(NSEG)

SUM=O

X3=0.

NSTEP=NSEG/4

DO 39 I-1,NSTEP

X1=X3

FX1=FX(X3)

X2=Xi+LSEG

X3=X2 +LSEG

SUM=SUM+FX1+4.•FX(X2)+FX(X3)

39 CONTINUE

AREA=2.•LSEG/3.•SUM

RETURN

END