Preprint
Article

Solution of Orifice Hollow Cathode Plasma Model Equations by Means of Particle Swarm Optimization

Altmetrics

Downloads

143

Views

45

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

27 June 2024

Posted:

27 June 2024

You are already at the latest version

Alerts
Abstract
Orificed Hollow Cathodes are electric devices necessary for the functioning of common plasma thrusters for space applications. From their reliability mainly depends the success of a spacecraft’s mission equipped by electric propulsion. The development of plasma models is crucial in the evaluation of plasma properties within the cathodes that are difficult to measure due to the small dimensions. Many models, based on non-linear system of plasma equations, have been proposed in the open literature. These are solved commonly by means of iterative procedures. This paper investigates the possibility to solve them by means of Particle Swarm Optimization method. The results of the validation tests confirm the expected trends for all the unknowns; the confidence bound of the discharge current as function of mass flow rate is very narrow (2÷5V), moreover the results match very well the experimental data except at the lowest mass flow rate (0.08mg/s) and discharge current (1A), where the computations underpredict the discharge current to the utmost by 40%; the highest data dispersion regards the plasma density in the emitter region (± 20 % of the average value) and the wall temperatures (± 50K respect to the average values) of orifice and insert; those of the others variables are very tiny.
Keywords: 
Subject: Engineering  -   Aerospace Engineering

1. Introduction

Orifice Hollow Cathodes (OHCs) [1] are used both as electronic source and/or neutralizer for Electric Propulsion (EP) devices e.g., Hall Effect Thrusters (HETs), Gridded Ion Thrusters (GITs), High Efficiency Multistage Plasma Thrusters (HEMPTs) [2]. A lot of computational tools, based on phenomenological models, also known as zero-dimensional (0D) or reduced-order models, aimed to predict OHCs’ plasma properties and performances and useful to support their design, have been developed since ‘80s’ [3,4,5,6,7,8,9,10,11,12]. The most recent contributions (last five years) in this field of research have particularly dealt with: the proposal of original model based on scaling relationship of the total pressure inside the OHCs [8]; the improvement of the existing models by means of the inclusion of such aspects neglected in past works [9]; the detailed plasma modeling within heater-less cathodes and coupling with thermal model [11,12].
Also the Italian Aerospace Research Centre (CIRA) has developed its own design and analysis tools for xenon fed OHCs (named PlasMCat) [13], that is based on Albertoni’s et al. model [6], whose plasma equations are corrected according to the suggestion by Wordingham et al. [14] and cathode’s internal pressure is computed by means of the empirical formula proposed by Taunay et al. [15]. Moreover a proper plasma model is included for the cathode-keeper gap.
The equations that constitutes the most of the models presented in the open literature [5,6,7,9,10,11,12] are commonly iteratively solved to face with their own non-linear characteristic. Each author adopts a different approach. PlasMCat follows a procedure similar to Albertoni et al. [6] and in agreement with them, assessed that the adopted solution technique has limited stability because abrupt changes in geometry, flow rate or discharge current can cause the solution to diverge. This has given motivation to explore different solution strategies.
In this work the possibility to solve the whole system of models’ equations by means of a Particle Swarm Optimization (PSO) algorithm is investigated. PSO is a population based stochastic optimization technique developed by Eberhart and Kennedy in 1995 [16], inspired by social behaviour of bird flocking or fish schooling. PSO has been successfully applied in many research and application areas [17] and recently it has also been successfully used in EP application to solve the complex systems of non-linear equations used to model Helicon Plasma Thrusters’ (HPTs) plasma [18].
Hence, because of the relevant number of involved variables and to the associated strong non-linearity, in order to benefit from this approach, a code that searches for all the possible solutions of the plasma equations for OHCs has been built and tested.
The article is structured as follows: the main definitions and assumptions of the plasma model (whose equations are recalled in Section 2.3 for completeness) are discussed in Section 2; the optimization technique and its detailed model are presented in Section 3; an application of the methodology is presented in Section 4 where the results of the simulation are discussed; finally, the conclusions are drawn in Section 5.

2. Plasma Model of Hollow Cathode

2.1. Cathode Architecture

In the typical OHCs’configuration (Figure 1) a conductive tube, made of refractory metal (e.g., tantalum or molybdenum) or graphite, drives the flow of a neutral gas (e.g., xenon, argon but also other non noble-gases). At the end of the main tube an orifice plate is usually placed. By changing the diameter of the orifice the internal pressure can be adjusted [1]. An hollow cylinder, namely insert, made of low-work function material, is placed inside the tube. The insert is the core part of the cathode because it emits electrons by means of thermionic effect. Barium-based and rare-earth emitters (e.g., LaB6, CeB6), which have respectively a work function of about 2 eV and 2.7 eV [1], allow to provide current densities as high as 105 A/m−2 at a surface temperature of about 1300 K and 1900 K, respectively. The main drawbacks for Barium-based emitters is that they are highly susceptible to oxygen contamination whereas rare-earth emitters require very high temperature to work and care must be taken when they are in contact with many refractory metals at high temperatures because they react. For this kind of coupling, a thin graphite sleeves is usually placed between the emitter and the main tube.
An heater is wrapped around the tube to trigger thermionic emission. If well design, a cathode can sustain the emitting temperature by means of plasma heating thus heater can be turned off. For this reason, thermal management is a key factor in OHCs (e.g., thermal shields around the heater are often used to reduce the radiation losses).
Another important element is the keeper electrode that encloses the heater-main tube assembly. It has the double function of to help the extraction of the electrons to the discharge plasma and to protect the cathode from the ion bombardment.

2.2. Definition and Assumptions

The hollow cathode internal plasma volume is divided into three regions, the keeper region, pedix "k" (which includes the keeper orifice region and the cathode-keeper gap), the insert or emitting region, pedix "e", and the orifice region, pedix "o", as depicted in Figure 1. A set of conservation equations is written for each region (Section 2.3) to compute averaged plasma properties (0-D model). The modelled plasma consists of neutrals, singly-charged ions and electrons. The unknown parameters are the electron temperatures ( T e , e , T e , o , T e , k ), the plasma densities ( n e , e , n e , o , n e , k ), the neutral densities ( n n , e , n n , o , n n , k ), the effective emission length ( L e f f ), the sheath potential ( V p ), the wall temperature of the emitter ( T s ), orifice ( T w , o ) and keeper ( T w , k ).
The main assumptions, in agreement with the norm for 0D models [5,6,9], are:
  • the heavy particles (ions and neutrals) are in a thermal equilibrium between each other and their temperature is assumed to be equal to the temperature of the wall, i.e., T i o n = T n = T w ;
  • a chocked gas flow model is used in the orifice regions of both cathode tube and keeper;
  • a double sheath potential drop is recorded at the orifice entrance region, thus a planar double sheath is modelled (Figure 1).
The neutral flow through the cathode’s insert, that is transitional, is not modelled as continuous and governed by the Poiseuille’s law, which is a well established assumption for 0D models [14]. Instead, the cathode internal pressure is computed by means of Taunay et al. formula [15].
For all the simulations described in this paper, xenon is used as propellant because it is widely used in EP devices. Particularly, the electron-neutral cross section has been updated in order to include only elastic collisions as suggested by Wordingham et al. [14] and to take into account the most recent experimental data, a fit of which is proposed in Ref. [13].

2.3. Plasma Model Equations

The following subsections recall the equations presented in Ref. [13] and used to model the plasma in the three regions in which cathode’s internal volume has been divided.

Keeper

In the keeper region the unknowns are the plasma number density ( n e , k ), electron temperature ( T e , k ) and neutral density ( n n , k ) and the equations are the current density balance, Equation 1(a), the balance of ion flux, Equation 1(b) and the equation of state, Equation 1(c).
n e , k = n e , o A o A k T e , o T e , k ( a ) q π r k 2 L k n n , k n e , k K i = j i , k 2 π r k L k + 2 j t h , k π r k 2 ( b ) n e , k + n n , k k B T w , k 1 + α T e , k T w , k = m ˙ π r k 2 R g T w , k γ 1 + α T e , k T w , k ( c )

Orifice

The system of equations for the cathode’s orifice region writes a balance for the ion flux Equation 2(a) and for the power Equation 2(b), and states that pressure calculated by means of kinetic theory is equal to that computed with continuum flow and chocked orifice, Equation 2(c). The plasma number density ( n e , o ), the electron temperature ( T e , o ), and the density of the neutrals ( n n , o ) in the orifice region are the unknowns.
q π r o 2 L o n n , o n e , o K i n ˙ i o n = j i , o 2 π r o L o + π r o 2 1 1 4 1 0.61 8 π + j t h , o π r o 2 ( a ) R I d 2 = n ˙ i o n ϵ i + 5 2 k B q I d T e , o T e , e ( b ) n e , o + n n , o k B T w , o 1 + α T e , o T w , o = m ˙ π r o 2 R g T w , o γ 1 + α T e , o T w , o ( c )
Here the plasma resistance is calculated as:
R = η L o π r o 2
The plasma resistivity is generally computed as:
η = ( ν e i + ν e n ) m e n e q 2
and collision frequencies are:
ν e i = 2.9 · 10 12 n e l n Λ T e k B q 3 2
v e n = σ e n n n k B T e n e
where the Coulomb logarithm is:
l n Λ = 23 1 2 log 10 6 n e T e k B q 3
and the fit for the Maxwellian average elastic electron-neutral cross section ( [ m 2 ] ) proposed in Ref. [13] for xenon is:
L o g 10 ( σ e n ) = p 1 T e * 5 + p 2 T e * 4 + p 3 T e * 3 + p 4 T e * 2 + p 5 T e * + p 6 T e * 2 + q 1 T e * + q 2 , T e * = L o g 10 T e k B q
p 1 p 2 p 3 p 4 p 5 p 6 q 1 q 2
Values 0.0350 0.3057 0.6698 18.5416 12.8821 9.9614 0.7093 0.5227

Insert

Four equations are written in order to get the unknown plasma number density ( n e , e ), electron temperature ( T e , e ), neutral density ( n n , e ) and emitter sheath voltage drop ( V p ) in the insert region. These are the balance for the ion flux, Equation 9(a), the current conservation, Equation 9(b) the power balance, Equation 9(c) and the continuity equations, Equation 9(d).
q π r e 2 L e f f n n , e n e , e K i + j i , o π r o 2 = j i , e A e f f + π r e 2 π r o 2 + j t h , e π r e 2 n ˙ o u t , e ( a ) I d A e f f = j i , e + j e m j e r ( b ) j i , o V d s + 2 k B T s q π r o 2 + R I d 2 + j e m V p + 3 2 k B T s q A e f f = n ˙ o u t , e ε i + 2 k B T s q + j e r 2 k B T e q A e f f + 5 2 k B T e , e q I d ( c ) n e , e + n n , e k B T s 1 + α T e , e T s = 1.36 · 10 9 I d 0.3 T n 0.95 m ˙ 0.57 L o 0.15 M 0.10 ε i 0.20 μ 0.35 d e 1.43 d o 1.71 ( d )
The plasma resistance in the insert region is computed as:
R = 3 η 4 π L e f f
A voltage drop is supposed to be between insert and orifice region; the intensity of this double sheath is computed according to literature [5,6,9] by using the following expression:
V d s = 9 I d k B T e , o 7.5 π n e , o r o 2 q 2 m e 2 q 2 / 3
The xenon viscosity in Ns/m2 is computed by means of the formula proposed by Goebel and Katz [1]:
μ = 2.3 · 10 5 · T r 0.71 + 0.29 T r , T r = T s / 289.7
True if T r > 1 . In the following, all the expressions for the current densities used in Equation 9 are recall for completeness:
j i = 0.61 q n e k B T e m i 1 / 2
j t h = 0.25 q n e 8 k B T e π m i 1 / 2
j e m = D T s 2 e x p q ϕ e f f k B T s
ϕ e f f = ϕ 0 q E c 4 π ϵ o
E c = n e k B T e ϵ 0 2 1 + 2 q V p k B T e 4 1 / 2
j e r = q 1 4 n e e x p q V p k B T s 8 k B T e π m e 1 / 2

3. Particle Swarm Optimization Methodology

3.1. Overview

The finding of a solution of a so complex system of equations would have had to be carried out by resorting to a very high number of manual attempts until the correct combination was reached. It was therefore decided to solve the equations system using a heuristic optimizer such as the Particle Swarm Optimization (PSO) [19,20,21].
PSO is a population based stochastic optimization technique developed by Eberhart and Kennedy in 1995 [16], inspired by the social behaviour of bird flocking or fish schooling. The particle swarm optimization concept consists in the generation of a swarm of individuals, each of them characterized by a proper set of variables’ value, and hence, a unique state variable vector, X ^ . During the generation of the swarm all these values are randomly assigned in such a way to cover the variables domain as it is defined in the initial problem setup. In fact, all the individuals try to find the problem’s solution searching inside the associated domain and their boundaries, using a performance index as main indicator. This latter is usually called fitness index or cost function. So, once the domain has been defined, it is necessary to determine the structure of the cost function to optimize. In this regard, the cost function can be written as a linear combination of target values and/or any constraints to be imposed. In this work it is required that all the residuals, as defined in the following, are minimized.
Now, deepening in the procedure, at the first iteration, after the swarm generation, the solver is applied to each individual in order to calculate the cost function that defines how much that solution would be good. The calculation of the cost function J is the starting point for the next iteration. In fact, to start the exploration of the variables’ domain, all the individuals need a proper velocity. This latter is calculated as the sum of many different contributions linked to their Personal Best and the Global Best index. The Personal Best (PB) and the Global Best (GB) are, respectively, the best solutions found from a single individual and from the whole swarm, in the exploration history. To understand what is best, each individual in the swarm must remember the set of variables associated with the best cost function it has previously encountered. In this way, the information coming from the next iteration q+1 can be compared with the previous one, q, and the exploration process is driven towards the most promising region of the domain. Another contribution to the velocity calculation comes from the velocity itself at the previous iteration. In fact, this term is also called Inertia. With the same criterion, the term due to the PB is called Nostalgia, and the other one, due to the GB, is called Social. The composition of these three terms defines the new velocity for the next iteration q+1, as showed in Figure 2.
So, the Personal Best and the Global Best of that iteration are used to update the individual velocity and hence their position (variables’ set). These positions, in the respective domains, are the starting point for the next iteration. At each time step, the domain can be explored, iteration after iteration, changing the velocity of all the individuals and hence accelerating toward the updated Personal Best and Global Best. The cycle stops with the criteria selected by the user, at a fixed number of iterations, or in dynamic way depending from the number of concluded iterations without significant optimization. In both cases the user obtain the final solution that should also be the optimal one in that domain.
The theory behind the PSO, as for all the others heuristic techniques, excludes that the solution obtained could be effectively the optimal one in absolute sense [21]. This situation could be relevant if this technique is used in a scientific field admitting high number of "local optima", such as astrodynamics. Unfortunately, the OHC model seems to be very similar to this situation because of the high number of variables and equations and their strong non-linearity. So, it is reasonable to start the optimization process many times, to overcome, or at least mitigate, the possibility that the result is not the real best solution, in absolute sense. Also for this reason, in the last years, PSO has been updated with many variants and advancements [22,23,24,25] in order to enhance the research capability and/or faster convergence.

3.2. PSO General Model

As aforementioned, the solution of the equation’s system has to be searched in a certain domain. The constrains are set to the ten independent variables i.e., neutral and plasma densities in the orifice and the emitter ( n n , e , n e , e , n n , o , n e , o ), the electron and wall temperatures ( T e , e , T e , o , T s , T w , o ), the emission length ( L e f f ) and the plasma potential ( V p ). In fact, by analysing Equations 1, it is simple to verify that the keeper plasma unknowns ( n n , k , n e , k , T e , k ) depend on the plasma density and the electron temperature of the orifice. Moreover the keeper wall temperature ( T w , k ) is set constant to typical value. The ten parameters govern the entire problem and therefore represent the 10-dimensional exploration domain of the PSO method. The boundaries of the exploration domain (min and max) for each variable can be arbitrarily defined; clearly, reasonable values must be used to foster the research and do not waste computational effort in regions representing non-physical operative conditions. Thus, a good way to proceed is to use values in agreement with the typical value reported for cathodes, well characterized in open literature.
Once the domain has been defined it is necessary to initialize the swarm assigning an initial state and velocity for each individual. Between the ten variables defined above, some of those, such the electronic temperature T e , have a range with a minimum and maximum in the same order of magnitude, and the value of the specific variable can be assigned with a simple random approach as in 19. Other variables, instead, are defined in a range with an order of magnitude much greater than one, such as neutral density n n , e , that spans from 10 21 up to 10 24 . For this second type of variable, a different equation for the state generation is required (20). Indeed, if the equation 19 were used for these variables, an in-homogeneous swarm distribution would be obtained, with a greater concentration near the upper boundary. This undesirable condition is overcome thanks to the equation 20. Using the notation of X m i n and X m a x , respectively, for the minimum and maximum of a general variable X of the state vector X ^ of the problem, the two different expressions can be considered:
X = X m i n + ( X m a x X m i n ) r a n d
with X = { V p , T e , T e , o , T s , L e f f , T w , o } and
X = 10 exp ( log ( X m i n ) + [ log ( X m a x ) log ( X m i n ) ] r a n d )
with X = { n n , e , n e , e , n n , o , n e , o } . Here r a n d is a random operator that randomly assumes values between 0 and 1.
The wider order of magnitude of certain variables is irrelevant for the initial velocity definition as it is quickly updated during the first iterations. Hence, in a simple manner, it also defines the first velocity of each variable as follow:
( V X ) z = ( X m i n ) z + [ ( X m a x ) z ( X m i n ) z ] r a n d
At this point, after the swarm initialization, each individual, with its own state vector, represents a possible solution of the problem, hence it is processed in the numerical model, in order to solve the associated equations (plasma density, electronic temperature, etc...), and obtain the residuals ϵ of that specific running setup. So, the cost function J is evaluated, for each individual, as the sum of the residuals ϵ ,
J z q = ϵ I + ϵ P o w , e + ϵ i o n , e + ϵ p r e s s , e + ϵ P o w , o + ϵ i o n , o + ϵ p r e s s , o
where:
  • ϵ I = RHS - LHS of the Equation 9b, (current)
  • ϵ P o w , e = RHS - LHS of the Equation 9c (emitter power)
  • ϵ i o n , e = RHS - LHS of the Equation 9a (emitter ions)
  • ϵ p r e s s , e = RHS - LHS of the Equation 9d (emitter pressure)
  • ϵ P o w , o = RHS - LHS of the Equation 2b (orifice power)
  • ϵ i o n , o = RHS - LHS of the Equation 2a (orifice ions)
  • ϵ p r e s s , o = RHS - LHS of the Equation 2c (orifice pressure)
where J z q is the cost function evaluated for the z t h individual at the q t h iteration. Furthermore, all the residuals, as defined above, have been divided to their reference variable (i.e., discharge current, power, and so on) in order to obtain dimensionless residuals with the same order of magnitude. This operation allows the optimization code to operate with different variables, but all with the same weight. Now, all these first cost functions J z q are assigned as the PB of each individual, so:
P B z q = J z q
while the GB has to be selected as the one with minimum value in the entire swarm, hence:
G B q = m i n J z q
With this latter step, the swarm initialization can be said concluded, and the real optimization process can start. So, to update the state of all the individuals, using the information coming from the exploration, and hence from the swarm, the new velocities can be calculated with the ultimate form:
V z q + 1 = C 1 V z q + C 2 r a n d ( P B z q X z q ) + C 3 r a n d ( G B q X z q )
where z is the z t h individual and q is the q t h iteration, C1-C2-C3 are constants in the range between 0 and 2, X refers to the variable whose velocity is being calculated, and r a n d is a random operator that assumes values between 0 and 1, as mentioned above.
Now, the individuals’ state can be updated with the velocity in order to start effectively the domain exploration:
X z q = X z q 1 + V z q
This updated state is considered as a temporary state since it has to be checked with respect to the original domain as defined at the beginning of this section. In fact, it could happen that the updated state results to be out of the domain bounds because of its movement due to the applied velocity. In this case, the state is corrected assigning the boundary value that it was trying to violate during the update.
With the updated state, also the cost function can be updated in order to compare the values with the ones calculated at the iteration before J z q 1 . If J z q results to be minor than the iteration before, for one or more individuals, than its own Personal Best would be updated and so:
P B z q = m i n J z q 1 , J z q
The same happens for the Global Best, where its own value is compared with all the PB of the entire swarm. If at least one individual has found a better solution P B z q , with respect to the GB of the iteration before G B q 1 , also the Global Best would be updated, so:
G B q = m i n G B q 1 , P B z q
This criterion is the fundamental of this optimization technique, since it allows to drive the individuals close to the promising region of the domain where better solutions are found. From this point ahead the cycle is continuously repeated until the stop criterion is satisfied.

3.3. PSO-Code Based Solver

It has been assessed that the developed PSO-code base solver can be used both to support the preliminary design of new hollow cathodes and for the rebuilding of plasma properties, if the electrical characteristics are known.
To design a new cathode, primarily the cathode geometry (i.e., insert, orifice and keeper diameters, orifice and keeper-orifice lengths), the propellant mass flow rate, the insert material and the current value have to be set. Secondly, a reasonable exploration domain must be defined. Then PSO-code can be run. It will give the evaluation of all the unknowns within a confidence bounds.
It has been verified (this become clear in the application example of Section 4) that the trend of the emitter wall temperature as function of mass flow rate can be better reproduced if a target total discharge voltage V t a r g e t is defined, which is calculated as the ratio of the absorbed power to the discharge current I d .
Thus it’s necessary to modify the formulation of the cost function J since it has to consider a new parameter, i.e., V t o t , as in Equation 29. In fact, in this case, the optimization process has to minimize also the difference between the total voltage V t o t and the target value V t a r g e t .
J z q = ϵ I + ϵ P o w , e + ϵ i o n , e + ϵ p r e s s , e + ϵ P o w , o + ϵ i o n , o + ϵ p r e s s , o + V t o t V t a r g e t
It is clear that this last approach can be used also for the rebuilding of plasma properties in an existing hollow cathode when the experimental electrical characteristics are known.

4. Results

This section reports the results of an application of the described methodology. The boundaries of the PSO exploration domain (Table 1) for each variable have been defined in agreement with the values estimated for the AlphCA Cathode, that is an OHC developed and tested by SITAEL [26,27]. A keeper surface temperature of 800K has been set for all the calculations according to the experimental measurements.
AlphCA has a Tantalum tube, 7.5 mm in diameter, which hosts a LaB6 insert that is 6 mm long and has an internal diameter of 3 mm. The tube is approximately thick 0.3 mm and 20.6 mm long and at its end there is an orifice that has a diameter of 0.4 mm and a length of 0.36 mm. The tube assembly is enclosed by means of a Molybdenum-alloy keeper that has an orifice (diameter: 0.6mm) at its end. Because the length of keeper orifice has not been found, a realistic value of 0.3 mm has been selected to perform the calculation. The propellant is pure xenon (contamination issues are not considered in the model). The cathode operates at a discharge current ranging from 1 to 3 A and at mass flow rates between 0.08 and 1 mg/s. The electrical characteristics in this operational envelope, measured in diode mode with keeper [27], are reported in Figure 3, as test2 and test3. The discharge voltage as function of the mass flow rate decreases reaching a plateau after a certain value; this value states the transition from the so called "plume" mode to the "spot" mode [1,5].
Table 1. Boundaries of the PSO exploration domain.
Table 1. Boundaries of the PSO exploration domain.
Variable min MAX Unit
n n , e 1 × 10 21 1 × 10 24 m 3
n e , e 1 × 10 18 1 × 10 21 m 3
V p 5 50 V
T e , e 1 4 e V
n n , o 1 × 10 21 1 × 10 23 m 3
n e , o 1 × 10 18 1 × 10 22 m 3
T e , o 1 4 e V
T s 1700 1950 K
L e f f 1 × 10 3 6 × 10 3 m
T w , o 1700 1950 K
The shaded zone in Figure 3 represents the confidence bound of the total discharge voltage, limited by the x ^ ± σ curves, in which the PSO solutions lay. More in detail, the mean x ^ , and the variance σ 2 , have been calculated for the best ten solutions obtained at each mass flow rate value. The mean and the variance are calculated with the standard equation, respectively as follow:
x ^ = 1 n z = 1 n x z
σ ^ 2 = 1 n z = 1 n x z 2 x ^ 2
Hence, the standard deviation σ = σ ^ 2 is used to generate the confidence bound ( x ^ ± σ ) around the mean values.
Several considerations arise from the observation of Figure 3:
  • the confidence bounds range between 2V and 5V for the current span examined;
  • the discharge current initially decreases with the mass flow rate, as expected, but it increases at the highest mass flow rates; this could be attributed to the absence of a second double sheath at the cathode orifice exit that could be negative because the plasma density decreases passing from the orifice region to the keeper region (see Figure 6a, Figure 7a);
  • the shaded zone perfectly covers the experimental measurements obtained for a current of 2A and it is slightly higher at 3A; at 1A and at lowest mass flow rates, the current is up to 40% lower, at its maximum, than the reference data, likely because the system of equations does not adequately describe the plasma physics in those conditions.
Figure 3. Present calculation’s confidence bound of the total discharge voltage as function of the propellant mass flow rate computed at 1A (a), 2A (b), 3A (c) by means of PSO; experimental measurements from Ref. [27] are reported for comparison as test2 and test3.
Figure 3. Present calculation’s confidence bound of the total discharge voltage as function of the propellant mass flow rate computed at 1A (a), 2A (b), 3A (c) by means of PSO; experimental measurements from Ref. [27] are reported for comparison as test2 and test3.
Preprints 110512 g003
Figure 5, Figure 6, Figure 7 report the best ten solutions obtained for the operative point at 1A and 2A in the whole mass flow rate span. These ten solutions have been selected from the PSO algorithm looking at the cost function. In fact, the cost function J is defined as the sum of the residuals ϵ evaluated for each equation, as reported in Section 3.2. Thus, with this assumption, the value of J can be used as a direct and global index of the solution quality and hence is presented in Figure 4.
Figure 4 presents the errors ϵ evaluated as the mean value of the ten solutions at 2A, as example, with respect to the mass flow rates. Similar figure would be drawn at 1A and 3A. The ten solutions presented at 1A and 2A, Figure 5, Figure 6 and Figure 7, are characterised by a J value that typically remains under the 5% from the minimum mass flow rate up to 0.6 mg/s. Beyond this value, and up to 1 mg/s, the residuals, ϵ , tend to increase as obviously do the cost function J, in consequence. This circumstance is also put in evidence from the greater dispersion of the solutions while moving toward the higher mass flow rate values. Anyway, also in this region of greater uncertainties, the J values are limited to 20% and 40%, respectively for the cases at 2A and 1A. Moreover, also in the worst case, the mean value of a single residual is not greater than 10% as in Figure 4, where the greater error is due to the calculation of the pressure in the orifice region.
All the trends of plasma properties (Figure 5, Figure 6 and Figure 7) reflect those reported in the open literature [5,6]: the neutral densities increases because of the internal pressure, particularly the trend is linear in the orifice and keeper regions (Equations 1, 2) and with 0.57 power in the emitter (Equation 9); the plasma densities trends have also an increasing trend because more neutrals means more collisions, thus ionization; the electron temperatures decrease with the mass flow rate in all regions because it decreases with the internal pressure as demonstrated in Ref. [1]; the orifice temperature slightly increases because the higher heat flux to the wall whereas no particular trend can be assumed for the emitter temperature where, as for the orifice region, an increasing trend is expected, too.
Figure 5. (a) Plasma density; (b) Neutral density; (c) Electronic temperature; (d) Emitter temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the emitter region, in case of operative points at 1A and 2A.
Figure 5. (a) Plasma density; (b) Neutral density; (c) Electronic temperature; (d) Emitter temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the emitter region, in case of operative points at 1A and 2A.
Preprints 110512 g005
The electron temperatures (Figure 5c, Figure 6c, Figure 7c), and the plasma (Figure 5a, Figure 6a, Figure 7a), and neutral (Figure 5b, Figure 6b, Figure 7b) densities found by PSO are very close to each other, generating a confidence region that seems to be very tight. The higher dispersion is that of emitter plasma density (i.e., ± 20 % of the average value). The wall temperatures of the emitter and orifice lay within less than 100K bound. Nevertheless it must be remarked that the PSO has been able to get the order of magnitude of neutral and plasma densities by searching within a range of three orders of magnitude.
Figure 6. (a) Plasma density; (b) Neutral density; (c) Electronic temperature; (d) Emitter temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the orifice region, in case of operative points at 1A and 2A.
Figure 6. (a) Plasma density; (b) Neutral density; (c) Electronic temperature; (d) Emitter temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the orifice region, in case of operative points at 1A and 2A.
Preprints 110512 g006
Figure 7. (a) Plasma density; (b) Neutral density; (c) Electronic temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the keeper region, in case of operative points at 1A and 2A.
Figure 7. (a) Plasma density; (b) Neutral density; (c) Electronic temperature as function of mass flow rate. All the figures refer to physical properties evaluated in the keeper region, in case of operative points at 1A and 2A.
Preprints 110512 g007
Figure 8 shows the PSO targeting the experimental measurements [27] of the discharge voltage as function of mass flow rate at 2A (i.e., PSO using Equation 29). It has been found that the results perfectly match the results calculated without the assumption made with Equation 29 (for this reason they are not shown). The only difference regards the emitter wall temperature. Figure 9 highlights that with this approach the expected increasing trend of the temperature can be detected. Moreover, the confidence bound is narrower. Here, in Figure 9, the two bounds have been calculated evaluating the mean and the standard deviation, as in equations 30 and 31, of the associated ten solutions provided from the two different simulations.

4.1. Computational Effort

The optimization code developed in this framework, namely the PSO, has been written in a Matlab environment. Simulations were conducted on a Personal Computer Intel-i7 @ 1.80GHz with 16.0GB RAM. For each research and optimization cycle an amount of 20mln possible solutions are evaluated in less than 7 minutes from which are extracted the best ten solutions. More in detail, a swarm of 2000 individuals for 1000 iterations has been used, for 10 cycles. As already mentioned, the strength of this method lies in its ability to be very fast and completely independent of the initial conditions. In fact, this independence allows to completely overcome the initial condition estimation problem, as reported in the open literature [6]. Anyway, faster or slower convergence can be selected changing the coefficients C1,C2 and C3 in the velocity equation (Equation 25) to balance the exploration and convergence dominance.

5. Conclusions

The paper has investigated the feasibility of using the Particle Swarm Optimization technique to solve a non-linear system of equations describing the plasma physics within orifice hollow cathodes.
The results of the assessment test cases, conducted on Sitael AlphCA, have shown that the coupling between PSO and plasma model reveals itself as fast, powerful and accurate preliminary design and rebuilding tool for OHCs.
The trends of all the unknowns of the problem are well predicted.
Particularly, the discharge voltage as function of mass flow rate is predicted with a confidence bound in the range 2 ÷ 5 V for the current span examined. The computations at 2A and 3A well match the experimental data whereas at 1A and low mass flow rates, the discharge current is underestimated (maximum by 40%), probably because the system of equation is inadequate to describe the plasma physics in those conditions. The discharge current trend as function of mass flow rate is of a parabolic type for all currents: it initially decreases, as expected, but it increases at high mass flow rates. This could be due to the absence of a double sheath at the cathode orifice exit in the model, which might dampen the discharge voltage increase at high mass flow rates.
All the others unknowns of the problem are computed with a very tiny confidence bound with the exception of the plasma density in the emitter region (± 20 % of the average value) and the wall temperature of the emitter and the orifice (± 50K respect to the average values). Regarding plasma densities, it is worth to remark that the PSO has been able to get within a range of three orders of magnitude the right one for neutral and plasma densities.
It has been verified that the confidence bound of the emitter temperature can be reduced and its trend can be improved by targeting the discharge voltage.
The difficulty to find solutions of the complex system of equation characterizing the plasma within the OHC has been fully overcame with the usage of the PSO technique. This enables the possibility to easily perform sensitivity studies on all the geometrical parameters of the cathode model, as well as quickly testing model modifications (i.e., the addition of the aforementioned double sheath at the cathode orifice exit).
As future development, a thermal model in FEMM [28] will be coupled to PSO-plasma model, according to Ref. [11]. This would give to the wall temperatures a new constraint thus more accurate plasma properties prediction are expected.

Author Contributions

Conceptualization M.P.; methodology, G.C. and M.P.; software, G.C.; validation, G.C.; formal analysis, G.C.; data curation, M.P. and G.C.; writing—original draft preparation, G.C. nd M.P.; writing—review and editing, F.B.; visualization, G.C.; supervision, F.B. and M.P.; project administration, F.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Italian research program in aerospace, PRORA DM662/20 (Programma Operativo Ricerche Aerospaziali) entrusted by MIUR (Ministero dell’Istruzione, Ministero dell’Università e della Ricerca).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors would like to thank A. Romano, C. Giaquinto and A. Smoraldi for the internal review of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
OHCs Orifice Hollow Cathodes
EP Electric Propulsion
HETs Hall Effect Thrusters
GITs Gridded Ion Thrusters
HEMPTs High Efficiency Multistage Plasma Thrusters
PSO Particle Swarm Optimization
HPTs Helicon Plasma Thrusters
CIRA Italian Aerospace Research Centre
PlasMCat Plasma Model for Cathodes
RHS Right Hand Side
LHS Left Hand Side
PB Personal Best
GB Global Best
Symbols
A = Area
L = Length
r = Radius
d = Diameter
q = Elementary charge
k B = Boltzmann’s constant
ϵ 0 = Vacuum permittivity
ϵ = Residual
E c = Electric field at the cathode sheath
I d = Discharge current
j = Current density
K i = Ionization rate coefficient
m e = Electron mass
m i = Ion mass
m ˙ = Mass flow rate
n = Density
n ˙ = Ion rate
p = Pressure
P = Power
R = Resistance
T = Temperature
V = Potential
ϕ = Work function
α = Degree of ionization
η = Plasma resistivity
ν = Collision frequency
σ = Cross section, standard deviation
γ = Specific heat ratio
M = Molecular Weight
μ = viscosity
Λ = Coulomb Logarithm
D = material-specific Richardson-Dushman constant
ϵ i = Average ionization energy (12.12 eV for xenon)
Subscripts
n = neutral
p = plasma
e = emitter, electron
o = orifice
k = keeper
e n = electron-neutral
e i = electron-ion
i = ion moving at Bohm speed
i o n = ion or ionization
e r = electron recombination
e m = thermionic emission
t h = thermal
e f f = effective
s = emitter surface
q = optimization iterations
z = swarm individual
w = wall
d s = double sheath

References

  1. Goebel, D.M.; Katz, I. Fundamentals of electric propulsion: Ion and Hall thrusters; John Wiley & Sons, 2008; Vol. 1. [Google Scholar]
  2. Lev, D.R.; Mikellides, I.G.; Pedrini, D.; Goebel, D.M.; Jorns, B.A.; McDonald, M.S. Recent progress in research and development of hollow cathodes for electric propulsion. Reviews of Modern Plasma Physics 2019, 3, 6. [Google Scholar] [CrossRef]
  3. Siegfried, D.; Wilbur, P. An investigation of mercury hollow cathode phenomena. In Proceedings of the 13th International Electric Propulsion Conference; 1978. [Google Scholar] [CrossRef]
  4. Salhi, A. Theoretical and experimental studies of orificed, hollow cathode operation. PhD thesis, Ohio State University, 1993. [Google Scholar]
  5. Domonkos, M.T. Evaluation of low current orificed hollow cathode. PhD thesis, Michigan University, 1999. [Google Scholar]
  6. Albertoni, R.; Pedrini, D.; Paganucci, F.; Andrenucci, M. A Reduced-Order Model for Thermionic Hollow Cathodes. IEEE Transactions on Plasma Science 2013, 41, 1731–1745. [Google Scholar] [CrossRef]
  7. Korkmaz, O. Global numerical model for the evaluation of the geometry and operation condition effects on hollow cathode insert and orfice region plasmas. MSc dissertation, Bogazici University, 2015. [Google Scholar]
  8. Taunay, P.Y.C.; Wordingham, C.J.; Choueiri, E. A 0-D model for orificed hollow cathodes with application to the scaling of total pressure. In Proceedings of the AIAA Propulsion and Energy, Indianapolis, IN; 2019. [Google Scholar]
  9. Gurciullo, A.; Fabris, A.L.; Potterton, T. Numerical study of a hollow cathode neutraliser by means of a zero-dimensional plasma model. Acta Astronautica 2020, 174, 219–235. [Google Scholar] [CrossRef]
  10. Gondol, N.; Tajmar, M. A volume-averaged plasma model for heaterless C12A7 electride hollow cathodes. CEAS Space Journal 2023, 15, 431–450. [Google Scholar] [CrossRef]
  11. Potrivitu, G.C.; Laterza, M.; Ridzuan, M.H.A.B.; Gui, S.Y.C.; Chaudhary, A.; Lim, J.W.M. Zero-Dimensional Plasma Model for Open-end Emitter Thermionic Hollow Cathodes with Integrated Thermal Model. In Proceedings of 37th International Electric Propulsion Conference Massachusetts Institute of Technology; Cambridge, MA, USA, 2022; Vol. IEPC-2022-116. [Google Scholar]
  12. Potrivitu, G.C.; Xu, S. Phenomenological plasma model for open-end emitter with orificed keeper hollow cathodes. Acta Astronautica 2022, 191, 293–316. [Google Scholar] [CrossRef]
  13. Panelli, M.; Giaquinto, C.; Smoraldi, A.; Battista, F. A Plasma Model for Orificed Hollow Cathode. In Proceedings of the 36th International Electric Propulsion Conference, Wien, Austria; 2019; Vol. IEPC-2019-452. [Google Scholar]
  14. Wordingham, C.J.; Taunay, P.Y.C.R.; Choueiri, E.Y. Critical Review of Orficed Hollow Cathode Modeling. In Proceedings of the 53rd AIAA/SAE/ASEE Joint Propulsion Conference, Atlanta, GA; 2017. [Google Scholar]
  15. Taunay, P.Y.C.; Wordingham, C.J.; Choueiri, E. An empirical scaling relationship for the total pressure in hollow cathodes. In Proceedings of the 2018 Joint Propulsion Conference. [CrossRef]
  16. Kennedy, J.; Eberhart, R. Particle swarm optimization. Proceedings of ICNN’95 - International Conference on Neural Networks; 1995; Vol. 4, pp. 1942–1948. [Google Scholar] [CrossRef]
  17. Kulkarni, M.N.K.; Patekar, M.S.; Bhoskar, M.T.; Kulkarni, M.O.; Kakandikar, G.; Nandedkar, V. Particle swarm optimization applications to mechanical engineering-A review. Materials Today: Proceedings 2015, 2, 2631–2639. [Google Scholar] [CrossRef]
  18. Coppola, G.; Panelli, M.; Battista, F. Preliminary design of helicon plasma thruster by means of particle swarm optimization. AIP Advances 2023, 13, 055209. Available online: https://pubs.aip.org/aip/adv/article-pdf/doi/10.1063/5.0149430/17349598/055209_1_5.0149430.pdf. [CrossRef]
  19. Zhang, Y.; Wang, S.; Ji, G.; et al. A comprehensive survey on particle swarm optimization algorithm and its applications. Mathematical problems in engineering 2015, 2015. [Google Scholar] [CrossRef]
  20. Houssein, E.H.; Gad, A.G.; Hussain, K.; Suganthan, P.N. Major advances in particle swarm optimization: Theory, analysis, and application. Swarm and Evolutionary Computation 2021, 63, 100868. [Google Scholar] [CrossRef]
  21. Bai, Q. Analysis of particle swarm optimization algorithm. Computer and information science 2010, 3, 180. [Google Scholar] [CrossRef]
  22. Jain, M.; Saihjpal, V.; Singh, N.; Singh, S.B. An overview of variants and advancements of PSO algorithm. Applied Sciences 2022, 12, 8392. [Google Scholar] [CrossRef]
  23. Kulkarni, M.N.K.; Patekar, M.S.; Bhoskar, M.T.; Kulkarni, M.O.; Kakandikar, G.; Nandedkar, V. Particle swarm optimization applications to mechanical engineering-A review. Materials Today: Proceedings 2015, 2, 2631–2639. [Google Scholar] [CrossRef]
  24. Ramírez-Ochoa, D.D.; Pérez-Domínguez, L.A.; Martínez-Gómez, E.A.; Luviano-Cruz, D. PSO, a swarm intelligence-based evolutionary algorithm as a decision-making strategy: A review. Symmetry 2022, 14, 455. [Google Scholar] [CrossRef]
  25. Moazen, H.; Molaei, S.; Farzinvash, L.; Sabaei, M. PSO-ELPM: PSO with elite learning, enhanced parameter updating, and exponential mutation operator. Information Sciences 2023, 628, 70–91. [Google Scholar] [CrossRef]
  26. Ducci, C.; Oslyak, K.; Dignani, D.; Albertoni, R.; Andrenucci, M. HT100D performance evaluation and endurance test results. In Proceedings of the Proceedings of 33rd International Electric Propulsion Conference, IEPC2013, Washington, DC, USA, 6-10 October 2013; Vol. IEPC-2013-140. [Google Scholar]
  27. Albertoni, R.; Pedrini, D.; Paganucci, F.; Andrenucci, M. Experimental Characterization of a LaB 6 Hollow Cathode for Low-Power Hall Effect Thrusters. Proceedings of Space Propulsion Conference, SPC2014, Cologne, Germany, 19-22 May 2014. [Google Scholar]
  28. Meeker, D. inite Element Method Magnetics–Version 4.0 User’s Manual. 2006. [Google Scholar]
Figure 1. (a) Schematic of typical Orifice Hollow Cathode and (b) sketch of plasma model regions.
Figure 1. (a) Schematic of typical Orifice Hollow Cathode and (b) sketch of plasma model regions.
Preprints 110512 g001
Figure 2. PSO: Schematic view of the velocity contributions. Full orange arrows are the single velocity contributions; The dotted orange arrows are the projections of the velocity terms that define the final velocity (blue/red arrow) and hence leads the updated state X z q + 1 .
Figure 2. PSO: Schematic view of the velocity contributions. Full orange arrows are the single velocity contributions; The dotted orange arrows are the projections of the velocity terms that define the final velocity (blue/red arrow) and hence leads the updated state X z q + 1 .
Preprints 110512 g002
Figure 4. Residual’s mean value evaluated at 2A. Each histogram represent the mean values of the single residuals obtained for the ten solutions provided by the PSO process.
Figure 4. Residual’s mean value evaluated at 2A. Each histogram represent the mean values of the single residuals obtained for the ten solutions provided by the PSO process.
Preprints 110512 g004
Figure 8. PSO targeting the experimental measurements (test3) [27] of the discharge voltage as function of mass flow rate at 2A.
Figure 8. PSO targeting the experimental measurements (test3) [27] of the discharge voltage as function of mass flow rate at 2A.
Preprints 110512 g008
Figure 9. Comparison of the expected emitter temperatures when PSO is targeting the experimental measurements [27] of the discharge voltage and when the discharge voltage is free. Both as functions of mass flow rate at 2A.
Figure 9. Comparison of the expected emitter temperatures when PSO is targeting the experimental measurements [27] of the discharge voltage and when the discharge voltage is free. Both as functions of mass flow rate at 2A.
Preprints 110512 g009

Short Biography of Authors

Preprints 110512 i001 Coppola Giovanni holds a Master degree in Space and Astronautical Engineering at the italian University of Rome, La Sapienza. He works in CIRA (Italian Aerospace Research Center), at the Space Division as researcher of the Space Propulsion Laboratory. He worked as researcher and system engineer in international projects related to hyper-sonic vehicles and electric propulsion systems. He is mainly involved in the CIRA electric propulsion laboratory projects and activity as researcher. The recent activities are associated with the development of numerical tools for the preliminary design of Gridded Ion Thrusters (GITs), Helicon Plasma Thrusters (HPTs) and HEMPT. Other activities are related to plasma modelling and optimization techniques.
Preprints 110512 i002 Mario Panelli holds a PhD in Aerospace Naval and Quality Management Engineering from the University of Naples, Federico II. He currently works in CIRA (Italian Aerospace Research Center), being affiliated to the Propulsion and Exploration Division as researcher of the Space Propulsion Laboratory. He worked as researcher, project and system engineer in different national and international project in the topics of methane-based propulsion systems, hybrid propulsion, solid propulsion and electric propulsion. He is involved in the CIRA electric propulsion projects as project engineer and researcher. His most recent activities are associated with design and testing of Orifice Hollow Cathodes (OHCs) and Hall Effect Thrusters (HETs), plasma and erosion modelling, and preliminary design of Gridded Ion Thrusters (GITs), Helicon Plasma Thrusters (HPTs) and HEMPT.
Preprints 110512 i003 Francesco Battista holds a PhD in Aerospace Engineering from University of Pisa. He currently works in CIRA (Italian Aerospace Research Center), being affiliated to the Propulsion and Exploration Division as responsible of the Space Propulsion Laboratory. He worked as system engineer and project manager in different national and international project mainly in the topics of methane-based propulsion systems, hybrid propulsion and electric propulsion. He is involved in the CIRA electric propulsion projects as project manager and technical responsible. He is author or co-author of more than 100 papers, published on international journals and conference proceedings.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated