Modelling an Argon Propellant Dual Stage Cylindrical Electrohydrodynamic Thruster at Low Pressure

In recent years, electric propulsion reached a high level of impact for space orbital maneuvers, planetary journeys, or deep space missions mostly because of its significant exhaust velocity, and its high specific impulse. In this study, we model a corona discharge to describe the behaviour of several key parameters in dual-stage electrohydrodynamic (EHD) thrusters investigating two different configurations. It was found that, between the two geometries, the three-electrode geometry was the better method since the four-electrode geometry would create a slowdown area, thus decreasing the output thrust. By setting the first cathode with a negative voltage and the second connected to the ground, the ions and electrons would be accelerated in the first cathode and would be neutralized, by virtue of secondary electron emission taking place in the second. Overall, the dual stage three electrodes EHD thruster spends 8mW of electric power and delivers a thrust of 157nN with an efficiency of 32mN/kW.


Introduction
Electric Propulsion (EP) is an alternative type of rocket propulsion in free space since it is more efficient at long operation manoeuvring in comparison to the old and conventional chemical propulsion. It has become a well-established technology for moving small satellites, such as CubeSats and spacecrafts in space (Levchenko et al. 2018). Electrically powered spacecraft propulsion systems make use of electrical energy to accelerate a propellant by different electrical and/or magnetic means. Plasma thrusters are usually classified into three main categories, according to their thrust generation processes: i) electrothermal; ii) electrostatic; and iii) electromagnetic. These three techniques, along with the associated plasma discharge and energy transfer mechanisms, include long-standing technologies, such as arc-jet thrusters, magnetohydrodynamics thrusters, ion engines, Hall thrusters and variants thereof (Keidar et al. 2015). Two-stage thrusters are an interesting concept and have been the subject of research as they can adjust thrust and specific impulse (two of the most important figures of merit in the subject of propulsion) and thus be potentially a more efficient alternative than the usual single-stage thrusters, although this has not yet been fully accomplished. The basic idea is to add an ionization stage, that is, an independent plasma source upstream of the acceleration stage. The underlying theory of Electrohydrodynamic (EHD) phenomena relies on the fact that charged species drift due to an applied electric field, as a consequence of a voltage drop between and anode and a cathode, and so doing they transfer momentum through collisions to the surrounding neutral species. Thus, a force, called the thrust, is exerted on the space vehicle through the ejection of, usually, a working gas (air, xenon, argon, nitrogen, or others) generally containing high kinetic energy. This force was first observed and studied by Thomas Townsend Brown and by Paul Alfred Biefeld when they observed, in 1929, the appearance of a net mechanical force in a Coolidge X-ray tube when its asymmetrical electrodes were applied to a high-voltage source, becoming known as the Biefeld-Brown effect (Fylladitakis, Theodoridis, and Moronis 2014). An electrical engine, such an EHD thruster, with a partially ionized plasma as the working gas, is fuel-efficient but at the expense of electrical power to sustain the plasma formation and ensuing gas acceleration. Usually, the design of a thruster is challenging, since its fundamental physics is very complex due to the combination of fluid mechanics, plasma physics, and chemistry. So, to produce an EHD thruster it is necessary to predict the flow pattern from the electrode geometry and the electric field profile, but it is important to have in mind that the discharge regime, and structure, are also closely connected to the electric field, and to the shape of the electrodes.

Materials and Methods
To develop the study of electrodynamics thruster, we had proposed a numerical model which comprises the coupling of several different modules. The complexity of the model derives from the underlying physics, which combines fluid mechanics, plasma physics and chemical reactions. Normally, to develop, and optimize this type of thruster, the programmer must predict the flow pattern inside its geometry, from the anode to the cathode, and the electric field profile. Yet, it is also important to redesign the discharge regime, which could lead to unwanted results, or even to errors during the simulation (Bedolla et al. 2017). To combine the appropriated physical environments, it was used a finite-element method (FEM), sub-dividing the computational domain into a mesh of points, where the differential equations were solved locally, considering the appropriate boundary conditions. The software utilized was the COMSOL Multiphysics® which provides the coupling of all these modules and solves all the differential equations with a time-dependent solver, and it can lead to a steadystate solution for each case study. In the Low-Earth-Orbit environment, with altitudes ranging from 20 to 100 km, pressure from roughly 40 to 0.25 Torr (1 Torr ≃ 133,3 Pa), and the ambient temperature from 190 to 270 K, the design of the electric propulsion system must be able to operate within these variations. Figure 1 shows the geometry of the thruster, with a needle-type anode and a hollow cylindrical cathode with variable inner radius and height with a few mm, and with 3 kV to 20 kV of the potential difference between the anode and the cathode (Granados 2018). The propellant gas used in this study was the noble gas argon ( ), with a pressure and temperature assumed to be 10 and 300 , respectively (Granados, Pinheiro, and Sá 2016). To control and stabilize the corona discharge, an external circuit with ballast resistor = 500 Ω and ballast capacitor = 1 , were connected. When modelling a plasma, usually space charge density and electric effects intervene with an atomic and molecular collisions set, using collisional cross sections, rate coefficients and swarm transport data from the literature. It is crucial to select a consistent set of reactions to attain a stabilized simulation of the real plasma.

Electrostatic model for electric field and potential distribution
In this model, the electric potential distribution, , under the presence of a space-charge distribution (r), is computed by solving the Poisson equation where is the plasma permittivity (ε = ε 0 ε ). Consistently in the same time step, the spacecharge density is computed considering the plasma chemistry and the species charged densities by means of where is the electron charge, is the ion charge number, the ion density for species and the electron number density. Finally, the calculation of the electric field, E is done under the relation between this physical quantity and the electric potential distribution

Plasma model for species spatial distributions
To simulate the plasma behavior under EHD body forces, we must consider the migration of all the particles species involved, such as the electrons, ions, and neutrals, since 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, ; the continuity equation for the electron energy density ε and the continuity equation for the gas particles, which are calculated, using the mass fraction of the -th species, ω (see COMSOL manual) (COMSOL, n.d.).
Here u is the main fluid velocity; ρ is the main gas mass density; Γ is the electron flux, Γ is the electron energy flux, Γ is the diffusive flux vector; finally, , and are the electron density, electron energy density of the -th specie density sources respectively. The second term in left-hand side of the set Equations (4-6) appears from the fact that COMSOL presents a generalized model for various plasma simulations. So, for example, when modelling fusion plasmas, this drift-diffusion approximation is not well suited, since it is required that the plasma must also be collisional, which means that the mean free path between electrons and the breakdown gas must be inferior to the characteristics length of the plasma (Granados, Pinheiro, and Sá 2017b).
To determine the , and source terms, it is necessary to know the chemical reactions that are generated inside the discharge reactor to determine their rate coefficients . Since we assumed that each gas owns several heavy species, the rate coefficients of their collision processes, such as inelastic collisions (excitation and ionization) or elastic collisions, are needed. To obtain the rate coefficients from these processes, we can compute the collision processes for electron collisions, using a set of the respective cross section data and the Maxwellian distribution function for electron energy ( ) through where is the electron energy in ( ); and ( ) represents the cross-section data for the considered reaction j. In Table 1, it is possible to observe how these coefficients were calculated. Note that for electron-impact collisions who have rate coefficients designed as ( ), their rate coefficients are computed from the available cross-section data. To reduce the computational time and increase efficiency in our simulation, the higher excited states of energy were assumed to decay instantly into the representative states of * (Fylladitakis, Theodoridis, and Moronis 2014).

Fluid model and coupling with electric field
To describe the main fluid motion of the EHD thruster, we must consider few assumptions.
The first one is that the fluid is incompressible, and the viscous flow will have, in a first approximation, a laminar behaviour, which can be described by a the Navier-Stokes fluid equation. Under an external force density term, this equation can be written in the wellknown form of u + (u ⋅ ∇)u = −∇ I + ∇ ⋅ [ (∇u + (∇u) )] + f where ρ is the fluid density, u is the main fluid velocity, is the absolute pressure and I is the identity of a symmetric tensor of rank two; is the dynamic viscosity, which are a specific property of the fluid, and finally f is the external force term, which as the form of: where is the Boltzmann constant and the variables and are, respectively, the ion and electron temperatures, given in electron-volt ( ) (Boeuf and Pitchford 2005;Boeuf et al. 2007).

Discussion
The design of the two-stage thruster separates the control of ionization and acceleration regions. Therefore, we propose and study two different concepts, one by introducing four electrodes with two different sets of a high voltage anode and a grounded cathode and another geometry that consists of three different electrodes: a high voltage anode, a negative voltage cathode, and a cathode ground. The design of the two-stage thruster separates the control of ionization and acceleration regions. Therefore, we propose and study two different concepts, one by introducing four electrodes with two different sets of a high voltage anode and a grounded cathode and another geometry that consists of three different electrodes: a high voltage anode, a negative voltage cathode, and a cathode ground.

Four-electrodes geometry
We started by introducing a new set of high voltage-anode and ground-cathode above the first set. The new electrode's group geometry was designed with the same dimensions as the previous one, with an anode to anode distance of 12 mm. Furthermore, we simulated several conditions for the applied voltages in each anode, either the first anode with higher input voltages than the second anode. The gathered results also showed that for the case where the first anode voltage was lower than the second one, a vortex would rise between the first cathode ground and the second anode, because of the electric potential distribution. On the other side, the results that converge without any type of vortices would be the cases where the first anode would have higher or equal voltage as the second anode. Then, we proceed to discover what should be the best potential difference set that would bring the best thruster's parameters. After several attempts, with thrust values in the range of the nano-newton, we found out that the best setting would be 20 kV / 0 kV / 20 kV / 0 kV with a secondary electron emission coefficient, , of 9 × 10−3 in each cathode. Although this was the best case, it is possible to conjecture from the electric potential distribution, displayed in Figure 2, that this type of geometry was not optimized mainly due to the potential gradient area between the first cathode and the second anode. Here we see that the ion particles suffer a great loss of momentum since they enter a deceleration region, where the electric field turns contrary to the direction of movement. Despite not having the optimal electric potential desired, this thruster compensates by delivering local ionization degrees of 10−7. In fact, both electrons and ions had their maximum peak concentration inside of the second hollow cathode. Electrons scattered in the middle region of this cylindrical electrode, reaching values of ≈ 1.1 × 10 16 m −3 , whereas the ions have its highest values near the second cathode's surface and distributed almost to the cathode's exit. They can reach concentrations above the 1.1 × 10 16 m −3 . The two-dimensional plasma velocity profile shows that the plasma reaches velocities near the first cathode of 30 / but this reduces slightly to the 20 cm/s at the entrance of the second cylindrical cathode. In this region, the velocity profile drops again to velocities near the 10 cm/s. Overall, the deceleration region can be observed since between the first hollow cathode to the second needle anode, the plasma only reaches speeds of 5 cm/s. The results of cathode's exit velocity profile do not show the quadratic profile, instead they present a maximum peak velocity in regions between the middle and the cathode's surface, where the middle area shows a reduction of the plasma's speed. After evaluating the thruster's parameters, the dual stage argon EHD thruster with a four electrode configuration, can deliver up to 53 nN with an efficiency of 9 mN/kW in the best simulation case, pointing out that these results are not very promising since the propelling force is very low (in the order of a few tens of nano-newton), because, presumably, we were not in its optimal configuration and geometry.

Three-electrodes geometry
The second part of this study was dedicated to the development of a new electrode geometry with a three-electrode configuration: one high voltage anode, one negative voltage cathode and one grounded cathode. This part was developed using two different configurations: connecting a negative voltage to the first cathode and set the second cylindrical hollow electrode at the ground voltage; then we simulate the opposite situation. However, we found out that the second case would result in instabilities and simulation errors, due to the fact that the region that neutralizes de gas was in the first cathode and few ions would be used to accelerate the plasma.
Then we tested several different sets of applied voltages and we found out the thruster would respond very well with tensions of i) 18 kV / (−2) V / 0 kV and ii) 15 kV / (−9) kV / 0 kV. Nevertheless, this would only happen if the cathode's inner radius would be increased up to = 18 , otherwise the results didn't converge or form instabilities within the plasma. After a set of numerical experiments, the obtained electric potential shows that this geometry set is the right choice, see Figure 3 and Figure 4. In direct comparison with the four-electrode configuration, we see that the three electrodes develop an electric potential distribution more favorable to the acceleration of the ion charged particles. We also found out that by reducing the negative voltage down to −9 kV, an acceleration region appeared, that would be neutralized in the second anode, as we can find in several grid ion engines configurations (Sangregorio et al. 2018). The particles distribution as different patterns for each voltage. In the first case (Figure 3, left), we saw electrons scatter along the middle of the thruster, and although it reached its maximum peak of 7 × 10 15 m −3 , it presents a uniform distribution from the anode to the second cathode. On the contrary, for the second case (Figure 4, left), the electrons form local peaks inside the middle of the first cathode, with concentration values of 1.9 × 10 17 −3 , two orders of magnitude above. Meanwhile in the second one (Figure 4, right), we observe the same local concentrations that appear in the electrons, reaching also values of 1 × 10 17 m −3 . From this point we can perceive that the ionization degree, in these regions is near 10 −6 , the highest value of this parameter, reached in all of these simulations. Since the particles had distinct profiles for each voltage set, it would be expected that the velocity profiles would be different. In fact, the set with a negative voltage −2 , presents a more uniform velocity profile along the three electrodes. After the anode, the plasma increases its speed up to 10 / within the first anode, decreases a little by the end of the first cathode and it accelerates again when the plasma reaches the second cathode, reaching is maximum velocity of 17 / . Then the plasma decreases its velocity, due to a stationary region inside the second cathode. The left side shows the electric potential distribution in , and the right side presents the two-dimensional velocity profile of the gas in / On the opposite side, when we applied a negative voltage of −9 kV (Figure 4), the twodimensional plasma velocity profile becomes more dynamic, with the appearance of two stationary regions, one at the exit of the first cathode and the other at the entry of the second cathode. Starting from the anode, the plasma speeds up to its maximum speed, with values of 20 cm/s, then it decreases when reaching the first electrode stage. Finally, at the exit of the dual-stage thruster the plasma's velocity speeds to 11 cm/s. The cathode's exit velocity profiles are also very distinct for each set of negative voltages. We obtained that for the first case (not shown in the figures), the radial component of the velocity has a considerable presence and the axial component did not reach its peak in a quadratic profile, but instead it created a plateau right in the surrounding zone of the center axis of the cathode. However, for the second set of negative voltages (case ii) and also not shown in the figures) it is observed the known behaviour of a quadratic profile with a maximum speed of 11 cm/s. After evaluating the thruster's parameters for each case, we discovered that the dual stage argon EHD thruster with a three-electrode configuration can deliver up to 158 with an efficiency of 20 mN/kW for the applied voltage set 18 kV / (−2) kV / 0 kV and 215 nN with an efficiency of 32 mN/kW for the applied voltage set 15 kV / (−9) kV / 0 kV. This allows us to assess that the three-electrode geometry produces better results than the four-electrode configuration. Yet, this part of the double stage thruster study is not fulling improved, lacking a deep investigation of the physics involved and the geometrical configuration.

Conclusions
An argon propellant dual-stage EHD thruster was modelled using finite element methods software. In this study, it was simulated two different configurations, one with four-electrodes geometry and the other with three-electrodes geometry. We concluded that the best one was the one with three-electrodes since it developed distinct regions of acceleration and neutralization, that are essential for the development of the thruster. Overall, by using a set of applied voltages of 15 kV / (−9) kV / 0 kV, the improved thruster delivered 215 nN with an efficiency of 32 mN/kW. This dual-stage EHD thruster is not as efficient as the others developed previously with only one stage (Granados, Pinheiro, and Sá 2017a, 2017b, 2016. We plan in the near future to pursue further investigation of the physical processes and the geometry of the configuration of these kind of dual-stage EHD thrusters. Also, in order to continue investigating on these topics, for future work we propose introducing a magnetic field across the cylindrical cathode to capture more electrons inside the cathode. This concept is already used in Hall-effect thrusters and potentially, the development of this thrust could be the future line of work as well.