Modeling of Electrohydrodynamic (EHD) Plasma Thrusters: Optimization of Physical and Geometrical Parameters

Abstract

This work aims to optimize a previous self-consistent model of a single stage electrohydrodynamic (EHD) thruster for space applications. The investigated parameters were the thruster performance (propulsion force T, the thrust to power ratio T/P, the electric potential distribution, the spatial distribution for the electrons and ions, and the laminar flow velocity) under several conditions, such as the design features related to the cathode’s cylindrical geometry (height and radius) and some electric parameters such as the ballast resistor, and the applied potential voltage. In addition, we examined the influence of the secondary electron emission coefficient on the plasma propellant parameters. The anode to cathode potential voltage ranges between 0.9 and 40 kV, and the ballast resistance varies between 500 and 2500 M. Argon and xenon are the working gases. We assumed the gas temperature and pressure constant, 300 K and 1.3 kPa (10 Torr), respectively. The optimal matching for Xe brings off a thrust of 3.80 μN and an efficiency T/P = 434 mN/kW, while for Ar, T = 2.75 μN, and thruster to the power of 295 mN/kW. To our knowledge, the missing data in technical literature does not allow the verification and validation (V&V) of our numerical model.

1. Introduction

Electric Propulsion (EP) is an alternative type of rocket propulsion with higher efficiency at long operation maneuvering when compared to conventional chemical propulsion. It has become a well-established technology in propulsion devices for small satellites, such as Cubesats [1]. The electrohydrodynamic effect (EHD) is related to the possibility of generating or altering a gas flow, also known as ionic wind, through the action of an electrical discharge, often a corona discharge. An EHD thruster, while working with partially ionized gas, is fuel-efficient, though at the expense of electrical power to sustain the plasma formation and ensuing gas acceleration.

The design of an EHD thruster is challenging as its fundamental physics is complex, combining fluid mechanics, plasma physics, and chemistry. The EHD thrust relies substantially on the flow pattern resulting from the electrode geometry and the electric potential morphology best suited for optimized performance. According to Masuyama and Barret [2], plasma propulsion is potentially the most energy-efficient mechanism, about 50 times more efficient than jet engines. However, this technology requires maturity for interplanetary-manned missions based on becoming practical.

This report introduces the model and accompanying study of an electrohydrodynamic thruster, which has an externally applied electric field that induces an ionic wind—the flow of an electrically charged fluid—that transfers momentum to the neutral heavy species in a gas, thus generating a flowing stream. Particle momentum is transmitted to the ions and the neutrals, but more efficiently to the later ones due to the low degree of ionization. In its principles, the EHD thruster is a non-equilibrium non-thermal, corona discharge, consisting basically of an emitter electrode (or anode) and a collector electrode (or cathode). With its very sharp endpoint, the anode promotes the gas ionization in the high-field ionization region when a high dc voltage is applied, thus producing the thrust due to the oriented migration of the charged cloud essentially in the drift region. The EHD device moves in the opposite direction of the ion swarm, obeying the action-to-reaction law. Numerical and experimental studies of EHD propulsion systems show that a thrust to power ratio of 100 mN/W can be achieved, despite an energy conversion efficiency of less than one. Other applications can be engineered, such as boundary layer enhancement, fluid pumping (including EHD micropumps), flow, heat transfer improvement, and drying and evaporation (see, e.g., Refs. [3,4,5]).

There are several theoretical, modulation, simulation, and experimental works that study the creation of an ionic wind resulting from the electrohydrodynamic force (EHD) created by a corona discharge, some with applications in electric propulsion (see, e.g., Refs. [6,7,8,9,10,11,12,13,14]). In 1967, Christenson and Moller [15] investigated ionic wind as a propulsion technique. A NASA technical note [16] examined how the thrust produced by the ionic wind can scale to values of interest in aeronautical propulsion—with no moving parts, virtually silent, and no fuel emission. With inconclusive results about possible advantages, the authors concluded that corona discharges are not practical in aeronautical propulsion. However, recently, Xu et al. [17] claim to have provided proof of aeronautical propulsion based on electro-aerodynamics, opening doors to avionics and aerodynamic devices that are quieter, mechanically simpler, and that do not produce fuel emission.

Bedolla et al. [18] simulated the EHD flux induced by a corona discharge, using an asymmetric high voltage capacitor to study the dependence of the thrust on parameters such as altitude, pressure, temperature, and humidity. Park et al. [19] used a source of weakly ionized gases generated in helium at atmospheric pressure to study the flow paths of different discharge parameters. The authors aimed to distinguish the effects of streamer propagation or drift of space charge that gives rise to the ionic wind and determine the role of electrons and (positive) ions in generating this wind. For a possible application to electric propulsion, Moreau et al. (see Refs. [20,21]), used and compared the EHD force experimentally produced, in air and at atmospheric pressure, by a corona discharge between a wire-shaped electrode and three cylindrical electrodes. The electrohydrodynamic force and the acceleration of the aerodynamic flow have also been the subject of theoretical and modulation studies in dielectric barrier discharges plasma actuators and surfaces (see Refs. [22,23]). However, none of these articles make simulations or experimental verifications that may be compared with those referred to in this article.

This paper aims to model and optimize the physics and geometric parameters of the EHD electrostatic propulsion devices for aerodynamic applications and is organized in the following parts: In Section 2, we briefly explain the model developed to analyze the behavior of an EHD thruster. Section 3 describes the chosen geometries and the results of the performed simulations, including the optimal thrust and thrust-to-power for the set of parameters. Finally, Section 4 comprises the findings of our work and conceivable guidelines for future research on the topic.

2. Electrohydrodynamic Model

A self-consistent model based on the corona discharge mechanism has been proposed in [24,25,26] to describe the behavior of several key parameters of a specific EHD thruster. Figure 1 shows the geometry of the thruster and the RC coupling circuit, with a needle-type anode and a hollow cylindrical cathode with variable height and internal radius. The applied potential difference between the anode and the cathode varies in the range 0.9–40 kV, the anode grounded and the cathode potential running between those values. The gas pressure and temperature assumed were 10 Torr (1.3 kPa) and 300 K, respectively.

To control and stabilize the corona discharge, an external RC series connects with a varying ballast resistor R = 500 − 2500 MΩ. This parameter is necessary to ensure the discharge does not reach the arcing region, while influencing the total current between electrodes and, consequently, the net thrust produced. A blocking capacitor of C = 1 pF avoids peaks between electrodes.

The input voltage, $V_{in}$, and the electrode voltage, $V_{plasma}$, are related by the circuit equation:

$$V_{plasma}=V_{in}-I_{P}R-RC\frac{d{V}_{plasma}}{dt}$$

(1)

where $I_P$ represents the current between the electrodes and includes the electron and ion current densities at the electrode walls.

The power $P$ dissipated to create and maintain the plasma is given by:

$$P=V_{plasma}{I}_P$$

(2)

The model used is well described in Granados et al. [24] but, in this paper, we make a synthesis of the model for the sake of completeness and clarify some physical aspects. Figure 2 is a schematic illustration of the pseudo-flow chart of the simulation, displaying the link between the several modules and entry parameters used in the EHD thruster. The commercial software COMSOL Multiphysics^®, version 5.3, based on finite element methods solved the equations included in the fluid and electrical modules [27]. The simulation region consists of a 2D axisymmetrical cylindrical with a variable-element mesh refined at the surface of the electrodes.

The model consists of three main modules (see Figure 2):

(i)

a heavy-species kinetic module including ground-states and ionized species;

(ii)

a plasma module to obtain the self-consistent electric field and space-charge density, considering the interaction with the electrode walls (including secondary electron emission), using an external electric circuit to power the discharge;

(iii)

a fluid module solving the Navier–Stokes equation under the assumption of laminar flow and with a force volume term due to electron-molecule and ion-molecule collisions acting on the neutral gas.

2.1. Electrostatics

In this model, the electric potential distribution, $V$, under the presence of a space-charge distribution ${\rho }_c\left(r\right)$, is obtained by solving the Poisson equation

$${\nabla }^2\mathrm{V}=-\frac{{\rho }_c\left(r\right)}{\epsilon },$$

(3)

where $\epsilon$ is the plasma permittivity. Consistently, in the same time step, the space-charge density is computed considering the plasma chemistry and the charged species densities by

$${\rho }_c\left(r\right)=e\left({{\displaystyle \sum }}_{j=1}^N{Z}_j{n}_j-n_e\right),$$

(4)

with $e$ the absolute value of the electron charge, $Z_j{n}_j$ is the ion charge number, $n_j\left(r\right)$ the ion density for species and $n_e\left(r\right)$ the electron number density. Finally, the calculation of the electric field, $E$, is performed under the relation between this physical quantity and the electric potential distribution

$$E=-\nabla V$$

(5)

The Poisson equation is solved self-consistently, giving the electric field at each point in space, showing the characteristics of the corona discharge produced, namely the rapid variation of the electric potential in the vicinity of the electrodes and the verification of quasi-neutrality in the remaining of the thruster volume:

$$\nabla ·J=0$$

(6)

where J is the total current density vector, including ions plus electrons current densities vectors, $J_i$ and $J_e$, respectively.

2.2. Plasma Kinetics

To simulate the plasma behavior under EHD body forces, we consider the migration of all the particles species involved, such as the electrons, ions, and neutrals, as they are fundamental to the calculation of the space charge density. The set of equations that are responsible for the evolution of each species are the continuity equation: for the electron density, $n_e$; for the electron energy density, $n_{\epsilon }$, and for the gas particles, which are calculated, using the mass fraction of the k-th species, ${\omega }_k$ [27]

$$\frac{\partial }{\partial t}{n}_e+\left(u·\nabla \right)n_e=-\nabla ·{\Gamma }_e+R_e$$

(7)

$$\frac{\partial }{\partial t}{n}_{\epsilon }+\left(u·\nabla \right)n_{\epsilon }=-\nabla ·{\Gamma }_{\epsilon }+R_{\epsilon }$$

(8)

$${\rho }_f\frac{\partial }{\partial t}{\omega }_k+\left(u·\nabla \right){\omega }_k=-\nabla ·{\Gamma }_k+R_k$$

(9)

Here, $u$ is the gas fluid velocity;${\rho }_f$is the gas mass density; ${\Gamma }_e$is the electron flux, ${\Gamma }_{\epsilon }$ is the electron energy flux, and ${\Gamma }_k$is the diffusive flux vector. Finally, the densities sources are$R_e$, $R_{\epsilon }$ and $R_k$ for the electron, electron energy, and the $k$-th species, respectively [24].

Table 1. Electron-impact and chemical reactions with rate coefficients for argon and xenon discharges. Reactions with a ($∎$) have units of ${\mathrm{m}}^3/\mathrm{s}$; (+) in units of ${\mathrm{m}}^6/\mathrm{s}$, and (*) in units of ${\mathrm{s}}^{-1}$.

Reaction	Rate Coefficient	References
Argon
$e^{-}+Ar\to e^{-}+Ar$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+Ar\to e^{-}+A{r}^{*}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+Ar\to 2e^{-}+A{r}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+A{r}^{*}\to 2e^{-}+A{r}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+A{r}_2^{+}\to 2e^{-}+A{r}^{+}$	$5.06\times {10}^{-15} T_e^{-0.67}(∎)$	[29]
$A{r}^{*}+A{r}^{*}\to e^{-}+Ar+A{r}^{+}$	$6.2\times {10}^{-16}(∎)$	[30]
$A{r}^{*}+A{r}^{*}\to e^{-}+2Ar$	$3.0\times {10}^{-16}(∎)$	[30]
$A{r}^{+}+2Ar\to A{r}_2^{+}+Ar$	$2.25\times {10}^{-43}$(+)	[31]
Xenon
$e^{-}+Xe\to e^{-}+Xe$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+Xe\to e^{-}+X{e}^{*}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+X{e}^{*}\to e^{-}+Xe$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+Xe\to 2e^{-}+X{e}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+X{e}^{*}\to 2e^{-}+X{e}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$X{e}^{+}+Xe\to Xe+X{e}^{+}$	$3.6\times {10}^{-10}(∎)$	[33]
$X{e}^{*}\to Xe$	$3.96\times {10}^{6 }$(*)	[33]
$X{e}^{*}+X{e}^{*}\to e^{-}+Xe+X{e}^{+}$	$5.0\times {10}^{-16}(∎)$	[33,34]

Table 1. Electron-impact and chemical reactions with rate coefficients for argon and xenon discharges. Reactions with a ($∎$) have units of ${\mathrm{m}}^3/\mathrm{s}$; (+) in units of ${\mathrm{m}}^6/\mathrm{s}$, and (*) in units of ${\mathrm{s}}^{-1}$.

Reaction	Rate Coefficient	References
Argon
$e^{-}+Ar\to e^{-}+Ar$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+Ar\to e^{-}+A{r}^{*}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+Ar\to 2e^{-}+A{r}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+A{r}^{*}\to 2e^{-}+A{r}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[28]
$e^{-}+A{r}_2^{+}\to 2e^{-}+A{r}^{+}$	$5.06\times {10}^{-15} T_e^{-0.67}(∎)$	[29]
$A{r}^{*}+A{r}^{*}\to e^{-}+Ar+A{r}^{+}$	$6.2\times {10}^{-16}(∎)$	[30]
$A{r}^{*}+A{r}^{*}\to e^{-}+2Ar$	$3.0\times {10}^{-16}(∎)$	[30]
$A{r}^{+}+2Ar\to A{r}_2^{+}+Ar$	$2.25\times {10}^{-43}$(+)	[31]
Xenon
$e^{-}+Xe\to e^{-}+Xe$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+Xe\to e^{-}+X{e}^{*}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+X{e}^{*}\to e^{-}+Xe$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+Xe\to 2e^{-}+X{e}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$e^{-}+X{e}^{*}\to 2e^{-}+X{e}^{+}$	$f\left[\sigma \left(\epsilon \right)\right]$	[32]
$X{e}^{+}+Xe\to Xe+X{e}^{+}$	$3.6\times {10}^{-10}(∎)$	[33]
$X{e}^{*}\to Xe$	$3.96\times {10}^{6 }$(*)	[33]
$X{e}^{*}+X{e}^{*}\to e^{-}+Xe+X{e}^{+}$	$5.0\times {10}^{-16}(∎)$	[33,34]

The source terms are obtained with the chemical reactions and associated rate coefficients. The set of kinetic reactions considered is shown in Table 1. The collisional processes of gain and loss of electron energy are adequately described using the Boltzmann equation [35], but as here we are only interested in the average influence of plasma parameters on thruster performance, we use a Maxwellian electron energy distribution function (EEDF), $f_M\left(\epsilon \right)$, for electron energy to find the coefficients for the electron–impact collision processes (elastic, inelastic, and ionization), with a set of cross-section data (for references, see Table 1), through

$$C_j={\left(\frac{2e}{m_e}\right)}^{1/2}{{\displaystyle \int }}_0^{\infty }\epsilon {\sigma }_j\left(\epsilon \right)f_M\left(\epsilon \right)d\epsilon$$

(10)

where $\epsilon$ is the electron energy (in eV units), $m_e$ is the electron mass, and ${\sigma }_j\left(\epsilon \right)$ represents the cross-section data for the considered reaction. As this information regarding electrons is necessary to predict the populations of the various species present in the corona discharge, the Maxwellian EEDF must be solved simultaneously with a system of rate balance equations for the various heavy species of particles (neutrals and ions) considered. Although the use in this work of a Maxwellian-type EEDF is a simplified form when compared to a more realistic EEDF obtained from solving the Boltzmann integro-differential equation, it allows us to obtain computational solutions faster in the framework of a two-dimensional finite element model with a dense mesh.

2.3. Governing Equations of Fluid Dynamics

We assume the fluid is incompressible, and the viscous flow is laminar, described by the Navier–Stokes fluid equation. As is well known, this equation is a balance of forces acting on the fluid. Therefore, with an appropriate external force density term, the equation is written in the form

$${\rho }_f\frac{Du}{Dt}=-\nabla pI+\nabla ·\left[\eta \left(\nabla u+{\left(\nabla u\right)}^T\right)\right]+f$$

(11)

where $p$is the absolute pressure, and $I$ is the identity symmetric rank two tensor; $\eta$ is the dynamic viscosity, which is a specific property of the fluid, and finally $f$ represents the external force density term. The left side of this equation is the total derivative of fluid momentum and represents a sum of inertial forces and convection forces; in addition, the first term on the right-hand side accounts for pressure forces (on the approach of pressure gradient), the second term represents the viscous forces, and the third term is the external electric force per unit volume acting on the charged particles that communicate momentum, per unit volume and per unit time, to the neutrals by the agency of electron–atom or ion–atom collisions. The electric current, represented by a drift-diffusion approximation form of the charged particle fluxes, appears in the Navier–Stokes equation through the force density, written here as in Ref. [22]

$$f=e\left(n_i-n_e\right)E-k_B{T}_i\nabla n_i-k_B{T}_e\nabla n_e$$

(12)

with $k_B$denoting the Boltzmann constant, and $T_i$ and $T_e$ the ion and electron temperatures, respectively. Equation (12) is a simplified form of the momentum–transfer equation (see, e.g., Ref. [36]). Eventually, other components of this force appear if a sophisticated form of the momentum transfer equation is used (see, e.g., Refs. [7,8]). These additional terms may be necessary under certain conditions but are not considered in this work for the sake of simplicity and because Equation (12) contains the more relevant terms (see, e.g., Ref. [37]). This equation for $f$ can be significant in some specific regions of the discharge and be practically irrelevant in others. Note that in a region where the plasma is perfectly neutral ($n_e=n_i$) and uniform ($\nabla n_e=\nabla n_i$), the force density term, Equation (12), is zero: the force exerted on the neutral particles due to the ions balances exactly that exerted by electrons. This is the result that, although the forces are proportional to the masses of the charged particles, the electron flux and electron-neutral collision frequency have a much higher value, respectively, than the ion flux and ion-neutral collision frequency.

2.4. EHD Thrust

The force provided to a thruster device, named the thrust, results from the ejection of propellant from the engine exhaust into space and is the time rate of change of the momentum of the propellant. The thrust $T$ is defined as: $T=\dot{m}{v}_{ex}$ where $\dot{m}$ is the mass flow rate and $v_{ex}$ the exhaust propellant velocity. Then, the thrust is created by the expulsion of a propellant mass.

The axial component of the steady-state mass flow velocity is $v_z\left(r\right)$. An infinitesimal element of mass crossing in a time $dt$ the circular area of the cathode outlet surface is $dm={\rho }_f\left[v_z\left(r\right)dt\right]dS$, where $dS=2\pi rdr$, and $dm/dt=2\pi {\rho }_f{v}_z\left(r\right)dr$. As the propulsion force through this surface element is $dT=v_z\left(r\right)dm/dt$, the total thrust value delivered by the fluid passing is estimated by means of the following equation

$$T=2\pi {\rho }_f\underset{0}{\overset{R}{{\displaystyle \int }}}r{v}_z^2\left(r\right)dr$$

(13)

where $v_z\left(r\right)$is the axial component of the steady-state velocity at the cathode’s exit; r is the distance from the center of the thrust symmetry axis (r = 0) to the cathode’s wall, located at r = R.

3. Influence of the Discharge Current

The starting point for increasing the cylindrical thruster was to adjust the degree of ionization, a physical quantity that expresses the proportion of charged particles in a neutral plasma, for it was verified that the degree of ionization, displayed in the study previously mentioned [24,25,26] present very low values, attaining the maximum of $2.4\times {10}^{-9}$. To increase the number of ionized particles, it is required to apply a stronger current between the electrodes. This is achieved by two methods, diminishing the ballast resistor with fixed applied voltage, or by raising the applied voltage at the fixed ballast resistor.

3.1. Variation of the Ballast Resistor

Maintaining the corona self-sustained discharge regime, we identified the influence of the ballast resistor in the argon single-stage EHD thruster, and we uncovered that, by decreasing the resistor value, we increased the number of charged particles, as was expected.

Both thrust and thrust-to-power ratios were computed directly by the software and shown in

Table 2. Simulation results for the study of the influence of discharge current in the single-stage argon EHD thruster by varying the ballast resistor for a fixed applied voltage of 3 kV.

Resistor MΩ	Current mA	Thrust nN	T/P Ratio mN/kW
2500	1.2	4.11	17.97
2000	1.5	5.46	19.44
1500	2.0	7.75	21.24
1000	3.0	12.27	23.31
750	4.0	16.75	24.55
500	6.0	25.86	26.27

. As can be observed, by changing the discharge current from 1 mA to 6 mA, the propellant force increases more than five times, from 4 nN to around 26 nN. Even with the increase in current, which consequently was accompanied by the power consumed on the plasma, the T/P efficiency ratio improved, highlighting the fact that there was a significant enhancement in the work regime of the thruster.

3.2. Variation of the Applied Voltage

Setting the most advantageous value of the ballast resistor in 500 MΩ, the impact of the applied voltage was examined. At this point, we were expecting that by increasing the voltage, the discharge current would also increase, producing a higher output thrust. We modulated our plasma to work in the range of low potential differences of 900 V between electrodes to a higher applied voltage of 40 kV.

One way that the applied voltage controls the operation of our thruster is the effect that it has on the electric potential distribution, which is a fundamental factor considering the electric field converts electrical to mechanical energy.

Table 3. Simulation results for the study of the influence of discharge current in the single-stage argon EHD thruster by varying the applied voltage for a fixed ballast resistor of 500 MΩ.

Voltage kV	Current μA	Thrust nN	T/P Ratio mN/kW
0.9	1.8	5.45	179.44
3	6.0	25.86	26.27
9	18	107.77	38.36
12	24	181.14	48.40
20	40	549.93	86.24
30	60	521.09	81.26
40	80	521.04	40.90

shows how the thruster parameters change with the applied voltage, and from them, we observe that the thrust reaches values of 550 nN, a 100-times increase from the previously developed thruster, reaching an efficiency of 86.24 mN/kW, the highest to date. However, we also found out that there was an optimal voltage of 20 kV, as higher voltages delivered the same thrust and reduce the thruster efficiency. For instance, at 30 kV, the thrust decreases more than 10 nN and at 40 kV, the thruster efficiency was halved.

This result is in agreement with that obtained by Pekker and Young [9] who showed, using a theoretical model, that an increase in the voltage between electrodes promotes an opposite behavior in terms of thrust and propeller efficiency: there is a decrease in efficiency of the thrust while there is an increase in the thrust; or, in other words, the greater the thruster current (greater propulsion force), the lower the thruster efficiency.

4. Influence of the Cathode’s Geometry

The conversion of electric to mechanical energy is made primarily by the electric field, making the potential distribution a key factor in the design of the thruster. In turn, the electric potential morphology relies on the electrode geometry. Therefore, the optimized performance of the EHD thrust, which includes the flow pattern, depends greatly on the thruster geometry.

In a first approach, we investigated the influence of the intrinsic cathode geometry on the thruster performance for two noble propellant gases, argon (Ar) and xenon (Xe). The first gas, because it was used in our previous studies [24,25,26] and so it serves well for comparison, between the results then obtained and the new ones; the second gas is here due to its widespread use in space applications. For example, xenon propellant is the most used in ion engines and Hall thrusters, despite its costs and availability, making other noble gases, such as argon and krypton (Kr), potential substitutes. However, the last gases, Ar and Kr, present more issues in terms of thruster performance specifications, as having a lesser molecular mass than xenon and need higher ionization energy. In general, a denser propellant will improve the output thrust and its efficiency as it has a higher energy per particle for the same exhaust velocity.

For the study of the influence of the cathode geometry on the thruster, the electric circuit developed to control and maintain the self-sustain plasma discharge (for both gases) were designed with 20 kV-applied DC voltage, a ballast resistor, and a blocking capacitor of 500 MΩ and 1 pF, respectively. As for the plasma fluid, the initial pressure was set at 10 Torr with an ambient temperature of 300 K.

4.1. Study of the Cylindrical Inner Radius

In terms of geometry, we began by changing the hollow radius of the cathode. Using Ar, we verified that the optimized width was 12 mm, as any other variation, upper or lower, would decrease the output thrust and the thruster efficiency. For xenon, we concluded that the optimal hollow radius was 20 mm, giving the best thrust and thruster efficiency (see Figure 3a,b).

4.2. Study of the Cylindrical Height

At 20 kV, with xenon as the working gas and starting from the initial cathode height value of 21 mm, we begin by reducing the cathode’s length down to 12.6 mm. These modifications can cause significant changes in the electric potential distribution that could improve our thruster. As we were reducing its length, we observed that values under 12.6 mm would not improve, but instead would unfold instabilities in the plasma. We conclude that at 20 kV, in the case of Ar, the cylinder height should be approximately 10.5 mm, and in the case of Xe, the cylinder height should be around 14.7 mm (see Figure 4a,b) to optimize the output. As we can see, the cathode’s height is not a determinant parameter in the efficiency of the EHD thrusters under study in the range considered.

With argon as the propellant gas, we verified that the width was already optimized in 12 mm, as any other variation, upper or lower, would decrease the output thrust and the thruster efficiency. On the other hand, when using xenon, we concluded that the optimal hollow radius was 20 mm. With the change of the cathode’s length, we found that the xenon propellant would attain its optimal state with a length of 14.7 mm. Overall, the xenon thruster delivered an output thrust of 3.80 μN and the thrust-to-power ratio of 434 mN/kW, values that when comparing with the current state-of-the-art of the modern electric propulsion systems reach the same order of magnitude.

5. Influence of Secondary Electron Emission Coefficient

The influence of the secondary electron emission coefficient, γ_i, in the cathode, was also investigated. This coefficient depends on the bombardment of the cathode surface by the kind and energy of particles (ions, neutrals, or electrons) from the gas and of the type of surface. We reported a decrease from the output thrust as γ_i increases, results consistent with the technical literature [38].

We plotted the thruster parameters such as the thrust versus γ_i and thrust-to-power ratio vs. γ_i in, respectively, Figure 5a,b. We observed that by reducing the electrons emitted by the cathode’s surface, the thruster reached an output thrust of 4 μN and efficiencies of nearly 460 mN/kW, whereas, for higher γ_i values, the thrust would only deliver a thrust of 0.24 μN and an efficiency of 80 mN/kW.

6. Gas Flow Velocity Distribution

Figure 5a,b present the flow optimal velocity distribution for the two gases utilized in our simulations. Notice that the scales are not the same to enhance images. In both figures, the voltage between the anode and cathode is 20 kV, the ballast resistance is 500 M$\mathsf{\Omega }$, and the block capacitor is 1 pF. For the argon thruster, Figure 5a, the cylindrical cathode has a height of 10.5 mm and an inside diameter of 24 mm, while in Figure 5b, for the xenon thruster, the dimensions of the cylindrical cathode are 14.7 mm in height, and the inside diameter is 40 mm.

The axial distance between the anode tip and the cathode inlet is about 30 mm on both thrusters in Figure 6a,b. A higher speed can be seen for argon (maximum value of 52.3 cm/s) than for xenon (maximum speed of 23 cm/s). Comparing the velocity distribution of gas flow for both thrusters, we observe for argon that the higher speeds can be found at points close to the entrance surface of the cathode, where the ions are more accelerated [24], and in the middle region of its output cross-section area.

For the xenon thruster, it is observed a similar behavior, although less pronounced, in the cathode inlet region, but with a more exuberant extension in the central region of cathode outlet cross-section. Another essential difference in these distributions at the cathodes operating with argon and xenon is that, in the latter case, the velocity distribution is higher only in a ring close to the entrance wall of the cathode, not forming the cylindrical configuration of high velocity existing when the working gas is argon.

Figure 7 shows the fluid velocity components (total, axial, and radial) as a function of the cathode diameter at the exit of the EHD thruster, where the similarity of the profiles for argon and xenon can be seen, as well as difference in the maximal values. We note that the values of the radial components are much lower than the axial.

7. Electric Field and Charged Particles Distributions

Figure 8 shows the electric potential morphology inside the thruster, with black arrows representing the electric field vectors at each grid point. For argon (Figure 8a), the electric potential attains the highest value in the central region with the field lines mostly confined inside the nozzle; the equipotential lines’ distribution is such that the electric field is radially predominant at the entrance of the cathode, which constitutes a less-than-ideal case for accelerating the ions in the axial direction.

For xenon (Figure 8b), the electric equipotential morphology extends further in the azimuthal direction, but a significant electric field radial component weakens the gas flow exhaust speed all along the height of the cylinder. Notice that the cross-sectional area of both cathodes is not equal, and although, for argon, the gas speed is higher than xenon, at the end the highest thrust is observed for xenon due to a higher mass and therefore a more effective momentum transfer.

The spatial distributions for electrons and ions for the two ionized gases are presented in Figure 9 and Figure 10, respectively. In the case of electrons, the spatial distributions in both working gases are different, but diffusion from anode to cathode is observed. In argon, electrons concentrate mostly near the cathode, while in xenon, electrons tend to the anode and are also constrained in a small central region inside the cylinder nozzle. In the case of the ions distributions, Ar⁺ concentrates near the cathode’s entry and for Xe⁺, they also tend to be constrained in the same small region inside the cylinder cathode as the electrons, and it appears a small concentration at the entrance of the cathode. This is a consequence of the electric potential morphology.

8. Conclusions

A parametric optimization of a single-stage EHD thruster was investigated. Based on an argon propellant EHD thruster developed by Granados [24], the parameters of the geometry of the cylindrical thruster, of the external control circuit, and the secondary electron emission at the cathode’s thruster were used for optimization. Adjusting the discharge current with the ballast resistor swapping from 2500 MΩ to 500 MΩ, and keeping the applied voltage constant at 3 kV, the thrust output increased 5 times; with a fixed ballast resistor of 500 MΩ, and altering the applied voltage from 0.9 kV to 40 kV, a maximum thrust of about 0.6 μN is attained for 20 kV, while in the range 20–40 kV, there is a stagnation of the thrust, whereas the T/P ratio decreased.

For argon and xenon, at fixed ballast resistor and applied voltage of 500 MΩ and 20 kV, respectively, the most performant geometrical dimensions were determined. The secondary electron emission coefficient appears to play a significant role, and for low values of γ_i, thrust and thrust-to-power ratios would significantly escalate. As it is expected, xenon is more efficient than argon, due to its higher atomic mass, and in the best cases of our work we obtained a thrust with xenon of 4 μN, while argon reaches 2 μN.

It is found that the flow velocity profile of the two gases, in the respective thrusters, as well as the distributions of charged particles (ions and electrons) are correlated with the behavior shown by the electric field/potential in the corona discharge region.

Although we have no information in the literature that permits us to verify and validate our model for the two noble gases studied, our results point to a promised future of the EHD thrusters in the astronautical industry, fundamentally for space missions. However, more fundamental research is necessary in the study of the geometry of the thruster, and in the understanding of its coupling with the physics of the electric field lines morphology inside the ion thruster discharge.