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

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., LaB₆, CeB₆), which have respectively a work function of about 2 eV and 2.7 eV [1], allow to provide current densities as high as 10⁵ A/m⁻² 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_{eff}$), 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_{ion}=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).

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{n}_{e,k}=n_{e,o}\frac{A_o}{A_k}\sqrt{\frac{T_{e,o}}{T_{e,k}}}\hfill & \left(a\right)\hfill \\ \\ q\pi r_k^2{L}_k{n}_{n,k}{n}_{e,k}{K}_i=j_{i,k}\left(2\pi r_k{L}_k\right)+2j_{th,k}\pi r_k^2\hfill & \left(b\right)\hfill \\ \\ \left({n_{e,k}+n}_{n,k}\right)k_B{T}_{w,k}\left(1+\frac{\alpha T_{e,k}}{T_{w,k}}\right)=\frac{\dot{m}}{\pi r_k^2}\sqrt{\frac{R_g{T}_{w,k}}{\gamma }\left(1+\frac{\alpha T_{e,k}}{T_{w,k}}\right)}\hfill & \left(c\right)\hfill \end{array}\right.\end{array}$$

(1)

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.

$$\begin{array}{c}\left\{\begin{array}{ll}\underset{{\dot{n}}_{ion}}{\underbrace{q\pi r_o^2{L}_o{n}_{n,o}{n}_{e,o}{K}_i}}=j_{i,o}\left[2\pi r_o{L}_o+\pi r_o^2\left(1-\frac{1}{4}\frac{1}{0.61}\sqrt{\frac{8}{\pi }}\right)\right]+j_{th,o}\pi r_o^2\hfill & \left(a\right)\hfill \\ R{I}_d^2={\dot{n}}_{ion}〈{ϵ}_i〉+\frac{5}{2}\frac{k_B}{q}{I}_d\left(T_{e,o}{-T}_{e,e}\right)& \left(b\right)\\ \left({n_{e,o}+n}_{n,o}\right)k_B{T}_{w,o}\left(1+\frac{\alpha T_{e,o}}{T_{w,o}}\right)=\frac{\dot{m}}{\pi r_o^2}\sqrt{\frac{R_g{T}_{w,o}}{\gamma }\left(1+\frac{\alpha T_{e,o}}{T_{w,o}}\right)}& \left(c\right)\end{array}\right.\hfill \end{array}$$

(2)

Here the plasma resistance is calculated as:

$$R=\frac{\eta L_o}{\pi r_o^2}$$

(3)

The plasma resistivity is generally computed as:

$$\eta =\frac{({\nu }_{ei}+{\nu }_{en})m_e}{n_e{q}^2}$$

(4)

and collision frequencies are:

$${\nu }_{ei}=2.9·{10}^{-12}\left(\frac{n_{e}ln\Lambda }{{\left(\frac{T_e{k}_B}{q}\right)}^{\frac{3}{2}}}\right)$$

(5)

$$v_{en}=〈{\sigma }_{en}〉n_n\left(\frac{k_B{T}_e}{n_e}\right)$$

(6)

where the Coulomb logarithm is:

$$ln\Lambda =23-\frac{1}{2}log\left(\frac{{10}^{-6}{n}_e}{{\left(\frac{T_e{k}_B}{q}\right)}^3}\right)$$

(7)

and the fit for the Maxwellian average elastic electron-neutral cross section ($\left[m^2\right]$) proposed in Ref. [13] for xenon is:

$$Lo{g}_{10}\left(〈{\sigma }_{en}〉\right)=\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^{*}=Lo{g}_{10}\left(\frac{T_e{k}_B}{q}\right)$$

(8)

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

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}q\pi r_e^2{L}_{eff}{n}_{n,e}{n}_{e,e}{K}_i+j_{i,o}\pi r_o^2=\underset{{\dot{n}}_{out,e}}{\underbrace{j_{i,e}\left(A_{eff}+\pi r_e^2-\pi r_o^2\right)+j_{th,e}\pi r_e^2}}\hfill & \left(a\right)\hfill \\ \frac{I_d}{A_{eff}}=j_{i,e}+j_{em}-j_{er}\hfill & \left(b\right)\hfill \\ \\ j_{i,o}\left(V_{ds}+2\frac{k_B{T}_s}{q}\right)\pi r_o^2+R{I}_d^2+j_{em}\left(V_p+\frac{3}{2}\frac{k_B{T}_s}{q}\right)A_{eff}=\hfill \\ \\ {\dot{n}}_{out,e}\left(〈{\epsilon }_i〉+2\frac{k_B{T}_s}{q}\right)+j_{er}2\frac{k_B{T}_e}{q}{A}_{eff}+\frac{5}{2}\frac{k_B{T}_{e,e}}{q}{I}_d\hfill & \left(c\right)\hfill \\ \\ \left(n_{e,e}+n_{n,e}\right)k_B{T}_s\left(1+\frac{\alpha T_{e,e}}{T_s}\right)=1.36·{10}^{-9}\frac{{I_d}^{0.3}{T_n}^{0.95}{\dot{m}}^{0.57}{L_o}^{0.15}}{M^{0.10}{\langle {\epsilon }_i\rangle }^{0.20}{\mu }^{0.35}{d_e}^{1.43}{d_o}^{1.71}}\hfill & \left(d\right)\hfill \end{array}\right.\end{array}$$

(9)

The plasma resistance in the insert region is computed as:

$$R=\frac{3\eta }{4\pi L_{eff}}$$

(10)

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_{ds}={\left(\frac{9I_d{k}_B{T}_{e,o}}{7.5\pi n_{e,o}{r}_o^2{q}^2}\sqrt{\frac{m_e}{2q}}\right)}^{2/3}$$

(11)

The xenon viscosity in Ns/m² is computed by means of the formula proposed by Goebel and Katz [1]:

$$\mu =2.3·{10}^{-5}·{T_r}^{0.71+\frac{0.29}{Tr}},T_r=T_s/289.7$$

(12)

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.61q{n}_e{\left(\frac{k_B{T}_e}{m_i}\right)}^{1/2}$$

(13)

$$j_{th}=0.25q{n}_e{\left(\frac{8k_B{T}_e}{{\pi m}_i}\right)}^{1/2}$$

(14)

$$j_{em}=\mathrm{D}{T}_s^{2}exp\left(-\frac{q{\varphi }_{eff}}{k_B{T}_s}\right)$$

(15)

$${\varphi }_{eff}={\varphi }_0-\sqrt{\frac{q{E}_c}{4\pi {ϵ}_o}}$$

(16)

$$E_c=\sqrt{\frac{n_e{k}_B{T}_e}{{ϵ}_0}}{\left[2\sqrt{1+2\frac{q{V}_p}{k_B{T}_e}}-4\right]}^{1/2}$$

(17)

$$j_{er}=\mathrm{q}\frac{1}{4}{\mathrm{n}}_{\mathrm{e}}exp\left(-\frac{q{V}_p}{k_B{T}_s}\right){\left(\frac{8k_B{T}_e}{{\pi m}_e}\right)}^{1/2}$$

(18)

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, $\widehat{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_{eff}$) 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 $Te$, 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_{min}$ and $X_{max}$, respectively, for the minimum and maximum of a general variable X of the state vector $\widehat{X}$ of the problem, the two different expressions can be considered:

$$X=X_{min}+(X_{max}-X_{min})\ast rand$$

(19)

with $X=\{V_p,T_e,T_{e,o},T_s,L_{eff},T_{w,o}\}$ and

$$X=10\exp (\log \left(X_{min}\right)+[\log \left(X_{max}\right)-\log \left(X_{min}\right)]\ast rand)$$

(20)

with $X=\{n_{n,e},n_{e,e},n_{n,o},n_{e,o}\}$. Here $rand$ 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:

$${\left(V_X\right)}_z={\left(X_{min}\right)}_z+[{\left(X_{max}\right)}_z-{\left(X_{min}\right)}_z]\ast rand$$

(21)

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=\parallel {ϵ}_I\parallel +\parallel {ϵ}_{Pow,e}\parallel +\parallel {ϵ}_{ion,e}\parallel +\parallel {ϵ}_{press,e}\parallel +\parallel {ϵ}_{Pow,o}\parallel +\parallel {ϵ}_{ion,o}\parallel +\parallel {ϵ}_{press,o}\parallel$$

(22)

where:

• ${ϵ}_I$ = RHS - LHS of the Equation 9b, (current)
• ${ϵ}_{Pow,e}$ = RHS - LHS of the Equation 9c (emitter power)
• ${ϵ}_{ion,e}$ = RHS - LHS of the Equation 9a (emitter ions)
• ${ϵ}_{press,e}$ = RHS - LHS of the Equation 9d (emitter pressure)
• ${ϵ}_{Pow,o}$ = RHS - LHS of the Equation 2b (orifice power)
• ${ϵ}_{ion,o}$ = RHS - LHS of the Equation 2a (orifice ions)
• ${ϵ}_{press,o}$ = RHS - LHS of the Equation 2c (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 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$$

(23)

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

$$G{B}^q=min\left(J_z^q\right)$$

(24)

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}=C1\ast V_z^q+C2\ast rand\ast (P{B}_z^q-X_z^q)+C3\ast rand\ast (G{B}^q-X_z^q)$$

(25)

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

(26)

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=min\left(J_z^{q-1},J_z^q\right)$$

(27)

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=min\left(G{B}^{q-1},P{B}_z^q\right)$$

(28)

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_{target}$ 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_{tot}$, as in Equation 29. In fact, in this case, the optimization process has to minimize also the difference between the total voltage $V_{tot}$ and the target value $V_{target}$.

$$J_z^q=\parallel {ϵ}_I\parallel +\parallel {ϵ}_{Pow,e}\parallel +\parallel {ϵ}_{ion,e}\parallel +\parallel {ϵ}_{press,e}\parallel +\parallel {ϵ}_{Pow,o}\parallel +\parallel {ϵ}_{ion,o}\parallel +\parallel {ϵ}_{press,o}\parallel +\parallel V_{tot}-V_{target}\parallel$$

(29)

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.

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	$eV$
$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	$eV$
$T_s$	1700	1950	K
$L_{eff}$	1 × ${10}^{-3}$	6 × ${10}^{-3}$	m
$T_{w,o}$	1700	1950	K

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	$eV$
$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	$eV$
$T_s$	1700	1950	K
$L_{eff}$	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 $\widehat{x}\pm \sigma$ curves, in which the PSO solutions lay. More in detail, the mean $\widehat{x}$, and the variance ${\sigma }^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:

$$\widehat{x}=\frac{1}{n}\sum _{z=1}^n{x}_z$$

(30)

$${\widehat{\sigma }}^2=\frac{1}{n}\sum _{z=1}^n\left(x_z^2-{\widehat{x}}^2\right)$$

(31)

Hence, the standard deviation $\sigma =\sqrt{{\widehat{\sigma }}^2}$ is used to generate the confidence bound ($\widehat{x}\pm \sigma$) 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.

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.

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.

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.

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.