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

1. Introduction

Orifice Hollow Cathodes (OHCs) [1] are used both as electronic sources and/or neutralizers for electric propulsion (EP) devices, e.g., Hall Effect Thrusters (HETs), Gridded Ion Thrusters (GITs), and High-Efficiency Multistage Plasma Thrusters (HEMPTs) [2]. Aot of computational tools, based on phenomenological models, also known as zero-dimensional (0D) or reduced-order models, aim to predict OHCs’ plasma properties and performances and are useful to support their design, and have been being developed since the 1980s [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 the original model based on the 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]; and the detailed plasma modeling within heaterless cathodes and coupling with the 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], which 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 constitute most of the models presented in the openiterature [5,6,7,9,10,11,12] are commonly iteratively solved to face their own non-linear characteristics. Each author adopts a different approach. PlasMCat follows a procedure similar to Albertoni et al. [6] and in agreement with them, and assesses that the adopted solution technique hasimited 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 of solving 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 the 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 applications 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 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 HollowCathode

2.1. CathodeArchitecture

In the typical OHC 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, and 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]. A hollow cylinder, namely an insert made of aow-work functional material, is placed inside the tube. The insert is the core part of the cathode because it emits electrons by means of a thermionic effect. Barium-based and rare-earth emitters (e.g., LaB₆, CeB₆), which have, respectively, a work function of about 2 eV and 2.7 eV [1] allow us to provide current densities as high as $10^{5}$ A/m² at a surface temperature of about 1300 K and 1900 K, respectively. The main drawback for Barium-based emitters is that they are highly susceptible to oxygen contamination, whereas rare-earth emitters require very high temperatures 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 sleeve is usually placed between the emitter and the main tube.

A heater is wrapped around the tube to trigger thermionic emission. If it is designed well, a cathode can sustain the emitting temperature by means of plasma heating; thus, the 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 radiationosses).

Another important element is the keeper electrode that encloses the heater–main tube assembly. It has the double function of helping the extraction of the electrons to the discharge plasma and protecting the cathode from the ion bombardment.

2.2. Definition andAssumptions

The hollow cathode internal plasma volume is divided into three regions: the keeper region, represented by “k” (which includes the keeper orifice region and the cathode–keeper gap), the insert or emitting region, represented by “e”, and the orifice region, represented by “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 emissionength ( $L_{e f f}$ ), the sheath potential ( $V_{p}$ ), and the wall temperature of the emitter ( $T_{s}$ ), the orifice ( $T_{w, o}$ ) and the keeper ( $T_{w, k}$ ).

The main assumptions, in agreement with the norm for 0D models [5,6,9], are as follows:

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 choked gas flow model is used in the orifice regions of both the cathode tube and the 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, which is transitional, is not modelled as continuous and governed by the Poiseuille’saw, 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 a 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 ModelEquations

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

2.3.1. 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), the balance of ion flux, Equation (2), and the equation of state, Equation (3).

${\begin{matrix} n_{e, k} = n_{e, o} \frac{A_{o}}{A_{k}} \sqrt{\frac{T_{e, o}}{T_{e, k}}} & (1) \\ 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} & (2) \\ ({n_{e, k} + n}_{n, k}) k_{B} T_{w, k} (1 + \frac{α T_{e, k}}{T_{w, k}}) = \frac{\dot{m}}{π r_{k}^{2}} \sqrt{\frac{R_{g} T_{w, k}}{γ} (1 + \frac{α T_{e, k}}{T_{w, k}})} & (3) \end{matrix}$

2.3.2. Orifice

The system of equations for the cathode’s orifice region writes a balance for the ion flux Equation (4) and for the power Equation (5), and states that pressure calculated by means of kinetic theory is equal to that computed with continuum flow and choked orifice, Equation (6). 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.

${\begin{matrix} \underset{{\dot{n}}_{i o n}}{\underset{︸}{q π r_{o}^{2} L_{o} n_{n, o} n_{e, o} K_{i}}} = j_{i, o} [2 π r_{o} L_{o} + π r_{o}^{2} (1 - \frac{1}{4} \frac{1}{0.61} \sqrt{\frac{8}{π}})] + j_{t h, o} π r_{o}^{2} & (4) \\ R I_{d}^{2} = {\dot{n}}_{i o n} 〈ϵ_{i}〉 + \frac{5}{2} \frac{k_{B}}{q} I_{d} (T_{e, o} {- T}_{e, e}) & (5) \\ ({n_{e, o} + n}_{n, o}) k_{B} T_{w, o} (1 + \frac{α T_{e, o}}{T_{w, o}}) = \frac{\dot{m}}{π r_{o}^{2}} \sqrt{\frac{R_{g} T_{w, o}}{γ} (1 + \frac{α T_{e, o}}{T_{w, o}})} & (6) \end{matrix}$

Here, the plasma resistance is calculated as:

$R = \frac{η L_{o}}{π r_{o}^{2}}$

(7)

The plasma resistivity is generally computed as:

$η = \frac{(ν_{e i} + ν_{e n}) m_{e}}{n_{e} q^{2}}$

(8)

and collision frequencies are:

$ν_{e i} = 2.9 \cdot 10^{- 12} (\frac{n_{e} n Λ}{{(\frac{T_{e} k_{B}}{q})}^{\frac{3}{2}}})$

(9)

$v_{e n} = 〈σ_{e n}〉 n_{n} (\frac{k_{B} T_{e}}{n_{e}})$

(10)

where the Coulombogarithm is:

$n Λ = 23 - \frac{1}{2} log (\frac{10^{- 6} n_{e}}{{(\frac{T_{e} k_{B}}{q})}^{3}})$

(11)

and the fit for the Maxwellian average elastic electron-neutral cross section ( $[m^{2}]$ ) proposed in ref. [13] (Table 1) for xenon is:

$L o g_{10} (〈σ_{e n}〉) = \frac{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} (\frac{T_{e} k_{B}}{q})$

(12)

2.3.3. Insert

Four equations are written in order to obtain 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 (13), the current conservation, Equation (14), the power balance, Equation (15), and the continuity equations, Equation (16).

${\begin{matrix} q π r_{e}^{2} L_{e f f} n_{n, e} n_{e, e} K_{i} + j_{i, o} π r_{o}^{2} = \underset{{\dot{n}}_{o u t, e}}{\underset{︸}{j_{i, e} (A_{e f f} + π r_{e}^{2} - π r_{o}^{2}) + j_{t h, e} π r_{e}^{2}}} & (13) \\ \frac{I_{d}}{A_{e f f}} = j_{i, e} + j_{e m} - j_{e r} & (14) \\ j_{i, o} (V_{d s} + 2 \frac{k_{B} T_{s}}{q}) π r_{o}^{2} + R I_{d}^{2} + j_{e m} (V_{p} + \frac{3}{2} \frac{k_{B} T_{s}}{q}) A_{e f f} = {\dot{n}}_{o u t, e} (〈ε_{i}〉 + 2 \frac{k_{B} T_{s}}{q}) + j_{e r} 2 \frac{k_{B} T_{e}}{q} A_{e f f} + \frac{5}{2} \frac{k_{B} T_{e, e}}{q} I_{d} & (15) \\ (n_{e, e} + n_{n, e}) k_{B} T_{s} (1 + \frac{α T_{e, e}}{T_{s}}) = 1.36 \cdot 10^{- 9} \frac{{I_{d}}^{0.3} {T_{n}}^{0.95} {\dot{m}}^{0.57} {L_{o}}^{0.15}}{M^{0.10} {〈 ε_{i} 〉}^{0.20} μ^{0.35} {d_{e}}^{1.43} {d_{o}}^{1.71}} & (16) \end{matrix}$

The plasma resistance in the insert region is computed as:

$R = \frac{3 η}{4 π L_{e f f}}$

(17)

A voltage drop is supposed to be between the insert and orifice region; the intensity of this double sheath is computed according to theiterature [5,6,9] by using the following expression:

$V_{d s} = {(\frac{9 I_{d} k_{B} T_{e, o}}{7.5 π n_{e, o} r_{o}^{2} q^{2}} \sqrt{\frac{m_{e}}{2 q}})}^{2 / 3}$

(18)

The xenon viscosity in $N s / m^{2}$ is computed by means of the formula proposed by Goebel and Katz [1]:

$μ = 2.3 \cdot 10^{- 5} \cdot {T_{r}}^{0.71 + \frac{0.29}{T r}}, T_{r} = T_{s} / 289.7$

(19)

which is true if $T_{r} > 1$ . In the following, all the expressions for the current densities used in Equations (13)–(16) are recalled for completeness:

$j_{i} = 0.61 q n_{e} {(\frac{k_{B} T_{e}}{m_{i}})}^{1 / 2}$

(20)

3. Particle Swarm OptimizationMethodology

3.1. Overview

The finding of a solution to 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 of 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, $\hat{X}$ . During the generation of the swarm, all these values are randomly assigned in such a way as to cover the variable’s domain as it is defined in the initial problem setup. In fact, all the individuals try to find the problem’s solution by searching inside the associated domain and their boundaries, using a performance index as the main indicator. Thisatter is usually called the 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 ainear 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. Thisatter is calculated as the sum of many different contributionsinked 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’s 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 a dynamic way depending on the number of concluded iterations without significant optimization. In both cases, the user obtains the final solution that should also be the optimal one in that domain.

The theory behind the PSO, as for all the other 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 a 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 ateast mitigate, the possibility that the result is not the real best solution, in absolute sense. Also for this reason, in recent years, PSO has been updated with many variants and advancements [17,22,23,24] in order to enhance the research capability and/or faster convergence.

3.2. PSO GeneralModel

As aforementioned, the solution of the equation’s system has to be searched in a certain domain. The constraints 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 emissionength ( $L_{e f f}$ ) and the plasma potential ( $V_{p}$ ). In fact, by analysing Equations (1)–(3), 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 the 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 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, which is well characterized in the 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 as 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 Equation (26). 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 (Equation (27)). Indeed, if Equation (26) were used for these variables, an in-hom*ogeneous swarm distribution would be obtained, with a greater concentration near the upper boundary. This undesirable condition is overcome thanks to Equation (27). 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 $\hat{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$

(26)

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)$

(27)

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 to 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 follows:

${(V_{X})}_{z} = {(X_{m i n})}_{z} + [{(X_{m a x})}_{z} - {(X_{m i n})}_{z}] * r a n d$

(28)

At this point, after the swarm initialization, each individual, with its own state vector, represents a possible solution to 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} ∥$

(29)

where:

$ϵ_{I}$ = RHS − LHS of Equation (14), (current);
$ϵ_{P o w, e}$ = RHS − LHS of Equation (15) (emitter power);
$ϵ_{i o n, e}$ = RHS − LHS of Equation (13) (emitter ions);
$ϵ_{p r e s s, e}$ = RHS − LHS of Equation (16) (emitter pressure);
$ϵ_{P o w, o}$ = RHS − LHS of Equation (5) (orifice power);
$ϵ_{i o n, o}$ = RHS − LHS of Equation (4) (orifice ions);
$ϵ_{p r e s s, o}$ = RHS − LHS of Equation (6) (orifice pressure).

where $J_{z}^{q}$ is the cost function evaluated for the z-th individual at the q-th iteration. Furthermore, all the residuals, as defined above, have been divided into 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}$

(30)

while the GB has to be selected as the one with the minimum value in the entire swarm; hence:

$G B^{q} = m i n (J_{z}^{q})$

(31)

With thisatter 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})$

(32)

where z is the z-th individual and q is the q-th 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}$

(33)

This updated state is considered 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 smaller than the iteration before for one or more individuals, then its own Personal Best would be updated, and so:

$P B_{z}^{q} = m i n (J_{z}^{q - 1}, J_{z}^{q})$

(34)

The same happens for the Global Best, where its own value is compared with all the PB of the entire swarm. If ateast one individual has found a better solution $P B_{z}^{q}$ , with respect to the GB of the iteration before $G B^{q - 1}$ , the Global Best would also be updated, so:

$G B^{q} = m i n (G B^{q - 1}, P B_{z}^{q})$

(35)

This criterion is fundamental to this optimization technique since it allows us 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 BasedSolver

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-orificeengths), 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 becomes clear in the application example of Section 4) that the trend of the emitter wall temperature as a 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 is 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 (36). 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} ∥$

(36)

It is clear that thisast 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 2) 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 [25,26]. A keeper surface temperature of 800 K 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 mmong and has an internal diameter of 3 mm. The tube is approximately thick 0.3 mm and 20.6 mmong and at its end there is an orifice that has a diameter of 0.4 mm and aength of 0.36 mm. The tube assembly is enclosed by means of a Molybdenum alloy keeper that has an orifice (diameter: 0.6 mm) at its end. Because theength 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 [26], are reported in Figure 3, as test2 and test3. The discharge voltage as a 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].

The shaded zone in Figure 3 represents the confidence bound of the total discharge voltage, imited by the $\hat{x} \pm σ$ curves, in which the PSO solutionsay. More in detail, the mean $\hat{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 follows:

$\hat{x} = \frac{1}{n} \sum_{z = 1}^{n} x_{z}$

(37)

${\hat{σ}}^{2} = \frac{1}{n} \sum_{z = 1}^{n} (x_{z}^{2} - {\hat{x}}^{2})$

(38)

Hence, the standard deviation $σ = \sqrt{{\hat{σ}}^{2}}$ is used to generate the confidence bound ( $\hat{x} \pm σ$ ) around the mean values.

Several considerations arise from the observation of Figure 3:

The confidence bounds range between 2 V and 5 V 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 4a and Figure 5a);
The shaded zone perfectly covers the experimental measurements obtained for a current of 2A and it is slightly higher at 3A; at 1A and atowest mass flow rates, the current is up to 40%ower, at its maximum, than the reference data, ikely because the system of equations does not adequately describe the plasma physics in those conditions.

Figure 4, Figure 5 and Figure 6 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 algorithmooking 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 7.

Figure 7 presents the errors $ϵ$ evaluated as the mean value of the ten solutions at 2A, for example, with respect to the mass flow rates. A similar figure would be drawn at 1A and 3A. The ten solutions presented at 1A and 2A, Figure 4, Figure 5 and Figure 6, 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. This circ*mstance is also evidenced by the greater dispersion of the solutions while moving toward the higher mass flow rate values. Also in this region of greater uncertainties, the J values areimited to 20% and 40%, respectively, for the cases at 2A and 1A. Moreover, in the worst case, the mean value of a single residual is not greater than 10% as in Figure 7, where the greater error is due to the calculation of the pressure in the orifice region.

All the trends of plasma properties (Figure 4, Figure 5 and Figure 6) reflect those reported in the openiterature [5,6]: the neutral densities increase because of the internal pressure, particularly the trend isinear in the orifice and keeper regions (Equations (3) and (6)) and with 0.57 power in the emitter (Equation (16)); the plasma densities trends have also an increasing trend because more neutrals means more collisions, and 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 of 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.

The electron temperatures (Figure 4c, Figure 5c and Figure 6c), and the plasma (Figure 4a, Figure 5a and Figure 6a), and neutral (Figure 4b, Figure 5b and Figure 6b) 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 orificeay withiness than 100 K bound. Nevertheless, it must be remarked that the PSO has been able to obtain the order of magnitude of neutral and plasma densities by searching within a range of three orders of magnitude.

Figure 8 shows the PSO targeting the experimental measurements [26] of the discharge voltage as a function of mass flow rate at 2A (i.e., PSO using Equation (36)). It has been found that the results perfectly match the results calculated without the assumption made with Equation (36) 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 (37) and (38), of the associated ten solutions provided from the two different simulations.

ComputationalEffort

The optimization code developed in this framework, namely the PSO, has been written in a Matlab environment (v.r2 2022). Simulations were conducted on an Intel-i7 Personal Computer @ 1.80 GHz with 16.0 GB RAM (Santa Clara, CA, USA). For each research and optimization cycle, an amount of 20 million possible solutions are evaluated iness than 7 min from which the best 10 solutions are extracted. In more detail, a swarm of 2000 individuals for 1000 iterations has been used for 10 cycles. As already mentioned, the strength of this methodies in its ability to be very fast and completely independent of the initial conditions. In fact, this independence allows us to completely overcome the initial condition estimation problem, as reported in the openiterature [6]. Anyway, faster or slower convergence can be selected by changing the coefficients C1, C2 and C3 in the velocity equation (Equation (32)) 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, Pisa, Italy, have shown that the coupling between PSO and plasma model reveals itself as a 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 a 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 andow 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 a 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 other 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 (±50 K respect to the average values). Regarding plasma densities, it is worth remarking that the PSO has been able to get within a range of three orders of magnitude of 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 in finding solutions to the complex system of equations characterizing the plasma within the OHC has been fully overcome with the usage of the PSO technique. This enables the possibility of easily performing 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 [27] will be coupled to the PSO-plasma model, according to ref. [11]. This would give the wall temperatures a new constraint; thus, more accurate plasma properties predictions 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. The data are not publicly available due to privacy.

Acknowledgments

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

Conflicts of Interest

The authors declare no conflicts 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
$\dot{m}$	Mass flow rate
n	Density
$\dot{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

Goebel, D.M.; Katz, I. Fundamentals of Electric Propulsion: Ion and Hall Thrusters; John Wiley & Sons: Hoboken, NJ, USA, 2008; Volume 1. [Google Scholar]
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. Rev. Mod. Plasma Phys. 2019, 3, 6. [Google Scholar] [CrossRef]
Siegfried, D.; Wilbur, P. An investigation of mercury hollow cathode phenomena. In Proceedings of the 13th International Electric Propulsion Conference, San Diego, CA, USA, 25–27 April 1978. [Google Scholar] [CrossRef]
Salhi, A. Theoretical and Experimental Studies of Orificed, Hollow Cathode Operation. Ph.D. Thesis, Ohio State University, Columbus, OH, USA, 1993. [Google Scholar]
Domonkos, M.T. Evaluation of Low Current Orificed Hollow Cathode. Ph.D. Thesis, Michigan University, Ann Arbor, MI, USA, 1999. [Google Scholar]
Albertoni, R.; Pedrini, D.; Paganucci, F.; Andrenucci, M. A Reduced-Order Model for Thermionic Hollow Cathodes. IEEE Trans. Plasma Sci. 2013, 41, 1731–1745. [Google Scholar] [CrossRef]
Korkmaz, O. Global Numerical Model for the Evaluation of the Geometry and Operation Condition Effects on Hollow Cathode Insert and Orfice Region Plasmas. Master’s Dissertation, Bogazici University, Istanbul, Turkey, 2015. [Google Scholar]
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, USA, 19–22 August 2019. [Google Scholar]
Gurciullo, A.; Fabris, A.L.; Potterton, T. Numerical study of a hollow cathode neutraliser by means of a zero-dimensional plasma model. Acta Astronaut. 2020, 174, 219–235. [Google Scholar] [CrossRef]
Gondol, N.; Tajmar, M. A volume-averaged plasma model for heaterless C12A7 electride hollow cathodes. CEAS Space J. 2023, 15, 431–450. [Google Scholar] [CrossRef]
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 the 37th International Electric Propulsion Conference Massachusetts Institute of Technology, Cambridge, MA, USA, 19–23 June 2022; Volume IEPC-2022-116. [Google Scholar]
Potrivitu, G.C.; Xu, S. Phenomenological plasma model for open-end emitter with orificed keeper hollow cathodes. Acta Astronaut. 2022, 191, 293–316. [Google Scholar] [CrossRef]
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, 15–20 September 2019; Volume IEPC-2019-452. [Google Scholar]
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, USA, 10–12 July 2017. [Google Scholar]
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, Cincinnati, OH, USA, 9–11 July 2018. [Google Scholar] [CrossRef]
Kennedy, J.; Eberhart, R. Particle swarm optimization. In Proceedings of the ICNN’95-International Conference on Neural Networks, Perth, WA, Australia, 27 November–1 December 1995; Volume 4, pp. 1942–1948. [Google Scholar] [CrossRef]
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. Mater. Today Proc. 2015, 2, 2631–2639. [Google Scholar] [CrossRef]
Coppola, G.; Panelli, M.; Battista, F. Preliminary design of helicon plasma thruster by means of particle swarm optimization. AIP Adv. 2023, 13, 055209. [Google Scholar] [CrossRef]
Zhang, Y.; Wang, S.; Ji, G. A comprehensive survey on particle swarm optimization algorithm and its applications. Math. Probl. Eng. 2015, 2015, 931256. [Google Scholar] [CrossRef]
Houssein, E.H.; Gad, A.G.; Hussain, K.; Suganthan, P.N. Major advances in particle swarm optimization: Theory, analysis, and application. Swarm Evol. Comput. 2021, 63, 100868. [Google Scholar] [CrossRef]
Bai, Q. Analysis of particle swarm optimization algorithm. Comput. Inf. Sci. 2010, 3, 180. [Google Scholar] [CrossRef]
Jain, M.; Saihjpal, V.; Singh, N.; Singh, S.B. An overview of variants and advancements of PSO algorithm. Appl. Sci. 2022, 12, 8392. [Google Scholar] [CrossRef]
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]
Moazen, H.; Molaei, S.; Farzinvash, L.; Sabaei, M. PSO-ELPM: PSO with eliteearning, enhanced parameter updating, and exponential mutation operator. Inf. Sci. 2023, 628, 70–91. [Google Scholar] [CrossRef]
Ducci, C.; Oslyak, S.; Dignani, D.; Albertoni, R.; Andrenucci, M. HT100D performance evaluation and endurance test results. In Proceedings of the 33rd International Electric Propulsion Conference, IEPC2013, Washington, DC, USA, 6–10 October 2013; Volume IEPC-2013-140. [Google Scholar]
Albertoni, R.; Pedrini, D.; Paganucci, F.; Andrenucci, M. Experimental Characterization of a LaB 6 Hollow Cathode for Low-Power Hall Effect Thrusters. In Proceedings of the Space Propulsion Conference, SPC2014, Cologne, Germany, 19–22 May 2014. [Google Scholar]
Meeker, D. Finite Element Method Magnetics–Version 4.0 User’s Manual. 2006. Available online: https://www.femm.info/wiki/HomePage (accessed on 5 June 2024).

Figure 1. (a) Schematic of typical Orifice Hollow Cathode and (b) sketch of plasma model regions.

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 henceead the updated state $X_{z}^{q + 1}$ . Circles represent figurative positions in the exploration domain.

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. [26] are reported for comparison as test2 and test3.

Figure 4. (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 the case of operative points at 1A and 2A.

Figure 5. (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 the 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 emitter region in the case of operative points at 1A and 2A.

Figure 7. Residual’s mean value evaluated at 2A. Each histogram represents the mean values of the single residuals obtained for the ten solutions provided by the PSO process.

Figure 8. PSO targeting the experimental measurements (test3) [26] of the discharge voltage as function of mass flow rate at 2A.

Figure 9. Comparison of the expected emitter temperatures when PSO is targeting the experimental measurements [26] of the discharge voltage and when the discharge voltage is free, both as functions of mass flow rate at 2A.

Table 1. Coefficients of xenon electron-neutral cross section numerical fit [13].

	$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$

Table 2. Boundaries of the PSO exploration domain.

Variable	Min	Max	Unit
$n_{n, e}$	$1 \times 10^{21}$	$1 \times 10^{24}$	m⁻³
$n_{e, e}$	$1 \times 10^{18}$	$1 \times 10^{21}$	m⁻³
$V_{p}$	5	50	V
$T_{e, e}$	1	4	eV
$n_{n, o}$	$1 \times 10^{21}$	$1 \times 10^{23}$	m⁻³
$n_{e, o}$	$1 \times 10^{18}$	$1 \times 10^{22}$	m⁻³
$T_{e, o}$	1	4	eV
$T_{s}$	1700	1950	K
$L_{e f f}$	$1 \times 10^{- 3}$	$6 \times 10^{- 3}$	m
$T_{w, o}$	1700	1950	K

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.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).