Dipartimento di Elettronica, Politecnico di Torino, Torino, I-10129, Italy

Advanced Concepts of Electromagnetic Generation, Confinement and Acceleration of High Density Plasma for Propulsion

Xxxxxxxx Xxxxxx, Xxxx Xxxxxxxxxxx, Xxxxxxxx Xxxxxxxx

Xxxxxxx Xxxxxxx, Xxxxxx Xxxxx

CISAS X. Xxxxxxx, Padova, I-35131 Italy

Xxxxxxxx Xxxxxxxx

Advanced Concepts team Researcher, ESA-ESTEC, Noordwijk, The Netherlands

ARIADNA id: 05/3202

Contract Number: 4919/05/NL/HE Final report

Version 1.02

November 2007

1 Introduction 5

1.1 Scientific rationale . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2 Study objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Description of the work accomplished . . . . . . . . . . . . . . . 7

1.3.1 RF systems modelling and design . . . . . . . . . . . . . 7

1.3.2 Modelling of the plasma devices and system-level model

of the thruster . . . . . . . . . . . . . . . . . . . . . . . . 9

2 ICRH unit RF modelling with TOPICA 11

2.1 The case for ICRH unit simulation 11

2.2 TOPICA overview 12

2.3 Antenna problem formulation 13

2.3.1 Application of the Equivalence Principle 15

2.3.2 Statement of the equations 18

2.4 The plasma Green’s function 19

2.4.1 Reduced Xxxxxxx’x equations in the plasma 21

2.4.2 FEM solution 23

2.5 The plasma model 25

2.6 Solution by hybrid Moment Method 30

2.6.1 Algebraic system 30

2.6.2 Spectral reaction integrals 31

2.6.3 Antenna loading calculation 32

3 Plasma device modelling 36

3.1 Global model of plasma discharge 36

3.1.1 Plasma reactions 37

3.1.2 Gas dynamic model 39

3.1.3 Magnetic mirror 41

3.1.4 Conditions on the particles reflection 42

3.1.5 Plasma parameters at the ICRH 42

3.1.6 ICRH effect 44

3.1.7 Magnetic nozzle 44

3.2 Model validation 46

3.3 Model optimisation 47

4 Numerical results 48

4.1 RF modelling of the ICRH antenna 48

4.2 Plasma device 52

4.3 Scalability criteria and results 55

5 Conclusions 59

List of Figures

1.1 An example of a RF plasma thruster: layout of the VASIMR ex- perimental engine (cf. [2]). . . . . . . . . . . . . . . . . . . . . . 6

1.2 System model block structure: m˙ is the mass flow rate, n stand for

the plasma components and neutral densities and T for the plasma components temperatures. . . . . . . . . . . . . . . . . . . . . . 9

2.1 A CAD (electromagnetic) model of a sample ICRH unit including the thruster walls, a counter-driven two-loop antenna, and a cylin- drical plasma flow: also shown is the 3D surface triangular-facet

mesh for simulations with TOPICA. 13

2.2 Cartoon of a typical ICRH unit featuring two loop antennas sur- rounding a plasma flow. Also sketched are the surfaces whereon the EP is applied and unknown surface current densities defined: in that regard, ST needs neither to be defined nor meshed, as the

plasma Green’s function already includes its effect. 14

2.3 First application of the EP to the geometry of Fig. 2.2: (a) equiv- alent surface current densities are introduced on the surfaces SA+ and ST and (b) SA+ is closed by a PEC cylinder, therefore only the

magnetic current M A+ density contributes to the fields in VC,0 16

2.4 EP applied to the antenna region of the geometry shown Fig. 2.2:

(a) equivalent current densities are introduced on a surface SC,1 partly wrapping all conductors, the aperture SA, the feeding aper- tures SP,k and partly constituting a fictitious boundary; (b) all con- ductors and the plasma are removed so as to obtain a classical prob-

lem of EM field propagation in free space. 17

2.5 Qualitative representation of the triangular pulses defined in Eq. 2.31. 24

2.6 Plasma tensor entry εzz: (top) imaginary part as a function of the longitudinal wavenumber kz [1/m]; (bottom) enlarged view show- ing the unwanted occurrence of a null and consequent sign change

at approximately kz = 48 1/m. 28

2.7 Derivation of the coupling resistance: (a) cartoon of the adopted intermediate model and (b) equivalent (spectral domain) circuit of

the intermediate model for each (kz, m) pairs. 35

3.1 Schematic configuration of the plasma generation stage. A neutral gas is injected from the left, it gets ionized by the helicon antenna and it is exhausted from the right. The coils are necessary to gen-

erate the magnetic field to confine the plasma. 37

3.2 Qualitative magnetic field axial profile. 41

3.3 Maxwellian distribution function with a drift velocity parallel to

v . The particles outside the dash line are reflected by the effect of

the peak of the magnetic field. 43

3.4 Example of distribution function at the ICRH. A fraction of the distribution is cut by the presence of a peak in the magnetic field. 45

3.5 Comparison between experimental data [48] and the numerical model: helium discharge, 3 kW RF power 47

4.1 Real part of Z˜/Z0 entries as a function of m and kz/k0, with Z0

(k0) the free space impedance (wavenumber). 49

4.2 Imaginary part of Z˜/Z0 entries as a function of m and kz/k0, with

Z0 (k0) the free space impedance (wavenumber). 49

4.3 Standard counter-driven two-loop antenna: sample electric current magnitude distribution on conducting bodies and at plasma/air in- terface. 50

4.4 Standard counter-driven two-loop antenna: sample magnetic cur-

rent magnitude distribution at plasma/air interface. 50

4.5 Standard counter-driven two-loop antenna: plasma loading as a function of the frequency computed with TOPICA for different loop-width to plasma-radius ratios. Also superimposed are mea-

sured data and simulations published in the work by Xxxx [11]. 52

4.6 Electron temperature and density of different species during the plasma discharge. Gas: helium, absorbed power on the helicon stage: 4,000 W, injected mass flow rate: 2 10−6 kg/s. 53

4.7 Electron temperature and density of different species during the plasma discharge. Gas: argon, absorbed power on the helicon stage: 1,000 W, injected mass flow rate: 4 10−6 kg/s. 54

4.8 Thrust efficiency as function of specific impulse using helium and argon. 55

4.9 Thrust parameters given by the optimization using helium as pro- xxxxxxx (power repartition: 0 means full power on the ICRH, 1 means full power to the helicon). 56

4.10 Thrust parameters given by the optimization using argon as propel- lant (power repartition: 0 means full power on the ICRH, 1 means

full power to the helicon). 57

4.11 Thruster mass versus input power using argon and helium. 58

4.12 Specific mass versus input power using argon and helium. 58

Chapter 1

Introduction

1.1 Scientific rationale

Plasma-based propulsion systems make it possible to attain a very high specific impulse combined with continuous thrust. This implies that even though the thrust is orders of magnitude lower than for chemical thrusters, the continuous accelera- tion gained by a spacecraft propelled by such engines allows accomplishing very ambitious missions, whilst requiring relatively higher transfer times. On the other hand these engines demand for less propellant mass, thanks to the higher specific impulse they afford. Nevertheless, they need a larger power amount in contrast to chemical thrusters, which may be attained through solar arrays with increased surface area or even nuclear sources, should the power level be quite high.

In the past, several studies have been conducted trying to convert technologies expressly developed for fusion applications into propulsion systems. The most interesting ones are focused on the possibility of transferring energy to the plasma via electromagnetic waves at radio frequencies (1-50 MHz, RF), exploiting the possibility of having very efficient devices to generate and heat the plasma. These studies lead to very interesting features as:

•

the possibility of building variable specific impulse (Isp) and thrust at maxi- mum power, offering a great mission flexibility;

•

the possibility of building electrode-less thrusters, which completely avoid the problem of electrode erosion that normally is a significant limitation for high power electric thrusters;

•

high power density: this issue becomes fundamental in space application where dry mass is very expensive.

The structure of the reference system assumed in this study (as in the SOW) comprises three stages, where plasma is respectively generated, heated and ex- panded in a magnetic nozzle. The first stage handles the main injection of pro- xxxxxxx gas and the ionization subsystem. In it plasma is generated by a so-called

Figure 1.1: An example of a RF plasma thruster: layout of the VASIMR experi- mental engine (cf. [2]).

helicon antenna and confined by a suitable magnetic field. The second stage acts as an amplifier to further energize ("heat") the plasma; here plasma is heated by radio frequency (RF) waves in the regime of ion cyclotron (IC) resonance. The third stage, called magnetic nozzle, converts the thermal energy of the plasma into directed flow, while protecting the nozzle walls and insuring efficient plasma de- tachment from the magnetic field. For ease of reference, a sketch of a representative concept, the NASA VASIMR project, is shown in Fig. 1.1. An engine of this type has the potential of effecting exhaust modulation at constant power. It means that this system will have the remarkable capability of "shifting gear" during normal operations. The development of systems of this kind is an interesting technology challenge with very high potentials. The interest on this kind of system could be summarized with the following considerations:

•

They could lead to a flexible and adaptable technology, scalable to high power, for human missions, large robotic missions, thus precluding the need to develop separate propulsion systems for each purpose.

•

The possibility to complete Mars type mission accomplished with a sin- gle type of propulsion system, reducing therefore significantly the costs for quick planetary escape and low propellant consumption for interplanetary cruise, appear to be feasible if such a system could be developed.

•

Technology growth is open-ended and it could lead to potentially very high power systems in the future.

1.2 Study objectives

The general objective of this activity is to examine the challenges of a variable- specific-impulse plasma-thruster based on fusion-technology, and assess its fea- sibility. The obvious reference is the VASIMR project of the NASA mentioned in the SOW, but of interest here is also the investigation on the feasibility of a lower-power thruster. In full compliance with the SOW, the specific objectives of the study concentrated on the modelling of the above-referenced system, with the aim of allowing the optimization of its design, and hence, the assessment of the required power and efficiency levels for given performances. This in turn should allow the design of an experiment with an estimate of the resources needed for its development and testing.

1.3 Description of the work accomplished

In order to present the activity, we first briefly list the main components of the system as envisioned above, and the requirements of the analysis. Next, we will describe the specific tasks performed to reach the stated objectives, and the method- ology employed to carry out these tasks.

The modelling has been further divided into two different and linked activi- ties: RF system modelling, and modelling of the plasma device. The RF systems modelling and design activity provides "interface" models between the RF power generation and the power deposition into the plasma; the plasma model of the var- ious constituent regions (helicon, ICRH, nozzle) affords the system-level model of the engine. This allows assessing its performance as a function of the chosen geometry, magnetic field structure, RF frequencies, etc.

1.3.1 RF systems modelling and design

The main RF components of a plasma-based propulsion system (i.e. plasma thruster) are the helicon antenna for plasma generation, and the double-loop antenna for plasma acceleration. The helicon source creates a (cold) plasma by ionizing the injected gas by an RF-sustained continuous discharge. The plasma flows into the accelerating section where ion cyclotron frequency heating (ICRH) by means of electromagnetic waves is the main mechanism for power deposition in the plasma.

In the ICRH section, the RF frequency antenna frequency should match the ion cyclotron frequency to ensure wave energy conversion into ion gyro-motion. The ICRH has two distinct features: first, each ion passes the resonance only once, gaining an energy that is much greater than the initial energy; second, the ion mo- tion is collisionless, i.e. the energy gain is limited not by collisions but by the time the ion spends at the resonance while moving along the field lines. The key fea- tures of the single-pass ICRH is a flow of cold ions in an equilibrium axisymmetric mirror magnetic field in the presence of a circularly polarized wave rotating in the ion direction, and launched nearly parallel to the magnetic field lines.

A good antenna design for the ICRH section will excite primarily the forward- resulting mode (called m = -1); the wave-plasma coupling should be high enough to allow a reasonable ratio between transferred power and (averaged) stored energy; seen from the circuit point of view, the ratio between the resistive "loading" and the reactive impedance should be enough to allow the transfer of the requested RF power available at the RF generators, yet with the aid of a matching network. The power absorbed by the plasma for a given antenna current determines the plasma loading resistance, which is a very important parameter for the antenna design. In order to efficiently couple RF power, the plasma loading resistance must be sub- stantially larger that the (negligible) loading resistance attained in vacuo, which (at these frequencies) is caused only by finite resistance effects throughout the entire circuit driving the antenna.

As recognized by the SOW, heating is a critical part of the system, and its efficiency is crucial in meeting the overall requirements, especially if one envisions lower powers than investigated so far (e.g. in the VASIMIR experiments). The goal of the modelling task is to arrive at a tool that allows predicting the RF power transfer to the plasma for a given antenna configuration and plasma parameters. The goal of the numerical simulations is to model the underlying physics processes and to design an antenna to maximize loading resistance, or more generally, reduce the impedance correction needed to match the power source. The quality factor (Q) of the antenna+matching network resonant circuit is likely be a relevant indicator of the overall antenna performance, but a study will be required to assess the most relevant parameter to be optimized by antenna design.

An innovative tool has been developed and used for the 3D simulation of ICRH antennae, i.e. accounting for antennas in a realistic 3D geometry and with an ac- curate plasma model. The tool is based on the TOPICA code [1]. The approach to the problem is based on an integral-equation formulation for the self-consistent evaluation of the current distribution on the conductors. The environment has been subdivided in two coupled region: the plasma region and the vacuum region. The two problems are linked by means of a magnetic current (electric field) distribution on the aperture between the two regions. In the vacuum region all the calcula- tions are executed in the spatial domain while in the plasma region an extraction in the spectral domain of some integrals is employed that permits to significantly reduce the integration support and to obtain a high numerical efficiency leading to the practical possibility of using a large number of sub-domain (rectangular or triangular) basis functions on each solid conductor of the system. The plasma en- ters the formalism of the plasma region via a surface impedance matrix; for this reason any plasma model can be used. The source term directly models the TEM mode of the coax feeding the antenna and the current in the coax is determined self-consistently, giving the input impedance/admittance of the antenna itself.

Magnetic field

Gas pressure

m& , n ,T

m& , n,T

m& , v

exhaust

,Thrust

Helicon

ICRH

Helicon

out i e

Coupled Power

Magnetic nozzle

Coupled Power

Magnetic nozzle

ICRH

i⊥

Figure 1.2: System model block structure: m˙ is the mass flow rate, n stand for

the plasma components and neutral densities and T for the plasma components temperatures.

1.3.2 Modelling of the plasma devices and system-level model of the thruster

A global time-dependent model is necessary to assess the system performance and optimize the thrusters design. "Global" means that it incorporates all the three main stages, i.e. the helicon plasma source, the ICRH acceleration region and the mag- netic nozzle. Each stage requires a different physics-based modelling, and there- fore three sub-models were developed and then connected. The plasma parameters at the exit of the helicon source are the inputs of the ICRH acceleration region, and the plasma parameters at the exit of the ICRH stage are the inputs for the magnetic nozzle stage. The output of the magnetic nozzle stage gives the characteristics of the exhaust and thus it permits to evaluate the propulsive parameters of the device (see Fig.1.2).

All this makes it possible to simulate the plasma behavior in the desired config- uration and to determine the efficiency and thrust performance. The model follows the plasma parameters history as a function of the axial coordinate. This allows evaluating the antennae/plasma coupling, and thus the characteristics of the match- ing circuit and the power coupling efficiency. The model described takes into ac- count all the relevant physical parameters like dimensions of the device, magnetic field configuration. This permits estimating masses, thermal and structural loads.

Model of the helicon plasma source In this stage the neutral gas is ionised by electron-particles collisions excited by RF helicon waves. Input of the sub- model are the inlet gas pressure, and the power coupled by the helicon an- xxxxx with plasmas. Outputs of the model are: the propellant mass flow rate, the density of each neutral and ionized species, the electron density and temperature, and the time history of each particle/energy loss channel. The model considers that the coupled power is absorbed by electron-impact reac- tions and by the increment of electron temperature. The model will combine

Magnetic field

a global plasma source simulation with a 0-dimensional gas-dynamic simula- tion. It will account for changes in the neutral density, ionization, excitation and dissociation.

Plasma modelling in the ICRH region The RF power transfer derives from the electromagnetic model of the antenna; the power deposited by the ICRH is absorbed by the plasma in the form of normal kinetic energy. The inputs of this stage are the output parameters of the helicon stage and the power coupled by the ICRH antenna with plasmas. The outputs are the propellant mass flow rate, the plasma density and the normal ion temperature.

Modelling of the Magnetic Nozzle Downstream the ICRH, a diverging magnetic field converts part of the normal kinetic energy of the plasma into parallel kinetic energy, and thus produces thrust. This conversion is effective until the plasma detaches from the magnetic field lines. The inputs of this stage are the output parameters of the ICRH stage. The outputs are the propellant mass flow rate, the exhaust velocity and thus the thrust.

It is assumed that the detachment happens when the parallel kinetic energy density becomes greater than the magnetic energy density. Applying the conservation of energy, magnetic momentum, particles flow, and magnetic flux, it is possible to evaluate the axial velocity of the exhaust particles, and thus the thrust and thrust efficiency.

Chapter 2

ICRH unit RF modelling with TOPICA

2.1 The case for ICRH unit simulation

In the ICRH section, the antenna frequency should match the ion cyclotron fre- quency to ensure wave energy conversion into ion gyro-motion. The ICRH has two distinct features: first, each ion passes the resonance only once, gaining an en- ergy that is much greater than its initial one; second, the ion motion is collisionless, in consequence the energy gain is limited not by collisions but by the time the ion spends at the resonance while moving along the field lines.

The key features of the single-pass ICRH is a flow of cold ions in an equilib- rium axis-symmetric mirror magnetic field in the presence of a circularly polarized wave rotating in the ion direction, and launched nearly parallel to the magnetic field lines. A good antenna design for the ICRH section will excite primarily the forward-resulting mode (labelled m = -1); the wave-plasma coupling should be high enough to allow a reasonable ratio between transferred power and (averaged) stored energy. Seen from the circuit viewpoint, the ratio between the resistive loading and the reactive impedance should be enough to allow the transfer of the requested RF power available at the generators, albeit with the aid of a matching network.

The power absorbed by the plasma for a given antenna current determines the plasma loading resistance, which is a very important parameter to assess the an- xxxxx performance with regard to its capability of conveying power to the plasma. In order to efficiently couple RF power, the plasma loading resistance must be sub- stantially larger than the (negligible) loading resistance attained in vacuo, which at the operating frequencies commonly met for this application is mostly caused by finite resistance effects throughout the entire circuit driving the antenna.

Since heating is a critical part of the system, and its efficiency is crucial in meeting the overall requirements, especially if one envisions lower powers than investigated so far (e.g. in the VASIMIR experiments), then our main goal is to

arrive at a tool that allows predicting the RF power transfer to the plasma for a given antenna configuration and plasma parameters.

The purpose of the numerical simulations is to model the underlying physics processes and to design an antenna to maximize loading resistance, or more gen- erally, reduce the impedance correction needed to match the power source. The quality factor (Q) of the antenna plus matching network resonant circuit is likely to be a relevant indicator of the overall antenna performance, but a study will be required to assess the most relevant parameter to be optimized by antenna design.

2.2 TOPICA overview

An innovative tool, based on the TOPICA code [1], has been developed for the simulation of ICRH antennas, i.e. allowing for antennas in a realistic 3D geometry and with an accurate plasma model.

The problem of determining the antenna input admittance matrix [Y ], in pres- ence of a cylindrical plasma flow, has been tackled by an integral-equation for- mulation for the self-consistent evaluation of suitable unknown electric (J ) and magnetic (M ) surface current densities, whence [Y ] can be derived.

A widely-used configuration for the RF booster consists of a counter-driven two-loop antenna encircling the plasma column, but we want to emphasize that the approach developed allows considering any shape of the antenna—and of the acceleration unit as well. On the contrary we do require the magnetically confined plasma to possess circular symmetry.

Upon invoking the Equivalence Principle (EP) [3, 4, 5], the environment has been formally subdivided into two coupled regions: the plasma region and the an- xxxxx region. The two equivalent problems are then linked by means of a magnetic current (electric field) distribution on the air-plasma separation surface (dubbed aperture). Two coupled integral equations (IEs) to be solved for J and M ensue by enforcing the boundary and continuity conditions the tangential fields must ful- fill over the surface whereon the EP was applied. The plasma enters the formalism via its Green’s function (a rank-2 tensor surface admittance), therefore any plasma model can be used as a matter of fact.

The source term directly models the TEM mode of the coax feeding the an- xxxxx, meaning that the ultimate forcing term is actually the voltage supplied by the microwave amplifier. Once the current within the coax has been self-consistently determined, the input admittance of the antenna itself can be calculated.

Eventually, the IEs are solved by applying the Moment Method [6], both in spatial and spectral domain, the latter being the one wherein the plasma is far more easily described. In particular, the unknown surface current densities are repre- sented by a linear combination of a finite set of subdomain basis functions, which are defined over pairs of triangular facets (cf. Fig. 2.1). The coefficients of these linear combinations represent the actual unknowns of the problem, which turns from a set of integral equations into an algebraic system.

Figure 2.1: A CAD (electromagnetic) model of a sample ICRH unit including the thruster walls, a counter-driven two-loop antenna, and a cylindrical plasma flow: also shown is the 3D surface triangular-facet mesh for simulations with TOPICA.

The entries of the system matrix are computed either in the spatial or in the spectral domain. To be specific, in the antenna region the calculations are executed in the spatial domain, allowing a large number of subdomain basis functions on each conducting surface and the accurate modelling of small geometrical details. On the contrary, in the plasma region the spectral (Fourier) domain is employed, for with that representation of fields and sources the plasma Green’s function is more easily obtained.

2.3 Antenna problem formulation

The main purpose is to determine the admittance matrix [Y ] of an arbitrarily shaped antenna operating within the acceleration unit of a plasma thruster. Assuming NP is the total number of antenna ports, then [Y ] comes to be an NP -order matrix. To help formulate the problem, schematically depicted in Fig. 2.2 is a cross-sectional view of a typical ICRH stage basically comprised of a cylindrical plasma flow, a multiport antenna and the thruster wall.

The chief idea underlying the formulation—which can be regarded as an ex- tension to cylindrical plasmas of the approach adopted in [1]—is that the antenna interacts with the plasma through the electromagnetic (EM) unknown fields exist- ing over the aperture SA, as is apparent from the schematic drawing of Fig. 2.2. Considering the antenna and the plasma as belonging to conceptually distinct EM regions is numerically convenient, for the various regions differ as per the compu- tational difficulties they offer.

ICRH ANTENNA

SA-

zˆ

SA-

FICTITIOUS BOUNDARIES

THRUSTER WALLS

SA-

SA+

PLASMA FLOW v|| , B0

SA-

Figure 2.2: Cartoon of a typical ICRH unit featuring two loop antennas surround- ing a plasma flow. Also sketched are the surfaces whereon the EP is applied and unknown surface current densities defined: in that regard, ST needs neither to be defined nor meshed, as the plasma Green’s function already includes its effect.

In effect, the antenna region (surrounding the plasma beam in Fig. 2.1) ex- hibits a greater geometrical complexity, due to the presence of feeding ports, an- xxxxx, curved thruster walls, and as such demand a fully 3D description in terms of triangular facets. Conversely, the plasma region (the innermost one in Fig. 2.1) possesses a rather simple geometrical structure, but nonetheless it comprises the plasma, which is a magnetized, anisotropic and inhomogeneous medium wherein the Xxxxxxx’x equations are to be solved as well. Therefore, if somehow the two regions are self-consistently separated in the formulation from the beginning, then they can be suitably treated according to the specific challenges they pose. The mathematical tool behind the separation alluded at above is the Equivalence Prin- ciple [3, 4, 5].

The proposed solution is full-wave and relies on the following rather weak restrictions, viz.

•

All the metals (antenna, walls) are considered perfect electric conductors (PEC), since finite-conductivity effects are usually negligible in the ICRF heating and for the materials typically employed. On the other hand, a finite conductivity can be easily taken into account as a perturbation (to evaluate power losses in metals) [7] or even retained as an impedance boundary con- dition [4].

•

The antenna region is supposed totally free of charged particles, i.e. in vac- uum, for the (possible) presence of charged particles around the antenna neg- ligibly affects the wave propagation, owing to the steep density gradient in the plasma region.

Eventually, we assume time-harmonic EM fields and imply a common factor

exp(jωt) throughout.

2.3.1 Application of the Equivalence Principle

As we will repeatedly invoke the EP, it is expedient to recall its two relevant asserts [4]:

1. the EM field within a given region of space V (even unbounded), enclosed by a surface S and containing arbitrary media and sources, is exactly equivalent to the EM field radiated by suitable equivalent magnetic and electric surface current densities placed on S;

2. the equivalent problem so defined departs from the original one as per the fields outside V , which are null.

While solving antenna problems, the latter remarkable result is commonly ex- ploited by filling the null-field region with a medium chosen so as to ease the computation of the EM fields within the volume of interest V [4]. In the present case, we need to employ the EP twice to achieve the desired self-consistent sepa- ration of the antenna from the plasma in the geometry shown in Fig. 2.2.

To begin with, we apply the EP to the plasma region, i.e. a suitable volume VC,0 enclosing the cylindrical plasma beam limited by the aperture surface SA and a fictitious boundary ST , whereon perfect electric boundary conditions (PEC) are assumed. The role of the fictitious boundary as well as the validity of this assumption will be better clarified later.

As a whole the surfaces SA and ST constitute an infinite circular cylinder, whereon the EP is applied. Namely, as shown in Fig. 2.3a, on SA and ST we introduce equivalent magnetic and electric surface current densities, to wit:

M C,0 = E0 × (−nˆ0), X X,0 = (−nˆ0) × H0, (2.1)

with nˆ0 being the unit vector normal to SA and pointing outward the volume1, E0 and H0 are the (yet to be determined) electric and magnetic fields on SA and ST .

As just reminded, the currents (2.1) radiate the right EM fields within VC,0 and no fields at all outside. Thus, we take advantage of this by closing the aperture SA with a PEC cylinder, in order to attain the desired separation from the antenna region. Then, we notice that, due to the boundary conditions at a PEC interface, the current X X,0 radiates no field [5] and can be discarded, whereas the magnetic current M C,0 is better written as:

M C,0 = M A+ = E0 × (−ρˆ), (2.2)

wherein the subscript A+ is to remind that the current is placed at an infinitesimal distance from SA inside the VC,0 and exist only over SA (see Fig. 2.3b).

At this point, it is quite apparent that the fictitious PEC wall ST serves to limit the computational domain. This assertion indeed has a double interpretation: on

1To say the truth, the unit normal vector is customarily taken as pointing inward the volume wherein the fields are computed; our choice, however, appears motivated for it unifies the notation in view of the next application of the EP.

PLASMA FLOW

SA+

nˆ

ICRH ANTENNA

SA+

zˆ

SA+

M C ,0

X X ,0

M C ,0

X X ,0

FICTITIOUS

BOUNDARIES

(a)

PLASMA FLOW

M A+

VC,0

M A+

SA+

nˆ0

SA+

zˆ

SA+

PEC WALL

(b)

Figure 2.3: First application of the EP to the geometry of Fig. 2.2: (a) equivalent surface current densities are introduced on the surfaces SA+ and ST and (b) SA+ is closed by a PEC cylinder, therefore only the magnetic current M A+ density contributes to the fields in VC,0.

the one hand, ST contributes to identify the plasma region as the volume inside a cylinder, on the other, it allows limiting the support of the unknown magnetic current M A+ to the surface SA of finite extent. In order to justify the introduction of ST , we note that the coupling of the RF waves to the plasma occurs mostly in a well defined region downstream the ICRH antenna and also that the fields decrease quite rapidly away from the antenna. Therefore, we do expect the presence of ST not to yield a major effect on the solution and ultimately on the antenna input parameters, which has been confirmed by numerical simulations.

To proceed further, we apply the EP in the antenna region, i.e. a suitable vol- ume VC,1 bounded by the surface SC,1. As can be seen in Fig. 2.4, part of SC,1 enfolds all conductors (antenna, thruster walls), the NP feeding apertures SP,k and the aperture SA, and part actually represents a fictitious PEC boundary, whose sole purpose again is to limit the computational domain. Concerning this, considera- tions quite similar to the ones discussed above for the surface ST still hold. In fact, the validity of our ICRH stage model can be backed as follows: for one thing, the whole structure size is usually several order in magnitude smaller than the vacuum wavelength, in view of the very low operating frequency (about 1 MHz is com- mon), then we have the distance (h) between the side fictitious walls large enough as compared to the antenna size. Under those circumstances, we do expect the fic- titious boundaries not to affect the antenna parameters significantly, although their

SP,k

XX-

X X,1

MC ,1

ANTENNA PORT

ANTENNA

SA-

nˆ1 zˆ

PLASMA FLOW

FICTITIOUS BOUNDARIES

nˆ1

X X,1

SC V

MC ,1

C,1

THRUSTER WALLS

(a)

SC X X,1 SA- M A−	M	P, k nˆ
			1 XX-


X X,1 nˆ1 J C,1			M A−

		SC VC,1

SP,k

zˆ

(b)

Figure 2.4: EP applied to the antenna region of the geometry shown Fig. 2.2: (a) equivalent current densities are introduced on a surface SC,1 partly wrapping all conductors, the aperture SA, the feeding apertures SP,k and partly constituting a fictitious boundary; (b) all conductors and the plasma are removed so as to obtain a classical problem of EM field propagation in free space.

sensitivity to h can be easily studied.

The equivalent magnetic and electric current densities to be placed on SC,1 are:

M C,1 = E1 × nˆ1, X X,1 = nˆ1 × H1, (2.3)

wherein nˆ1 is the unit vector normal to SC,1 and pointing inwards VC,1, and E1, H1 are the unknown fields on SC,1. To ease calculation of the EM fields in VC,1 and to achieve the separation from the plasma—and from the feeding network as well—we end the application of the EP by filling the volume outside VC,1 with vacuum, i.e. the same medium inside the cavity, which entails formally removing all conductors, feeding lines and also the plasma, as shown in Fig. 2.4b. Since this way now the currents (2.3) radiate in vacuum (free space), the problem has become far simpler than the original, as desired. Furthermore, while X X,1 exists all over SC,1, M C,1 is nonzero only over SA and SP,k, due to the pristine boundary condition at a PEC surface, thus it is better specified as:

Σ=

M C,1

k=1

e# × nˆ

, on ∪ S

= M A− = E1 × ρˆ, on SA,

(2.4)

P,k

−

wherein the subscript A means that the current lies infinitesimally close to SA but inside the antenna cavity and it exists only over SA. Finally, e# (Vk) is the dominant TEM modal eigenfunction (voltage) [8] of the k-th feeding coax.

As argued in [1], the second of (2.4) yields the most accurate description of the coax feeding; hence, upon setting the NP voltages Vk, the resulting magnetic currents M P,k come to represent the independent sources of our EM problem.

2.3.2 Statement of the equations

The currents M A−, M A+ and X X,1 introduced so far constitute three unknowns to be determined, thus we need to state as many equations. To this aim, we note that inside VC,1 the fields E1 and H1 can be split into primary (labelled by p), gener- ated by the sources M P,k, and scattered (or secondary, labelled by s) radiated by the equivalent currents X X,1 and M A−; within VC,0, instead, only the secondary fields generated by M A+ exist.

\ \

The equations that M A−, M A+ and X X,1 satisfy ensue upon enforcing the boundary conditions that the total tangential fields must fulfill over the surfaces whereon the EP has been applied. To be concrete, the tangential electric field must either vanish on SC,1 SA SP,k or equal the relevant magnetic current densities on SA and SP,k, and finally both tangential magnetic and electric fields must be continuous across SA. With the aid of the characteristic functions:

χ rα

( ) = 1, r ∈ Sα

0, r /∈ Sα

, α = A, (C, 1), (P, k), (2.5)

the preceding statements can be written in symbols as:

tan

χC,1 (Ep + Es)| = nˆ1 × (χP M P + χAM A−) (2.6)

tan

0 tan

χA (Hp + Hs)| = χAHs| , (2.7)

M A− = −M A+, (2.8)

where equation (2.8) is a trivial consequence of (2.2) and the first of (2.4), and merely means a reduction of the number of unknowns from three to two.

To complete the derivation of (2.6), (2.7), we have to link the fields to their sources.

In the antenna region the pertinent relationships are surface integral operators involving the free-space tensor Green’s function [8] that is known analytically. For the sake of brevity, we will not report the resulting formulas for the scattered (Es,

Hs) and primary (Ep , Hp

) fields, as they follow immediately from Eqs. (17),

1 k1 k1

(18), (60), (61) of [1], upon formally appending the proper subscripts to fields, sources, unit vectors and integration domains.

As for the scattered magnetic field Hs over the surface SA, its dependence on the unknown magnetic current density M A+ can still be stated through a surface integral operator, namely:

Hs(z, θ, a) × ρˆ = ∫ d2r'YP (z − z', θ − θ', a, a') · M A−(z', θ') × ρˆ, (2.9)

wherein (2.8) has tacitly been used, primed (unprimed) quantities denote observa- tion (source) points on SA, a is the plasma radius (or more generally the radius of the cylindrical volume VC,0) and the kernel YP is the spatial plasma tensor Green’s function, whose calculation is addressed in Section 2.4.

− −

That the kernel has to be convolutional, i.e. depend on θ θ', z z', is a conse- quence of the speculated rotational (along θ) and longitudinal (along z) invariance of the plasma. From a physical standpoint, and in accordance with Fig. 2.3b, the entries of YP represent the component of the tangential magnetic field excited in (z, θ, a) by an infinitesimal magnetic dipole of unit intensity located at (z', θ'), namely, on the inner wall of a plasma-filled PEC cylinder.

Due to (2.4) or equivalently (2.2), we are as well permitted to deem (2.9) a quite general form of boundary condition that links the tangential magnetic and electric fields on the separation surface SA.

Finally, it is worthwhile noticing that the formal separation between antenna and plasma has been maintained in the equations, for the plasma enters only (2.7) via YP .

2.4 The plasma Green’s function

The main difficulty associated with the direct use of (2.9) is that for a magnetized anisotropic and inhomogeneous cylindrical plasma YP cannot be given in closed form and as such has to be computed numerically.

To this end, we have to solve the time-harmonic Xxxxxxx’x curl equations for the EM fields within the volume VC,0, i.e. the plasma, which in turn enters through

−

a current density X X . The solution is further complicated, for the plasma is also spatially dispersive, meaning that in the spatial domain the constitutive relation- ship J P = J P (E) is an integral one [9], whose kernel is the conductivity tensor σ(ω; r, r'). In the important case of a spatially homogeneous equilibrium, σ de- pends not on r and r' separately but rather on r r' [9]. Accordingly, fields and currents can be represented in the spectral (Fourier) domain, where calcula- tions become simpler, as the plasma constitutive relationship—and (2.9) for that matter—becomes algebraic.

The mathematical tool behind the transformation is the following Ansatz for fields and sources:

A (z, θ, ρ) = 1

∫ dk Σ e−jkzz−jmθ A˜

(k , m; ρ), (2.10)

−∞ ∞

4π2a

with ν = z, θ, ρ and m (kz) the azimuthal (longitudinal) modal index; both the integral and the summation extend from to + , but kz spans a continuous spectrum of modes, whereas m takes on discrete values.

As alluded at before, in the spectral (kz, m) domain, the plasma current reads:

J˜P (kz, m; ρ) = σ˜(kz; ρ) · E˜(kz, m; ρ), (2.11)

wherein σ˜ constitutes the plasma conductivity tensor. Upon pairing the displace- ment and the plasma currents, we formally obtain the plasma permittivity tensor as

ε˜ = ε0I − jσ˜/ω, (2.12)

wherein ε0 is the vacuum permittivity and I = ˆzˆz + θˆθˆ the identity tensor; the explicit form of ε˜is a somewhat delicate issue and is addressed in Section 2.5.

Thanks to (2.10), the link between YP and its spectral representation can be written as:

Y (z − z', θ − θ', a, a') = 1

∫ dkzΣ e−jkz(z−z')−jm(θ−θ')Y˜P

(kz, m; a, a'),

(2.13)

4 2π a

with

˜ '  Yθθ Yθz 

Yzθ Yzz

˜−1

YP (kz, m; a, a ) =

= ZP . (2.14)

It can be anticipated that, despite the presence of the magnetizing field that makes the plasma gyrotropic (i.e. non-reciprocal), the tensor Y˜P proves to be symmetric, as symmetry of the spectral Green’s function is not a sufficient condition for the background medium to be reciprocal.

Although in principle YP could be obtained by carrying out the integral and the summation indicated in (2.13), however, in Section 2.6 we will show that Y˜P is all we need to effect the numerical solution of (2.6) and (2.7).

Eventually, inserting (2.13) in (2.9) and performing a little algebra yields:

A−

H˜ s(k , m; a) × ρˆ = Y˜ (k , m; a, a') · M˜ (k , m; a') × ρˆ,

= −Y˜ (k , m; a, a') · E˜s (k , m; a), (2.15)

which will be used to compute Y˜P in what follows. In this form, the meaning of (2.15) as a boundary condition—which assumes different values for any pair (kz, m) of spectral variables—is far more apparent than it was in (2.9).

2.4.1 Reduced Xxxxxxx’x equations in the plasma

To compute the EM fields within the plasma beam, we start from the source-free Xxxxxxx’x curl equations and recast them in cylindrical coordinates ρ, θ, z. For the sake of simplicity the symmetry axis of the ICRH unit is taken coincident with the z axis of the reference frame. Under this assumption, the plasma permittivity tensor in cylindrical coordinates is represented by the same matrix as in cartesian coordinates, viz.

 ερρ ερθ 0   εxx εxy 0 

ε˜(kz; ρ) = 

εθρ εθθ 0

0 0 εzz

 = 

εyx εyy 0

0 0 εzz

 , (2.16)

with the well known properties [13]:

εxx = εyy, εyx = −εxy. (2.17)

Then, upon substituting each field component with its spectral representation as given by (2.10), making use of (2.17) and the constitutive relationships within the plasma:

D˜ = ε0ε˜· E˜, (2.18)

B˜ = µ0H˜ , (2.19)

with ε0 (µ0) the vacuum permittivity (permeability), we arrive at a set of six first order differential equations, where the fields still depend on the radial coordinate ρ. After a great deal of tedious but straightforward algebra omitted for brevity, one

obtains four equations involving only the transverse-to-ρˆfield components E˜θ , E˜z ,

Hz, (2.20)

H˜θ and H˜z , that read:

ρ (ρEθ) = jm

∂ρ ε

∂ ˜

εxy ˜

jmkz ˜

2 m2 ˜

ρEθ − ωε ε

ρHθ −jωµ0

− k2ε

∂ ˜

ρ (ρHθ) =

∂ρ

ρEθ + jωε0εzz

jmkz ˜

ρEθ + jωµ0

ωµ0

0 xx

m2 ˜

k2εxx

ρHθ +

Hz,(2.22)

0 xx

− k2ε

Ez, (2.21)

Ez = jkz

∂ ˜

εxy ˜

k2 ˜

jmkz ˜

∂ρ

z 0

εxx

1−

ρ ∂ H˜ = −jωε

2 + ε2

−

k2 !

ωε0εxx

ρE˜θ

∂ρ εxx 2

+jk

εxy ρH˜

−

z εxx θ

jmkz E˜ ωµ0 z

jmεxy H˜

−

εxx z

, (2.23)

where k0 = (ε0µ0)1/2 is the vacuum wavenumber. The remaining two equations provide the radial components E˜ρ and H˜ρ as a function of the transverse fields only, namely:

H˜ = m E˜ − kz E˜ , (2.24)

ρ ωµ0ρ z ωµ0 θ

E˜ = − m H˜

+ kz H˜

− εxy E˜ . (2.25)

ρ ωε0εxxρ z

ωε0εxx θ

εxx θ

Equations (2.20)-(2.25) have to be supplemented with proper boundary and regu- larity conditions in order for the solution to be unique and represent an EM field.

In that regard, since we need the admittance tensor, the logic would suggest to force the electric field component at the air-plasma interface ρ = a and to deter- mine the magnetic field at the same position. However, proceeding that way, the numerical solution proves to be stiff and the convergence poor, for, as will be seen,

the electric field E˜θ exhibits a steep variation just inside the plasma and close to

the air-plasma boundary. Accordingly, the admittance Xxxxx’s function obtained through this procedure happens to be less accurate and possibly non-physical. To circumvent this hurdle, since the relationship between tangential fields at the air- plasma interface (2.15) can also be written as:

E˜s(k , m; a) = −Z˜ (k , m; a, a') · H˜ s(k , m; a') × ρˆ, (2.26)

it clearly suffices to enforce the magnetic field components and determine the elec- tric ones, which intrinsically points at the impedance Xxxxx’s function Z˜ as a result.

Now, in order to evaluate the four components of Z˜, thanks to linearity, it is convenient to impose the following sets of boundary conditions at ρ = a:

H˜θ H˜z

= 1 , (2.27)

= 0 , (2.28)

alternatively, whence in the light of (2.26), the entries of the impedance Xxxxx’s function follow as:

  

E˜θ (a)

−

E˜θ (a)



Zθθ Zθz

 Zzθ Zzz

 = 

H˜ (a)

z ˜



Hθ(a)=0

Ez(a)

(

)

−

H˜ (a)

θ ˜

Hz (a)=0

Ez(a)

H˜

(a)

(2.29)

Hz a H˜θ (a)=0 θ H˜z (a)=0

wherein apparently the first and the second column have been obtained applying (2.27) and (2.28), respectively.

To say the truth, there would be no need to determine the two off-diagonal entries separately, since the tensor is expected to be symmetric in theory. However, a number of factors (e.g. machine precision, numerical round-offs, convergence of the numerical algorithm) can actually prevent Zθz and Zzθ from being perfectly coincident, thus an independent computation is rather useful, in that it helps check the accuracy of the numerical algorithm.

Eventually, according to (2.14), a simple inversion provide the sought for Y˜P ,

which is not critical from a numerical standpoint.

→

Concerning the regularity issue, since ρ = 0 constitutes a singular point of the equations (2.20)-(2.25), we must exclude solutions that are unbounded or discon- tinuous on the z-axis. Stated another way, acceptable fields must possess a finite and unique (i.e. independent of θ) limit as ρ 0. This requirement translates into a condition for any spectral component, to wit:

A˜ν (m, kz ; 0) = constant m = 0

, (2.30)

0 m 0

in view of (2.10). As discussed in Section 2.4.2 below, (2.30) has to be directly included in the numerical solution of (2.20)-(2.23).

2.4.2 FEM solution

A brand new module of TOPICA has been purposefully coded to solve (2.20)- (2.23) via a Finite Element Method (FEM) numerical scheme using a set of piece- wise linear interpolating functions (triangular-pulses) defined as:

f (ρ) = ρ1 − ρ , ρ ≤ ρ ≤ ρ ,

ρ1 − ρ0

− −



ρ − ρn−1 , ρ

≤ ρ ≤ ρ

fn(ρ) =

ρn ρn 1

ρn+1

− ρn

 ρn+1 − ρ

n−1

, ρn ≤ ρ ≤ ρn+1,

n = 1, . . . , N, (2.31)

fN +1

(ρ) = ρ − ρN , ρ

ρN +1 − ρN

≤ ρ ≤ ρN+1,

f 0

f n

f N +1

ρ0 ρ1

ρn−1 ρn

ρn+1

ρN ρN +1

Figure 2.5: Qualitative representation of the triangular pulses defined in Eq. 2.31.

wherein ρ0 = 0, ρ1, . . . , ρN+1 = a represent N + 2 nodes distributed on a non- uniform grid (see Fig. 2.5).

The FEM solution begins with expressing the transverse field components as linear superposition of functions fn, namely:

ρH˜θ Z0 = fn(ρ)I1,n + fN+1I1,N+1, (2.32)

n=1

H˜z Z0 = fn(ρ)I2,n + fN+1I2,N+1, (2.33)

n=0

N +1

ρE˜θ = fn(ρ)V1,n, (2.34)

n=1

N +1

ρE˜z = fn(ρ)V2,n, (2.35)

n=0

wherein Z0 = (µ0/ε0)1/2 stands for the vacuum intrinsic impedance and serves as a normalizing constant. Equations (2.32)-(2.33) also provide a simple mean to incorporate the boundary conditions (2.27)-(2.28), for:

aZ0H˜θ (a) = I1,N+1, (2.36)

aZ0H˜z (a) = I2,N+1, (2.37)

thus I1,N+1, I2,N+1 are to be considered known quantities, as a matter of fact.

Furthermore, it is worth explaining why ρH˜θ and ρE˜θ rather than the bare

θ-components were expanded in terms of fn. For one thing, that choice implies I1,0 = V1,0 = 0, which evidently allows reducing the number of overall unknown coefficients to be calculated. Even more to the point, however, the assumption appears most useful in that it naturally xxxxx with the regularity issue as per the

θ-components: in fact, since ρH˜θ and ρE˜θ intrinsically vanish on the z-axis, we

rest assured that the numerical algorithm yields solutions that possess a unique and finite limit as ρ → 0.

By contrast, regularity of E˜z and H˜z on the z-axis has to be imposed directly, which, in the light of (2.30), implies

I2,0 = V2,0 = 0, m 0, (2.38)

whereas the case m = 0 requires just leaving I2,0, V2,0 as unknowns in (2.33), (2.35). In consequence, the total number of unknown coefficients comes to be 4N + 4 for m = 0 and 4N + 2 for every other m.

To proceed further with the FEM, we substitute (2.32)-(2.35) into (2.20)-(2.23), and then by means of the symmetric inner product:

∫ a

< a, b >=

dρ a(ρ)b(ρ). (2.39)

we project the differential equations onto the expansion functions fn(ρ). The over- all procedure yields a very sparse algebraic system, whose matrix is complex and not symmetric.

System inversion is carried out by standard LU factorization, which is preferred to an iterative method, since we need to solve the system twice with the boundary conditions (2.27)-(2.28). The solution does not pose particular challenges until the spectral variable kz is not large as compared to k0, otherwise the first term on the right hand side of (2.23) can become dominant over the other contributions and make the whole system poorly conditioned. When necessary, we cure the problem by means of a Xxxxxx’x preconditioning procedure before applying the LU factorization.

After the unknown expansion coefficients in (2.32)-(2.35) have been obtained, we can compute the plasma impedance Xxxxx’s function through (2.29), that now reads:

  

Z0V1,N +1

—

Z0V1,N +1 

Zθθ Zθz



 = 

aI2,N +1 Z V

I1,N +1=0

I1,N +1



Z aV

I2,N +1=0

 . (2.40)

Zzθ Zzz

 — 0

2,N +1

0 2,N +1

I2,N +1

2.5 The plasma model

I1,N +1=0

I1,N +1

I2,N +1=0

In this Section we discuss the explicit form of the constitutive relation of the plasma (2.18) in the spectral domain, which is used by the numerical module that eventu- ally provides the plasma surface Green’s function Y˜P through (2.40).

A review of the literature was carried out to select an accurate plasma model: in the end we adopted the same plasma description as in [11, 12, 14]. Specifically, the model, which assumes a linearized warm collisionless plasma, allows for radially inhomogeneous density ne,i and temperature Te,i profiles both for electrons and ions. More importantly, we also account for the macroscopic plasma flow velocity v, which manifests its effect by shifting the ion cyclotron frequency [12]. As a

further assumption, we consider the high speed plasma flow, occurring in the RF thrusters, to rapidly and wholly absorb the ion cyclotron waves launched by the antennas [10].

It seems not superfluous to remind at this point that the velocity v will not be chosen arbitrarily but rather will come out from the global plasma model developed in this very project.

Zxx x−1α

, (2.41)

Basically, the derivation of the permittivity tensor entries start with the solution to the linearized Vlasov kinetic equation coupled with the Xxxxxxx’x equations. As the procedure can be found in any book of plasma physics, such as [9], we only report the result, for the sake of brevity. To be concrete, when the antenna operating frequency is the order of the fundamental ion cyclotron resonance Ωcα, as in the case of interest here, the plasma tensor entries take on the form:

εxx ≈

—

pα

2ω2

Zxx x+1α

1 Σ ω2 [ (

) + ( )]

εxy ≈ —j

2 2ω

Σ ω [ ( ) ( )]

pα Zxy x+1α — Zxy x−1α , (2.42)

εzz ≈ 1 —

Σ ω

pα Zzz

ω2

(x0α), (2.43)

wherein the summation is over electron and all ion species, ωpα is the plasma frequency for the species α. Furthermore:

Z (x

) = — T⊥α — 1 — Ωcα + u0α + T⊥α ζ

Z(ζ

), (2.44)

xx +1α

T||α

kzvth||,α 2

T||α

+1α

Z (x

) = — T⊥α — 1 — Ωcα + u0α + T⊥α ζ

Z(ζ

), (2.45)

xx −1α

T||α

kzvth||,α 2

T||α

−1α

Zxy(x+1α) = Zxx(x+1α), (2.46)

Zxy(x−1α) = Zxx(x−1α), (2.47)

Zzz(x0α) = —2

ω2 k2v2

! [1 + ζ0αZ(ζ0α)] , (2.48)

wherein

ζ = x

z th||,α

0α

— u = ω — Ωcα — kzv||α , (2.49)

+1α

kzvth||,α ω

ζ = x

+ u = ω — Ωcα — kzv||α , (2.50)

0α

−1α −1α

kzvth||,α ω

0α

ζ = ω — kzv||α , (2.51)

kzvth||,α ω

x0α =

ω kzvth||,α

, (2.52)

vth||,α

= 2κT||α , (2.53)

mα

with Ωcα the cyclotron frequency, vth||,α the parallel-to-ˆz thermal velocity, v||α, the parallel-to-ˆz flow velocity, mα the mass, T⊥α (T||α) the perpendicular-to-ˆz (parallel-to-ˆz) temperature, each for the species α with evident notation, and κ the Boltzmann’s constant. Finally, Z stands for the plasma dispersion function [9] defined by:

∫

1 ∞

Z(ζ) = √π

−∞

e−t2

t — ζ

— j√

πSign(k )e−ζ2 , Imζ = 0, (2.54)

when ζ is either ζ+1α or ζ−1α or ζ0α.

The resulting tensor entries (2.41)-(2.43) are quite the same as in [11], save for few apparent typos that mar the formulas displayed therein and perhaps a bit different notation. However, the approach we follow to include the plasma effects within the formulation departs substantially from the procedure outlined in [11].

In fact, we notice that in [11, 12], after the plasma tensor is formally derived in the spectral (Fourier) domain, i.e. as a function of k|| = kz, it is then evaluated only for specific values of k||(r, z), which reportedly are the slow wave solution to the kinetic dispersion relation for parallel propagation. Accordingly, ε˜ is interpreted as function of (r, z) rather than of kz and as such it is substituted in the electric field wave equation in the spatial domain, for subsequent solution via a 2D finite difference scheme [12].

By contrast, our approach hinging on a integro-differential formulation, we include the plasma effects through a suitable Green’s function YP or better, as we will show in Section 2.6.2, via its spectral counterpart Y˜P (kz, m), whose numerical calculation was addressed in Section 2.4. The latter entails solving Eqs. (2.20)- (2.23) and correspondingly knowing ε˜(kz) over a suitable range of kz. In that regard, the largest value of kz—and m for that matter, but ε˜(kz) is independent of the azimuthal wavenumber—is determined by the convergence properties of the spectral integrals arising in the numerical solution of (2.7) via the Method of Moments, as discussed in Section 2.6. We anticipate that, in consequence, kz can commonly take on values as large as 104 1/m.

Now, numerical evaluation of (2.41)-(2.43) is straightforward even for large values of kz, the main issue being the calculation of the plasma dispersion function (2.54). Nevertheless, it can be proven that for certain combinations of legit plasma parameters (2.41)-(2.43) can yield non-physical results, even for moderate values of kz. To be concrete, in Fig. 2.6 we report the imaginary part of εzz as a function of kz and computed for a typical set of relevant parameters. As can be seen, although the top plot suggests that Imεzz ought to be monotonically negative, however the bottom enlarged view clearly highlights the occurrence of a null and a consequent sign change.

Figure 2.6: Plasma tensor entry εzz: (top) imaginary part as a function of the longitudinal wavenumber kz [1/m]; (bottom) enlarged view showing the unwanted occurrence of a null and consequent sign change at approximately kz = 48 1/m.

Thus, more generally, non-physical results may take place if the plasma con- ductivity tensor—essentially the anti-hermitian part of ε˜—exhibits entries that do not possess the proper sign for the plasma to be a passive medium, meaning

∂W = 1 E∗ · σ˜ · E ≥ 0, (2.55)

∂t

with ∂W/∂t the average power density. Otherwise stated, (2.55) is not fulfilled, if σ˜ does not happen to be a positive definite tensor; this circumstance ultimately manifests itself as a negative real part of the Green’s function diagonal elements

Y˜θθ , Y˜zz , or as a negative antenna loading as well. On the other hand, a plasma that

does not absorb but rather emits energy seems at odds with what is usually recorded in laboratory experiments on thrusters; therefore, we are induced to believe that non-physical entries of σ˜ do not correspond to some overlooked exotic plasma behavior.

More simply, the possible inadequacy of (2.41)-(2.43) can be traced back to the very plasma flow or rather to the way we inserted it in the model. In fact, accord- ing to [9], we accounted for the parallel velocity of ions and electrons by means of a displaced maxwellian distribution centered around the average (macroscopic) velocity v. The latter approach, which is neat and well suited for energetic parti- cle beams, works pretty fine for small or moderate values of kz/k0, but it breaks down, as argued above, for high kz/k0 ratios, such as the one needed to carry out the spectral integrals (cf. (2.66) farther below) required by TOPICA.

—

To circumvent this drawback, we also tried a less accurate version of the model consisting of a cold plasma tensor (i.e. independent of kz), wherein a heuristic fac- tor in the form of a non-null collision frequency takes into account the absorption, which boils down to a complex frequency ω jν. The underlying idea is pretty much the same as in [12], wherein, though, the cold plasma tensor is modified upon inserting a complex particle mass. Seemingly, the ε˜ obtained this way does not suffer from the inconsistencies outlined previously, but nonetheless it only rep- resents a fairly realistic situation, to say nothing of the uncertainty associated with the choice of the collision frequency.

Then, in order to determine the antenna loading while at the same time main- taining an accurate plasma model, we provisionally resorted to an intermediate approach, which proved to be less sensitive to the flaws of ε˜ as given by (2.41)- (2.43): the details of the method will be discussed in Section 2.6.3.

Meanwhile, since the implementation of the module computing the plasma Green’s function (Section 2.4) as well as of the module filling the matrix [GP ] (Section 2.6.2) is largely independent of the contingent form of ε˜, we continue our quest for other plasma models that can provide well-behaved tensors in wider kz ranges.

2.6 Solution by hybrid Moment Method

In this Section we address the solution of (2.6), (2.7) by reducing them to an al- gebraic system in a two-step weighted-residual finite-element procedure; in the context of computational electromagnetics this approach is ubiquitously known as “Method of Moments” (MoM) [6], and we will preserve this name in what fol- lows. Moreover, since some entries of the resulting system matrix are computed in the spatial domain and some other in the spectral domain (see (2.66)), we call this solution scheme a “Hybrid Spatial-Spectral Method of Moments”.

The procedure here goes along the same lines extensively described in Section 4 of [1] and therefore we will only summarize the main issues.

2.6.1 Algebraic system

To begin with, the unknown vector functions X X and M A− are approximated by a

linear combination of NC and NA vector basis functions {f n}NC , {gm}NA with

unknown coefficients In, Mm, to wit

n=1

m=1

X X(r) = Inf n(r), (2.56)

n=1

ΣM r M r M g r .mA− m

( ) = NA ( ) = ( ) (2.57)

A−

m=1

Secondly, (2.56) and (2.57) are inserted in (2.6), (2.7) and an algebraic system is obtained upon projecting the IEs onto a set of weighting (test) functions, namely

⟨f k, (2.6)⟩ = 0, k = 1, . . . , NC, (2.58)

⟨gl, (2.7)⟩ = 0, l = 1, . . . , NA, (2.59)

wherein

⟨a(r), b(r)⟩ ≡ ∫

d2r a(r) · b(r) (2.60)

is a non-hermitian inner product defined via integration on a suitable surface Σ. In the scheme to be used (also known as Galerkin’s method) the set of test functions is identical to the sets of basis functions.

The final system can be written succinctly as:

" [ ] [ ] #

1 [ ]  

− [ ]

— G11

— [G

] [G

G12

] — [G ] 

—Z0 2 [M ]

− =

Z0 2

Z 2 [H]

 , (2.61)

∫

where evidently [I], [M ] are column arrays that collect all the unknown expansion coefficients appearing in (2.56)-(2.57) and the system matrix is size (NC + NA)- by-(NC + NA). The elements of the matrix [GP ] are given by:

P = Z0

d2ρgl(ρ) × ρˆ· ∫

d2ρ'YP (ρ — ρ') · gk(ρ') × ρˆ, (2.62)

whereas the entries of the other blocks are displayed in Equations (45)-(48) of [1]. Finally, the forcing term blocks, that depend on the primary field, have elements given by (50)-(51) again in [1]. The system matrix in (2.61) is usually called inter- action matrix, and its entries are said interaction- or reaction integrals, because of their meaning in the context of Xxxxxxx equations with respect to the Reciprocity Theorem [5]. As a result of this theorem, in a reciprocal medium (as in the an- xxxxx region) the overall matrix can be proved to be symmetrical, in particular this translates into the notable identities:

[G11]T = [G11] , [G22]T = [G22] , — [G21]T = [G12] , (2.63)

so that only half of the entries need to be computed. No such property holds true for a magnetized plasma, hence the block matrix [GP ] is not symmetric and all its N 2 entries are to be evaluated.

2.6.2 Spectral reaction integrals

As it was shown in Section 2.4, the kernel YP may be naturally expressed in the spectral domain, due to the translational invariance of the Green’s function over the cylindrical air-plasma interface. Thus, it is necessary to express the re-

action integrals GP

given by (2.62) employing Y˜P provided by (2.14). This is

accomplished by inserting (2.13) into the reaction integrals (2.62) and substituting

gl(θ, z), gk(θ', z') by their Fourier transform, to wit

∫ ∫

g˜l(—m, —kz) = a ∫ ∫ dzdθ gl(θ, z)e−jkzz−jmθ, (2.64)

g˜k(m, kz) = a dz'dθ' gk(θ', z')ejkzz'+jmθ' , (2.65)

with Tl (Tk) denoting the domain whereon the basis function gl (gk) is non-zero. Then performing a change of order of integration and a little algebra yields the result:

GP = Z0

∫ dk

Σ g˜ (—m, —k )×ρˆ·Y˜

(m, k ; a, a')·g˜ (m, k )×ρˆ. (2.66)

4π2a

Now, transforming the reaction integrals GP from the spatial to the spectral

domain has turned the two-fold double integration along θ and z into an integration over kz extending over the whole real axis and an infinite summation over m, which demands for a discussion on the convergence. It can be shown, however, that the

asymptotic convergence, i.e. for large values of m and kz, is guaranteed by the asymptotic behavior of Y˜P and of the Fourier transform of the basis functions, if these latter are correctly chosen, as described in Section 5 of [1].

In practice, to save computational time the integral and the summation are re- arranged in order to involve only positive values of kz and m. Integration over

∞ ∞

[0, + ] is carried out by trapezoidal rule and stopped at a suitable kz,max < + , which may depend strongly on the distance and the support dimension of the two basis functions involved in (2.66) and obviously on Y˜P . Thus kz,max is adap- tively chosen so as the corresponding value of the integrand function is smaller than an externally provided relative threshold δ (suitable values range from 0.01 and 0.0001). Convergence of the summation over m is addressed in the same man- ner but is far faster and does not pose particular concerns.

As regards the time required to perform the spectral calculations, this has been considerably reduced with respect to the past, thanks to the new parallelized ver- sion of TOPICA, that now can be run on clusters of workstations. Although this was not an effortless task at all, it should be noted, however, that the process of filling a matrix like [GP ] is intrinsically parallel, as each integral-summation ap- pearing in (2.66) can be performed independently of all the others. Hence, the overall time required to fill a given matrix decreases with the number of available processors and ultimately shrinks to the very time needed to effect the a single entry calculation.

2.6.3 Antenna loading calculation

After the system (2.61) has been solved for the coefficients In and Mk, the surface current densities on the conductors and the aperture can be computed via (2.56)- (2.57). Incidentally, as a remarkable consequence of the adopted formulation, the tangential electric field at the air-plasma interface ensues immediately from M A− via (2.4). Knowing X X and M A− enables us to compute a number of parameters, such as the antenna input admittance (and quantities derived thereof), the radiated fields around the antenna as well as the field and power distribution inside the plasma. In this Section, however, we only focus on the derivation of the antenna admittance and the plasma loading resistance.

For a start, we assume for the time being that the antenna is a single-port device feeded by a coaxial cable. Then, we recall that the true forcing term, hidden in M P , is the modal voltage V # germane to the TEM electric mode eigenfunction propagating in the feeding coax. Therefore, once the modal current I# has been computed or, better, related to the fields in the cavity, the antenna admittance can be evaluated via:

Ya = I#/V #, (2.67)

where voltage and current are taken in the end section of the truncated coaxial cable. To determine I#, we have to consider the transverse magnetic field over the coax aperture SP , which, owing to continuity, is also equal to the transverse magnetic field of TEM mode just inside the coax, to wit

H1 = H# = h#I#, on SF , (2.68)

where h# stands for the TEM transverse magnetic eigenfunction [8]. Upon adopt- ing the customary normalization condition for the TEM transverse eigenfunctions

[8], namely:

⟨e#, h# × nˆ⟩ = 1, (2.69)

the modal current can be deduced from (2.68) by taking the inner product of both sides to nˆ × e#, viz.

V #

s p

I# = ⟨nˆ × e#, H

⟩ = —⟨M P , H1⟩ = —⟨M P , H1⟩ — ⟨M P , H1⟩, (2.70)

where use has been made of (2.7). Now, in view of definition (2.67), the antenna input admittance becomes:

Ya = —⟨M P , Hs⟩/V #2 — ⟨M P , Hp⟩/V #2, (2.71)

wherein it can be verified that generally the first term is dominant and need be calculated. By contrast, the second contribution in (2.71), which represents the self-admittance of the truncated coax radiating in free-space, is usually negligible. Besides, it can be proved that (2.71) is a variational formula [5], which makes it numerically robust and somewhat insensitive to the approximation error that affects both M A— and X X.

Although (2.71) already constitutes the sought for result, we can arrive at an

exactly equivalent expression more suited for numerical evaluation. In fact, one notes that the inner products involved in (2.71) represent so-called reactions of fields on sources [5]. In particular the first one is the reaction of the field Hs,

scattered by X X and M A—, on the source M P and it can be rewritten as

s p p

—⟨M P , H1⟩ = ⟨X X, E1⟩ — ⟨M A—, H1⟩, (2.72)

due to the strong form of Reciprocity Theorem [5]. The latter finds application here, since the equivalent problem obtained via the EP in the antenna region con- sists of surface current densities that radiate in free-space.

Equation (2.72) provides a more efficient way of computing the dominant re- action appearing in (2.71), in that involves directly the calculated current densities and the primary fields, which were already computed to get the forcing term in the RHS of the algebraic system (2.61).

Finally, to obtain the admittance matrix [Ya] in the more common instance of multiport antennas, one energizes the k-th coax with V #, while shorting out the remaining coax apertures, and computes the modal currents I# in each coax. Then, owing to linearity, the admittance matrix entries ensues through:

Ymk = m

V #=0, l/=k

, (2.73)

which generalizes equation (2.71). Even in this case Ymk can be linked to the sur- face current densities, upon expressing the modal currents in the form of reactions and invoking again the reciprocity.

The admittance matrix [Ya]—or equivalently any other set of circuit parameters derivable thereof—can be used to design the antenna and its feeding network and also to compute the radiated power PRF. On the other hand, to assess the antenna capability of conveying power to the plasma, it is customary to introduce a figure of merit, namely, the so-called coupling resistance,

Rc = 2PRF /|I|2, (2.74)

wherein I denotes the current entering the antenna terminals. Even though the coupling resistance, strictly speaking, is not an output of TOPICA code, it can be evaluated through (2.74), after PRF has been determined via standard circuit techniques applied to the antenna admittance matrix.

We notice that [Ya] appears quite sensitive to changes in the plasma density and temperature profiles, whereas the converse holds true of Rc, due to its definition, whence we do expect Rc to be meaningful, as it were, even when the adopted plasma model may become instable, as was pointed out in Section 2.5.

Therefore, as far as we are concerned mainly with the power transferred to the plasma by the antenna, it is desirable to obtain Rc through an alternate route that does not rely on the admittance [Ya]. To accomplish this, we resorted to an ap- proximate approach that works pretty well, providing the antenna does not exhibit too much geometrical complexity, in which case the full-wave solution becomes definitely necessary.

To be specific, in our simplified model we solely retain the plasma column and those parts of the surrounding ICRH antenna that strictly face the flow. Outside the plasma beam—which still comes into play via Y˜P —we account for an infinite free space region by means if the relevant spectral Green’s function, which is known

in closed form. The resulting electromagnetic problem is schematically depicted in Fig. 2.7a, wherein, for illustration purposes only, the antenna part facing the plasma has been assumed as comprised of two circular loops.

Once turned in the spectral domain, the problem of determining the power conveyed to the plasma and hence the coupling resistance is reduced to the solution of the equivalent modal circuit shown in Fig. 2.7b: the dyadic admittances Y˜P and

Y˜F S represent the spectral plasma and free space Green’s functions, respectively,

whereas the current generator is the spectral counterpart of the electric current densities directly facing the plasma beam.

Now, we can effect the calculation of PRF from the circuit with standard net- work analysis and then calculate Rc from (2.74). To this aim, a key step is the proper choice of the current density over the conductors and ultimately of the cur- rent generator in Fig. 2.7b. In this respect, we do not need to resort to any educated guess, but rather we are in a position of knowing the right current: as a matter of fact, we can solve the original EM problem with TOPICA, i.e. in a self-consistent and full-wave manner, and thus obtain, among others, the electric current distri- bution over the antenna. In consequence, using this result in our intermediate ap- proach affords a pretty high degree of accuracy.

FREE SPACE

CURRENTS

ρˆ

J2

PLASMA FLOW J1

YFS (kz , m)

Jfacing

ρˆ

YP (kz , m)

(a) (b)

Figure 2.7: Derivation of the coupling resistance: (a) cartoon of the adopted in- termediate model and (b) equivalent (spectral domain) circuit of the intermediate model for each (kz, m) pairs.

Chapter 3

Plasma device modelling

A numerical model that describes the physics of the thruster is useful for a number of reasons: it permits to understand the physics of the mechanisms involved, it gives direction on the optimization and design of the experiment, it permits to predict the performance of the thruster, and it contributes to reduce the number of experiments required.

The model resembles the three physical stages of the thruster: the plasma dis- charge, the RF plasma acceleration, and the magnetic nozzle.

3.1 Global model of plasma discharge

Aim of this section is to provide a global time dependent model of a helicon source developed for producing plasma density above 1018 m-3.

The model combines a global plasma source simulation with a 0-dimensional gas-dynamic simulation. It accounts for changes in the neutral density, ionization, excitation and dissociation. Outputs of the model are: the density of each neutral and ionized specie, the electron density and temperature and the time history of each particle/energy loss channel [19]-[31].

The device consist of helicon antenna wrapped around a quartz tube which serves as the vacuum boundary. It is surrounded by a solenoidal magnetic coil producing an axial magnetic. Another magnetic coil surrounds the exhaust area. The magnetic field configuration is adjusted by changing the currents applied to each magnetic coil. Figure 3.1 provides a general description of the system.

The neutral gas flow is provided by a tank through a duct connected by a dis- penser to one side of the source tube. The gas flows along the source chamber, through the exhaust, into the vacuum chamber, where it is evacuated by a turbo- molecular pump.

Magnetic coils

Helicon antenna

Quartz tube

Neutral gas inflow

plasma exhaust

Figure 3.1: Schematic configuration of the plasma generation stage. A neutral gas is injected from the left, it gets ionized by the helicon antenna and it is exhausted from the right. The coils are necessary to generate the magnetic field to confine the plasma.

3.1.1 Plasma reactions

A global model has been developed to better understand experimental observa- tions and to lead the experiment design. This approach is similar to other global models previously developed for simulating process plasma sources [19]-[31]. The plasma balance equations, for particles and energy, have been written for describ- ing a uniform distributed plasma inside of a region determined by the magnetic field configuration.

Many studies have been done about interactions between plasma and neutrals, including the effect of neutral losses to ionization [19]-[22], [40, 41] and neutral heating [32, 33]. Several models take into account the neutrals density inserting a source term and a sink term into the neutrals balance equation [19]-[23], [26]-[31]. These terms are related respectively with the feeding flow from the reservoir and the flow to the vacuum pump. In other models the plasma neutral interactions are not considered at all and no equations are written to follow the neutrals density behaviour [24, 25].

Due to the specific gas-dynamic configuration of the device, the neutral inter- action with plasma has been considered in this work by coupling a 0-dimensional gas-dynamic model of the entire system, with a global plasma model of the source. This model provides an estimate for the pop-off feeding-valve operation, ef- ficiency of neutral pumping by the vacuum pump, efficiency of a gas trap in the source to increase the ionization efficiency. The interactions that are taken into

account in the model are:

• neutral density reduction due to ionization;

• neutral dissociation (molecular species-atom species);

• zero dimensional gas dynamic analysis behavior in the plasma source and in

the vacuum chamber;

• wall recombination and volume recombination in the main vacuum chamber.

Plasma is generated in the source chamber. A preliminary investigation shows that in specific magnetic field configurations or in a specific operation mode (he- xxxxx mode), the plasmas could be confined inside of a volume smaller than the source chamber volume [15]. Plasma has been considered confined in a cylindrical volume Ve defined by a radius, plasma radius rp, and having the same length than the source chamber, L. Inside this volume different species are considered for ev- ery gas. The model follows the density of all of these species. Plasmas also flows and diffuse through the external surfaces of the volume Ve. These surfaces will be named in different way to highlight the different process involved. The back axial surface is the surface in front of the feeding orifice, plasma in this zone is electro- statically confined and the mass loss is calculated using Godyac and Maximov [39] solution of diffusion equation. Plasma also diffuse through the radial surface, but in this zone the magnetic field generated by the solenoid coil improves the confine- ment. The particle loss in this area has been calculated using again the Godyac and Maximov solution modified by Cheetham [24] to take into account the magnetic field contribute in the confinement. The axial surface toward the vacuum chamber is named exhaust surface. Plasma flows in this zone with a speed that is a fraction of the ion sound velocity. The speed strongly depend on the shape of the plasma potential in this area. Being this calculation beyond the purpose of this model, a parameter has been introduced into the numerical analysis named cs. Therefore the exhaust velocity is the ion sound velocity multiplied by the cs coefficient that has been considered as a free parameter.

Particles diffusing through the lateral surface and through the back axial surface are neutralized. As will be explained later, plasma equations are coupled with neutral equations since in the source chamber the neutrals density is not constant but free to change in relation to the neutral flow, the dissociation processes and the plasma-neutral interaction. The reactions involving ionized species and electrons are found in literature [37, 42]. The particle balance equations for the ionized particles and electrons are written in a particle flux form,(i.e. particles/(s m3)). The general form for the balance equations of charged particles is:

W,i

dni = Γs — Γl — Γ — Γ

EXH,i

, (3.1)

— 2

4πν2dν, (3.2)

wherein Γs represents the source term due to plasma processes for the i-species, ΓL is the loss term due to plasma processes, ΓW i is for the i-species the loss term due to particle recombination at the wall (the particle diffuses through the wall sheath before reaching the wall), and finally ΓEXH,i is the loss term due to the particle flow through the exhaust. The reaction rates were obtained averaging the cross section for the specific reaction over a maxwellian distribution [42]:

Kiz =

σ(ν)ν exp

m 3/2 ∫ ∞

πkTe

mν2

where Te is the electron temperature in eV, and m is the electron mass, and σ the cross section. Wall losses are calculated as done in other works in literature[19]- [21], [31]. Xxxx lost at the exhaust are computed as:

√

LEXH = niuBAEXH (3.3)

uB = kTe/mi (3.4)

with uB the ion Xxxx’x velocity and AEXH is the geometrical exhaust area.

Another parameter affects the exhaust flow: in fact, at the exit of the plasma generation section the magnetic field increases and then decreases. This peak acts as a magnetic mirror that partly reflects the plasma flow. Therefore the net flow is given by the difference between the incident flow and the reflected flow. The reflected flow depends on the configuration of the magnetic field and on the plasma parameters as explained later in the section "magnetic mirror".

To calculate the electron temperature, the power balance equation has been written as follows (units W/m3):

PABS

= d 3 en T + P + P

+ Σ P , (3.5)

EXH

PABS is the deposited power into the plasma, which is assumed to be known, e is the electron charge, ne the electron density, Te is the electron temperature, Ve again the plasma volume, PW is the power lost at the wall due to the electron-ions flow, PEXH is the power loss associated with the electron and the ion flux at the exhaust, assuming that the escaping velocity is the ion-Xxxx velocity. Eventually, the Pi terms constitute the power that is lost in the i-th reaction; their general expression is:

Pi = XxXXX,inenj, (3.6)

where Ki is the rate constant for the specific reaction, ETH—i the threshold energy for the ith-reaction [37], nj the density of the species involved in the ith-reaction. Experimental results [47] indicate the presence of a hot tail in the electron pop- ulation in hydrogen and helium discharge. This distribution has been modelled summing two maxwellian distributions: one with the temperature of the bulk of the plasma and one with the temperature of the hot tail.

3.1.2 Gas dynamic model

The gas-dynamic configuration of the engine is the following: a pressure reservoir provides, through a feeding line and a feeding valve, a mass flow of H2, He, or Ar. The gas expands in the source chamber and flows through the exhaust into the vacuum chamber where it is evacuated by a turbomolecular pump. The system was modelled considering three different volume connected by valves. Valve equation provides the mass exchanged by volumes at different pressure. The problem of how characterize the mass flow through the orifice and the source involve a deep

understanding of the gas-dynamic insight of the engine. Being this problem beyond the points of this work, only a preliminary investigation has been performed.

The mass flow through a valve is simulated by means the following equation [34]:

m˙ = sign · CdAνg2 √



1 —

× 1/γr 1—1/γr

1/2

P1 ≤ Pcrit

Tr  g2 Pr pr

P1 ≤ Pcrit

(3.7)

Pcrit = Pr



γr + 1

γr γr−1

t ≤ tresp

(3.8)



Aν = 4 D ×

g1 =

tresp

t > tresp

2γr

Rr(γr — 1)

1/2

(3.9)

(3.10)

γr−1

(3.11)

" γ 2 γr +1 #1/2

γr + 1

where Cd is the discharge coefficient, Aν is the area of the valve orifice, D its diameter. In case of pop-off valve the area changes with time. tresp is the net valve opening-time, different response time can be used for opening and closing operations. g1 and g2 are constants calculated using the physics characteristics γr and Rr of the gas flowing through the valve. γr is the adiabatic constant and Rr the gas constant, for the considered species.

Pa, Pb, Ta, Tb are respectively the pressure and the temperature of the gas on each side of the valve and, basically are the drivers of the mass flow through the orifice. The units are pascal for the pressure and xxxxxx for the temperatures. PR is maximum pressure between Pa and Pb, Tr its corresponding temperature, P1 the lower pressure between Pa and Pb. Pcrit is a critic pressure, defined as an edge parameter for two different gas-dynamic behavior through the valve, sonic for P1 < Pcrit, subsonic for P1 > Pcrit.

Sign is 1 depending on how the valve-mass-flow term is used in the equa- tions. It is +1 wherever this term is a source term, for example when the volume is connected with another volume at higher pressure behaving as a source. On the other side sign is -1 when this term is a sink term, as in the case of a volume connected with a volume at lower pressure than therefore is receiving mass flow. The source chamber and the vacuum chamber have been modelled assuming an isothermal gas behaviour, using particle balance equations with sink source terms and flow terms. The general form for the neutral particle balance equation is:

dni = —ΓL + Γs + Γ

valve

rec

+ Γ , (3.12)

ICRH

Helicon

Bmax

B2 B1

z2 z

Figure 3.2: Qualitative magnetic field axial profile.

where Γs is for the specie i the source term due to plasma processes, ΓL is the loss

i i

term due to plasma processes, Γvalve is the flow term through valves, Γrec is the volume recombination.

3.1.3 Magnetic mirror

In the second stage of the thruster the plasma is heated by an ICRH antenna. The magnetic field axial profile between the helicon antenna and the ICRH antenna is shaped as in Fig. 3.2. The peak of the profile acts as a magnetic mirror reflecting part of the ion flow. Furthermore this configuration changes the distribution func- tion of the plasma from z1 to z2. In order to obtain the correct plasma parameters at the immediate upstream of the ICRH, these two effects have to be considered.

The model considers the plasma parameters variations in the axial directions and neglects the variations in radial direction (1-dimensional model). This means that the plasma radial drift and diffusion are neglected. The area of the plasma flow channel changes along the axial coordinate z. From the conservation of the magnetic flux:

A(z) = A

(3.13)

1 B(z)

where A1 is the area and B1 is the magnetic field in z1.

The distribution function is expressed as a function of the parallel and normal components of the velocity. The integration of the distribution function f (v , v⊥) gives the plasma density n, the axial drift velocity u, the mean parallel and normal

temperatures T and T⊥:

∫

n = 2π ∞ f (v , v⊥)v⊥dv v⊥ (3.14)

—∞

u = 2π ∫ ∞ f (v , v

)v v dv v

(3.15)

n —∞

⊥ ⊥ ⊥

T = 2π ∫ ∞

1 m(v — u)2f (v , v )v

dv v

(3.16)

n —∞ 2

⊥ ⊥ ⊥

T = 2π ∫ ∞

1 mv2 f (v , v )v

dv v

(3.17)

⊥ n —∞ 2

⊥ ⊥ ⊥ ⊥

Since the distribution functions are generally not maxwellian, the temperatures are intended as mean energies of the particles motion.

3.1.4 Conditions on the particles reflection

The presence of a peak in the magnetic field profile creates the conditions of a partial particles reflection as in a magnetic mirror. For the conservation of the magnetic momentum µ = mv2 /(2B) and of the kinetic energy, only the particles

which satisfy v /v > (B ⊥/B — 1)1/2 pass beyond the maximum, the others

⊥

are reflected back [50].

max 1

Considering the distribution function of the plasma exhausted from the helicon as a maxwellian with a drift velocity, the particles reflected are the ones outside the cone in Fig. 3.3. The ratio between the "passing" flow and the total incident flow constitutes the chocking coefficient used in the modelling of the plasma discharge stage.

3.1.5 Plasma parameters at the ICRH

The plasma flow exhausted by the helicon is:

∫

Γ = A12π ∞ v

—∞

1f (v

1, v⊥1)v⊥1dv

1v⊥1 (3.18)

which can be written as

Γ = πA1 ∫ ∞ f (v , v

)dv2 v2

(3.19)

1 2 —∞ 1

⊥1 1 ⊥1

applying the conservation of the magnetic momentum and of the kinetic energy, we can write:

2 = v2

+ v2

⊥1

= v2

(1 — r) (3.20)

r (3.21)

where r = B2/B1.

⊥2 ⊥1

v⊥

reflected

passing

v//

Figure 3.3: Maxwellian distribution function with a drift velocity parallel to v . The particles outside the dash line are reflected by the effect of the peak of the magnetic field.

Changing the coordinates from the system (v2 , v2 ) to the system (v2 , v2 ),

1 ⊥1 2 ⊥2

whose Jacobian is 1/r (i.e. dv2 dv2 = (1/r)dv2 dv2 1/r), we have:

1 ⊥1

Γ = πA1 ∫ ∞ f (v (v , v

2 ⊥2

), v (v , v

))dv2 v2

(3.22)

1 2 —∞ 1 2 ⊥2

⊥1 2

⊥2 2 ⊥2

If g(v 2, v⊥2) is the distribution function in the section just before the ICRH, the particles flow is:

Γ = πA1 ∫ ∞ g(v , v

)dv2 v2

(3.23)

2 2 —∞ 2

⊥2 2 ⊥2

Under the hypothesis that no particle is reflected by the magnetic mirror, the flows in the two sections are equals, namely Γ 2 = Γ 1. Considering that A1/A2 = r, from the two previous equations we have:

g(v 2, v⊥2) = f (v 1(v 2, v⊥2), v⊥1(v 2, v⊥2)) (3.24)

(v 1 — u1)2

—

⊥1

(3.25)

Considering the distribution function at the exit of the helicon as a maxwellian with a drift, to wit:

f (v 1, v⊥1) = (2

)

exp

—

n1 "

π 3/2v3

v2 #

2v2

Therefore at the ICRH the distribution function is:

 q(v2

— v2

(1 — r)/r — u1)2 2 

g(v

2, v⊥2

) = n1

(2π)3/2v3

exp —

2 ⊥2

—

2v2

v⊥2

2rv2

(3.26)

Some particles are reflected due to the relation v 1/v⊥1 > (R — 1)1/2. If we apply

the change of coordinates from (v2 , v2 ) to (v2 , v2 ), we obtain:

2 ⊥2 1 ⊥1

2 = v2

+ v2

⊥2

1 — r r

(3.27)

v2 = v2

⊥1 ⊥2

r (3.28)

and thus we have to integrate with the condition:

g(v

2, v⊥2

) = 0 when v 2

< v⊥2

R — r . (3.29)

Integrating the distribution function g(v 2, v⊥2) (cf. Fig.3.4), calculated in this way, we can obtain n2 and u2.

3.1.6 ICRH effect

The ICRH is expected to heat the plasma to ion temperatures the order of 100 eV. Since the helicon heats the plasma to ion temperatures less than 1 eV, it seems rea- sonable to consider cold (Ti = 0 eV) the plasma arriving at the ICRH. Therefore, it has been assumed that the plasma flow entering the ICRH region is monoenergetic with an axial velocity equal to u2.

The power deposited by the ICRH is absorbed by the plasma in the form of

normal kinetic energy:

1 mv2

= PABS

. (3.30)

3.1.7 Magnetic nozzle

2 ⊥2

n2u2A2

Finally, downstream the ICRH, the magnetic nozzle converts part of the normal kinetic energy of the plasma into parallel kinetic energy, and thus produces thrust. This conversion is effective until the plasma detaches from the magnetic field lines. We assume that the detachment happens when the parallel kinetic energy density (mnv2/2) becomes greater than the magnetic energy density (B2/2µ0).

Applying the conservation laws of energy, magnetic momentum, particles flow,

and magnetic flux, we have:

v2 + v2

⊥

= const. (3.31)

v2 /B = const. (3.32)

⊥

Figure 3.4: Example of distribution function at the ICRH. A fraction of the distri- bution is cut by the presence of a peak in the magnetic field.

nv A = const. (3.33)

BA = const. (3.34)

from the last equation we have the area of the flow channel along the axial coordi- nate z:

A(z) = A

(3.35)

from the second, the normal velocity:

v2 (z) = v2

2 B(z)

= B(z)

(3.36)

⊥ ⊥2 B2

thus, from the first equation we have the parallel velocity:

v2(z) = u2 + v2

1 — B(z) (3.37)

and from the third the density:

n(z) =

⊥2 B2

n2u2B(z) v (z)B2

(3.38)

the axial position where the plasma detaches comes from the solution of the fol- lowing equation:

B2(z∗) = 1

( ∗)

2( ∗) (3.39)

2µ0

2 mn z v z

it is possible now to compute the thrust:

T = mionn(z∗)v2(z∗)A(z∗) (3.40)

and rearranging:

ion

T = m

n u A su2 + v2

1 — B(z∗) (3.41)

3.2 Model validation

Isp =

⊥2

m˙ g0

(3.42)

The numerical outputs of the model described above were compared with experi- mental data found in literature [48] about plasma parameters in the helicon stage. This operation made it possible to validate the model and to determine the value of few parameters. Figure 3.5 shows the comparison between the model results and experimental data found in literature.

Figure 3.5: Comparison between experimental data [48] and the numerical model: helium discharge, 3 kW RF power.

3.3 Model optimisation

An optimisation technique is required in order to find the best performance attain- able by the thruster. The model takes as inputs many parameters, including: total power, power repartition between the two RF antennas, propellant mass flow rate, magnetic field configuration and thruster dimensions. An evolutionary algorithm was used in order to find the combination of parameters that gives the best perfor- xxxxx in terms of thrust efficiency and specific impulse.

Chapter 4

Numerical results

4.1 RF modelling of the ICRH antenna

The numerical approach outlined in Section 2 has been implemented in an ex- tended version of TOPICA code, so that, once the current densities X X, M A— have been computed, the antenna admittance matrix [Y ] along with other common parameters [1] can be obtained. The code has been thoroughly tested in vacuum, while preliminary results with plasma show good agreement with data available in literature.

To give an example of TOPICA’s capabilities, we dealt with a neutral electron- helium (4He+) plasma and a two-loop ICRH antenna as in Fig. 2.1: listed in Table

4.1 are the main parameters involved in the simulation. The density profile has been assumed parabolic and given analytically, even though it may be read from file, if desired.

An important intermediate step towards the solution is the calculation of the Green’s function Y˜P (m, kz), which accounts for all of the plasma effects. As men- tioned in Section 2, this is accomplished upon solving the Xxxxxxx’x equations within the plasma; the finite-element numerical scheme appears more stable if the

magnetic field at the air-plasma interface is enforced and then the electric field is computed therein. In consequence, a tensor impedance Z˜P being the natural out- put, Y˜P (m, kz) is obtained by inversion, in accordance with (2.14).

Extensive numerical experiments have shown that, owing to the average size

—

(the order of cm) of the triangular facets forming the aperture mesh, it suffices to compute Y˜P (m, kz) for m within the range from 60 to 60 (i.e. a grand total of 121 azimuthal modes) or even less, in order for the spectral integrals to converge [1]. For the sake of completeness, Fig. 4.1 and 4.2 show the real and imaginary part of the entries of Z˜P germane to the case under study. The tensor is symmetric, as theory predicts, but each of the components lacks any symmetry along m, because of the confining magnetic field aligned with ˆz—which is indeed suggestive of the non-reciprocal (gyrotropic) nature of the plasma.

We considered the standard counter-driven two-loop antenna shown in Figs. 4.3

−5

−10

−5

−10

Re Zθ,θ

−2000 0 2000

k /k

z 0

Re Z

−2000 0 2000

k /k

z 0

0.035

0.03

0.025

0.02

0.015

0.01

0.005

x 10−4

−2

−4

−6

−5

−10

−5

−10

Re Zθ,z

−2000 0 2000

k /k

z 0

Re Z

z,z

−2000 0 2000

k /k

z 0

x 10−4

−2

−4

−6

x 10−5

Figure 4.1: Real part of Z˜/Z0 entries as a function of m and kz/k0, with Z0 (k0) the free space impedance (wavenumber).

−5

−10

−5

−10

Im Zθ,θ

−2000 0 2000

k /k

z 0

Im Zz,θ

−2000 0 2000

k /k

z 0

−0.01

−0.02

−0.03

−0.04

−0.05

−0.06

x 10−3

−1

−2

−3

−5

−10

−5

−10

Im Zθ,z

−2000 0 2000

k /k

z 0

Im Zz,z

−2000 0 2000

k /k

z 0

x 10−3

−1

−2

−3

x 10−4

−1

−2

−3

Figure 4.2: Imaginary part of Z˜/Z0 entries as a function of m and kz/k0, with Z0

(k0) the free space impedance (wavenumber).

Figure 4.3: Standard counter-driven two-loop antenna: sample electric current magnitude distribution on conducting bodies and at plasma/air interface.

Figure 4.4: Standard counter-driven two-loop antenna: sample magnetic current magnitude distribution at plasma/air interface.

Table 4.1: Geometrical and physical data, single species (4He+) plasma

Loop radius (m)	0.021
Loop distance (m)	0.02
Plasma radius (m)	0.02
Density at center (m—3)	5.55 1019
Density at edge (m—3)	0
Axial magnetic field (T)	0.25
Parallel velocity (m/s)	8 103
Ion perp./par. temp. (eV)	2
Electron perp./par. temp. (eV)	4

—

— —

and Fig. 4.4, wherein colors refer to the magnitude of the electric and magnetic sur- face current densities over the conducting bodies at the air-plasma interface, when the antenna terminals are driven by {+V , V , +jV , jV ,}, respectively, since we want the whole antenna to launch an electric field with the circular polarization that maximizes the wave coupling to the m = 1 mode [11]. The electric cur- rent appears almost constant over the antenna loops, for the vacuum wavelength is pretty larger than the structure size. Besides, the magnetic current magnitude (that is the tangential electric field) is quite strong under the loops, which agrees with the expectation.

The main TOPICA’s output is the antenna 4-by-4 admittance matrix [Y ], which in this instance may be further reduced to a 2-by-2 matrix, since each loop is driven with opposite voltages. In particular, upon enforcing the voltages at the antenna ports, we can compute the radiated power via [Y ] and use it to evaluate the coupling resistance, according to (2.74).

However, due to the intrinsic limitations posed by the plasma model adopted so far, namely, the occurrence of possibly non physical entries of the plasma con- ductivity tensor (Section 2.4), for the time being we resorted to the intermediate approach outlined in Section 2.6.3 as per the calculation Rc.

As an example, with the data reported in Table 4.1, we obtained Rc for the two-

loop antenna under investigation as a function of the frequency; results are plotted in Fig. 4.5, where the parameter of the curve is the ratio of the loop width to the plasma radius. For the sake of comparison, the same graph also shows measured and simulated results after [11]. As can be seen the plasma loading predicted by TOPICA not only is well in the ballpark but also compares favorably with the measurements, and especially so at frequencies above the ion cyclotron resonance. In addition, the dashed lines in Fig. 4.5 reportedly pertain to different plasma flow profiles and were obtained through the approach described in [11], using just the same expression for the plasma tensor entries as given in (2.41)-(2.43). Nonetheless, as outlined in Section 2.5, the two numerical solutions differ from each other as for the way they use the spectral plasma tensor and this help explain

why the simulated results are slightly different.

Simulated (after [2])

Measured (after [2])

w/a = 0.6

w/a = 0.8

w/a = 1

w/a = 1.2

w/a = 1.4

0.25

0.2

Plasma loading [Ω]

0.15

0.1

0.05

0.9 1

ω/ω

1.1

Figure 4.5: Standard counter-driven two-loop antenna: plasma loading as a func- tion of the frequency computed with TOPICA for different loop-width to plasma- radius ratios. Also superimposed are measured data and simulations published in the work by Xxxx [11].

In the end, we cursorily mention that to carry out our simulations we resorted to velocity values typical for application to RF thrusters, for in [11] we could not find any clue on the adopted velocities of ions and electrons. In particular, as can be seen in Table 4.1, we chose coincident ion and electron flow velocities, in accordance with [12], which is a work companion to [11], even though, to say the truth, this assumption does not appear enough justified in general, the masses of ions and electrons being significantly different.

4.2 Plasma device

Figures 4.6 and 4.7 show the results obtained from the model respectively with helium and argon in operative conditions. The plots show the density of the species considered and the electron temperature time evolution in the helicon stage. These values represent the estimate of the plasma parameters that will occur during the experiment. Therefore they give an indication of the requirements for the plasma probes needed.

electr on temperature [eV]

0 0.005 0.01 0.015 0.02 0.025 0.03

He e-

density [m-3]

0 0.005 0.01 0.015 0.02 0.025 0.03

Time [s]

Figure 4.6: Electron temperature and density of different species during the plasma discharge. Gas: helium, absorbed power on the helicon stage: 4,000 W, injected mass flow rate: 2 10—6 kg/s.

Plasma densities from 1017 m—3 to 1020 m—3 are obtained using the two gases. Electron temperatures up to 20 eV are obtained. The time evolution shows that the plasma probes should have a frequency pass-band from 0 to 1 kHz. The stationary conditions are reached after few cents of second. This may be a useful parameter for the design of the experiment that could last less than one second.

As expected, at constant mass flow rate, heavy gases are easier to ionize than lighter gases, simply because their flow has less particles per second be to ionized. Therefore, less power is used on the helicon antenna when argon is used than when hydrogen is used. Considering the entire thruster, which shares a constant total amount of power between the ionizing stage and the heating stage, heavy gases seems appropriate for low specific impulse ranges and light gases for high specific impulse ranges where the acceleration of a heavy gas may be difficult.

The plasma discharge model provides the density of the neutral and ionized species, and the electron temperature at the exit of the helicon antenna stage. The equations written so far make it possible to follow the evolution of the plasma along the axis of the thruster and to obtain the fundamental parameters of an elec- tric thruster: specific impulse, thrust efficiency and thrust. Fig. 4.8 shows the thrust performance of a 10 kW thruster obtained from the optimisation by using the evo-

electr on temperature [eV]

0 0.02 0.04 0.06 0.08 0.1

density [m-3]

0 0.02 0.04 0.06 0.08 0.1

Time [s]

Figure 4.7: Electron temperature and density of different species during the plasma discharge. Gas: argon, absorbed power on the helicon stage: 1,000 W, injected mass flow rate: 4 10—6 kg/s.

lutionary algorithms. Analogous simulations made for input powers ranging from 1 kW to 100 kW gave results very similar to those shown. This means that the thruster’s efficiency seems not to depend on the input power.

The thrust performance increases with the specific impulse basically because less power is needed to ionise the plasma and hence more power is diverted on the ICRH antenna that boosts the plasma. Figs. 4.9 and 4.10 show how control parameters like power repartition and propellant mass flow vary in function of the specific impulse. It is interesting to note that the high specific impulse is achieved by directing more power to the ICRH (taken from the helicon) and by reducing the propellant flow rate.

Although the results obtained are susceptible of uncertainty, the relative per- formance indicates a preference in the use of heavy gases and a behaviour of the thruster which privileges high specific impulse. Light gases become interesting only at very high specific impulses (> 15,000 s), where it may be difficult to ac- celerate a heavy gas. A drawback in the use of heavy gases is the necessity of producing more intense magnetic field in order to confine the plasma. Indeed, the radius of gyration of a magnetized plasma (Larmor radius) is proportional to the

Argon Helium

0.7

0.6

0.5

thrust efficiency

0.4

0.3

0.2

0.1

0 2000 4000 6000 8000 10000 12000

specific impulse [s]

Figure 4.8: Thrust efficiency as function of specific impulse using helium and argon.

atomic/molecular mass of the element considered and inverse proportional to the magnetic field applied. Therefore, the use of argon may need magnetic fields of more than 2 T. This is a serious challenge for the superconductor technology and for the thermal control. On the other side, tanks for light gases are very heavy and may force the use of cryogenic storage systems. This solution is not necessary with heavy gases which can be easily contained in gaseous form.

The subsystems analysis makes it possible to evaluate the mass and dimensions of the thruster. With a power of 10 kW, the mass of the whole thruster would be 25 kg using helium and 37 kg using argon, including structure (He: 4 kg, Ar: 10 kg), superconducting magnets (He: 4 kg, Ar: 10 kg), thermal insulators (1 kg), antennae and equipment (16 kg). Besides that, the masses of propellant and tanks should be considered. The thruster would have a cylindrical shape, 800 mm length and 300 mm diameter.

4.3 Scalability criteria and results

Two factors induce to believe that the use of this technology at high power is very promising:

thrust efficiency

0.5

0 2000 4000 6000 8000 10000 12000

power repartition

0.5

propellant flow [kg/s]

0 2000 4000 6000 8000 10000 12000

-5

x 10

0 2000 4000 6000 8000 10000 12000

specific impulse [s]

Figure 4.9: Thrust parameters given by the optimization using helium as propellant (power repartition: 0 means full power on the ICRH, 1 means full power to the helicon).

•

the absence of electrodes excludes any problem related to cathode erosion which in other electric thrusters is a limit in terms of life and attainable power;

•

the absence of space-charge distortions in the applied electric field, which in electrostatic thrusters constitutes a limit in the thrust to nozzle area ratio.

Basically, two strategies can be followed when scaling up the power: the pro- xxxxxxx flow can be increased proportionally with power, or it may remain constant. In the former case the specific impulse stays roughly constant because the energy imparted to each ion does not change. The thrust increases proportionally with power because the propellant flow increases. The confining magnetic field remains unchanged because the ion temperature is constant. If the diameter of the plasma tube is kept constant, plasma density scales up. This affects the antennas’ effi- ciency most likely in a negative way, therefore it may be useful to increase the radial dimension of the thruster.

Otherwise, if the propellant flow is kept constant, the specific impulse, and thus the thrust, scales up because more energy is imparted to each ion. Detachment efficiency scales favorably with specific impulse. The confining magnetic field needs to be increased in order to limit the ion’s Larmor radius because the ion temperature increases.

It is possible to analyze how the mass of the thruster scales with power. The mass of the superconducting magnets is independent from the power if the mag-

thrust efficiency

0.5

2000 2500 3000 3500 4000 4500 5000 5500 6000

power repartition

0.4

0.3

0.2

propellant flow [kg/s]

2000 2500 3000 3500 4000 4500 5000 5500 6000

-5

x 10

2000 2500 3000 3500 4000 4500 5000 5500 6000

specific impulse [s]

Figure 4.10: Thrust parameters given by the optimization using argon as propellant (power repartition: 0 means full power on the ICRH, 1 means full power to the helicon).

netic field intensity does not change, and so it is the mass of the structure that has to support them. The mass of the antennae and equipment can instead be considered proportional to the power. These considerations indicate that there is an inferior limit in power, below which the ratio thruster mass over power is inconveniently high (see Figs. 4.11-4.12).

Ar He

200

180

160

140

Thruster mass [kg]

120

100

0 2 4 6 8 10

Power [W ]

x 10

Figure 4.11: Thruster mass versus input power using argon and helium.

Ar He

0.02

0.018

0.016

Thruster mass/power [kg/W]

0.014

0.012

0.01

0.008

0.006

0.004

0.002

0 2 4 6 8 10

Power [W ]

x 10

Figure 4.12: Specific mass versus input power using argon and helium.

Chapter 5

Conclusions

According to the objectives of this activity, we have developed and described a methodology to address the design of RF-based plasma thrusters.

We consider a three-stage thruster with magnetic transversal confinement. It is composed of a generation unit, a helicon plasma source; a booster section employ- ing radio frequency (RF) power at the ion cyclotron (IC) frequency (ion-cyclotron resonant "heating", ICRH); and finally a magnetic nozzle.

Our approach relies on both the self-consistent TOPICA code for modeling the RF behavior of the ICRH antenna, and a fast zero-dimensional numerical model purposefully developed for the global description of the plasma-device. TOPICA has been modified to handle the specific geometry and plasma of the thruster; it allows the evaluation of the ICRH antenna parameters (strictly related to the an- xxxxx shape and plasma composition) and the ensuing radiated power. The zero- dimensional global model provides the relevant global system parameters, neces- sary to assess the engine performance.

TOPICA results for the parameters of the ICRH antenna of the booster unit were validated against published measured data (for the so-called coupling resis- tance). The predictions of the global modelling were likewise validated against published experimental values.

The developed tool set as a whole can be used, for instance, to investigate the feasibility of a low-power (a few kW) system and then to address its scaling to high power. Because of the reduced computational load of the zero-dimensional model, evolutionary algorithms can be employed to optimize the model.

With this approach, optimal values of power partition, mass flow rate, type of propellant, magnetic field and dimensions of the thruster have been be computed as an example. The results of the model give the indication of the values of plasma pa- rameters that should be expected during the experiment, and also the requirements for the needed diagnostics. In this respect, the carried out prototypical optimization exercise confirms that the use of Argon as propellant is preferable to other lighter gases, mainly because it can be ionized with significantly lower power.

Bibliography

[1] X. Xxxxxxxxxxx, X. Xxxxxxxxx, X. Xxxxxxxx, X. Xxxxxx, X. Kyrytsya 2006

Nucl. Fusion 46 S476

[2] F. R. Xxxxx-Xxxx 2000 Sci. Am. 283 90

[3] I. V. Xxxxxxx 1992 Methods for Electromagnetic Field Analysis Oxford Engi- neering Science Series (Oxford: Clarendon Press)

[4] X. X. Xxxxxxxx, X. X. Xxx, X. Xxxxxx 1998 Computational Methods for Elec- tromagnetics (New York: IEEE Press)

[5] X. Xxxxxxxxxx 1961 Time-harmonic Electromagnetic fields (New York: Wi- ley) 2001

[6] X. Xxxxxxxxxx 1993 Field Computation by Moment Methods (New York: Ox- ford)

[7] X. X. Xxxxxx 2001 Foundations for Microwave Engineering (New York: IEEE Press)

[8] X. Xxxxxxxxx 1951 Waveguide Handbook (New York: XxXxxx-Xxxx)

[9] X. Xxxxxxxxx 1998 Kinetic theory of plasma waves (Oxford: Clarendon Press)

[10] E. Bearing, X. Xxxxx-Xxxx, X. Xxxxxx 2004 The Radio Science Bulletin, No. 310

[11] X. X. Xxxx et al. 2004 Proc. 42 AIAA, No. 0151

[12] X. Xxxx et al. 2005 Proc. 43 AIAA, No. 0949

[13] X. X. Xxxx 1992 Waves in Plasmas, (New York: American Institute of Physics)

[14] X. Xxxxxxxxx, X. X. Xxxxxxx 1983 Nuclear Fusion 23

[15] Pavarin D., “Diagnostic and optimization procedures for gas and plasma propulsion system”, Ph.D. thesis, University of Padova, 2001

[16] Xxxxxxxx X. X., Xxxxxxx D., Xxxxxxxx F., Xxxxxx X. X., Xxxxxx M. D., Xxxxxxxx X., Xxxxxx D. O., “Helicon plasma source configuration analysis by means of density measurements”, Proceedings of the ICEAA, Torino, Italy, 1999

[17] Xxxxx S., “Design and modeling of an experiment for the performance evalu- ation of a variable specific impulse magnetoplasmadynamic thruster”, Ph.D. thesis in Mechanical Measurements for Engineering, University of Padova, 2005

[18] Xxxxxxxx X. X., Xxxxx X. X., Xxxxxx X. X., Xxxxxx M. D., Xxxxx-Xxxx F. R., Xxxxxxx D., Xxxxxx D. O., Xxxxxx X. X., “Helicon plasma source optimization studies for VASIMR”, Meeting APS 1999

[19] Suwon C., “A self consistent global model of neutral gas depletion in pulsed helicon plasmas”, Physics of plasmas, Vol. 6, January 1996, pp. 359-365

[20] Xxxxxxx J., “Neutral Pumping in a Helicon Discharge”, Plasma Sources Sci.

Techn. No. 7, 1998, pp. 416-422

[21] Xxxxx X., “Global model of a radio frequency H2 plasma in XXXXXX”, Plasma Sources Sci. Techn. No. 9, 2000, pp. 161-168

[22] Xxxxx X. X., “Global model and scaling laws for inductively coupled oxygen discharge plasmas”, J. Appl. Phys. Vol. 86, No. 7, October 1999

[23] Xxx X., X. Electrochem. Soc. Vol. 14-1, No. 6, June 1994

[24] Xxxxxxxx A. D., “Characterization and modelling of a helicon plasma source”, X. Vac. Sci. Techn. A Vol. 16, No. 5, September-October 1998

[25] Xxxxxxxxx M. A., “Global Model of pulse-power-modulate high-density, low-pressure discharges”, Plasma Sources Sci. Techn. Vol. 5, 1996, pp. 145- 158

[26] Xxx X., “Global Model of Plasma Chemistry in a High Density Oxygen Dis- charge”, J. Electrochem. Soc., Vol. 141, No. 6, June 1994

[27] Ashida S. , “Spatially averaged (Global) model of time modulated high den- sity chlorine Plasmas, J. Appl. Phys. Vol. 36, 1997, pp. 854-861

[28] Xxxxxx X. , “Spatially averaged (global) model of time modulated high den- sity argon plasmas”, X. Vac. Sci. Techn. A Vol. 13, No. 5, September-October 1995

[29] Xxxxxxxxx X., “A spatially-averaged model for high density discharges”,

Vacuum, Vol. 47, No. 3, pp. 215-220, 1996

[30] Gordies B., “Self-consistent kinetic model of low-pressure N2-H2 flowing discharges I-II”, Plasma sources Sci. Techn., Vol. 7, 1998, pp. 378-388

[31] Xxx X., “Global model of Ar, O2, Cl2, and Ar/O2 high-high density plasma discharges”, J. Vac. Sci. Techn. A Vol. 13, No. 2, March-April 1995

[32] Xxxxxx T.,“Ion and neutral temperature in electron cyclotron resonance plasma reactors”, Appl. Phys. Lett. Vol. 58, No. 5, 1991, pp. 458-460

[33] Xxxxxxx J. ,“Neutral gas temperature in a multipolar electron cyclotron res- onance plasmas”, Appl. Phys. Lett. VOl. 58, No. 22, p. 3, June 1991

[34] Xxxxxx X. X.,“Quickgun: an algorithm for estimating the performance of two- stage light gas guns”, ORNL/TM-11561 report

[35] Xxxxxxxx X., Partially Ionized Xxxxx, X. Xxxxx, 1973, New York.

[36] XxXxxxxx X. X., Collision phenomena in ionized gases, Xxxx Xxxxx & Sons,

New York, 1964

[37] Xxxxx X. X., Elementary processing hydrogen helium plasmas, Springer- Verlag, 1987.

[38] Barnet C. F., “Collision of H, H2, He and Li atoms and ions with atoms and molecules”, ORNL report

[39] Godyak V. A., 1986

[40] Mozetic M., “Atomic hydrogen density along continuously pumped glass tube”, Vacuum, Vol. 5, No. 3-4 pp. 319-322, 1998

[41] Xxxxx X. X. “Neutral depletion and transport mechanism in large-area high density plasma sources”, J. Appl. Phys., Vol 86, No. 10, November 1999

[42] Xxxxxxxxx M. A., Principles of Plasma Discharges and Material Processing,

New York Wiley, 1994

[43] Xxx X. X. Electroch. Soc, Vol 141, No. 6 June 1994

[44] Dushmann, Scientific foundation of vacuum technique, Xxxx Xxxxx & Sons,

New York 1964

[45] Holland L., Vacuum Manual, F. N. Spon, London

[46] Xxxxxxx, The mathematical theory of non uniform gas, Cambridge At the University Press, 1970

[47] Xxxxxxxx M. I., “Characterization of the Resonant Electromagnetic Mode in Helicon Discharges” Ph. D. Thesis, University of Texas at Austin, 2003

[48] Xxxxxx X.X. et al. “Experimental research progress toward the VASIMR en- gine”, Proc. of the International Electric Propulsion conference Toulouse, 2003

[49] Hochstim A. R., Kinetic processes in gases and plasmas, Academic Press,

London, NY, 1969

[50] Xxxxxxxx X. X., Xxxxxxxxxx P. H., Introduction to Plasma Physics, pp. 39-40

Document Metadata

Table of Contents

Dipartimento di Elettronica, Politecnico di Torino, Torino, I-10129, Italy

Document Metadata