1. Introduction
Since the Wright Brothers achieved the first controlled human flight over a hundred years ago, aircraft within the atmosphere have been propelled by propellers and gas turbine engines, mostly powered by fossil fuels. The current trend is to replace chemical fuel propulsion systems with electric propulsion systems to reduce carbon emissions. In the aerospace field, electric propulsion technology is considered an alternative to fossil fuel technology [
1]. Ion thrusters have been widely and successfully used in space applications [
2]. In recent years, there has been a growing interest in Electrohydrodynamic (EHD) propulsion in this field, with initial applications in atmospheric flight. Ion propulsion, as an alternative to traditional engines, offers many advantages, such as no moving parts, a high thrust-to-power ratio, low noise, and high sustainability of the propulsion system [
3].
In 1928, Brown constructed a pair of asymmetric electrodes that generated directional propulsion in the air using high voltages (i.e., tens of kilovolts) [
4]. Two electrode rods with different radii of curvature constituted the propulsion structure; the emitter electrode had the lower radius of curvature, and the collector electrode rod had the larger radius of curvature. The asymmetric electrode structure produces thrust based on the principle that a high-voltage electric field generates ions via the ionization of neutral air molecules, and the field also accelerates the ions. Because of the isotropic motion, the accelerated ions collide with neutral molecules, exchanging momentum with them to drive the ions forward [
5]. The study of these dynamics is known as electroaerodynamics (EAD). EAD has been studied and utilized in the areas of gas pumping [
6,
7], flow control [
8,
9], and air purification [
10,
11].
Because EAD devices do not have moving elements, such as turbines and propellers, they have progressively garnered more interest; in fact, they may be able to replace other propellers used in aircraft propulsion. In 2009, researchers from NASA analyzed and demonstrated the use of plasma propulsion technology (EAD propulsion) for unmanned aerial vehicles (UAVs). They concluded that because of the technology’s low thrust, it was unlikely to be used as the primary thrust generator in UAVs [
12]. Pekker and Young [
13] concluded that the thrust density of plasma propulsion was insufficient for aircraft propulsion. However, research on plasma propulsion technology is still ongoing, and a number of researchers have pioneered the use of plasma propulsion in hot-air balloons and airships. For example, in 2014, Johnson et al. [
14] proposed a pulsed EAD thruster design for high-altitude, long-endurance airships. Similarly, in 2016, Wynsberghe et al. [
15] suggested using an EAD thruster to drive a stratospheric hot-air balloon. In 2017, Khomich and Rebrov [
16] incorporated a transmitter and receiver into a wireless power transmission device for use in an EAD thruster designed for vertical takeoffs. In 2018, Drew and Lambert [
17] developed and successfully launched a small (2 cm × 2 cm) and lightweight (30 mg) microrobot with a thrust-to-weight ratio of 10. Also in 2018, researchers [
2] from the Massachusetts Institute of Technology successfully tested the world’s first plasma-propelled drone (as shown in
Figure 1), demonstrating that plasma propulsion technology is suitable for drone flight. In 2019, Chirita and Ieta [
5] completed their first unmanned flight with a rotary ion engine (RIE), which used ionized wind to blow the rotor blades.
The fixed-wing plasma aircraft flown by the Massachusetts Institute of Technology, which featured a high-voltage power converter and an EAD thruster, is shown in
Figure 2. Compared with traditional aviation propellers (piston and jet engines), plasma propulsion vehicles exhibit low fatigue, are ultra-quiet, and have long service lives owing to their simple structure and lack of moving parts. In addition, electrical propulsion does not contaminate the environment. Compared to the typical thrust-to-power ratio of conventional aerospace propulsion technology (3 N/kW), the predicted thrust-to-power ratio of the EAD aircraft (50 N/kW) was substantially higher in the laboratory [
18]. However, the low propulsion and thrust density of EAD thrusters have constrained their application in small UAVs. Accordingly, in recent years, many researchers have worked on designing lightweight and high-voltage power converters (HVPCs) [
19,
20,
21] that will improve thruster propulsion performance.
The fundamental component of an EAD thruster, also called a thruster unit, is an asymmetric pair of electrode rods. One electrode is called the emitter and the other is called the collector; they are separated by an electrode gap d and connected to a high-voltage source. The collector typically has a larger frontal area, such as that of a cylinder or airfoil, and the emitter typically has a very small radius of curvature. In early studies, such as those carried out by Moreau and Benard [
22], the effects of the electrode gap d, air humidity, and cylinder diameter of the grounded electrodes on the thruster were investigated. Subsequent experiments with EAD thrusters have shown that the thrust-to-power ratio and thrust density of EAD thrusters in small unmanned aircraft can be increased by using series and parallel arrangements [
23,
24]. In this context, a series connection refers to a plurality of thruster units connected head to tail and placed in the same plane, which form a multistage thruster. A parallel connection refers to a plurality of thrusters, with emitters in one plane and collectors in the other; the two planes are parallel. Examples of thruster unit geometries include wire to cylinder [
25], wire to airfoil [
17], and pin to mesh [
16]. Early research primarily focused on the wire-to-cylinder geometry, but switching to an airfoil significantly minimizes the drag caused by the plasma flow. This finding was confirmed by Belan et al. [
26], who also investigated the effect of using collectors with different airfoil shapes on the propulsion performance. Their findings demonstrated the necessity for further systematic studies on electrode geometries. In addition, Belan et al. [
27] tested the effect of emitter density on the plasma thruster performance and determined that an optimal collector–emitter spacing ratio exists.
In previous performance optimization experiments, two primary electrode distribution patterns were examined for their effects on the electrical parameters (e.g., voltage and polarity) and geometrical parameters (e.g., the thruster unit spacing s between two parallel emitters or collectors and the electrode gap d between the emitters and the collectors). As shown in
Figure 3, when the emitter is located above the collector, it is called the collinear distribution, which is used in the majority of studies [
18]. When the emitter is placed in the middle of succeeding collectors, it is called the staggered distribution, which operates according to the same physical principles [
17]. Moreau’s experimental data indicated that the current-to-thrust conversion is more efficient for the staggered electrode distribution [
3]. This phenomenon was confirmed in Monrolin’s study [
25]; however, he did not further theorize about its cause.
The effects of a single-variable thruster unit geometry and the coupling of two to three variables on the propulsion performance have been well-studied. However, when the collectors are airfoils, the plasma thruster design requires accounting for additional variables and is characterized by nonlinearity and principal unknowns. As mentioned above, researchers have not yet conducted a comprehensive and detailed investigation into the impact of the geometric configuration of the thruster on the performance of the plasma thruster due to the various challenges involved in conducting physical experiments on plasma thrusters, which is what we are considering implementing. To simplify the thruster design, a proxy model is required. A proxy model is a mathematical or statistical model used to represent, simulate, or predict a system, phenomenon, or process. It is typically a simplification or approximation of a real system that captures its key factors and interactions for simulation purposes. It can contribute to understanding problems, making predictions, and optimizing decisions [
28].
The main purpose of this article is to conduct a detailed study of the five most significant parameters that affect the performance of plasma thrusters, aiming to achieve performance optimization. This is beneficial for improving the thrust density and thrust-to-power ratio issues currently faced by plasma thruster drones, thereby enhancing their utility. A thruster unit configuration consisting of a staggered distribution of emitters and collectors was tested because it produces less parasitic drag and yields superior aerodynamic effects. We aimed to select a reasonable sampling method to collect the sample points and establish an approximate proxy-based model, thereby ensuring that the output response values of the model would achieve the desired accuracy. Finally, we aimed to optimize the performance of the thruster unit through a multi-objective genetic algorithm.
The remainder of this paper is organized as follows.
Section 2 briefly elaborates on the one-dimensional EAD theory, explaining why electrodes in a staggered distribution can decrease drag and provide a useful parametric design space.
Section 3 introduces the customized experimental setup.
Section 4 explains the optimization design method, including the sampling method, proxy model, and multi-objective genetic algorithm.
Section 5 presents the predicted results of the proxy model and the optimized parameters, followed by a comparison and validation. Finally,
Section 6 contains the conclusions of the study.
3. Optimization of the Design Methods
3.1. Experimental Design
The optimal design process for the target parameters in this study was divided into three steps: experimental design, proxy model building, and genetic algorithm optimization. The experimental design employed Latin hypercube sampling (LHS), which has the significant advantage of being able to construct a high-precision analytical model that is very similar to the actual model but with a smaller number of experimental samples, thereby substantially reducing the computational load. This method can be used to select the sample points more evenly within a sampling range. The choice of sampling strategy has a certain impact on the performance and generalization ability of the model, but if the ranges of parameter values are too large or too many samples are collected, the efficiency of the experiment will decrease and the experiment will be prolonged. However, too few samples may cause the proxy model to overfit or underfit, increase the instability during the sample training process, and reduce the generalization ability.
LHS is an equal-probability stratified sampling method that randomly generates samples within a range of values for each parameter [
30,
31]. For this study, LHS indicated that the sample size should be more than three times that of the number of design parameters. There were five design parameters for the thruster; therefore, a minimum of 15 tests should be performed. However, we chose to perform three times this number of tests to obtain more accurate results. Thus, 45 samples were collected. According to the sampling principle, the independent space of each parameter was evenly divided into 45 parts, and then a parameter point was randomly selected in each part of the independent space. A random combination of five sets of design parameters yielded 45 sample points. When a comparatively limited number of samples was used for probing, these 45 sample points performed better because they were highly homogenized, covered all regions of the parameter space, and were well-covered.
Table A1 lists the 45 sample points selected for this study.
3.2. Proxy Model Building
Proxy models based on artificial neural networks are frequently used to establish relationships between inputs and outputs [
32]. Neural networks attempt to (at least partially) simulate the structure and function of the brain and nervous system in living organisms [
33,
34]. Accordingly, they consist of interconnected neurons that can simulate the signal transmission that occurs between neurons in the human brain, thereby automating the learning and decision-making process. Thus far, the error backpropagation algorithm combines several optimization algorithms, such as gradient descent, to continuously adjust weights and biases based on training data, thereby enhancing the model’s accuracy and generalization capability [
33]. Therefore, it is considered one of the most successful learning algorithms for feedforward neural networks. BP neural networks consist of a forward propagation process for the working signal and a backward propagation process for the error signal. BP neural networks process
m inputs and produce
n outputs, and they typically consist of three layers: input, hidden, and output. However, several hidden layers can be placed between the input and output layers. A three-layer BP network is capable of mapping from
m-dimensions to
n-dimensions. The numbers of neurons in the input and output layers are fixed. The number of neurons in the hidden layer can be determined according to the empirical formula
where
q,
m, and
n are the numbers of neurons in the hidden, input, and output layers, respectively, and
a is a conditioning constant between 1 and 10. In this study, the input layer defined five input neurons (
d,
s,
t,
c, and
θ), and the output layer defined two output neurons (
Te and
P). Based on Equation (22), the preliminary determination of the number of neurons in the hidden layer was between 4 and 13. The best structure was determined via trial and error. Testing the performance of the neural network in each of these cases revealed that the performance parameters of the plasma thruster can be effectively predicted using six hidden-layer neurons.
The dataset, consisting of 45 sample points, was divided into three sets: training, validation, and test. The training, validation, and test sets contained 31, 7, and 7 sample points, respectively. To speed up the training process, we used the purelin linear transfer function and the Levenberg–Marquardt algorithm as training methods for the neural network. Compared to traditional gradient descent algorithms, the approach typically converges more quickly to local optimal solutions. It is beneficial for both the training speed and performance of neural networks, while mitigating issues, such as gradient vanishing and exploding. The dataset was normalized to linearly scale the input values within the range (0, 1]. The network training error was the average relative error between the predicted and actual values of the sample points.
Because the selection of the neural network structure, initial connection weights, and thresholds has a significant impact on network training but is difficult to obtain accurately, we used a genetic algorithm coupled with a neural network to initialize the initial weights and thresholds and perform the optimization. The basic idea was to use an individual to represent the initial weights and thresholds of the network, calculate the average relative error between the predicted and actual values of the test set as the output of the objective function, and then calculate the fitness value of the individual. Through the selection, crossover, and mutation operations, we searched for the individual that represented the optimal initial weights and thresholds of the BP neural network. This enabled the optimized BP neural network to make better sample predictions.
The optimization steps were as follows. First, the population was initialized, and the individuals were encoded using binary coding. Each individual was represented by a binary string consisting of four parts: the connection weights between the input and hidden layers, the threshold of the hidden layer, the connection weights between the hidden and output layers, and the threshold of the output layer. Each weight and threshold was encoded using
M bits of binary code, and the coding of all of the weights and thresholds was concatenated to form the encoding of each individual. The fitness function used a ranking fitness assignment function, the selection operator used random traversal sampling, and the crossover operator used the simplest single-point crossover operator. Mutations occurred with a certain probability and resulted in a certain number of mutated genes. The mutated genes were randomly selected. If the encoding of the selected gene was 1, it was changed to 0; otherwise, it was changed to 1. The running parameters for the genetic algorithm are listed in
Table 2.
After establishing the appropriate initial weights and thresholds of the neural network using the genetic algorithm, the Bayesian regularization method was used to train the BP neural network. The Bayesian regularization approach can achieve a strong generalization performance by automatically choosing the optimal regularization parameter, thus avoiding the issues of local minima and overfitting, which are typical of BP neural network algorithms.
3.3. Multi-Objective Genetic Algorithm NSGA-II Optimization
In this study, the thrust, thrust-to-power ratio, and thrust density were selected as the thrust performance indicators. Because of this selection, there were multiple optimization objectives that needed to be satisfied simultaneously. Therefore, single-objective optimization methods were not suitable for solving this problem. Traditional multi-objective solution methods primarily include the weighted sum, ε-constraint, and equality constraint methods. These methods essentially transform a multi-objective problem into a single-objective problem; however, they can only obtain one optimal solution at a time, which is inefficient. NSGA-II, a multi-objective genetic algorithm, is a non-dominated sorting algorithm well-suited for circumventing this issue. The algorithm consists of the following steps. First, the initial population (parent population) of the design variables is randomly generated. Second, after evaluating and determining the ranks of individuals in the population based on the non-dominance criteria and the crowding distance, an intermediate population is generated via the selection, crossover, and mutation operators. Third, by combining the intermediate population with the initial population, the ranks of the members are determined, and the next generation is selected from this population based on the fitness criteria. This process comprises a certain number of iterations that eventually generate a Pareto frontier. The goal is not to find the optimal solution for each sub-objective, because guaranteeing that all parameters are optimized under ideal conditions in the presence of multiple conflicting objectives is difficult. There is not simply one such solution but rather a continuous set of infinite solutions, which form the Pareto frontier.
In engineering practice, it is essential to seek a unique solution. Merely calculating the Pareto frontier is insufficient. Therefore, establishing an individual evaluation model for the Pareto frontier to assess the individuals within it is necessary. To address differences in the physical meaning and magnitude of the objective functions, this study normalizes each set of non-dominated solutions in the Pareto frontier. Subsequently, the Analytic Hierarchy Process is utilized to narrow down the Pareto solution space, thus determining the optimal compromise solution for the thrust, the thrust-to-power ratio, and the thrust density. The process involves a thorough analysis of the essence of complex decision problems, influencing factors, and their internal relationships. It utilizes limited quantitative information to mathematize the decision-making process, offering a straightforward approach to decision making for complex problems with multiple objectives, criteria, or other unstructured characteristics.
The analysis steps are as follows. The first involves categorizing the parameter objectives, performance indicators, and Pareto frontiers from high to low in terms of their relationships as the objective layer, the criterion layer, and the scheme layer. For adjacent layers, the higher one is called the objective layer, and the lower one is called the factor layer. Starting from the second layer of the hierarchy model, judgment matrices are constructed for each set of factors belonging to the same layer as the one above until the lowest layer. The influencing factors from the criterion layer to the scheme layer have already been identified in the Pareto front, so only the consistent matrix method is needed to construct judgment matrices from the objective layer to the criterion layer. The consistent matrix method does not compare all factors together but rather pairwise, minimizing the difficulty of comparing various factors with different properties to enhance accuracy. Next, consistency indices are used to check the judgment matrices; when the consistency ratio is <0.1, the inconsistency level of the judgment matrix is considered within an acceptable range, indicating satisfactory consistency for it to be used as a judgment matrix. Finally, the judgment matrix is transformed into weights for performance compromise indicators and allocated to the normalized coordinates of each solution group in a Cartesian coordinate system. The distances from the normalized Pareto front for each individual to the coordinate origin are calculated, and then the parameter combination corresponding to the solution group with the greatest distance is selected as the optimal parameter scheme.
4. Experimental Setup
Figure 6 shows the customized experimental device and its description. Based on the outcomes of the LHS procedure described above, 45 airfoils were 3D-printed. The outer surface of each airfoil was then covered with aluminum foil, which served as the collector electrode. The geometric parameters of all wing shapes complied with the NACA symmetric airfoil standards. The length of the airfoil collector
b was 150 mm. A tungsten wire with a diameter of 50 μm was used as the emitter. It passed through a preset circular hole on the emitter slider and was connected in parallel to a high-voltage power supply, which provided a maximum voltage of 50 kV with an adjustment accuracy of 0.01 kV. The power supply was also equipped with an overcurrent-protection mechanism. The output cable of the power supply was fixed to a stable platform and carefully connected to the emitter to ensure that the weight and stiffness of the cable did not affect the measurements. The emitters were held straight. Once the wire was slightly bent, it created higher localized electric fields, which changed the uniformity of the discharge and may have led to premature arc discharge.
As shown in
Figure 6, the customized experimental setup was utilized to perform direct thrust and electric field measurements. The device consists of a support framework, a two-stage thruster support structure, and two-stage thrusters. Each thruster stage was composed of five emitters and six collectors arranged in a staggered distribution. The emitters of the first stage were installed at the top of the device. The collectors of the first stage and the electrode assembly of the second stage could slide together within the pillars of the device structure, allowing for adjustments in the electrode gap
d and thruster stage spacing
θ. To adjust the thruster unit spacing
s, the structures supporting both the emitters and the collectors were equipped with grooves that had smooth surfaces, which enabled the sliding blocks used for installing the electrodes to move smoothly within the grooves. Scale rulers were placed near the grooves of the pillars and structures supporting the test device, which facilitated the adjustment of the distance parameters during each sampling. The airfoil collectors and supporting airfoil sliding blocks were connected using quick connectors, allowing for a convenient method to replace the airfoil electrodes after each sampling. The plastic clip was used to fix the position of the three-stage electrode below after adjusting the distance. All components of the experimental setup were fabricated using 3D-printing technology with insulating materials, such as polylactic acid (PLA). The adjustment accuracy for the electrode gap d, the thruster unit spacing s, and the thruster stage spacing θ was 1 mm.
In this study, two precision measurement balances with an accuracy of 0.01 g were symmetrically placed on both sides of the device to prevent downward ion wind from interfering with the measuring instruments. The weight range of each balance was 0–10 kg. To measure the current passing through the thruster, an ammeter that had an accuracy of 1 μA was connected in series between the collector electrode and the ground wire. The experimental procedure was as follows. Firstly, select the airfoil with a corresponding thickness and chord length before the experiment starts, then adjust the three distance parameters for installing the airfoil, and start the experiment. During measurement, fix the power supply voltage at 20 kV. Use the sum force from two balances minus their initial readings as the thrust indication at the current moment. The experiment involves certain uncertainties, which are typically greater when the ion wind is stronger, resulting in more significant fluctuations in balance readings. This uncertainty cannot be precisely measured. So, the experimental procedure involves setting the high voltage to the desired value for 5 s and then recording average data within 30 s to reduce experimental errors. All experiments were conducted at a voltage at 20 kV, a temperature between 19 and 21 °C, an atmospheric pressure between 0.97 and 1.01 bar, and a relative humidity between 28 and 32%. The specific experimental parameters are summarized as shown in
Table 3.
6. Conclusions
This study predicted and optimized the influence of the main geometric parameters of plasma thrusters on their performance. The main geometric parameters included the electrode gap d, the electrode unit spacing s, the thruster stage spacing θ, the maximum thickness of the collector airfoil t, and the chord length of the collector airfoil c. The staggered distribution was adopted for the thruster’s configuration. The advantages of the staggered distribution in reducing drag were briefly analyzed from an aerodynamic perspective based on the one-dimensional EAD theory. To reduce the experimental and computational costs, a proxy model was established using a BP neural network that was improved via a genetic algorithm. By customizing the measurement of thrust and current, a more detailed analysis of the impact of the EAD thruster performance was carried out based on 45 sets of experimental data. Using the NSGA-II multi-objective genetic algorithm for optimizing the prediction results, we applied the Analytic Hierarchy Process to enhance the Pareto solution set through balancing and selection. This process helped us identify seven geometric parameter combinations with superior overall performance, each reflecting distinct performance preferences. However, this study suffers from a small sample size, insufficient data, and potential biases in the predictions of the experimental setup and proxy models, which have implications and limitations for the conclusions of this paper. In addition, the optimized performance may still not meet the performance requirements of plasma aircrafts.
In the staggered distribution of 5E/6C (five emitters/six collectors), an increase in thrust led to a decrease in the thrust-to-power ratio, which was consistent with the experimental and theoretical results. Within the parameter range investigated in this study, reducing the electrode gap d and the electrode unit spacing s resulted in a higher thrust density; however, increasing both values improved the thrust-to-power ratio. Increasing the thruster stage spacing θ effectively suppressed the generation of a reverse-ion wind, and its impact on the thrust density and thrust-to-power ratio was the opposite of that on the thrust. Thinner and smaller airfoil collector electrodes were beneficial for increasing the thrust. Changing the thickness t and chord length c had a negligible impact on the working volume of the plasma thruster; therefore, the variation in the thrust density followed a trend similar to that of the thrust.
To determine the predicted values of the maximum thrust, the maximum thrust density, and the maximum thrust-to-power ratio, geometric parameter combinations were obtained using a single-objective genetic algorithm. Compared to the average values of the original samples, the thrust increased by 225.23%, the thrust density increased by 357.14%, and the thrust-to-power ratio increased by 195.31%.
The normalized Pareto solution set provided a balanced performance across multiple objectives. Therefore, it discovered a range of solutions that performed well across different objectives rather than optimizing only for a single objective. Using the NSGA-II genetic algorithm to acquire the Pareto optimal solution set, a decision was made utilizing the Analytic Hierarchy Process to pinpoint seven representative optimal parameter combinations. Subsequent researchers could select an appropriate solution based on the drone’s payload capacity and endurance requirements. For instance, the first solution set was suggested for prioritizing payload capacity improvement while ensuring decent endurance. It was beneficial to address the current performance shortcomings of plasma propulsion drones and provide greater flexibility and choice for subsequent researchers.
Follow-up research could focus on the use of more accurate proxy models and experimental data, as well as the use of higher-thrust thruster configurations, such as the multi-staged ducted thrusters proposed by Gomez-Vega [
35].